#  Inferring the drivers of language change using spatial models

James Burridge  School of Mathematics and Physics, University of Portsmouth,
Portsmouth PO1 3HF, United Kingdom  Tamsin Blaxter  Gonville and Caius
College, Cambridge CB2 1TA, United Kingdom

###### Abstract

Discovering and quantifying the drivers of language change is a major
challenge. Hypotheses about causal factors proliferate, but are difficult to
rigorously test. Here we ask a simple question: can 20th Century changes in
English English be explained as a consequence of spatial diffusion, or have
other processes created bias in favour of certain linguistic forms? Using two
of the most comprehensive spatial datasets available, which measure the state
of English at the beginning and end of the 20th century, we calibrate a simple
spatial model so that, initialised with the early state, it evolves into the
later. Our calibrations reveal that while some changes can be explained by
diffusion alone, others are clearly the result of substantial asymmetries
between variants. We discuss the origins of these asymmetries and, as a by-
product, we generate a full spatio-temporal prediction for the spatial
evolution of English features over the 20th Century, and a prediction of the
future.

##  1 Modelling language evolution

Modelling the collective behaviour of systems of sentient agents, from flocks
of birds [1, 2], to economies [3], cities [4] and languages [5, 6, 7], is
attractive, but not easy. To quote Emanuel Derman [8], an early adopter of
financial modelling,

In physics there may one day be a Theory of Everything; in finance and the
social sciences, you’re lucky if there is a usable theory of anything.

The problems which arise when trying to model social systems depend to some
extent on the system in question, but share a good deal in common. In many
cases we have only one realisation of the change process which we want to
understand, and the data available about the change may be sparse. Because we
cannot observe the thought processes and motivations of the agents, we
caricature them with a simple model. However, many different models may fit
observations similarly well. This problem recedes if we have more data and
more realisations, which allow us to apply Statistical Learning methods [9],
such as cross validation, to rigorously select the ``best model'' by
maximizing its ability to predict unseen data. If this is not possible,
because we have only one realisation of a particular change process, then
another approach is to start with some simple default theory, or null model
[10], and then seek to determine whether it is sufficient to explain
observations. We may also consider an alternative model with incrementally
increased complexity and minimal additional assumptions, the aim being to
understand what is missing from the null model, while avoiding unjustified
assumptions about what should replace it.

The broad question which interests us is: what are the processes which drive
the evolution of languages, and how should they be represented in a coarse
grained spatial model? A longer term aim is to be in possession of sets
stochastic equations, derived from simple assumptions about speaker-level
behaviour, which describe how language evolves over space and time, and
amongst different social groups. Such equations should calibrate to available
data on the language state up until the present, then allow us to understand
historic changes, and make future predictions. Here we take a step toward
addressing these aims, by fitting a coarse grained spatial model to two large
scale spatial datasets which capture the states of various English language
features at the beginning and end of the 20th Century.

The mathematical modelling of language evolution, often using methods of
statistical physics, is now an established field. Many aspects of language
change have been studied using non-spatial models [11, 5, 12, 13, 14], and an
important focus has been to understand the replacement of one linguistic
feature by another, typically following the ubiquitous temporal S-curve [15,
16]. Such changes can be driven by numerous mechanisms, including those which
favour particular variants (e.g. regularization [17], social biases [18], or
exogenous factors [16]), by stochastic effects [19, 20], and by biases driven
by linguistic variation over the age spectrum [12, 13, 21, 22, 23]. In this
work we allow speakers a learning bias toward individual linguistic features,
without explicitly modelling any particular mechanisms which might generate
this. We seek only to determine whether some form of bias existed, given the
observed changes.

Spatial modelling of language change is a topic of growing interest [24, 25,
26, 27, 28, 7, 29, 30, 31]. Whereas non-spatial linguistic data is plentiful,
detailed spatial datasets are rarer. Their construction typically requires
either large scale collaborative effort [32, 33, 34], or viral success online
[35, 36, 37]. Spatial modelling presents additional challenges in the form of
greater analytical and computational complexity, particularly when matching
models to real world spatial language distributions. There is also long
standing interest within the quantitative linguistics community in the
statistical analysis of spatial, social and temporal variations within real
linguistic domains [34, 38, 39].

Recent spatial models [7, 29, 30, 31], which accurately model geography and
population distributions, suggest that the shapes of linguistic domains and
the locations of cities and towns can have a partially predictable influence
on the evolution of language. In fact, the importance of geography has been
well known to dialectologists for a long time [18, 40]. From a mathematical
perspective, this influence arises from the impact of geography on the shapes
of boundaries (isoglosses) between alternative linguistic features, which
behave like two dimensional bubbles. The existence of this surface tension
effect has recently been tested using historical English dialect data [31]. We
therefore believe that space and spatial processes must be part of models
which seek to understand language change processes. Here we provide, to our
knowledge, the first detailed spatial model which matches the evolution of
real linguistic features within an extended spatial domain (England in the
20th Century), accounting for realistic migration patterns, and plausible
learning behaviour.

##  2 Data and existing theories

###  2.1 Linguistic variables and variants

For the purpose of academic study, languages may be broken down into distinct
components: single units of sounds (phones), rules for combining sounds
(phonology), words (the lexicon), rules for word construction (morphology),
and rules for sentence construction (syntax). A language may be viewed as a
complete specification of all its components, and language change as the
process of progressively exchanging components for alternatives. We refer to a
language feature for which there are a set of alternatives as a variable, and
to the alternatives as variants. These might be different words for the same
object or idea, different sounds playing the same role in certain groups of
words, or alternative syntactic rules. As an example, the word for the season
after summer historically had variants autumn, backend and fall in England
[32]. The fundamental quantities which we build models of are the relative
usage frequencies of variants at different locations (or regions) in space.
Simultaneous spatial variations in many linguistic variables can create
distinctive regional dialects.

###  2.2 The Survey of English Dialects and the English Dialect App

We model the spatial evolution of a set of linguistic variables over the 20th
Century in England, the largest and most densely populated nation within the
British mainland. As initial and final conditions, we use data for the same
variables from two surveys: the Survey of English Dialects (SED) [41] and the
English Dialects App (EDA) [35, 36].

The SED was carried out in the 1950s at 313 localities. With 986 respondents,
it is a small survey by modern standards. The localities were selected to
provide a relatively evenly-spaced sample across England, but included almost
no speakers from urban areas. This was a deliberate choice; the survey sought
out the most conservative demographic (typically farm-labourers born in the
1870s and 80s) to capture the network of ``traditional'' dialects before they
disappeared. Methodologically, the SED was rigorous, making few allowances for
speed or convenience at the expense of data volume and quality: fieldworkers
visited respondents in person and recorded their responses to over a thousand
questions in narrow phonetic transcription (and later also on tape). Despite
some criticisms arising from systematic differences in how fieldworkers asked
questions or transcribed responses [42, 43], the SED is a rich and relatively
trustworthy source of data and has provided the material for more than half a
century's worth of secondary analysis.

The EDA represents a new generation of dialect data collection methods, using
digital technologies to reach a large number of respondents very cheaply. Over
50,000 speakers answered 26 questions about their usage through a smartphone
app in 2016, and all but one question duplicated variables surveyed in the
SED. The EDA exemplifies a different set of trade-offs between data quality,
quantity and speed. It is possible to mine the SED data for patterns across
many responses and examine variables which were not specifically targeted,
whereas the EDA required speakers to decide between predefined categories,
limiting analysis to 26 variables and predetermined sets of variants. Instead,
the EDA's major advantage is its 45,287 speakers in 39,590 locations within
the region covered by the SED. It is also more representative, covering urban
areas and all demographics, although, because it was carried out via
smartphones, it is somewhat skewed towards younger and more affluent speakers.

In spite of methodological differences, we argue that comparison of the SED
and EDA is a viable way to explore language change in 20th century English
English. Laboratory perception tests [36] suggest that speakers are largely
able to discriminate between the variants of the phonetic variables elicited
in the EDA. Spatial patterns in the results of the EDA, when compared to the
SED, show that levelling (loss of variants) and isogloss (linguistic boundary)
movements are in line with our prior understanding of the mechanisms of
language change. However, because typical SED respondents were old, while EDA
respondent were younger, we are not simply mapping the 60 or so years of
language change between the median dates of the two language surveys. If we
assume that speakers' linguistic norms are set early in life and change little
in adulthood (the underlying assumption of much historical linguistic [44] and
sociolinguistic [45] work on language change), the relevant comparison is
instead speakers' dates of birth: the median date of birth is the 3rd of May
1881 for SED respondents, and the 20th of July 1983 for EDA respondents,
implying that we are looking at around 100 years of change.

Because urban areas are under-represented in the SED, we must consider whether
city varieties are likely to have differed substantially from the varieties in
the regions surrounding them. In that case SED maps would be missing islands
of highly divergent usage in places where spatial models [7, 29, 30] suggest
that population density gradients should preserve those distinctions or
encourage the spread of urban variants. However, the rapid urbanisation of the
industrial revolution had started only a couple of generations before the SED
speakers were born; less than the three generations required for new dialect
formation, according to standard models [46], and certainly not long enough
for urban varieties to have substantially diverged from their regional inputs.
This is borne out by evidence. For example, Coates [47], surveying research on
traditional Bristol English, notes: ``what can be shown to distinguish
Bristolian at all linguistic levels from other dialects of the region is
relatively little''. Clearly city varieties at the time of the SED will have
had distinguishing features, just as varieties at any point in space differ in
some ways from those nearby; but we can assume that they tended to agree in
most respects, and so their omission is not likely to be any more problematic
than a situation where any other location happened to be lacking in samples.

###  2.3 Current hypotheses about linguistic changes

The overwhelming story of change in English English dialects over the 20th and
early 21st centuries is one of loss of diversity. In the words of Britain
[48]:

There has been such considerable and ongoing dialect attrition that the
language use reported across the country by Ellis's survey of 1889 seems, in
many cases and in many places, quite alien to that spoken just over one
hundred years later.

In place of smaller traditional dialects we find distinctions at larger,
regional levels [49, 50, 51], and much research has focused on these
regiolects and the processes of regional dialect levelling [52, 53, 54, 55,
49] which generated them. Reductions in geographical linguistic differences
are not only due to convergence to one of several local variants but also to
geographically widespread adoption of common innovations [50]. These changes
are part of a near-universal pattern across the traditional varieties of
Europe [56, 57, 58].

It is believed that interactions between individuals can lead to a decline in
linguistic variation via accommodation, where conversation partners adjust
their speech to better match each other, and by child learners acquiring
accommodated forms [59, 60, 61, 62, 63]. However, in certain social contexts,
children may also learn variants directly from mobile outsiders [64]. The fact
that accommodation is mediated through interactions between speakers has lead
some linguists to conclude that the primary drivers in the decline of
linguistic diversity are travel, commuting and migration [59, 65, 49, 55, 66,
50, 67]. Other potential drivers discussed include changes in social network
structure [59, 68, 65, 49, 50, 55, 66], the age structure of the community
[69], the influence of mass media [58, 70], normative attitudes and education
[71, 72, 58] and relatedly the salience and stereotype status of particular
variants [73, 55], identity factors [74], the informalisation of public life
[58] and socio-economic forces [55]. Purely linguistic internal factors such
as structural regularity, functional economy, or naturalness may determine
which variants win out in the levelling process [59, 67, 49, 55] or they may
not be relevant at all [71, 50]. Arguments have also been made for the central
importance of idealogical factors, in particular strong normative attitudes
towards the standard, alongside mobility and contact [74, 72, 75]; parallel to
these arguments, it has been suggested that strong alignment of speaker
identities with the local community may be enough to check the levelling
process and so preserve a distinct local variety [58, 51].

There are clearly numerous potential drivers of linguistic change, and it is
beyond the scope of this work to quantify their relative importance. It may be
that driving processes which operate in different ways at small scales (at the
level of motivations, interactions and contacts between individuals) yield
similar or identical terms in evolution equations for coarse grained
population averages, making the task of inferring the importance of individual
effects impossible using such models. However, macroscopic models have been
used to infer different classes of driving mechanism from non-spatial
linguistic time series, specifically exogenous (population level) vs
endogenous (individual level) drivers [16]. We take a similar approach here,
but our two classes are spatial processes (migration, movement), and processes
which introduce asymmetry with respect to variants (ideology, internal
linguistic effects, social prestige, normative bias toward a standard etc.).
We do this by defining coarse grained evolution equations which account for
both movement and biased copying. We then fit our model to the initial and
final conditions supplied by the SED and EDA data, to determine the relative
importance of each process. Since movement is often viewed as the primary
driver of change, we view variant symmetric dynamics as our null model, and
determine the extent to which it can explain the observed changes, before
breaking variant symmetry by allowing for biasing factors.

##  3 The model

We consider a spatial domain divided into $L$ cells, each containing
(approximately) $N$ speakers. The centroids of these cells are written
$\boldsymbol{\mathbf{r}}_{1},\boldsymbol{\mathbf{r}}_{2},\ldots,\boldsymbol{\mathbf{r}}_{L}$.
Consider a linguistic variable with $q\in\\{1,2,\ldots\\}$ variants, and let
the relative frequencies with which these variants are used within cell
$\boldsymbol{\mathbf{r}}$ define the frequency vector
$\boldsymbol{\mathbf{f}}(\boldsymbol{\mathbf{r}})=(f_{1}(\boldsymbol{\mathbf{r}}),\ldots,f_{q}(\boldsymbol{\mathbf{r}}))^{T}\in\Delta^{q}$,
where $\Delta^{q}$ is the $q$-dimensional simplex and T denotes transpose. We
have suppressed time dependence for brevity. The cells we use in our analysis
are Middle Layer Super Output Areas (MSOAs), a set of $L=7,201$ geographical
polygons, each with a similar number of residents, used for reporting census
data in England and Wales [76]. The mean MSOA population in 2011 was $7,787$,
and the modal area of English MSOAs is approximately $0.29\text{km}^{2}$,
corresponding to a few tens of streets within a densely populated city.
Interactions with areas outside the English borders, for which we lack survey
data, are neglected. We justify this simplification on the basis that
population densities are low at the Scottish and Welsh borders, the English
make up the great majority ($\approx 87\%$) of the British mainland
population, and Wales and Scotland possess distinctive cultural and linguistic
identities. We assume that much longer range interactions, for example from
American TV, cinema, and music, will be subsumed into biases in the learning
function (see section 3.3) which determines the adoption probabilities of
different variants by young speakers.

###  3.1 The language community

We refer to the language that a speaker is exposed to as their linguistic
environment, and assume this environment consists predominantly of voices from
their own and nearby cells. We allow for biasing factors by weighting the
influence of variants using a set of biases which affect the learning process
(see section 3.3). We capture the local environment of a speaker in cell
$\boldsymbol{\mathbf{r}}$ using the community influence matrix, $W$, where
$W(\boldsymbol{\mathbf{r}},\boldsymbol{\mathbf{r}}^{\prime})$ gives the
influence that speakers from $\boldsymbol{\mathbf{r}}^{\prime}$ have on
speakers from $\boldsymbol{\mathbf{r}}$. The matrix is stochastic, meaning
that it is square with unit row sums. We define the community frequency vector
of a variable as follows

$\widehat{\boldsymbol{\mathbf{f}}}(\boldsymbol{\mathbf{r}})\triangleq\sum_{\boldsymbol{\mathbf{r}}^{\prime}}W(\boldsymbol{\mathbf{r}},\boldsymbol{\mathbf{r}}^{\prime})\boldsymbol{\mathbf{f}}(\boldsymbol{\mathbf{r}}^{\prime})$
(1)

where $\triangleq$ denotes a definition. We assume that local influence
strengths are mediated by two factors: physical separation, and population
distribution. A definition which incorporates both is

$W(\boldsymbol{\mathbf{r}},\boldsymbol{\mathbf{r}}^{\prime})=\frac{\exp\left(-\frac{|\boldsymbol{\mathbf{r}}-\boldsymbol{\mathbf{r}}^{\prime}|^{2}}{2\sigma^{2}}\right)}{\sum_{\boldsymbol{\mathbf{r}}^{\prime\prime}}\exp\left(-\frac{|\boldsymbol{\mathbf{r}}-\boldsymbol{\mathbf{r}}^{\prime\prime}|^{2}}{2\sigma^{2}}\right)},$
(2)

where we call $\sigma$ the interaction range. To illustrate the competing
effects of interaction range and population distribution we consider the one
dimensional small-cell limit of the influence matrix, which takes the form of
a kernel giving the influence of location $x$ on location $x_{0}$

$\tilde{W}(x_{0},x)=\frac{\rho(x)\exp\left(-\frac{(x-x_{0})^{2}}{2\sigma^{2}}\right)}{\int\rho(y)\exp\left(-\frac{(y-x_{0})^{2}}{2\sigma^{2}}\right)dy}$
(3)

where $\rho(x)$ is population density at $x$, inversely proportional to cell
size at this position. Population density is automatically accounted for in
definition (2) because higher density areas have more cells. Figure 1 shows
the spatial distribution of influences on three speakers who live in or
between two cities with radii and separation similar to Birmingham (population
$\approx 10^{6}$) and Coventry (population $\approx 3\times 10^{5}$).

Figure 1:  Effects of towns and cities on the spatial distribution of the 1D
continuum influence kernel $\tilde{W}(x_{0},x)$, given by (3). Brown curve is
proportional to population density. Coloured curves show the influences on
speakers in the city centres and at the midpoint, when the interaction range
is $\sigma=10$km.

Consider a speaker in central Coventry ($x_{0}=15$). Most of the influence on
them comes from within the city, because there are relatively few people to
interact with outside. This makes their effective interaction range smaller. A
speaker in central Birmingham ($x_{0}=-15$) will have a wider range of
influence because their city is larger. For a speaker in between the two
cities ($x_{0}=0$), influence comes from both, with the larger city dominating
(see blue curve in Figure 1). In the limit of very large cities (or other
regions of approximately uniform population density) the influence
distribution becomes normal so that $\approx 63\%$ of influence comes from
within $\sigma$km.

We assume that speakers' linguistic states are fixed in their youth (see
section 3.3), so the interaction range should be chosen to capture the
community state observed by a younger speaker. An empirical measure of such
separations is provided by school travel distances; secondary education has
been compulsory since 1918, so we would expect children to regularly interact
with others from their school catchment. In the early 21st century secondary
pupils travelled, on average, $\approx 5.5$km to school [77]. We will call
this the mean catchment radius, denoted $\mu_{R}$. The number of state schools
has declined somewhat over the 20th century; in 1951 there were 5,900 vs.
4,072 in 2010 [78]. This implies that average catchment radii have increased,
with 1951 radii being $\approx 83\%$ early 20th century radii. Scaling 2010
travel distances by this factor yields an average travel distance of $\approx
4.5$km. Suppose we model the distribution of displacements between home and
school as a bivariate normal distribution

$f(x,y)=\frac{e^{-\frac{x^{2}+y^{2}}{2\omega^{2}}}}{2\pi\omega^{2}},$ (4)

then the mean home-school distance is $\mu_{R}=\sqrt{\tfrac{\pi}{2}}\omega$
and the mean distance between two homes is
$2\omega=\sqrt{\tfrac{8}{\pi}}\mu_{R}\approx 1.6\mu_{R}$. This suggests a
typical distance between interacting speakers is between 5 and 10km. Since
$\sigma$ functions as an upper bound in cities, we take our interaction range
as the upper limit of this interval, so $\sigma=10\text{km}\approx
6\text{miles}$.

We assume that speakers from a given cell have similar environments, so that
$\widehat{\boldsymbol{\mathbf{f}}}(\boldsymbol{\mathbf{r}})$ approximates the
perceptions of all speakers in zone $\boldsymbol{\mathbf{r}}$. The assumption
that all the speakers in a cell access the same primary linguistic data is a
simplification of reality. We know that within many geographical areas there
are subgroups whose language differences may depend on class, sex or ethnicity
[18]. We simplify in this way because we are interested in national level
spatial variations for which we have detailed geographical data, but no
consistent breakdown by social subgroup, so cannot infer any effects of
subgroups. We implicitly assume that social subgroups can induce biases on
variants.

###  3.2 Migration

Migratory processes may be internal, where individuals move between locations
in a single language area (typically a country), or external, where they
arrive from another country or language area. During the 20th Century,
foreign-born people have made up a small and slowly changing percentage of the
English population; 4.2% in 1951, for example, rising to 6.7% in 1991 [79].
The cumulative total number of migrants to Britain from the early 1800s to
1945 is estimated as $2.34$ million [80]. Since numbers are small, we do not
explicitly model external migration, instead viewing it as one of many
possible factors which may induce biases. To model internal migration, we need
to know how often people move, and how far. The English Longitudinal Survey of
Ageing (2007), showed that for individuals aged 50-89, the average number of
different residences occupied for at least six months during their lifetime
was 5.6 for those born 1918-1927, rising to 6.42 for those born 1948-1957
[81], giving an average of approximately one move every ten years. Inter-
county migration distances may be extracted from birthplace and place of
residence recorded in censuses, providing data back to the mid nineteenth
century [82]. Such data is, however, spatially coarse; it provides no
information about short range migrations. Modern lifestyle survey data can
help fill the gap. Analysis [83] of large scale research polls by Acxion Ltd
(2005-2007), which recorded current and previous postcode, generated $\approx
1.25\times 10^{5}$ migration distances, having mean 25.77km, median 2.89km and
standard deviation 63.91km. This implies a migration distribution concentrated
at short distances, with sub-exponential decline at large distances. This
large distance migration behaviour is often modelled with a power, or
"gravity" law [83, 84].

Models of migration must be consistent with migration distance statistics, and
constrained so that total flows in and out of cells match observed values.
This can be achieved [85] by writing the expected number of people who leave
their address in cell $\boldsymbol{\mathbf{r}}$ and move to a new address in
$\boldsymbol{\mathbf{r}}^{\prime}$ as

$m(\boldsymbol{\mathbf{r}},\boldsymbol{\mathbf{r}}^{\prime})\triangleq
A(\boldsymbol{\mathbf{r}})B(\boldsymbol{\mathbf{r}}^{\prime})O(\boldsymbol{\mathbf{r}})D(\boldsymbol{\mathbf{r}}^{\prime})h(|\boldsymbol{\mathbf{r}}-\boldsymbol{\mathbf{r}}^{\prime}|)$
(5)

where $O(\boldsymbol{\mathbf{r}})$ is the total number of movers whose
previous address was in $\boldsymbol{\mathbf{r}}$,
$D(\boldsymbol{\mathbf{r}}^{\prime})$ is the total number of people whose new
address is in $\boldsymbol{\mathbf{r}}^{\prime}$, and
$h(|\boldsymbol{\mathbf{r}}-\boldsymbol{\mathbf{r}}^{\prime}|)$ is a distance
function which, typically, decreases for longer moves. The conditions

$\displaystyle\sum_{\boldsymbol{\mathbf{r}}^{\prime}}m(\boldsymbol{\mathbf{r}},\boldsymbol{\mathbf{r}}^{\prime})$
$\displaystyle=O(\boldsymbol{\mathbf{r}})$ (6)
$\displaystyle\sum_{\boldsymbol{\mathbf{r}}}m(\boldsymbol{\mathbf{r}},\boldsymbol{\mathbf{r}}^{\prime})$
$\displaystyle=D(\boldsymbol{\mathbf{r}}^{\prime}),$ (7)

imply that

$\displaystyle A(\boldsymbol{\mathbf{r}})$
$\displaystyle=\left(\sum_{\boldsymbol{\mathbf{r}}^{\prime}}B(\boldsymbol{\mathbf{r}}^{\prime})D(\boldsymbol{\mathbf{r}}^{\prime})h(|\boldsymbol{\mathbf{r}}-\boldsymbol{\mathbf{r}}^{\prime}|)\right)^{-1}$
(8) $\displaystyle B(\boldsymbol{\mathbf{r}}^{\prime})$
$\displaystyle=\left(\sum_{\boldsymbol{\mathbf{r}}}A(\boldsymbol{\mathbf{r}})O(\boldsymbol{\mathbf{r}})h(|\boldsymbol{\mathbf{r}}-\boldsymbol{\mathbf{r}}^{\prime}|)\right)^{-1}.$
(9)

A mover might shift to a new address in the same cell, and the quantity
$m(\boldsymbol{\mathbf{r}},\boldsymbol{\mathbf{r}})$ gives the expected number
of people who do this within cell $\boldsymbol{\mathbf{r}}$. The sets of
constants $A\triangleq\\{A(\boldsymbol{\mathbf{r}}_{k})\\}_{k=1}^{L}$ and
$B\triangleq\\{B(\boldsymbol{\mathbf{r}}_{k})\\}_{k=1}^{L}$ are found by
generating a sequence of sets $A_{0},B_{0},A_{1},B_{1},\ldots$ where $A_{0}$
is the set for which $A(\boldsymbol{\mathbf{r}})=1$ for all
$\boldsymbol{\mathbf{r}}$, $B_{k}$ is computed from $A_{k}$ using (9), and
$A_{k+1}$ is computed from $B_{k}$ using (8). The sequence is then iterated to
convergence.

Figure 2: Cumulative distribution of migration distances generated using
distance function (10) with $\omega=0.75$, $a=3$km and $\rho=1.5$. The mean,
median and standard deviation are 31.3km , 3.6km and 63.4km respectively. For
comparison, the corresponding statistics computed from Acxion survey data [83]
are 26.3km, 3.6km, and 63.7km. Horizontal dashed lines show the fraction of
moves less then ten miles (71%) and fifty miles (87%).

We assume that our cell populations are in equilibrium and without loss of
generality assume there is one mover per cell, so
$O(\boldsymbol{\mathbf{r}})=D(\boldsymbol{\mathbf{r}})=1$, and
$m(\boldsymbol{\mathbf{r}},\boldsymbol{\mathbf{r}}^{\prime})$ is a stochastic
matrix with elements which represent transition probabilities for a single
mover. The matrix can be rescaled to capture realistic migration rates. We
define the distance function

$h(r)=\omega e^{-(r/a)^{2}}+\mathbf{1}_{\\{r>0\\}}(1-\omega)r^{-\rho}$ (10)

where the first term models short range moves with typical range $a$, and the
second term is the long distance fat-tail of the distribution, with gravity
exponent $\rho$. The parameter $\omega\in[0,1]$ interpolates between entirely
short range, and pure power law distributions. Having calculated
$m(\boldsymbol{\mathbf{r}},\boldsymbol{\mathbf{r}}^{\prime})$ for our system,
the cumulative distribution of the migration distance $X$ is calculated as

$\mathbb{P}(X<r)=\frac{1}{L}\left[\sum_{\underset{\boldsymbol{\mathbf{r}}\neq\boldsymbol{\mathbf{r}}^{\prime}}{\boldsymbol{\mathbf{r}},\boldsymbol{\mathbf{r}}^{\prime}}}H\left(r-|\boldsymbol{\mathbf{r}}-\boldsymbol{\mathbf{r}}^{\prime}|\right)m(\boldsymbol{\mathbf{r}},\boldsymbol{\mathbf{r}}^{\prime})+\sum_{\boldsymbol{\mathbf{r}}}H\left(r-R(\boldsymbol{\mathbf{r}})/\sqrt{2}\right)m(\boldsymbol{\mathbf{r}},\boldsymbol{\mathbf{r}})\right]$
(11)

where $H$ is the Heaviside step function and $R(\boldsymbol{\mathbf{r}})$ is
the radius of a circle with identical area to cell $\boldsymbol{\mathbf{r}}$.
The second term in the square brackets accounts for all intra-cell moves,
assumed to be of distance $R(\boldsymbol{\mathbf{r}})/\sqrt{2}$, following
[86]. We select distance function parameters to match the short range (median)
and long range (standard deviation) statistics of the inter- and intra-cell
distance distribution computed in [83] for MSOAs, using the Acxion data. This
was achieved by parameter grid search to minimize the sum of the relative
errors in the median and standard deviation between the model and the Acxion
data. The resulting distribution (Figure 2) is our estimate for migration
distances at the turn of the 21st century. Its shape reflects two kinds of
migratory events: frequent short moves, perhaps to change accommodation, and
rarer long distance moves, leaving behind friends, work and other social
contacts, where the distance travelled becomes relatively less significant.

Figure 3:  Yearly percentage changes in population for English and Welsh
counties, averaged over 1851-1951, calculated from the census based estimates
of Friedlander [82]. Counties with the highest and lowest average rates of
change are labelled.

We now consider whether this distribution may be used to approximate migration
during in earlier periods. Increased move frequency over the 20th century does
not imply that the distribution of distances has also changed. Although we
lack detailed information about historical short range moves, we can compare
our distribution to statistics of longer migrations extracted [87] from the
National Marriage and Fertility survey (1959-1960), which sampled more than
2000 couples aged 16-60. The survey showed that 16% of moves were of 10-49
miles and 15% were 50 miles or more. Our corresponding percentages are 16% and
13% respectively (Figure 3). The similarity of these estimates, and the sample
error associated with historical percentages (estimated at $\pm 1\%$),
suggests that the statistics of long distance internal migrations have not
changed dramatically. Because people are most likely to migrate in their mid
twenties [81], many of the moves in the 1959-1960 marriage survey may have
been made much earlier in the century.

We also consider whether our assumption that each cell is in equilibrium is
reasonable. Using estimates of net inter-census population changes calculated
by Friedlander [82] for counties in England and Wales, we have computed the
yearly percentage population change by county, averaged over 1851-1951 (Figure
3). The results show south east England experienced the greatest inflow of
migrants, and that the root mean squared net flow rate over all counties is
$\approx 0.45\%$. By comparing this figure with the overall yearly migration
rate of $\approx 10\%$ we see that, on average, in cells that are gaining, for
every 20 people who leave, 21 will arrive. Since net changes in population are
therefore typically at least an order of magnitude smaller than than absolute
population flows, we assume that non-equilibrium effects can be ignored.

Long range migration in principle allows variants to jump between locations
within the spatial domain. This jumping process is the subject of the
linguistic gravity model [88, 18], which attempts to quantify the spread of
variants between a finite set of population centres, viewed as points in
space. In contrast, the cells in our model cover the entire spatial domain,
and population centres appear as densely packed clusters. Other than cell
size, there is no intrinsic difference between city and rural cells. The
effects of allowing city speakers to behave differently to their rural
counterparts, in a spatial model similar to ours, was explored in [29], along
with connections to other classical models of spatial linguistic spread
including the wave model and hierarchical diffusion [89, 40].

###  3.3 Language learning and evolution

We assume an iterative model of language evolution, in which dying adults are
replaced by new speakers who select a variant using the probability mass
function

$\boldsymbol{\mathbf{p}}(\boldsymbol{\mathbf{r}})\triangleq\boldsymbol{\mathbf{g}}\left(\widehat{\boldsymbol{\mathbf{f}}}(\boldsymbol{\mathbf{r}})\right)$
(12)

where $\boldsymbol{\mathbf{g}}:\Delta^{q}\rightarrow\Delta^{q}$ is the
learning function. An alternative interpretation, which yields dynamics which
have identical deterministic component, but no stochasticity, is that the
components of $\boldsymbol{\mathbf{p}}(\boldsymbol{\mathbf{r}})$ give the
relative frequencies with which new speakers use different variants. We could,
if we wished, interpolate between these two interpretations by assuming a
mixture of the two behaviours. In either case equation (12) captures the
learning process which converts infants into adult speakers. It is an
approximation of the cumulative effect of countless small iterations made
within a changing learning environment. Underlying this model is the
assumption that linguistic behaviour is mainly acquired in childhood, and
changes little in later life. Learning to speak requires the assimilation of a
variety of structures, from sets of individual sounds and sound patterns,
through words and morphological rules, to syntax. Learning processes differ
between structures, and we do not aim to capture these processes in detail.
Rather we propose a simple but plausible learning function which maps current
speech patterns to learned forms. We derive evolution equations by considering
the case which generates maximum stochasticity: where new speakers select a
single variant. We will then consider the importance of this stochasticity.

A simple learning model, analogous to neutral evolution in genetics [90, 91,
92, 5, 14], is that variants are selected with a probability equal to their
current frequency. This has two problems as a model of language evolution.
First, as we will show below, linguistic communities would take an
extraordinarily long time to settle on one feature. Second, the geographical
boundaries between language features (isoglosses) revealed by linguistic
surveys, appear to require some form of social conformity [7, 29, 31] meaning
that speakers preferentially select the most common variants with a
probability which exceeds their relative frequency. Such a selection strategy
is optimal if advantage can be gained by matching the speech of others. We can
incorporate both conformity and asymmetric bias toward different variants by
defining the learning function

$\left[\boldsymbol{\mathbf{g}}(\boldsymbol{\mathbf{f}})\right]_{k}\triangleq\frac{h_{k}f_{k}^{\beta}}{\sum_{i=1}^{q}h_{i}f_{i}^{\beta}}$
(13)

where $[\boldsymbol{\mathbf{f}}]_{k}$ denotes the $k$th component of the
vector $\boldsymbol{\mathbf{f}}\in\Delta^{q}$. The variables $h_{k}>0$ are
biases, and $\beta\geq 1$ is the conformity number. If $\beta>1$ then
$\boldsymbol{\mathbf{g}}(\boldsymbol{\mathbf{f}})$ increases the frequencies
of already popular variants, and reduces the frequencies of less popular ones.
The biases $h_{1},\ldots,h_{q}$, which we write as a vector
$\boldsymbol{\mathbf{h}}$, introduce variant asymmetry, with the selection
probability of the $k$th variant an increasing function of $h_{k}$. However,
$\boldsymbol{\mathbf{g}}(\boldsymbol{\mathbf{f}})$ is invariant under a
constant re-scaling of all biases so it is their ratios which determine the
extent to which variants are favoured.

To derive equations for the evolution of frequencies, we denote the counts of
speakers using the variants of a given linguistic variable within cell
$\boldsymbol{\mathbf{r}}$ at discrete time step $t$ as

$\boldsymbol{\mathbf{X}}(\boldsymbol{\mathbf{r}},t)=(X_{1}(\boldsymbol{\mathbf{r}},t),\ldots,X_{q}(\boldsymbol{\mathbf{r}},t))^{T}$
(14)

where

$\sum_{i=1}^{q}X_{i}(\boldsymbol{\mathbf{r}},t)=N.$ (15)

We model renewal (birth and death), and migration. At each time step, $n_{R}$
adults die and are replaced with young speakers, and $n_{M}$ migrate. The
renewal and migration rates are

$\displaystyle\lambda_{R}$ $\displaystyle\triangleq\frac{n_{R}}{N}$ (16)
$\displaystyle\lambda_{M}$ $\displaystyle\triangleq\frac{n_{M}}{N}.$ (17)

We set $\lambda_{M}=0.1$, corresponding to one migration per decade, and
$\lambda_{R}=0.02$ corresponding to a typical period of 50 years between
reaching linguistic maturity, and leaving the language community. We justify
this latter estimate on the basis that linguistic maturity occurs in the late
teens [93], and that UK life expectancy in the mid 20th century was in the
late sixties [94]. The total number of speakers who die or migrate in each
cell is $n\triangleq n_{R}+n_{M}$. These speakers are selected uniformly at
random from the existing population, implying that both lifetimes and times
between migrations are geometrically distributed with parameters $\lambda_{R}$
and $\lambda_{M}$, respectively.

If each migrator were to select their new cell using the symmetric transition
probability matrix
$m(\boldsymbol{\mathbf{r}},\boldsymbol{\mathbf{r}}^{\prime})$, then the
expected number of migrators received by each cell would be $n_{M}$, and the
expectation of the mean linguistic state of the migrators incident on cell
$\boldsymbol{\mathbf{r}}$ would be

$\sum_{\boldsymbol{\mathbf{r}}^{\prime}}m(\boldsymbol{\mathbf{r}}^{\prime},\boldsymbol{\mathbf{r}})\boldsymbol{\mathbf{f}}(\boldsymbol{\mathbf{r}}^{\prime},t)\triangleq\bar{\boldsymbol{\mathbf{f}}}(\boldsymbol{\mathbf{r}},t).$
(18)

We model migration by replacing our $n_{M}$ migrators with $n_{M}$ speakers
with states selected according to the probability vector
$\bar{\boldsymbol{\mathbf{f}}}(\boldsymbol{\mathbf{r}},t)$. We write the
variant counts of these speakers
$\boldsymbol{\mathbf{X}}_{M}(\boldsymbol{\mathbf{r}},t)\sim\text{multinomial}(n_{M},\bar{\boldsymbol{\mathbf{f}}}(\boldsymbol{\mathbf{r}},t))$.
The variant counts of the speakers who replace the dead is a random variable,
$\boldsymbol{\mathbf{X}}_{R}(\boldsymbol{\mathbf{r}},t)\sim\text{multinomial}(n_{R},\boldsymbol{\mathbf{p}}(\boldsymbol{\mathbf{r}}))$.
The variant counts of the speakers who die or migrate is a random variable,
$\boldsymbol{\mathbf{Y}}(\boldsymbol{\mathbf{r}},t)$, drawn from the multi-
hypergeometric distribution [95] with parameters
$(n,\boldsymbol{\mathbf{X}}(\boldsymbol{\mathbf{r}},t))$. The change in
variant counts between $t$ and $t+1$ is then

$\Delta\boldsymbol{\mathbf{X}}(\boldsymbol{\mathbf{r}},t)=\boldsymbol{\mathbf{X}}_{R}(\boldsymbol{\mathbf{r}},t)+\boldsymbol{\mathbf{X}}_{M}(\boldsymbol{\mathbf{r}},t)-\boldsymbol{\mathbf{Y}}(\boldsymbol{\mathbf{r}},t)$
(19)

where
$\Delta\boldsymbol{\mathbf{X}}(\boldsymbol{\mathbf{r}},t)\triangleq\boldsymbol{\mathbf{X}}(\boldsymbol{\mathbf{r}},t+1)-\boldsymbol{\mathbf{X}}(\boldsymbol{\mathbf{r}},t)$.
Dividing by $N$, we have

$\displaystyle\Delta\boldsymbol{\mathbf{f}}(\boldsymbol{\mathbf{r}},t)$
$\displaystyle=\mathbb{E}[\Delta\boldsymbol{\mathbf{f}}(\boldsymbol{\mathbf{r}},t)]+(\Delta\boldsymbol{\mathbf{f}}(\boldsymbol{\mathbf{r}},t)-\mathbb{E}[\Delta\boldsymbol{\mathbf{f}}(\boldsymbol{\mathbf{r}},t)])$
(20)
$\displaystyle=\lambda_{R}\left(\boldsymbol{\mathbf{g}}\left(\widehat{\boldsymbol{\mathbf{f}}}(\boldsymbol{\mathbf{r}},t)\right)-\boldsymbol{\mathbf{f}}(\boldsymbol{\mathbf{r}},t)\right)+\lambda_{M}\left(\bar{\boldsymbol{\mathbf{f}}}(\boldsymbol{\mathbf{r}},t)-\boldsymbol{\mathbf{f}}(\boldsymbol{\mathbf{r}},t)\right)+\Delta\boldsymbol{\mathbf{Z}}(\boldsymbol{\mathbf{r}})$
(21)

where
$\Delta\boldsymbol{\mathbf{Z}}(\boldsymbol{\mathbf{r}})\triangleq\Delta\boldsymbol{\mathbf{f}}(\boldsymbol{\mathbf{r}},t)-\mathbb{E}[\Delta\boldsymbol{\mathbf{f}}(\boldsymbol{\mathbf{r}},t)]$
is a zero mean noise term. The variance of the $k$th component of this term,
using the properties of the multinomial and multi-hypergeometric
distributions, and suppressing $\boldsymbol{\mathbf{r}}$ dependence for
brevity, is

$\mathbb{V}[\Delta
Z_{k}]=\frac{1}{N}\left(\lambda_{R}p_{k}(1-p_{k})+\lambda_{M}\bar{f}_{k}(1-\bar{f}_{k})+\frac{\lambda(1-\lambda)}{1-N^{-1}}f_{k}(1-f_{k})\right)$
(22)

where $\lambda\triangleq\lambda_{R}+\lambda_{M}$. The magnitude of the
stochastic element of cell dynamics therefore scales as $N^{-1/2}$. For
reasons set out below, we assume
$\Delta\boldsymbol{\mathbf{Z}}(\boldsymbol{\mathbf{r}})\approx 0$, in the
current paper.

###  3.4 The importance of noise

To explore the significance of the noise term in (21), we consider a small,
well connected community where linguistic evolution is neutral [96, 90, 91],
meaning that new speakers adopt variants with probabilities equal to their
relative frequencies, so
$\boldsymbol{\mathbf{p}}(\boldsymbol{\mathbf{r}})=\boldsymbol{\mathbf{f}}(\boldsymbol{\mathbf{r}})$.
In the case of a binary variable, the relative frequency of variant 1 admits
the diffusion approximation [97]

$df=\sqrt{\frac{2\lambda_{R}}{N}}\sqrt{f(1-f)}dW_{t}$ (23)

where $W_{t}$ is a standard Brownian motion [98]. Equation (23) is known as
the Wright-Fisher diffusion [97] and describes the evolution of the relative
frequency of one of two variants (alleles) of a gene in a haploid organism.
Notice that the evolution is driven entirely by noise in this case. Let $T$ be
the time taken for the cell population to settle on one variant ($f=0$ or
$1$). This is the fixation time. Starting from the initial condition
$f=\tfrac{1}{2}$, we have

$\mathbb{E}[T]=\frac{N}{\lambda_{R}}\ln 2.$ (24)

To be concrete, let us consider the case where $N$ is the average population
of a single cell ($N=7,787$) in our model, equivalent to a 0.25km2 block of
city streets. Using $\lambda_{R}=0.02$ (as above) we have
$\mathbb{E}[T]=269,877$ years. From this we see that even in a relatively tiny
population, the time for one variant to become dominant by noise alone is
comparable to the time for which humans have existed as a species. We
therefore neglect the stochastic element of our dynamics on the basis that it
is inadequate to explain changes over a plausible time scale. Many of the
linguistic variants which interest us have disappeared over the course of a
single century and only the deterministic terms in (21) are capable of
producing such a change.

Having decided to neglect the stochastic term in our model, we now consider
whether other forms of stochasticity could explain the observed changes. We
have assumed our small community is ``well connected'', meaning that all
learners are exposed the same inputs. Suppose instead that children learn
variants from only a small subset of the community. This effectively
subdivides the population into smaller groups (e.g. families) within which
variants are reproduced by neutral selection, and between which individuals
are exchanged (e.g. by marriage). Such processes are studied in models of
neutral genetic evolution in subdivided populations, such Wright's island
model [99], and Kimura's stepping stone model [96, 91]. At low exchange rates
between groups, the time to fixation increases relative to the fully mixed
(single group) case. At sufficiently high mixing rate the population becomes
effectively panmictic, and fixation time statistics match the fully mixed
case. Subdividing the population therefore does not change our decision to
neglect the stochastic term in our model. One way for unbiased random copying
to be a more powerful driver of change is via the topology and dynamics of the
social contact network [92]. For example, if it is strongly clustered then
certain individuals or groups become particularly influential. Another
possibility is that variants spontaneously gain linguistic momentum [21, 12,
100] starting from small stochastic fluctuations. We do not rule out either
possibility, but emphasize that stochasticity as we have modelled it is
inconsequential in large populations. This does not mean that we think
language change is deterministic. Rather, unpredictable driving forces must
result from factors which are not part of our model. Whether their origin is
stochastic momentum-like effects, social network dynamics, internal linguistic
biases or societal changes outside the linguistic system, they all exert an
effective bias (real or apparent) on language change, but in a direction that
cannot be known in advance. We can view these factors as random variables
whose values are realised as history unfolds. Having observed our linguistic
system's historical state, we can then ask: how strong must these biases have
been, acting in combination with migration, mixing and social conformity, to
produce today's language distributions? Having established which variants
require bias in our simple model, we can then turn to understanding their
origin.

###  3.5 Interpreting model parameters

Our deterministic evolution equation

$\Delta\boldsymbol{\mathbf{f}}(\boldsymbol{\mathbf{r}},t)=\lambda_{R}\left(\boldsymbol{\mathbf{g}}\left(\widehat{\boldsymbol{\mathbf{f}}}(\boldsymbol{\mathbf{r}},t)\right)-\boldsymbol{\mathbf{f}}(\boldsymbol{\mathbf{r}},t)\right)+\lambda_{M}\left(\bar{\boldsymbol{\mathbf{f}}}(\boldsymbol{\mathbf{r}},t)-\boldsymbol{\mathbf{f}}(\boldsymbol{\mathbf{r}},t)\right)$
(25)

models four processes, two spatial, two not. The two right hand terms are the
replacement term, with rate constant $\lambda_{R}$, and the migration term
having rate $\lambda_{M}$. Mundane mobility is represented by community matrix
$W(\boldsymbol{\mathbf{r}},\boldsymbol{\mathbf{r}}^{\prime})$ in our model,
which is used to compute
$\widehat{\boldsymbol{\mathbf{f}}}(\boldsymbol{\mathbf{r}},t)$ and migration
by the matrix $m(\boldsymbol{\mathbf{r}},\boldsymbol{\mathbf{r}}^{\prime})$,
used to compute $\bar{\boldsymbol{\mathbf{f}}}(\boldsymbol{\mathbf{r}},t)$;
the effects of these are modulated by conformity $\beta$. Together, these
represent a purely spatial, accommodation- and diffusion-driven model which
does not differentiate between variants. Factors which might create bias (see
section 2.3) are modelled with the bias vector $\boldsymbol{\mathbf{h}}$. If
dialect change has mostly been driven by mobility, as many scholars have
suggested, we would expect
$\boldsymbol{\mathbf{h}}\approx\boldsymbol{\mathbf{1}}$. However, if
linguistic factors have been crucial in selecting variants, if the standard
ideology has driven levelling, or if local identity factors are important in
explaining differences between conservative and innovative regions, then
$\boldsymbol{\mathbf{h}}$ will play a more substantial role. One possibility
is that bias drives the selection of the levelling target, but that migration
and other spatial factors determine the speed and course of levelling. This
would be reflected in the model by high but spatially uniform
$\boldsymbol{\mathbf{h}}$. In such a case, bias might reflect a linguistic
constraint derived from a property common to the linguistic system at all
localities, or a near-universal social pressure (such as normative pressure
from universal schooling). Alternatively, observed regional differences in
levelling in the EDA data might not be explained by migration and spatial
factors alone and must be modelled with a strongly spatially varying bias.
This would imply that local differences in language ideologies play an
important role, or that interactions between dialect features create internal
pressures varying between dialects.

##  4 Interaction between conformity and migration

In the absence of bias, when $\beta>1$, then the replacement and migration
terms in (25) represent two competing processes. The replacement term moves
the community toward a single local variant, whereas the migration term pulls
local forms toward the average linguistic state in a much wider linguistic
area. Before comparing our model to data, we illustrate how it operates.

###  4.1 Survival of an isolated domain

Figure 4:  Phase diagram showing regions of $(\beta,r)$ space for which
$\mathfrak{D}$ can retain variant 1 when $\phi_{1}=0$, for migration-
replacement ratios, $\rho=\lambda_{M}/\lambda_{R}$, in the set
$\rho\in\\{1,5,20\\}$. If a region of the phase diagram is stable for some
$\rho^{\ast}$, then it is stable for all $\rho<\rho^{\ast}$.

We consider a binary variable, and let
$\boldsymbol{\mathbf{f}}(t)=(f_{1}(t),f_{2}(t))^{T}$ be its relative frequency
vector within a dialect area $\mathfrak{D}$. Let $r$ be the escape distance:
the expectation of the distance from a randomly chosen point in
$\mathfrak{D}$, to the boundary of $\mathfrak{D}$, in a randomly selected
direction. The fraction of migrators who escape $\mathfrak{D}$ is then
approximately $1-\mathbb{P}(X<r)\triangleq\epsilon(r)$ where $\mathbb{P}$ is
the migration measure defined in (11). By the symmetry of the migration
matrix, this is equal to the fraction of speakers who enter $\mathfrak{D}$.
Writing $\phi_{1}$ for the relative frequency of variant 1 outside
$\mathfrak{D}$, then the average state of the migrators incident on
$\mathfrak{D}$ (including internal movements) will be

$\bar{f}_{1}=\underbrace{(1-\epsilon(r))f_{1}}_{\text{moves within
$\mathfrak{D}$}}+\underbrace{\epsilon(r)\phi_{1}}_{\text{migration from
outside $\mathfrak{D}$}}.$ (26)

The language state of $\mathfrak{D}$ then evolves (assuming
$\boldsymbol{\mathbf{h}}=\boldsymbol{\mathbf{1}}$) as follows

$\Delta
f_{1}=\underbrace{\lambda_{R}\left[\frac{f_{1}^{\beta}}{f_{1}^{\beta}+(1-f_{1})^{\beta}}-f_{1}\right]}_{\text{local
conformity}}+\underbrace{\epsilon(r)\lambda_{M}\left[\phi_{1}-f_{1}\right]}_{\text{long
range levelling}}.$ (27)

The first term generates conformity to the majority variant within
$\mathfrak{D}$. The second generates levelling toward the external language
state. Suppose that most speakers in $\mathfrak{D}$ use variant 1, which
distinguishes them from the nation as a whole, where variant 2 is used (so
$\phi_{1}\approx 0$). The equilibria ($\Delta f_{1}=0$) of (27) determine the
conditions under which $\mathfrak{D}$ will retain its distinctive feature.
Informally, we look for the values of $f_{1}$ for which the rates of
conformity and levelling are in balance

$|\text{local conformity}|=|\text{national levelling}|.$ (28)

When $\phi_{1}=0$, this condition may be written

$\frac{f_{1}^{\beta}}{f_{1}^{\beta}+(1-f_{1})^{\beta}}=\left(1+\frac{\epsilon(r)\lambda_{M}}{\lambda_{R}}\right)f_{1}.$
(29)

An upper bound on the migration rate which allows variant 1 to persist within
$\mathfrak{D}$ is obtained by considering solutions in the infinite conformity
limit $\beta\rightarrow\infty$ where learners always select the majority local
variant. In this case the left hand side of (29) is a unit step at
$f_{1}=\tfrac{1}{2}$, and there is a stable equilibrium
$f_{1}^{\ast}>\tfrac{1}{2}$ provided

$\frac{\epsilon(r)\lambda_{M}}{\lambda_{R}}<1.$ (30)

This is an intuitively reasonable condition; it says that a local distinctive
feature can in principle survive provided the rate of inward migration does
not exceed the rate at which young indigenous speakers replace dying adults.
We estimate that during the 20th century, $\lambda_{M}\approx 0.1$. There is
some flexibility in the way we interpret $\lambda_{R}$. If we consider
speakers to be members of the speech community from birth, then
$\lambda_{R}^{-1}$ is simply life expectancy. However, if speakers are only
influential once older, and cease to be so when very old, then
$\lambda_{R}^{-1}$ is the average of this time interval. As noted in section
3.3, we have set $\lambda_{R}=0.02$. What determines domain stability is not
however the absolute values of the migration and replacement rates. It is
their ratio

$\rho\triangleq\frac{\lambda_{M}}{\lambda_{R}}.$ (31)

For given $\rho$, the stability of $\mathfrak{D}$ can be computed for any
value of $\beta$ and $r$ by checking that the dynamics (27) has a stable fixed
point $f_{1}^{\ast}>\phi_{1}$. The phase diagram in Figure 4 shows that for
lower migration-replacement ratios smaller domains can persist for lower
values of conformity. In systems where the migration-replacement ratio grows
slowly, such as in the English dialect area in previous centuries, we would
expect smaller dialect areas to disappear as the phase boundary shifted.

###  4.2 Competition between migration and conformity in the English dialect
area

We consider a hypothetical period of 500 years, beginning with a domain in a
completely randomised state. This state is constructed by setting the initial
frequency $f_{1}(\boldsymbol{\mathbf{r}},0)$ of each cell to be an independent
uniform random variable on $[0,1]$. We allow the migration rate to be a
linearly increasing function of time, starting with zero migration and
reaching late 20th century levels by end of the period

$\lambda_{M}(t)=0.1\times\frac{t}{500}.$ (32)

Within a few decades, recognisable dialect domains emerge (Figure 5).

Figure 5:  A 500 year simulation starting from randomized initial conditions
with linearly increasing migration rate $\lambda_{M}(t)=0.1\times t/500$.
Other parameters $\lambda_{R}=0.02,\beta=2.5,\sigma=10$km.

These are analogous to domains in ferromagnetic materials which form because
neighbouring atoms tend to align their magnetic directions. Domain boundaries
evolve as if they feel a form of surface tension [101] which penalises
curvature. In physical systems, domains gradually expand as their boundaries
straighten. The analogy between these magnetic domain boundaries, and
isoglosses, was explored in detail in [7, 29], including the effects of non-
uniform population distribution. It was shown in [7] that concentrated human
settlements (villages, towns, cities) balance the surface tension effect by
inducing isogloss curvature, implying that stable isogloss shapes are
partially predictable from population distributions. At low migration levels,
competition between isogloss smoothing, induced by surface tension, and
curvature, induced by population distribution (as well as other natural
barriers) leads to stable geographical distributions. In Figure 5 the map
obtained after 200 years would be stable in the absence of increasing
migration. This stability is revealed in Figure 6, which shows how the system
evolves starting from an initial condition which is an exact copy of that in
Figure 5, but under conditions of consistently low migration
($\rho=\tfrac{1}{2}$).

Figure 6:  A 500 year simulation starting from randomized initial conditions
with constant low migration rate $\lambda_{M}=0.01$. Other parameters
$\lambda_{R}=0.02,\beta=2.5,\sigma=10$km.

This shows that the pattern reached after 200 years remains stable for the
next 300 years, and includes a number of isolated dialect areas with small
escape radii. Returning to Figure 5 we see that as migration increases these
small dialect areas disappear, being replaced with the variant which is in the
majority further afield. Eventually the system reaches a state observed in
both the SED and EDA surveys, which consists of a major isogloss bisecting
England. The stable north east dialect centred on Tyneside is also seen in
survey data. The major Wash-Severn isogloss is particularly stable because
neither variant is in the majority in the country as a whole, and because it
bisects the country via one of the shortest routes. A detailed mathematical
explanation of these effects, which are related to continuum percolation, is
given in [7, 29]. The north east dialect is stable in this simulation because
of its status as a dense but geographically isolated population centre.

##  5 Inference for English dialects

Table 1:  Questions, variables and variants for the EDA (and SED). All phonetic variants are written using the International Phonetic Alphabet [102]. ME = Middle English  Question in the EDA |  Variable |  Variants
---|---|---
1\. The word ``tongue'' can end in two different ways. Which do you use? |  Coalescence of [Ng] |  1\. [Ng], 2. [N]
2\. Which pronunciation of ``new'' is the most similar to your own? |  yod-dropping |  1\. no [j], 2. [j]
3\. How do you pronounce the word ``last''? |  TRAP-BATH split |  1\. [A:], 2. [a:], [æ:] 3. [a], [æ]
4\. In the word ``butter'', I pronounce the letter ``u'' as… |  FOOT-STRUT split |  1\. [U], 2. [2]
5\. Do you pronounce the ``r'' in ``arm''? |  rhoticity |  1\. no /r/, 2. /r/
6\. How do you pronounce the word ``three''? |  realisation of word-initial /T/ |  1\. [T], 2. [d], [ d], 3. [f], 4. [t]
7\. Do you pronounce the ``r'' before the ``-ing'' in ``thawing''? |  intrusive /r/ |  1\. /r/ 2. no /r/
8\. What is the season that follows summer? |  autumn lexical item |  1\. autumn, 2. backend, 3. fall
9\. How do you pronounce the ``l'' in ``shelf''? |  l-vocalisation |  1\. [], 2. [u], [U], 3. [l]
10\. What do you call a piece of wood stuck under your skin? |  splinter lexical item |  1\. shiver, 2. sliver, 3. speel, 4. spelk, 5. spell, 6. spile, 7. spill, 8. splint, 9. splinter, 10. spool
11\. How do you pronounce the word ``room''? |  realisation of ME /o:/ in ``room'' |  1\. [Y:], [Y], 2. [u:], 3. [U]
12\. How do you pronounce the word ``chicken''? |  weak vowel merger |  1\. [@], [2], 2. [I]
13\. How pronounce the word ``night''? |  realisation of ME /ix/ |  1\. [AI], 2. [A:], 3. [O], [6], 4. [Ei], 5. [2i], [@i], 6. [aI], [æI], 7. [a:], [æ:], 8. [i:]
14\. He wasn’t careful with the knife and managed to cut… |  initial element in 3sg m. reflexive pronoun |  1\. himself, 2. hisself
15\. If it belongs to a woman its… |  3sg f. possessive pronoun |  1\. hern, 2. hers
16\. Do you pronounce the ``h'' in ``hands''? |  h-dropping |  1\. /h/, 2. no /h/
17\. If you’d like someone to pass you something, would you say: |  dative alternation |  1\. give it me, 2. give it to me, 3. give me it
18\. How do you pronounce the word ``off''? |  LOT-CLOTH split |  1\. [O:], 2. [A:], 3. [6]
19\. How do you pronounce the word ``house''? |  realisation of ME /u:/ |  1\. [@u], [2u], [œ7], 2. [E:], 3. [Eu], 4. [aI], 5. [a:], 6. [æu], [æ7], 7. [au], 8. [u:]
20\. How do you pronounces the word ``bacon''? |  realisation of ME /a:/ |  1\. [I@], [Ia], 2. [æI], 3. [e:], [e@], 4. [ei], [Ei], [E:], [E;@]
21\. How do you pronounce the word ``five''? |  realisation of ME /i:/ |  1\. [AI], 2. [A:], 3. [OI], 4. [EI], 5. [2I], [@I], 6. [aI], [æi], 7. [a:], [æ:]
22\. How do you pronounce the last sound of ``bit'' when saying ``a bit of''? |  coda /t/ |  1.[P] , 2. [d], [R] , 3. [t]
23\. Fill in the gap: Every day on her walk past the lake, she … the ducks. |  habitual present 3sg. |  1\. do feed, 2. feed, 3. feeds
24\. How do you pronounce the word ``happy''? |  happY tensing |  1\. [9], 2. [I], 3. [e], 4. [i]
25\. An animal that carries its house on its back is a… |  snail lexical item |  1\. dod-man, 2. hodmedod, hoddy-dod, hoddy-doddy, 3. snail

###  5.1 Methodology

We now explore how our model predicts the evolution of linguistic variables
which have been recorded in both the SED and the EDA. These are listed in
Table 1. For each variable we evolve our model $T=100$ years forward in time
using its SED state as our initial condition. To generate this initial
condition from raw survey data, we set variant frequencies in each model cell
(MSOA) to match the variant frequencies from the nearest SED survey location.
If multiple variants were recorded in this location, they are assumed to all
have equal frequency. This process yields a tessellation of England into
single and shared-variant regions. Relative frequencies in each cell are then
replaced with average frequencies over its 20 nearest neighbours (mean
separation 5.4 km). This amounts to a form of smooth interpolation between
relative frequencies observed at SED survey points, which have typical
separation $\approx 20$km (calculated as square root of mean land area per
location). The resulting distribution has smooth transitions with widths
typically 10-30km, consistent with English transition zones discussed in Ref.
[18]. While we would not expect our results to be sensitive to the number of
nearest neighbours used in the smoothing (interpolation) step, we have avoided
smoothing to a level where frequency values at SED locations would be changed
substantially in transition regions (e.g. if the smoothing range approached
the typical separation of survey locations).

We consider both bias free ($\boldsymbol{\mathbf{h}}=\boldsymbol{\mathbf{1}}$)
and biased evolution. We allow for the possibility that biases may be
different in different parts of the country by introducing a spatially varying
field

$\boldsymbol{\mathbf{h}}(\boldsymbol{\mathbf{r}})=(h_{1}(\boldsymbol{\mathbf{r}}),\ldots,h_{q}(\boldsymbol{\mathbf{r}}))^{T}$
(33)

where $h_{k}(\boldsymbol{\mathbf{r}})$ is the bias on variant $k$ in cell
$\boldsymbol{\mathbf{r}}$. We do not allow the bias field to depend on time
because we have only initial and final conditions, and cannot therefore
calibrate time dependence. We wish to fit
$\boldsymbol{\mathbf{h}}(\boldsymbol{\mathbf{r}})$ so that the final state of
our simulation, initialised with an SED variable, closely matches the
corresponding EDA data. This is achieved by iterating the following learning
step, which increases or decreases local biases on variants which are
respectively under- or over-represented in the model, when compared to the EDA

$\widehat{\boldsymbol{\mathbf{h}}}_{n+1}(\boldsymbol{\mathbf{r}})=\eta(\boldsymbol{\mathbf{f}}_{\text{EDA}}(\boldsymbol{\mathbf{r}})-\boldsymbol{\mathbf{f}}_{n}(\boldsymbol{\mathbf{r}},T))+\theta(\boldsymbol{\mathbf{1}}-\boldsymbol{\mathbf{h}}_{n}(\boldsymbol{\mathbf{r}})).$
(34)

Here $\boldsymbol{\mathbf{h}}_{n}(\boldsymbol{\mathbf{r}})$ is the $n$th
iteration of our bias estimate and
$\boldsymbol{\mathbf{f}}_{n}(\boldsymbol{\mathbf{r}},T)$ is the frequency
distribution obtained using this estimate. The learning rate, $\eta$, controls
how rapidly adjustments are made, and the reversion rate, $\theta$, is a
regularization parameter that prevents runaway bias increases, and maintains
$\boldsymbol{\mathbf{h}}=\boldsymbol{\mathbf{1}}$ as the ``no-bias'' reference
point so that calibrations for different variables can be compared. We use the
initial condition
$\boldsymbol{\mathbf{h}}_{0}(\boldsymbol{\mathbf{r}})=\boldsymbol{\mathbf{1}}$.

To avoid over-fitting we introduce a smoothing factor, $\sigma_{s}$, which
interpolates between independent variation in bias between cells
($\sigma_{s}=0$) and constant bias in all cells
($\sigma_{x}\rightarrow\infty$). The quantity
$\widehat{\boldsymbol{\mathbf{h}}}_{n+1}(\boldsymbol{\mathbf{r}})$ is the new
bias estimate before spatial smoothing (indicated by the hat symbol). Having
applied the learning step (34), we apply a smoothing step using the
interaction matrix (2) with interaction range $\sigma_{s}$,

$\boldsymbol{\mathbf{h}}_{n+1}(\boldsymbol{\mathbf{r}})=\sum_{\boldsymbol{\mathbf{r}}^{\prime}}W_{\sigma_{s}}(\boldsymbol{\mathbf{r}},\boldsymbol{\mathbf{r}}^{\prime})\widehat{\boldsymbol{\mathbf{h}}}_{n+1}(\boldsymbol{\mathbf{r}}^{\prime}).$
(35)

This two-step iterative estimation method is continued until the bias field
stabilizes. The value of $\sigma_{s}$ is determined by optimising the trade
off between model complexity and model error. We measure complexity, $K_{N}$,
as the average $\ell_{1}$ norm of the difference between the bias in each
cell, and the average,
$\overline{\boldsymbol{\mathbf{h}}}_{N}(\boldsymbol{\mathbf{r}})$, of its $N$
nearest neighbours

$K_{N}=\sum_{\boldsymbol{\mathbf{r}}}\lVert\boldsymbol{\mathbf{h}}(\boldsymbol{\mathbf{r}})-\overline{\boldsymbol{\mathbf{h}}}_{N}(\boldsymbol{\mathbf{r}})\rVert_{1}.$
(36)

This measure approximates the magnitude of the second spatial derivative of
the bias field, averaged over the whole system. The number of neighbour cells
($N$) is chosen to be large enough to smooth out inhomogeneity in the pattern
of cell centroids, but small enough so that shorter range fluctuations in the
bias field are not ignored. We take $N=20$, corresponding to an average
neighbour distance of 5.4 km, on the assumption that repeated application of
the smoothing step will ensure that fluctuations on shorter scales are small.
Complexity will be high for fields with large variations in bias over short
distances throughout the system, and low for fields which are approximately
constant over large geographical areas. Higher learning rates will generate
bigger local changes at each step, but these are counteracted by increased
smoothing range. Reducing the learning rate to $\eta^{\prime}=\eta/k$,
$k\in\mathbb{N}$, is approximately equivalent to applying $k$ smoothing steps
between each learning step, or increasing the smoothing parameter to
$\sigma_{s}^{\prime}=\sqrt{k}\sigma_{s}$, so we can explore the behaviour of
our calibration method by fixing $\eta$ and varying $\sigma_{s}$.

Figure 7:  Model complexity versus model error for variable 7, using learning
parameters $\eta=0.15,\theta=0.03$ and smoothing parameters in the range
$\sigma_{s}\in\\{1,2,\ldots,10,12,\ldots,20,25,30\\}$. Model parameters
$\lambda_{M}=0.1,\lambda_{R}=0.02,\beta=2.0,\sigma=10$km.

Figure 7 shows the relationship between model complexity and model error using
variable 7 as an example. For high smoothing parameter ($\sigma_{s}=30$km) the
error approaches that of the maximally parsimonious model, where bias is
constant in space. Reducing $\sigma_{s}$ improves the model error, until we
reach an optimal region ($\sigma_{s}\approx 6$km), the ``knee'' point, with
low model errors and minimal complexity. Further increasing complexity has
minimal effect on error. We take this knee as our heuristic for determining
the optimal level of complexity because this is the simplest explanation that
fits the data (Occam's Razor [103]). Exploration of the error vs. complexity
trade-off for other questions reveals similarly located knee points. For
consistency we use $\sigma_{s}=6$km for all questions (results with
$\sigma_{s}=10$km are shown in supplementary material [104]). For comparison,
we also calibrate with very high smoothing $\sigma_{s}=150$km in the case
$\beta=2$, which produces an approximately constant bias field. We select
learning parameters $\eta=0.15,\theta=0.03$ to allow convergence within a
feasible time frame, while maintaining stability in the iteration.

We measure the difference between modelled and observed variant distributions
at time $T$ using the total variation distance [105] between the variant
frequency distributions in each cell, averaged over all cells

$\text{error}\triangleq\frac{1}{2L}\sum_{\boldsymbol{\mathbf{r}}}\lVert\boldsymbol{\mathbf{f}}_{\text{model}}(\boldsymbol{\mathbf{r}})-\boldsymbol{\mathbf{f}}_{\text{EDA}}(\boldsymbol{\mathbf{r}})\rVert_{1}\in[0,1]$
(37)

where $\lVert\cdot\rVert_{1}$ is the $\ell_{1}$ norm. The factor $(2L)^{-1}$
ensures that the maximum error is 1. In this extreme case no speaker in any
cell uses a variant observed in that cell in the survey.

The model is implemented in vectorized Python using sparse arrays which
accelerate the calculation of the community frequency vector. The number of
bias field iterations is set between 25 and 35, sufficient for the error (37)
to stabilize. To calibrate a single question using a single set of model
parameters takes approximately 5-10 minutes on a modest laptop (1.7 GHz, 4
cores). For much larger systems than ours, there is potential for the
migration matrix to create a bottleneck, requiring further optimization.

###  5.2 Bias free results

We begin by considering the bias free case
($\boldsymbol{\mathbf{h}}=\boldsymbol{\mathbf{1}}$) where the only factors
driving evolution are migration, daily movement and conformity. Figure 8 (a)
shows the range of model errors in this case. In approximately $25\%$ of cases
the model generates a final distribution with error ($<0.2$) comparable with
errors obtained by fitting a bias field (Figure 8 (b)). These are variants for
which spatial process and conformity play a substantial role in explaining
their evolution. Many exhibit one of two phenomena seen in two dimensional
coarsening systems studied by physicists [101]: stripe states, and shrinking
droplets.

Figure 8: Stacked histograms of model errors over all variables for conformity
$\beta\in\\{1,2,3\\}$, $\lambda_{R}=0.02,\lambda_{M}=0.1$. (a) Errors in the
bias free case $\boldsymbol{\mathbf{h}}=\boldsymbol{\mathbf{1}}$. (b) Errors
after calibrating bias using (34) and (35). Very high smoothing parameter
$\sigma_{s}=150$km leads to approximately constant spatial bias.

####  5.2.1 Stripe states

In the bias free case, provided $\beta>1$, our model behaves like a
subcritical coarse grained two-dimensional ferromagnet [101, 106]. The
condition $\beta>1$ allows interfaces to form between single variant domains,
analogous to domains of uniform spin alignment in magnetic systems. These
interfaces feel a form of surface tension, and in finite systems, stable
configurations often emerge where interfaces connect opposite system
boundaries [106], creating stripes. Variant [a] of variable 3 in Table 1
illustrates this phenomenon (Figure 9 (a)).

Figure 9:  SED (left), no-bias model results at $T=100$ (middle) and EDA
(right) distributions for (a) Variant [a] of variable 3, the trap-bath split,
(b) Variant [U] of variable 4, the foot-strut split. Model parameters
$\lambda_{M}=0.1,\lambda_{R}=0.02,\beta=2.0$.

Here, we see that the geographical boundary between those who pronounce the
vowel in the bath lexical set [107] (last, staff, brass, …) in the same way as
vowel as the trap set (tap, back, badge, …), from those who do not, has been
stable over the 20th Century. Note that lexical sets (intoduced by John C.
Wells [107]) are groups of words defined by some common phonological feature.
The variant [a] (the vowel in northern ``bath'') dominates in the northern
region, where the vowels in trap and bath do not differ. This same stability
is predicted by the model, and derives from the surface tension effect which
causes interfaces to find the shortest routes across the system, which often
begin in boundary indentations [7] (the Wash in Figures 9 (a) and (b)). One
important difference between our model and simple magnetic systems is
migration, which induces effective long-range interactions equivalent to local
biases in favour of the national majority variant. Because the trap-bath split
divides the country into approximately equal parts, these biases are close to
zero, stabilizing the boundary.

Another example of this effect, shown in Figure 9 (b), is the distribution of
the variant [U] of variable 4. This is the northern English form of the
stressed vowel in ``butter'' which belongs to the Wells' lexical set strut
(cup, suck, …). Speakers from the north typically use the same vowel sound for
the lexical set foot (put, bush). Although the no-bias model captures the
broad overall pattern of foot-strut, we can see from the EDA distribution that
the northern form is in fact receding, which indicates a bias in favour of the
southern pronunciation, [2]. Further evidence of such a bias, is that the
model predicts levelling of both variants, when in fact the southern form
remains dominant in the EDA.

The effect of increasing the conformity number is to make single variant
regions increasingly pure, boundaries sharper, and more strongly affected by
variations in population distribution. In this paper, for every variable, we
have run the model with $\beta\in\\{1,2,3\\}$. Mean errors in the biased and
no-bias model are summarized in Table 2. For the no-bias model we compute
average errors excluding cases where these exceed a given threshold, $C$,
allowing us to consider only cases where the model provides a plausible
description of changes. In the biased model we find that conformity generates
a better fit than no conformity (a paired t-test with null hypothesis of equal
mean error between $\beta=1$ and $\beta=2$ yields p-value 0.0016). There are
small differences ($\lessapprox 0.01$) between the cases $\beta=2$ and
$\beta=3$. In the no-bias model, when errors are within one standard deviation
of the biased model error, the conformity model performs marginally better
than the conformity free case. If we raise the error threshold to two standard
deviations the results are equivocal. Table 2 therefore provides only limited
evidence that $\beta>1$ in the no-bias case. We suggest that $\beta>1$ is more
plausible because in the long term, linguistic boundaries could not
spontaneously emerge or be maintained if $\beta=1$. Moreover, without
conformity the unbiased dynamics simply mixes variants, dissolving boundaries,
but leaving total variant frequencies over the system as a whole unchanged,
meaning that levelling is impossible. Previous modelling work [31] also
indicates that the SED distributions are more likely to have been generated by
a model with conformity. Finally, we note that the substantial linguistic
changes observed during our period of interest may be driven by a wide range
of factors, of which conformity is just one, makings its effects harder to
discern. A complete set of maps for all variants of all variables for the SED,
EDA, no-bias, and bias models (bias fields included) for $\beta=2$, is given
in the supplemental material for this paper. For the remainder of the paper we
focus on the case $\beta=2$, and consider spatial varying bias
($\sigma_{s}=6$km), approximately constant bias ($\sigma_{s}=150$km) and no
bias.

Table 2:  Mean errors for $\beta\in\\{1,2,3\\}$. Here $\varepsilon_{\text{bias}}$ and $\varepsilon$ denote errors in the bias and no-bias models respectively. The shorthand $\langle\varepsilon|\varepsilon<C\rangle$ denotes mean error conditional on $\varepsilon<C$. Counts give the numbers of questions satisfying the condition. Note that $C=0.17$ is one standard deviation above the mean error in the $\beta=2$ bias model, and $C=0.22$ is two standard deviations above. For comparison, the mean error with $\beta=2$ and spatially constant bias is $0.15\pm 0.09$. $\beta$ | $\langle\varepsilon_{\text{bias}}\rangle\pm\text{std. dev.}$ | $\langle\varepsilon|\varepsilon<0.17\rangle$ | Count($\varepsilon<0.17$) | $\langle\varepsilon|\varepsilon<0.22\rangle$ | Count($\varepsilon<0.22$)
---|---|---|---|---|---
1 | $0.133\pm 0.070$ | 0.105 | 6 | 0.121 | 7
2 | $0.109\pm 0.058$ | 0.087 | 5 | 0.123 | 8
3 | $0.112\pm 0.061$ | 0.089 | 5 | 0.126 | 8

####  5.2.2 Shrinking droplets

Another phenomenon seen in ferromagnetic interfaces, and in our model, is the
shrinking droplet. In two dimensions a droplet is simply a closed curve,
evolving under surface tension, which acts to shorten its length. As a result,
droplets become more circular and eventually shrink to nothing. In our model,
droplets of variants which are in the minority are further reduced by
migration. As an example, consider the form of the 3sg feminine possessive
pronoun, hers or hern (variable 15, Figure 10).

Figure 10:  SED (left), no-bias model results (middle) and EDA (right)
distributions for spatial distribution of variant 1 (hern) of SED variable 15.
Model parameters $\lambda_{M}=0.1,\lambda_{R}=0.02,\beta=2.0,\sigma=10$km.
Model error = 0.13.

Hern is a minority variant in the SED (29% vs. 71% nationally), surrounded by
regions in which hers is dominant. Without conformity, the mixing effect of
migration would move the system as a whole to a spatially uniform state with
all cells having frequency vector $\boldsymbol{\mathbf{f}}=(0.29,0.71)$. With
conformity, surface tension acting at the boundary of the hern region, plus
the effects of migration, almost extinguish the hern variant within 100 years.
Examples of droplet shrinkage from a different class of linguistic feature are
found in the reflexes of Middle English /u:/ (variable 19), where minority
variants such as [E:] and [aI] are predicted to undergo near-complete
levelling (disappearance) without bias.

###  5.3 Bias driven evolution

Figure 11: Variants with the maximum average bias per cell (given by bar
height) in the bias model with parameters
$\lambda_{M}=0.1,\lambda_{R}=0.02,\beta=2.0$. Each bar is labelled with the
variant that experiences the maximum bias. Starred bars indicate variables for
which the no-bias model with $\beta=2$ generates an error $<0.2$.

The evolution of many variables in Table 1 may be understood as the result of
combined spatial and bias effects. The primary effect of bias is to lend
certain variants a special status which may have either linguistic (internal)
or social (external) origin. In some cases, a bias appears to accelerate
primarily spatial processes. Figure 11 shows the average bias per cell applied
to the most biased variant of each variable, using $\beta=2$ (similar results
are obtained with $\beta=3$). Variables for which the no-bias model produces
an error $<0.2$ are starred and, with two exceptions, these variables have the
smallest maximum bias per cell. Bias magnitude therefore provides an
indication of the extent to which the pure spatial model fails to explain
observed changes between the SED and EDA. For the more extreme of the two
exceptions (variable 6, Figure 12), the no-bias model correctly predicts the
main 20th century change: disappearance of the stopped pronunciation of
``three'' from the West-Country, but fails to predict that the highly
localised Essex pronunciation ``free'', although it remained a minority
variant, would spread across the country, requiring a substantial bias in its
favour, to prevent complete levelling.

Figure 12: SED (left), no-bias model results (middle) and EDA (right)
distributions for variant 2 ([di:]) and variant 3 ([fi:]) of variable 6:
pronunciation of ``three''. Model parameters
$\lambda_{M}=0.1,\lambda_{R}=0.02,\beta=2.0$.

####  5.3.1 Special status variants

Consider variable 10, the word for a small piece of wood stuck under the skin.
Minority variants like shiver and sliver (see Figure 13) are correctly
predicted by the no-bias model to undergo near-total levelling.

Figure 13: SED (left), no-bias model results (middle) and EDA (right)
distributions for variants shiver and sliver of variable 10, splinter. Model
parameters $\lambda_{M}=0.1,\lambda_{R}=0.02,\beta=2.0$.

Such patterns are most typical of variables with numerous variants, many of
which were historically restricted to small regions. In contrast, the present
day distribution of splinter cannot be explained without bias (Figure 14). In
the SED this was the dominant southern variant, and the calibrated bias field
indicates that it must have attained a special status as a standard variant in
order to spread north. However, this progress was halted in far north east
England, where spelk persists. Our spatially invariant bias field (Figure 14
(b)) indicates that this effect can be partially, but not completely explained
by the region's geographical isolation. According to the spatially varying
bias field, spelk not splinter is the variant with special status in the north
east (Figure 14 (a)), without which it would be dominated by splinter, except
in the major population centre of the region (Newcastle). A similar effect is
observed in variable 12, with the preservation of weak-vowel merger in the
north east, when the rest of England has moved to the standard variant (no
merger).

Figure 14: (a) SED (far left), bias model results (middle left) and EDA
(middle right) distributions for variant splinter of variable 10. Far right
map shows calibrated bias field using $\sigma_{s}=6$km. Model parameters
$\lambda_{M}=0.1,\lambda_{R}=0.02,\beta=2.0$. (b) As for (a), but using
$\sigma_{s}=150$km to yield constant bias field.

The no bias model predicts the initial expansion of splinter due to its
dominance in the densely populated south. However, its progress north is
arrested by repulsion effects [7] created by the major urban conurbations of
Manchester, Sheffield and Leeds (Figure 15). We find that these cities play an
important role in the pronunciation of ``thawing'' (variable 7, Figure 18).

An important question is: how should the bias fields be interpreted? In the
strictest sense the fields are merely model fitting parameters which measure
how much asymmetry needs to be injected into the model to explain the observed
changes. Where the variant symmetric model adequately explains the observed
changes without the need for additional bias, we expect to see a calibrated
bias field close to one. This does not mean that in reality there is no
linguistic or social bias, only that conformity and migration alone are
sufficient to explain the changes. For example, in Figures 14 and 15 because
the no bias model predicts the first part of the northerly expansion of
splinter, minimal bias is needed in the south. The major obstacles to further
expansion are the cities of Manchester and Sheffield, which consequently have
high bias field. This does not imply that surrounding areas necessarily have
lower bias, only that it is not needed to explain observations. The bias field
is therefore most informative when viewed in conjunction with the predictions
of the no-bias model (a complete catalogue is in the supplemental material
[104]).

Figure 15: SED (left), no-bias model results (middle) and EDA (right)
distributions for variant splinter of variable 10, with model parameters
$\lambda_{M}=0.1,\lambda_{R}=0.02,\beta=2.0$.  Figure 16: All variants of
variable 19: reflexes of ME /u:/ in the no-bias model. For each variant we
have SED (left), no-bias model results (middle) and EDA (right) distributions.
For reference, variants are 1. [@u], [2u], [œ7], 2. [E:], 3. [Eu], 4. [aI], 5.
[a:], 6. [æu], [æ7], 7. [au], 8. [u:]. Model parameters
$\lambda_{M}=0.1,\lambda_{R}=0.02,\beta=2.0$.

In several cases the no-bias model would fail to predict levelling at all, or
would even predict levelling in the opposite direction to that observed. Here
the direction of levelling is driven by bias: some variant-specific external
or internal factor such as normative bias or markedness must be responsible
for the rise in the successful variant. An example of the first type is
variable 19 (Figure 16): the reflexes of ME /u:/ (generally Wells' mouth).
Although the no-bias model correctly predicts that small variants ([@U], [E:],
[aI]) are largely levelled out, since the major variants at the time of the
SED ([u:], [EU], [æU], [a:] and to a lesser extent [aU]) were found in complex
domains with none in the majority, the no-bias model incorrectly predicts a
great deal of mixing without much change in relative proportions (i.e. without
levelling). A powerful bias term (Figure 11) is needed to predict the observed
near-complete levelling to [aU] (see Figure 17).

Figure 17: (a) SED (left), bias model results (middle left) and EDA (middle
right) distributions for variant 7, [aU] of variable 19, reflex of ME /u:/.
Model parameters $\lambda_{M}=0.1,\lambda_{R}=0.02,\beta=2.0$. Far right map
shows calibrated bias field. (b) As for (a), but using $\sigma_{s}=150$km to
yield constant bias field.

In a different linguistic domain, we consider the initial element of the
3sg.m. reflexive pronoun (variable 14). Users of the standard variant,
him(self), represented a small minority in the SED and the no-bias model
predicts that this variant would disappear; instead, the opposite happened,
with substantial levelling towards him(self) in the EDA data. In the model,
this requires strong bias in favour of him(self) (see Figure 11).

####  5.3.2 Isogloss positioning

In certain cases, although the direction of change is driven by positive bias
throughout the spatial domain, the model correctly predicts the new location
of an isogloss as a product of coastline and population density. In Figure 18,
the model captures the movement northward of the intrusive r isogloss so that
its current position divides the country in a line stretching from the Humber
estuary to the Ribble estuary (the result of coastline indentations) with
irregularly shaped differentiation within the northern area (the result of
population density patterns). Here the spatially varying and spatially
constant bias field produce very similar results (Figure 18 (a) and (b))
indicating that isogloss positioning is influenced in a partially predictable
way by geography and population distribution, as predicted in refs. [7] and
[29]. The same isogloss appears (spuriously) in the bias free evolution of
variable 10 (Figure 13).

Figure 18: (a) SED (left), bias model results (middle left) and EDA (middle
right) distributions for variant /r/ of variable 7, intrusive r, with model
parameters $\lambda_{M}=0.1,\lambda_{R}=0.02,\beta=2.0$. Far right map shows
calibrated bias field. (b) As for (a), but using $\sigma_{s}=150$km to yield
constant bias field.

####  5.3.3 Accelerating spatially driven changes

For some variables, a pure spatial model would predict the direction of change
correctly, but at a slower pace than observed: here the bias term is needed to
increase the rate of change. An example is the changes affecting the
circumstances in which the rhotic consonant /r/ is pronounced (variable 5):
the no-bias model correctly predicts that the rhotic-non-rhotic isogloss in
the south is retreating westwards and that rhoticity is becoming less
categorical within the rhotic regions, but greatly underestimates this trend;
bias in favour of non-rhoticity is required to explain how advanced the change
is (see Figure 19). Another example is the levelling of the 3sg.f. possessive
pronoun (variable 15). The no-bias model correctly predicts that the
traditional southern and Midlands variant hern is levelled in favour of
standard hers to the extent that no clear isogloss is left, but incorrectly
predicts remaining hotspots of traditional usage in high population density
areas (Birmingham, Portsmouth/Southampton, Bristol). Bias towards the standard
variant accelerates the change in the model, so that these hotspots too are
levelled out. A more striking example is variable 8 (the word for autumn)
where the no bias model predicts that the autumn-backend isogloss will move
$\approx 100$km north, whereas in reality backend has all but disappeared,
being used by only $\approx 20\%$ of the population in the far north.

Figure 19: SED (left), no-bias model results (middle) and EDA (right)
distributions for variant /r/ of variable 5, rhoticity, with model parameters
$\lambda_{M}=0.1,\lambda_{R}=0.02,\beta=2.0$

####  5.3.4 What does bias represent?

In some cases bias plausibly reflects the existence of a phonetically natural
direction of change (the pronunciation of ``l'' in shelf shifting from clear,
through dark, to vocalised [l], for example, but not the reverse). In
principle, interaction between different dialect features could even result in
such internal bias being spatially variable, although there are no absolutely
convincing examples of this in these results. In rather more cases, however,
it seems most likely that bias represents external (non-linguistic) factors
and in particular the standard ideology. We can see this wherever bias favours
a long-standing standard variant. Examples include the 3sg.m. reflexive
pronoun (variable 14), rhoticity (variable 5), and the 3sg.f. possessive
pronoun (variable 15). Another example is the decrease in ``h dropping''
(variable 16) between the SED and EDA (see Figure 20); because this is a
phonologically unnatural change which undoes a simplification of the
pronunciation system, the bias in the model in favour of [h] may reflect the
normative pressure of the standard ([h] has gained the status of standard
variant). These results are inconsistent with an account which centres on
mobility and convergence as necessary and sufficient explanatory factors for
levelling. Instead, other, variant-specific factors—primarily external factors
such as the standard ideology—are required.

Figure 20: (a) SED (left), bias model results (middle left) and EDA (middle
right) distributions for variant [h] of variable 16, h-dropping, with model
parameters $\lambda_{M}=0.1,\lambda_{R}=0.02,\beta=2.0$. Far right map shows
calibrated bias field. (b) As for (a), but using $\sigma_{s}=150$km to yield
constant bias field.

A final pattern worth discussing is that of innovation diffusion. In four
variables, we can plausibly label the changes we see between the SED and EDA
as the expansion of a recent innovation: the [f] variant of /T/ (variable 6);
the [P] variant of coda /t/ (variable 22); the [] and [U] variants of coda /l/
(variable 9); and intrusive r (variable 7). Any account of the dynamics of
language change—including non-spatial accounts—must offer an explanation for
how innovations take hold. An innovation, by definition, begins as an extreme
minority variant, and so we are faced with a challenge to explain why such
variants are not immediately levelled in favour of the overwhelming majority.
Possible explanations include internal linguistic reasons why particular
variants might be favoured (although such explanations must also answer the
question of why the change had not already taken place), and social
explanations such as age vectors and incrementation or divergence for the
purpose of group identity formation [21, 100, 108]. This is a larger question
and not one that this paper seeks to answer; in our models, this too is simply
subsumed under the bias term.

###  5.4 Time series

Figure 21: One hundred year time series (in 25 year steps) for variant /r/ of
variable 7, intrusive r in ``thawing'', with model parameters
$\lambda_{M}=0.1,\lambda_{R}=0.02,\beta=2.0$. Figure 22: Predicted future of
foot-strut. Maps show predicted future of variant 1 [U] of variable 4 (Table
1); the vowel in ``butter''. Bias field calibrated to SED/EDA using
$\sigma_{s}=150$km, and model initialized with EDA. Mean calibrated bias
parameters $\boldsymbol{\mathbf{h}}=(0.71,1.29)$, and model parameters
$\lambda_{M}=0.1,\lambda_{R}=0.02,\beta=2.0$. Start year = 2016 (date of EDA
survey).

So far we have presented only the initial and final states generated by our
model, but we can also provide a complete reconstruction of the history of
each variant over the 20th Century. For example, Figure 21 shows our
predictions for the evolution of intrusive r in ``thawing''. In principle such
time series can be compared to mid-century data to validate our reconstruction
of the historical changes which have taken place. Such reconstructions may
have utility in historical linguistics, which aims to understand the sequences
of linguistic changes which have led to the current state of the world's
languages. If we were to make an initial guess at a historical language state,
a model similar to ours, adapted to allow for sequential changes in variants,
could be calibrated to later states to provide a predicted history. This could
then be compared to predictions made by more traditional linguistic methods.
In principle it would be possible to impose additional constraints on the
model, set by further data, at intervening points in its evolution. Where
changes are still in progress, the model might also be used to predict the
future. As an example, we have calibrated a constant bias field for variable 4
(foot-strut), giving model error 0.106. We then used this field to evolve the
model into the future, using the EDA as initial condition. The results are
shown in Figure 22, up to the year 2066. Constant bias was used to reduce the
chance of over fitting.

##  6 Conclusion

We have developed an explicit spatial model of language evolution, accounting
for mundane daily mobility, realistic migration patterns, spatial distribution
of population, and social conformity. One motivation is that many scholars
examining diffusion and levelling in English English over the last century
have argued for the primacy of movement (both ``routine'' mobility and
migration) in driving change (cf. section 2.3). This—tacitly or
explicitly—implies that ideological and normative factors and internal
linguistic pressures play only secondary roles: that widespread levelling
would have taken place even without the strong standard ideology that
characterises English language attitudes, and that this ideological context
served only to hasten or to determine the precise direction of processes which
were the product of more fundamental social-technological changes. In
scientific terms, we can think of the extreme form of this view, where all
change is driven by purely spatial processes (diffusion, long range mixing,
interface motion), as a null model \- the variant symmetric form of our model.
We have shown that although the evolution of some variants can be explained in
this way, the majority of historical changes could not have occurred without
certain variants gaining special status. In some cases this bias merely
accelerates or further drives changes which would have partially occurred by
spatial effects alone, requiring a relatively low level of asymmetry. In other
cases, strong biases are required either to drive minority variants such as
him(self) to national dominance, or to preserve distinctive local forms such
as spelk. In most cases the origins of these biases appear more likely to be
external social factors such as normative standards, rather than internal
linguistic factors. We have also seen that, even when bias is needed to
explain observations, spatial effects often play an important role. Our
conclusion is therefore that spatial process are important, but insufficient
to explain the majority of the changes we have seen. A complete catalogue of
results for all variants, including bias fields may be found in the
supplementary material [104].

Beyond assessing the importance of spatial processes to language change, this
modelling exercise represents a step towards the construction of minimalist
spatial models which can be used to reconstruct more general historical
linguistic changes, and potentially predict future change. In providing a
complete hypothetical time series for a long period of time, we are offering
explicit predictions which can be compared to data collected in England at any
point through the century. Since the bias field is not an observable quantity,
it is only these intermediate states that are directly falsifiable. We hope
that such comparisons may be used to test, improve or refute our model in
future.

### Code and data availability

No new data were created or analysed in this study. Burridge is happy to be
contacted with reasonable requests for additional maps, or source code.

### Acknowledgements

The authors thank the referees for their diligent and insightful reviews,
which raised a number of important questions. They are also grateful for the
support of a Royal Society APEX award APX\R1\180117 (funded by the Leverhulme
Trust).

## References

  * [1] C. Castellano, S. Fortunato, and V. Loreto.  Statistical physics of social dynamics.  Rev. Mod. Phys., 81:591, 2009.
  * [2] W. Bialek, A. Cavagna, I. Giardina, T. Mora, E. Silvestri, M. Viale, and A. Walczak.  Statistical mechanics for natural flocks of birds.  Proceedings of the National Academy of Sciences, 109(13):4786–4791, 2012.
  * [3] Yakovenko V. Dragulescu, A.  Statistical mechanics of money.  Eur. Phys. J. B., 17:723–729, 2000.
  * [4] M Barthelemy.  The statistical physics of cities.  Nature Reviews Physics, 1:406–415, 2019.
  * [5] G. J. Baxter, R. A. Blythe, W. Croft, and A. J. McKane.  Utterance selection model of language change.  Phys. Rev. E, 73(4):046118, 2006.
  * [6] J. Burridge and S. Kenney.  Birdsong dialect patterns explained using magnetic domains.  Phys. Rev. E, 93(6):062402, 2016.
  * [7] J. Burridge.  Spatial evolution of human dialects.  Phys. Rev. X, 7:031008, 2017.
  * [8] E Derman.  My Life as a Quant P: Reflections on Physics and Finance.  John Wiley and Sons, Hoboken, New Jersey, 2007.
  * [9] T. Hastie, R. Tibshirani, and J. Friedman.  The Elements of Statistical Learning.  Springer, New York, 2009.
  * [10] R.A. Blythe.  Neutral evolution: a null model for language dynamics.  Advances in Complex Systems, 15(3):1150015, 2012.
  * [11] M. A. Nowak, N. L. Komarova, and P. Niyogi.  Evolution of universal grammar.  Nature, 291:114–117, 2001.
  * [12] W. G. Mitchener.  A mathematical model of prediction-driven instability: How social structure can drive language change.  Journal of Logic, Language and Information, 20:385–396, 2011.
  * [13] W. G. Mitchener.  A stochastic model of language change through social structure and prediction driven instability.  SIAM Journal of Applied Mathematics, 77:2272–2293–396, 2017.
  * [14] R. A. Blythe and A. J. McKane.  Stochastic models of evolution in genetics, ecology and linguistics.  Journal of Statistical Mechanics: Theory and Experiment, P07018, 2007.
  * [15] R. A. Blythe and W. Croft.  S-curves and the mechanisms of propagation in language change.  Language, 88(2):269–304, 2012.
  * [16] F. Ghanbarnejad, M. Gerlach, J.M. Miotto, and EG Altmann.  Extracting information from s-curves of language change.  J. R. Soc. Interface, 11:20141044, 2014.
  * [17] E. Lieberman, Jean-Baptiste M., J. Jackson, T. Tang, and M. Nowak.  Quantifying the evolutionary dynamics of language.  Nature, 449:713–716, 2007.
  * [18] J.K. Chambers and Peter Trudgill.  Dialectology.  Cambridge University Press, 1998.
  * [19] A. Karjus, R. A. Blythe, S. Kirby, and K. Smith.  Challenges in detecting evolutionary forces in language change using diachronic corpora.  Glossa, 5:45, 2020.
  * [20] M. Newberry, C. Ahern, R. Clark, and J. Plotkin.  Detecting evolutionary forces in language change.  Nature, 551:223–226, 2017.
  * [21] K. Stadler, R. A. Blythe, K. Smith, and S. Kirby.  Momentum in language change.  Language dynamics and change, 6:171–198, 2016.
  * [22] W. Labov.  Principles of Linguistic Change.  Blackwell, Malden, Massachusets, 2001.
  * [23] G. Baxter and W. Croft.  Modeling language change across the lifespan: Individual trajectories in community change.  Language Variation and Change, 28:129–173, 2016.
  * [24] H. Kauhanen, D. Gopal, T. Galla, and R. Bermudez-Otero.  Geospatial distributions reflect temperatures of linguistic features.  Science Advances, 7:eabe6540, 2021.
  * [25] J. Stanford and L. Kenney.  Revisiting transmission and diffusion: An agent-based model of vowel chain shifts across large communities.  Language Variation and Change, 25:119–153, 2013.
  * [26] J. Nerbonne.  Measuring the diffusion of linguistic change.  Phil. Trans. R. Soc. B, 365:3821–3828, 2010.
  * [27] T. Tesileanu and H. Meyer-Ortmanns.  Competition of languages and their hamming distance.  Int. J Modern Physics C, 17:259, 2006.
  * [28] A. Baronchelli, L. Dall’Asta, A. Barrat, and Loreto V.  Topology-induced coarsening in language games.  Phys. Rev. E, 73:015102, 2006.
  * [29] J. Burridge.  Unifying models of dialect spread and extinction using surface tension dynamics.  R.Soc.opensci., 5:171446, 2018.
  * [30] J. Burridge, B. Vaux, M. Gnacik, and Y. Grudeva.  Statistical physics of language maps in the usa.  Physical Review E, 99:032305, 2019.
  * [31] J. Burridge and T. Blaxter.  Using spatial patterns of english folk speech to infer the universality class of linguistic copying.  Physical Review Research, 2:043053, 2020.
  * [32] H. Orton and E. Dieth.  Survey of English dialects.  E. J. Arnold, Leeds, 1962.
  * [33] A. Ellis.  On Early English Pronunciation, with Especial Reference to Shakspere and Chaucer. Part V.  Trübner and Co., London, 1889.  p226.
  * [34] Wilbert Heeringa and John Nerbonne.  Dialect areas and dialect continua.  Language Variation and Change, 13:375–400, 2001.
  * [35] A. Leemann, M-J. Kolly, R. Purves, D. Britain, and E. Glaser.  Crowdsourcing language change with smartphone applications.  PLoS ONE, 11:1–25, 2016.
  * [36] Adrian Leemann, Marie-José Kolly, and David Britain.  The English Dialects App: The creation of a crowdsourced dialect corpus.  Ampersand, November 2017.
  * [37] B. Vaux and M. L. Jøhndal.  Cambridge survey of world englishes, 2017.
  * [38] J. Grieve, D. Speelman, and D. Geeraerts.  A statistical method for the identification and aggregation of regional linguistic variation.  Language Variation and Change, 23:193–221, 2011.
  * [39] Y. Huang, D. Guo, A. Kasakoff, and J. Grieve.  Understanding u.s. regional linguistic variation with twitter data analysis.  Computers, Environment and Urban Systems, 59:244–255, 2016.
  * [40] L. Bloomfield.  Language.  Holt, Rinehart and Winston, 1933.
  * [41] Harold Orton and Eugen Dieth.  Survey of English Dialects.  Published for the University of Leeds by E.J. Arnold, Leeds, 1962\.
  * [42] Andrea Mathussek.  On the problem of field worker isoglosses.  In Marie-Hélène Côté, Remco Knooihuizen, and John Nerbonne, editors, The Future of Dialects, pages 99–116. Language Science Press, Berlin, 2016.
  * [43] K.J. Payne.  The Survey of English Dialects notebooks.  Tradition Today: The Journal of the Centre for English Traditional Heritage, 6:45–68, 2017.
  * [44] Donald A. Ringe and Joseph F. Eska.  Historical Linguistics: Toward a Twenty-First Century Reintegration.  Cambridge Univ. Press, Cambridge, 2013.
  * [45] William Labov.  Principles of Linguistic Change, volume 1: Internal Factors of Language in Society.  Blackwell, Oxford, UK ; Cambridge [Mass.], 1994.
  * [46] Peter Trudgill, Elizabeth Gordon, Gillian Lewis, and Margaret Maclagan.  Determinism in new-dialect formation and the genesis of New Zealand English.  Journal of Linguistics, 36(2):299–318, July 2000.
  * [47] Richard Coates.  Steps towards characterizing Bristolian.  In Laura Wright, editor, Southern English Varieties Then and Now, pages 188–226. de Gruyter Mouton, Berlin, July 2018.
  * [48] David Britain.  One foot in the grave? Dialect death, dialect contact, and dialect birth in England.  International Journal of the Sociology of Language, 196/197:121–155, 2009.
  * [49] David Britain.  Phoenix from the ashes?: The death, contact, and birth of dialects in England.  Essex Research Reports in Linguistics, 41:42–73, 2002.
  * [50] Eivind Torgersen and Paul Kerswill.  Internal and external motivation in phonetic change: Dialect levelling outcomes for an English vowel shift.  Journal of Sociolinguistics, 8(1):23–53, February 2004.
  * [51] Sali A. Tagliamonte.  Changing Places: Tracking Innovation and Obsolescence across Generations.  In Chris Montgomery and Emma Moore, editors, Language and a Sense of Place: Studies in Language and Region, pages 15–37. University of Cambridge, Cambridge, United Kingdom ; New York, NY, 2017.
  * [52] Paul Kerswill.  Koineization and Accommodation.  In J. K. Chambers, Peter Trudgill, and Natalie Schilling-Estes, editors, The Handbook of Language Variation and Change, Blackwell Handbooks in Linguistics, pages 669–702. Blackwell Publishers, Malden, MA, 2002.
  * [53] Paul Kerswill.  Models of linguistic change and diffusion: New evidence from dialect levelling in British English.  Reading working papers in linguistics, 6(02):187–216, 2002.
  * [54] Paul Kerswill.  Dialect levelling and geographical diffusion in British English.  In David Britain and Jenny Cheshire, editors, Social Dialectology. In Honour of Peter Trudgill, pages 223–243. Benjamins, Amsterdam, 2003.
  * [55] David Britain.  Supralocal Regional Dialect Levelling.  In Carmen Llamas and Dominic James Landon Watt, editors, Language and Identities, pages 194–204. Edinburgh University Press, Edinburgh, 2009.
  * [56] Peter Auer.  Europe's sociolinguistic unity, or: A typology of European dialect/standard constellations.  Perspectives on variation: Sociolinguistic, historical, comparative, 7:7–42, 2005.
  * [57] Unn Røyneland.  Vertical convergence of linguistic varieties in a language space.  In Peter Auer and Alfred Lameli, editors, Language and Space: An International Handbook of Linguistic Variation, number Bd. 30.1-Bd.30.3 in Handbücher Zur Sprach- Und Kommunikationswissenschaft = Handbooks of Linguistics and Communication Science, pages 259–274. De Gruyter Mouton, Berlin ; New York, 2010.
  * [58] Johan Taeldeman.  Linguistic stability in a language space.  In Peter Auer and Alfred Lameli, editors, Language and Space: An International Handbook of Linguistic Variation, number Bd. 30.1-Bd.30.3 in Handbücher Zur Sprach- Und Kommunikationswissenschaft = Handbooks of Linguistics and Communication Science, pages 355–374. De Gruyter Mouton, Berlin ; New York, 2010.
  * [59] Peter Trudgill.  Dialects in Contact.  Number 10 in Language in Society. B. Blackwell, Oxford, UK ; New York, NY, USA, 1986.
  * [60] H. Giles and T. Ogay.  Communication accommodation theory.  In B. B. Whaley and W. Samter, editors, Explaining communication: contemporary theories and exemplars.
  * [61] M. Babel.  Dialect divergence and convergence in new zealand english.  Language in Society, 39:437–456, 2010.
  * [62] W. Doise, A. Sinclair, and R. Y. Bourhis.  Evaluation of accent convergence and divergence in cooperative and competitive intergroup situations.  British Journal of Social and Clinical Psychology, 15:247–252, 1976\.
  * [63] J. S. Pardo, R. Gibbons, A. Suppes, and R. M. Krauss.  Phonetic convergence in college roommates.  Journal of Phonetics, 40:190–197, 2012.
  * [64] Peter Trudgill.  Sociolinguistic Typology: Social Determinants of Linguistic Complexity.  Oxford University Press, Oxford, 2011.
  * [65] Lesley Milroy.  Introduction: Mobility, contact and language change - Working with contemporary speech communities.  Journal of Sociolinguistics, 6(1):3–15, February 2002.
  * [66] David Britain.  The role of mundane mobility and contact in dialect death and dialect birth.  In Daniel Schreier and Marianne Hundt, editors, English as a Contact Language, Studies in English Language, pages 165–181. Cambridge University Press, Cambridge ; New York, 2013.
  * [67] Frans Hinskens, Peter Auer, and Paul Kerswill.  The study of dialect convergence and divergence: Conceptual and methodological considerations.  In Peter Auer, Frans Hinskens, and Paul Kerswill, editors, Dialect Change: Convergence and Divergence in European Languages, pages 1–48. Cambridge University Press, Cambridge, UK ; New York, 2005.
  * [68] Salikoko S. Mufwene.  The Ecology of Language Evolution.  Cambridge Approaches to Language Contact. Cambridge University Press, Cambridge, UK ; New York, 2001.
  * [69] Jack Chambers.  Dialect acquisition.  Language, 68:673–705, 1992.
  * [70] Jane Stuart-Smith, Gwilym Pryce, Claire Timmins, and Barrie Gunter.  Television can also be a factor in language change: Evidence from an urban dialect.  Language, 89(3):501–536, 2013.
  * [71] Dominic J. L. Watt.  Phonetic parallels between the close-mid vowels of Tyneside English: Are they internally or externally motivated?  Language Variation and Change, 12(1):69–101, March 2000.
  * [72] Inge Lise Pedersen.  Processes of standardisation in Scandinavia.  In Peter Auer, Frans Hinskens, and Paul Kerswill, editors, Dialect Change: Convergence and Divergence in European Languages, pages 171–195. Cambridge University Press, Cambridge, UK ; New York, 2005.
  * [73] Paul Kerswill and Ann Williams.  "Salience" as an explanatory factor in language change: Evidence from dialect levelling in urban England.  In Mari C. Jones and Edith Esch, editors, Language Change: The Interplay of Internal, External and Extra-Linguistic Factors, pages 81–110. Mouton de Gruyter, Berlin, 2002.
  * [74] Paul Foulkes and Gerard Docherty.  Urban Voices: Accent Studies in the British Isles.  Edward Arnold, London, 1999.
  * [75] Nanna Haug Hilton.  Regional Dialect Levelling and Language Standards: Changes in the Hønefoss Dialect.  Ph.D., University of York, York, 2010.
  * [76] 2011 census: Population and household estimates for small areas in england and wales.  Technical report, Office for national statistics, 2011.
  * [77] National travel survey 2014: Travel to school.  Technical report, Department for Transport, 2014.
  * [78] Paul Bolton.  Education: Historical statistics.  Technical report, House of Common Library, 2012.
  * [79] F Ambrose et al.  Focus on People and Migration: 2005 edition.  Office for national statistics, 2005.
  * [80] P Panayi.  An Immigration History of Britain.  Routledge, 2010.
  * [81] A. Bernard.  Cohort measures of internal migration: Understanding long-term trends.  Demography, 5:2201–2221, 2017.
  * [82] D. Friedlander and R. J. Roshier.  A study of internal migration in england and wales: Part i.  Population Studies, 19:239–279, 1966.
  * [83] J. Stillwell and M. Thomas.  How far do internal migrants really move? demonstrating a new method for the estimation of intra-zonal distance.  Regional Studies, Regional Science, 3:28–47, 2016.
  * [84] J. Stillwell.  Interzonal migration: some historical tests of spatial interaction models.  Environment and Planning A, 10:1187–1200, 1978.
  * [85] A. G. Wilson.  A family of spatial interaction models, and associated developments.  Environment and Planning, 3:1–32, 1971.
  * [86] M Batty.  Urban Modelling: Algorithms, Calibrations, Predictions.  Cambridge University Press, 1976.
  * [87] D. Friedlander and R. J. Roshier.  A study of internal migration in england and wales: Part ii.  Population Studies, 20:45–59, 1966.
  * [88] P. Trudgill.  Linguistic change and diffusion: description and explanation in sociolinguistic dialect geography.  Language and Society, 3:215–246, 1974.
  * [89] J. Schmidt.  The relations of the Indoeuropean languages.  Weimar: Hermann Böhlau, 1872.
  * [90] M. Kimura and Maruyama.  Pattern of neutral polymorphism in a geographically structured population.  Genetical Research, 18(2):121–131, 1971.
  * [91] J. T. Cox and R. Durrett.  The stepping stone model: New formulas expose old myths.  The Annals of Applied Probability, 12(4):1348–1377, 2002.
  * [92] H. Kauhanen.  Neutral change.  Journal of Linguistics, 53:327–358, 2017.
  * [93] William Labov.  Principles of Linguistic Change, volume 2: Social Factors of Language in Society.  Blackwell, Oxford, UK ; Cambridge [Mass.], 2001.
  * [94] Decennial Life Tables.  Office for National Statistics, 2021.
  * [95] W. Feller.  An Introduction to Probability Theory and Its Applications, Volume 1.  Wiley, 1968.
  * [96] M. Kimura.  Stepping stone model of population.  Ann. Rept. Nat. Inst. Genetics, Japan, 3:62, 1953.
  * [97] A. Etheridge.  Some Mathematical Models from Population Genetics.  Springer, London, 2009.
  * [98] B. Oksendal.  Stochastic Differential Equations.  Springer, 2010.
  * [99] S. Wright.  Isolation by distance.  Genetics, 28:114–138, 1943.
  * [100] J. Michaud.  Dynamic preferences and self-actuation of changes in language dynamics.  Language Dynamics and Change, 9:61–103, 2019.
  * [101] A.J. Bray.  Theory of phase-ordering kinetics.  Adv. Phys., 43:357–459, 1994.
  * [102] International Phonetic Association, editor.  Handbook of the International Phonetic Association: A Guide to the Use of the International Phonetic Alphabet.  Cambridge University Press, Cambridge, U.K. ; New York, NY, 1999.
  * [103] D. J. C. MacKay.  Information Theory, Inference and Learning Algorithms.  Cambridge University Press, 2003.
  * [104] J. Burridge and T. Blaxter.  Supplemental material, 2021.
  * [105] I. Csiszar and P. Shields.  Information Theory and Statistics: A Tutorial.  Now, Boston, MA, 2004.
  * [106] K. Barros, P. L. Krapivsky, and S. Redner.  Freezing into stripe states in two-dimensional ferromagnets and crossing probabilities in critical percolation.  Physical Review E, 80:040101, 2009.
  * [107] John C. Wells.  Accents of English.  Cambridge University Press, Cambridge [Cambridgeshire] ; New York, 1982.
  * [108] H. Kauhanen.  Replicator–mutator dynamics of linguistic convergence and divergence.  Royal Society Open Science, 7:201682, 2020.