Matchmaker, Matchmaker, Make Me a Match: Migration of Populations via Marriages in the Past

The study of human mobility is both of fundamental importance and of great potential value. For example, it can be leveraged to facilitate efficient city planning and improve prevention strategies when faced with epidemics. The newfound wealth of rich sources of data—including banknote flows, mobile phone records, and transportation data—has led to an explosion of attempts to characterize modern human mobility. Unfortunately, the dearth of comparable historical data makes it much more difficult to study human mobility patterns from the past. In this paper, we present an analysis of long-term human migration, which is important for processes such as urbanization and the spread of ideas. We demonstrate that the data record from Korean family books (called “jokbo”) can be used to estimate migration patterns via marriages from the past 750 years. We apply two generative models of long-term human mobility to quantify the relevance of geographical information to human marriage records in the data, and we find that the wide variety in the geographical distributions of the clans poses interesting challenges for the direct application of these models. Using the different geographical distributions of clans, we quantify the “ergodicity” of clans in terms of how widely and uniformly they have spread across Korea, and we compare these results to those obtained using surname data from the Czech Republic. To examine population flow in more detail, we also construct and examine a population-flow network between regions. Based on the correlation between ergodicity and migration in Korea, we identify two different types of migration patterns: diffusive and convective. We expect the analysis of diffusive versus convective effects in population flows to be widely applicable to the study of mobility and migration patterns across different cultures.


I. INTRODUCTION
Since Quetelet's advocacy of "social physics" in the 1830s [1] and Ravenstein's seminal work later in the nineteenth century [2], quantitative studies of human mobility have suggested that human movements follow statistically predictable patterns [3][4][5][6][7][8][9][10]. Such systems-level studies are an important complement to individual-based approaches, as they can reveal population-level phenomena that are difficult to deduce by focusing on the characteristics of isolated members [11].
Research that takes a physics-based approach has focused predominantly on modern mobility-rather than historical mobility and migration-due to the disproportionate availability of large, rich data sets from modern life [12][13][14][15][16].By contrast, historical data tends to be sparse, incomplete, and noisy.These constraints limit the scope of conclusions that one can draw about how humans mingled, mixed, and migrated over long time scales [17,18].In this paper, we investigate historical human mobility and associated human migration by studying the matchmaking process for traditional marriages in Korea combined with modern census data in South Korea.We obtain our data from Korean "family books" called jokbo (ᄌ ᅩ ᆨᄇ ᅩ in Korean).Such a confluence of historical and modern data is rare, and it allows a novel test of generative models for human mobility.the clan "Lee from Hakseong" is distinct from the clan "Lee from Jeonju (ᄌ ᅥ ᆫᄌ ᅮ ᄋ ᅵ)" [the royal clan of the Joseon dynasty and the Great Korean Empire ].When two Koreans marry, the bride's clan is customarily recorded in the jokbo owned by the groom's family.These jokbo are kept in the groom's family and passed down through the generations; they serve primarily as a record of the names and birth years of all male descendants [19,20].In previous work, researchers used the marriage data contained in these books to estimate the population sizes and distributions of clans in Korea as far as 750 years in the past [21][22][23].Such distributions are useful for understanding quantitative aspects of human culture, and we proceed even further by conducting a systematic investigation of the geographical information embedded in jokbo.
We examine a set of ten jokbo to try to understand how geographical separation affected human interaction in the past in Korea.Specifically, we examine how inter-clan marriage rates can be predicted by physical distance and how clans themselves have spread across the country during the past several hundred years.To do this, we apply two generative models for describing human mobility patterns to jokbo records of past marriages between two clans.Note that the identification of clans with specific geographical origins is not unique to Korea: for example, the origins of British and Czech surnames were also the subject of recent investigations [24,25].
Our analysis consists of two parallel approaches: we use marriages recorded in jokbo to obtain snapshots of migrations (mainly of individual women) for a "marriage-flux analysis."We apply two generative models for population flow, discuss the results of applying these models, and explain the limitations that arise from the wide variety in the geographicaldistribution patterns of the clans.To consider the geographical spread of clans in more detail, we conduct an "ergodicity analysis."We use the modern geographical distribution of clans from census data to infer "ergodicity" of clans (mainly caused by past movement of male descent lines).To provide an additional perspective, we also use this data to construct a network model of population flows.To the best of our knowledge, the notion of diffusive versus convective population flow is new for data-driven studies of human mobility and migration, and we believe that this kind of approach can provide valuable insights for many problems in population mobility and migration.In the present paper, we focus on long-term migration, which has significant effects on many processes over a variety spatial and temporal scales.Such processes include urban population growth and the demographic structure of cities [26]; city infrastructure and planning [27]; unemployment [28]; and the spread of culture, religion, and other ideas [29].Most early studies of migration emphasized so-called "internal migration" (i.e., movement within a country) [2,3,30,31] like the phenomena that we study, though international migration is also a prominent field of study [32,33].International migration has been a more popular field of study than internal migration during the past few decades, but present-day urbanization processes in Asia, Africa, and Latin America have led to renewed interest in internal migration [26,31].We hope that our work provides useful ideas to help solve some of the fundamental questions in the migration literature: who migrates, why people migrate, and the consequences of migration (e.g., rural depopulation).
The remainder of our paper is organized as follows.In Sec.II, we introduce the jokbo and census data that we use in our investigation.In Sec.III, we present our primary methodology for data analysis: the gravity and radiation models for marriage-flux analysis, a special case of the gravity model that we call the population-product model, and a diffusion model for ergodicity analysis.We present our main results in Sec.IV, and we conclude in Sec.V. We include detailed information on the data sets, data cleaning, additional results and practical considerations for our analysis, an investigation of a network model for population flow, and various other results in Appendices A-I.

II. DATA SETS A. Jokbo Data Sets
For our marriage-flux analysis, we use the same ten jokbo data sets that were employed in [21][22][23].An individual book contains between 1 873 and 104 356 marriage entries, and there are a total of 221 598 entries across all books.(See Ta-ble I and Figs. 8, 9 in Appendix A for details.)Each entry contains the bride's clan and year of birth [34].The oldest book has entries that date back to the 13th century.
Previous studies of this data set [21][22][23] did not use any of the information that is encoded implicitly in the geographical origins of each clan.Such information, together with the modern geographical distribution of clans, comprises a key ingredient of our analysis.We convert location names to geographical coordinates using the Google Maps Application Programming Interface (GMAPI) [35].Because of the much sparser coverage of North Korean regions by Google Maps (see Fig. 12 in Appendix D), this geolocation data is a biased sample of the full data.However, data for the southern half of the Korean peninsula is rich [36], and it is sufficient to draw interesting and robust conclusions.For example, the effect of a change in the legality of intra-clan marriage in 1997 is clearly observable in the data.

B. Modern Name Distributions
In addition to the jokbo data sets that we employ for marriage-flux analysis, we also use data from two Korean census reports (1985 and 2000) to evaluate the current spatial distribution of clans in Korea [37,38].As illustrated in Fig. 1, some clans have dispersed rather broadly but others remain localized (usually near their place of origin).Drawing on ideas from statistical mechanics [39,40], we use the term "ergodic" as an analogy to describe clans that have spread broadly throughout Korea.We suppose that such clans have reached a dynamic equilibrium: an ergodic clan is "spread equally" throughout Korea in the sense that one expects it to have roughly the same geographical distribution as the population as a whole.Note that we do not expect an ergodic clan to reach a spatially uniform state for the same reasons that the full population is not spatially uniform (e.g., inhomogeneities in natural resources, advantages to congregating in cities, etc.).
Non-ergodic clans should have rather different distributions from those that we dub ergodic, because their distribution must differ significantly from that of the full population.One can construe the notion of ergodicity as a natural extension of other physical analogies that were used in previous quantitative studies (including the original ones) on human migration [1][2][3][4][5][6][7][8][9][10].As we discuss later, we can quantify the extent of clan ergodicity.

A. Generative Models for Marriage-Flux Analysis
We compute a "marriage flux"-the rate of marriage of women from clan i into clan j-for all clan pairs (i, j) in our data [41].Historically, professional matchmakers were employed to travel between families to arrange marriages [42], so we posit that physical distance plays a significant role in determining marriage flux.We examine this hypothesis using two generative models: a conventional gravity model with adjustable parameters that incorporates the distance between regions and the effects (or lack thereof) of each regions's population [4,5]; and a recently developed, parameter-free radiation model [43,44].
The gravity model has been used to explain phenomena such as commuting patterns and disease spread [45][46][47][48].In this model, the flux of population G i j from a site i to a site j is where α, β, and γ are adjustable exponents.For our purposes, G i j is proportional to the flux of women from clan i to clan j through marriage.The total population of clan i is given by m i , and the variable r i j is the distance between the centroids of clans i and j.We employ census data from 2000 to calculate centroids using the spatial population distribution for each clan [37].Importantly, note that choosing γ = 0 in the gravity model yields a special case in which flux is independent of distance.As we will see in Sec.IV A, this situation arises when large uncertainties in geographical locations (due to clan ergodicity) hinder the accuracy of estimations of distances.
Determining the centroid locations of clans from modern census data is more accurate than attempting to locate the origins' names themselves [49] for two reasons.First, for many clans, origin-place names have differed from geographical clan centres since the beginning of recorded Korean history, in particular, during the period spanning our jokbo data sets [20,50].Second, the origin-place names for many clans have become outdated and cannot be located accurately via the names of modern administrative regions.For instance, the clan origin 'Hakseong' of the first author is an old name for the city Ulsan in South Korea, but the name 'Hakseong' is currently only used to describe the small administrative region 'Hakseong-dong' in Ulsan.However, as we demonstrate in Fig. 1, using the centroid location of 'Lee from Hakseong' correctly gives the modern city Ulsan.This procedure works in part because Lee from Hakseong is a non-ergodic clan; for ergodic clans such as Kim from Gimhae, the spatial precision is much worse.This is an important observation that we will discuss in detail later.
We use a version of the radiation model that takes finite-size effects into account [44].The population flux R i j from clan i to clan j is where R i = j R i j is proportional to the total population that marries from clan i into any other clan, N is the total population, and s i j is the exclusive population within a circle of radius r i j centered on the centroid of clan i.Note that members of clans i and j are not included in computing s i j [43].As before, m i is the population of clan i, members of clan i marry into clan j, and clan j keeps the marriage records.In contrast to the gravity model, the radiation model does not include any external parameters.Importantly, this renders it unable to describe the geographically-independent situation that we need to consider in our study (and which we can obtain by setting γ = 0 in the gravity model).
For both the gravity and radiation models, we use census data from the year 2000 [37] as a proxy for past populations.This allows us compute the quantities r i j , m i , and s i j .Our approximation is supported by previously reported estimates of stability in Korean society: historically, most clans have grown in parallel with the total population, so we assume that the relative sizes of clans have remained roughly constant [23].In both Eqs. ( 1) and ( 2), only the relative sizes m i /N and s i j /N matter for calculating the flux (up to a constant of proportionality).

B. Human Diffusion and the Ergodicity Analysis
One way to quantify the notion of clan ergodicity is to examine what we call the "clan density anomaly," which describes the local deviation in density of members of a given clan.The clan density anomaly is φ i (r, t) = c i (r, t) − [m i (t)/N(t)]ρ(r, t) at position r = (x, y) and time t, where c i (r, t) is the (spatially and temporally varying) local clan concentration (i.e., the clan population density), m i (t) is the total clan population, ρ(r, t) is the local population density (i.e., the total population of all clans at point r and time t, divided by the differential area), and N(t) is the total population of all of the clans at time t.If a clan were to occupy a constant fraction of the population everywhere in the country, then φ i = 0 everywhere because its local concentration would be c i = (m i /N)ρ.(This situation corresponds to perfect ergodicity.)The range of typical values for the clan density anomaly depends a clan's aggregate concentration in the country.Examining the anomaly relative to clan concentration, the year-2000 numbers for φ i /(m i ρ/N) range from −1700 to 7400 for Kim from Gimhae and from −19000 to 87000 for Lee from Hakseong.Clearly, the distribution of the latter is much more heterogeneous (see Fig. 17 in Appendix I).
Combining the notion of clan density anomaly with traditional arguments-flow ideas based on Ohm's law and "molecular weights for population" are mentioned explicitly in [6,10]-about migration from population gradients [2][3][4][5][6][7][8][9][10] suggests a simple Fickian law [51] for human transport on long time scales: we propose that the flux of clan members is J i ∝ ∇φ i , so individuals move preferentially away from high concentrations of their clans.This implies that ∂c i /∂t = ∇ • J i ∝ ∇ 2 φ i (where we have assumed that the constant of proportionality is independent of space), which yields the diffusion equation We thereby identify the constant of proportionality as an average diffusion constant D i with dimensions [length 2 /time].This prediction of diffusion of clan members is consistent with past theories that posited human diffusion (e.g., cultural [52] and demic [53] diffusion).An important distinction is that we are proposing a process of diffusive mixing of clans rather than diffusive expansion of an idea or group.If this theory is correct, then one should expect clan density anomalies to simply diffuse over time.One should also be able to estimate diffusion constants by comparing the spatial variance at two points in time.
One can gain insights into the above diffusion process by calculating the radius of gyration (a second moment) of the clan density anomaly as a proxy for measuring ergodicity.Suppose that clan i's concentration c i (r, t) is known on a set of discrete regions {S k } with areas {A k }.We define the centroid coordinates for the k-th region as where |S k | is the total number of coordinate points r in S k for normalization, and we henceforth use The centroid of the clan's anomaly has coordinates where r(k) = [x(k), y(k)] gives the coordinates of the centroid of region k and the normalization constant is where φ i (k, t) is the anomaly of clan i in region k at time t.Note that we calculate the centroid of population for the ith clan (as opposed to the centroid of its anomaly) using analogous formulas to Eqs. ( 5) and (6) in which φ i is replaced by the concentration c i .The radius of gyration (i.e., the spatial second moment) r g i (t) of clan i at time t is then defined by where • is the Euclidean norm.We can use the set of radii of gyration {r g (t)} from Eq. ( 7) as a proxy for ergodicity, because (by construction) r g i (t) quantifies how widely the clan density anomaly of clan i has spread across Korea [54].
We simulate Eq. (3) between the known anomaly distributions from census data at t 1 = 1985 and t 2 = 2000 to estimate a best-fit diffusion constant D i for each clan.We compare our results to a null model in which movement is diffusive but driven by the aggregate population density in each region rather than by clan population anomaly.Our clan-based diffusive model performs better than the null model for approximately 84% of the clans.

A. Marriage-Flux Analysis Based on Jokbo and Modern Census Data
We apply a least-squares fit on a doubly logarithmic scale to determine the coefficients α and γ from Eq. (1) (along with the proportional coefficient a G , which is essentially a normalization constant, for the total number of marriages).The parameter β is irrelevant for the aggregated entries in a single jokbo because m j is constant (and is equal to the total number of grooms in that jokbo).The strongest correlation between the gravity-model flux and the number of entries for each clan appearing in jokbo 1 occurs for α ≈ 1.0749 and γ ≈ −0.0349, which suggests that the frequency of marriage between two families is proportional to the product of the populations of the two clans and, in particular, that there is little or no geographical dependence.The likely explanation is that the clan in jokbo 1 is ergodic, so the grooms could have been almost anywhere in the country, which would indeed make geographical factors irrelevant.(In the context of population genetics, this corresponds to "full mixing" [55][56][57][58].)In other words, as we discussed in Sec.III A, this special case of gravity model (for which we use γ = 0 in our analysis) corresponds to having geographical independence.Consequently, we will henceforth use the term population-product model for the gravity model with γ = 0.For our analysis of other jokbo and additional details, see Appendix A (and Tables I and II therein).
With little loss of accuracy for the fit, we take γ = 0 (i.e., we use the population-product model) to avoid divergence in the rare cases in which a bride comes from the same clan as the groom (for which the distance is r i j = 0).We also take α = 1 with little loss of accuracy.Using γ = 0 allows us to include data from the approximately 22% of clans for which geographical origin information is not available.In Fig. 2, we show the fit for jokbo 1, where we have used linear regression to quantify the correlation between the populationproduct-model flux and the number of entries for each clan in the jokbo.The noticeably lower outlier to the right of the line is the data point that corresponds to the clan of jokbo 1, and we remark that this deviation results from a cultural taboo against marrying into one's own clan.Women from the same clan as the owners of a jokbo have traditionally been strongly discouraged from marrying men listed in the jokbo (it is possible that they were even recorded under false names in the book), and it was illegal until 1997 [59].For the other jokbo, see Fig. 6 in Appendix A. In the bottom panel of Fig. 2, we illustrate that the radiation model does not give a good fit to the data.Recall from our discussion in Sec.III A that the lack of parameters in the radiation model does not allow us to explicitly consider a geographically-independent special case when using it.We emphasize, however, that this does not imply that the gravity model is "better" than the radiation model, as a direct comparison between the two models is hampered by the ergodicity of clans.In other words, the current formulations of the gravity and radiation models do not provide a solution for how to estimate fluxes between the clan centroids.Consequently, to investigate population fluxes, we incorporate modern census data.See our discussions in the next subsection and in Appendix H.

B. Ergodicity Analysis Based on Modern Census Data and a Simple Diffusion Model
We use census data from the year 2000 [37] to examine the ergodicity of clans in three different ways: (1) the number of administrative regions quantifies how "widely" each clan is distributed; (2) the radius of gyration, which we calculate from the clan density anomaly using Eq. ( 7), quantifies how "uniformly" each clan is distributed; and (3) the standard deviation of anomaly values, which measures how much the anomaly varies across regions.For instance, using data from the 2000 census and considering all of the clans and the 199 standardized regions, we find that 3.04% of the clans have a member in every region but that 22.1% of the clans have members in 10 or fewer regions.
We illustrate the dichotomy of ergodic versus non-ergodic clans with the bimodal distribution in Fig. 3(a).However, from the perspective of individual clan members [see Fig. 3(b)], such a dichotomy is not apparent.We show the radii of gyration that we calculate from the 2000 census data in Figs.3(c) and (d).We can again see the bimodality in Fig. 3(c).In Fig. 17 in Appendix I, we illustrate the dichotomy for Kim from Gimhae and Lee from Hakseong.I in Appendix A, all ten of the clans for which we have jokbo are fairly ergodic, so the variables associated with the j indices (i.e., the grooms) in Eqs. ( 1) and ( 2) have already lost much of their geographical precision, which is consistent both with γ = 0 (i.e., with using the populationproduct model) and with α = 0. Again see the scatter plots in Fig. 2, in which we color each clan according to the number of different administrative regions that it occupies.Note that the three different ergodicity diagnostics are only weakly correlated (see Fig. 18 in Appendix I).

As we indicate in Table
Our observations of clan bimodality for Korea contrast sharply with our observations for family names in the Czech republic, where most family names appear to be nonergodic [25] (see Fig. 19).One possible explanation of the ubiquity of ergodic Korean names is the historical fact that many families from the lower social classes adopted (or even purchased) names of noble clans from the upper classes near the end of the Joseon dynasty (19th-20th centuries) [20,60].At the time, Korean society was very unstable, and this process might have, in essence, introduced a preferential growth of ergodic names.
In Fig. 4, we show the distribution of the diffusion constants that we computed by fitting to Eq. (3).Some of the values are negative, which presumably arises from finite-size effects in ergodic clans as well as basic limitations in estimating diffusion constants using only a pair of nearby years.In Fig. 20 in Appendix I, we show the correlations between the diffusion constants and other measures.2. Flux predictions from the population-product model (i.e., the special case of the gravity model with γ = 0) with α = 1 and the radiation models for jokbo 1.(a) Scatter plot of the number of clan entries in jokbo 1 versus the corresponding centroid in 2000 using the population-product-model flux with α = 1.We compute the line using a linear regression to find the fitting parameter a G ≈ 6.55(4) × 10 −11 (with a 95% confidence interval) to satisfy the expression N i = a G G i j , where G i j is the population-product-model flux and N i is the total number of entries from clan i in the jokbo.(b) The same clan entries compared to the radiation model.We compute the line using a linear regression to find the fitting parameter a R ≈ 0.049(2) to satisfy the expression N i = a R R i j , where R i j is the radiation-model flux and N i is the total number of entries from clan i in the jokbo.In both panels, we color the points using the number of administrative regions that are occupied by the corresponding clans [see Figs.3(a) and (b)].The red markers (outliers) in both panels correspond to the clan of jokbo 1 (i.e., the case i = j).

C. Convection in Addition to Diffusion as Another Mechanism for Migration
The assumption that human populations simply diffuse is a gross oversimplification of reality.We will thus consider the intriguing (but still grossly oversimplified) possibility of simultaneous diffusive and convective (bulk) transport.In the past century, a dramatic movement from rural to urban areas has caused Seoul's population to increase by a factor of more than 50, tremendously outpacing Korea's population growth as a whole [61].This suggests the presence of a strong attractor or "sink" for the bulk flow of population into Seoul, as has been discussed in rural-urban labor migration studies [28].The density-equalizing population cartogram [62] in Fig. 21 in Appendix I clearly demonstrates the rapid growth of Seoul and its surroundings between 1970 and 2010.
If convection (i.e., bulk flow) directed towards Seoul has indeed occurred throughout Korea while clans were simultaneously diffusing from their points of origin, then one ought to be able to detect a signature of such a flow.In Fig. 5(a), we show what we believe is such a signature: we observe that the fraction of ergodic clans increases with the distance between Seoul and a clan's place of origin.This would be unexpected for a purely diffusive system or, indeed, in any other simple model that excludes convective transport.By allowing for bulk flow, we expect to observe that a clan's members preferentially occupy territory in the flow path that is located geographically between the clan's starting point and Seoul.For clans that start closer to Seoul, this path is short; for those that start farther away, the longer flow path ought to contribute to an increased number of administrative regions occupied and hence to a greater aggregate ergodicity.We plot the fraction of ergodic clans versus the distance a clan has moved (which we estimate by calculating distances between clan origin locations and the corresponding modern clan centroids) in Fig. 5(b).This further supports our claim that both convective and diffusive transport have occurred.To further examine clan ergodicity, we also compare each clan's radii of gyration r g to the distances of their origin location to (1) Seoul and (2) its present-day centroid (see Fig. 22 in Appendix I).The latter shows the same general tendency as in Fig. 5.We speculate that the absence of statistical significance in the correlation between r g and the distances from between clan origin locations and Seoul is a sampling issue, as we could not include many of the small clans in this calculation because we cannot estimate the locations of their centroids from our data (see Appendix B).We assume that clans that have moved a larger distance have also existed for a longer time and hence have undergone diffusion longer; we thus also expect such clans to be more ergodic.This is consistent with our observations in Fig. 5(b) for distances less than about 150 km, but it is difficult to use the same logic to explain our observations for distances greater than 150 km.However, if one assumes that long-distance moves are more likely to arise from convective effects than from diffusive ones, then our observations for both short and long distances become understandable: the fraction of moves from bulk-flow effects like resettlement or transplantation is larger for long-distance moves, and they become increasingly dominant as the distance approaches 325 km (roughly the size of the Korean peninsula).We speculate that the clans that moved farther than 150 km are likely to be ones that originated in the most remote areas of Korea, or even outside of Korea, and that they have only relatively recently been transplanted to major Korean population centers, from which they have had little time to spread.This observation is necessarily speculative because the age of a clan is not easy to determine:  the first entry in a jokbo (see Table I in Appendix A for our ten jokbo) could have resulted from the invention of characters or printing devices rather than from the true birth of a clan [20].Ultimately, our data are insufficient to definitively accept or reject the hypothesis of human diffusion.However, as our analysis demonstrates, our data are consistent with the theory of simultaneous human "diffusion" and "convection."Furthermore, our analysis suggests that if the hypothesis of pure diffusion is correct, then our estimated diffu-sion constants indicate a possible time scale for relaxation to a dynamic equilibrium and thus for mixing in human societies.In mainland South Korea, it would take approximately (100000 km 2 )/(1.5 km 2 /year) ≈ 67000 years for purely diffusive mixing to produce a well-mixed society.A convective process thus appears to be playing the important role of promoting human interaction by accelerating mixing in the population.Despite the limitations imposed by our data, we try to estimate and quantify the centrality of Seoul using a networkflow model for population, and we find suggestive differences between the flow patterns of ergodic and non-ergodic clans.For details, see Appendix H.

V. CONCLUSIONS
The long history of detailed record-keeping in Korean culture provides an unusual opportunity for quantitative research on historical human mobility and migration, and our investigation strongly suggests that both "diffusive" and "convective" patterns have played important roles in establishing the current distribution of clans in Korea.By studying the geographical locations of clan origins in jokbo (Korean family books), we have quantified the extent of "ergodicity" of Korean clans as reflected in time series of marriage snapshots.This underscores the utility of investigating the location distributions of individual clans.Additionally, by comparing our results from Korean clans to those from Czech families, we have also demonstrated that our approach can give insightful indications of different mobility and migration patterns in different cultures.Our ergodicity analysis using modern census data clearly illustrates that there are both ergodic and non-ergodic clans, and we have used these results to suggest two different mechanisms for human migration on long time scales.Many mobility processes involve a balance between diffusive spreading and attractiveness of a central location (and between more general diffusive and convective fluxes), so we believe that our approach in the present paper will be valuable for many situations.
A noteworthy feature of our analysis is that we used both data with high temporal resolution but low spatial resolution (jokbo data) and data with high spatial resolution but low temporal resolution (census data).This allowed us to consider both the patterns of human movement on short time scales (mobility via individual marriage processes) and its consequences for human locations on long time scales (human migration via clan ergodicity).An interesting further wrinkle would be to compare such mobility-derived time scales for human mixing patterns to genetically-derived patterns [55][56][57][58].
From a more general perspective, our research has allowed us to test the idea of using a physical analogy for modeling human migration-an idea put forth (but not quantified) as early as the 19th century [1][2][3][4][5][6][7][8][9][10].Physics-inspired ideas have been very successful for the study of human mobility, which occurs on shorter time scales than human migration, and we propose that Ravenstein was correct when he posited that such ideas are also useful for human migration.

ACKNOWLEDGMENTS
We thank Hawoong Jeong (ᄌ ᅥ ᆼᄒ ᅡᄋ ᅮ ᆼ) for providing data from Korean family books and Josef Novotný for providing data on surnames in the Czech Republic.We thank Tim Evans for introducing us to helpful references and Erik Bollt, Valentin Danchev, Sandra González-Bailón, James Irish, Philip Kreager, Michael Murphy, and Tommy Murphy for helpful comments and discussions.We thank Marc Barthelemy and Richard Morris for details about their work on constructing flow networks [83], and we thank the anonymous referees for their helpful comments and suggestions.M.A.P. and S.H.L. acknowledge a grant (EP/J001759/1) from the EPSRC.B.J.  In our investigation, we examine ten digitized jokbo that were first studied in Ref. [21].In Table I, we give basic information about the ten jokbo and here we summarise the results of some of our computations.
First, we applied the same gravity-model fit that we used for jokbo 1 to all of the jokbo data, and the results do not deviate much from those for jokbo.That is, γ ≈ 0 and α ≈ 1, so we can apply the population-product model with α = 1.The largest deviations in the two parameter values are α ≈ 1.4930 (for jokbo 7) and γ ≈ 0.5377 (for jokbo 6).Interestingly, we could not find any empirical value of γ < 0.6 reported in literature [4,5,[44][45][46][47][48], and it seems to be extremely rare to report any empirical values at all for gravity-model parameters.As one can see in Fig. 6, the choice of α = 1 and γ = 0 fits the data reasonably well.Note that the suppressed case of bride and groom being from the same clan is apparent in Fig. 6.This is indicated by the red markers, which are significantly below the other points in some of the panels and do not exist at all in other panels.We show the radiation-model results for the other jokbo in Fig. 7.
Additionally, we can see that all of the clans in the jokbo data that we study (i.e., the grooms' side of marriages) are "ergodic" in the sense that they are widespread across the nation in 2000.This is not surprising, as the availability of digitized jokbo data itself reflects clan popularity.We present the gravity-model fitting results for temporally divided jokbo 1 data in Table II, and we give results that use clan origin locations instead of population centroid in 2000 in Table III.(We also temporally divided the data from jokbo 6-because, as shown in Table I, it has the largest γ value among the 10 clans-and we found that it too does not exhibit systematic changes over time.)With these calculations, we again find that α ≈ 1 and γ ≈ 0 appear to be reasonable.The general trend of population change in Korea is also reflected in the jokbo data.In Figs. 8 and 9, we show the time series of entries for each jokbo normalized by the total entries.These plots suggest that jokbo of the different sizes at different times tend to follow TABLE I. Number of entries and other information available in each jokbo, values that we determined by using additional data that we obtained from other sources, and a summary of some of our computational results for the clan corresponding to each jokbo.For each jokbo, we indicate the ID (1-10), the year t 0 of its earliest entry, its number of entries N e , and the number of distinct clans (including at least one bride for each clan) N c among those entries [21].The quantity N γ=0 gives the number of clans from the 2000 census (which is 4 303) plus the number of clans in each jokbo that are not already in the census.We can use the latter set of clans in the gravity model when γ = 0 (i.e., for the population-product model, which is applicable without geographical information) and α = 1.(See the discussion in Appendix B.) We also indicate the best values for the fitting parameters α and γ of the gravity model in Eq. ( 1).We apply this fit to the brides' side of a marriage, and we calculated these values by minimizing the sum of squared differences using the scipy.optimizepackage in Python [63] (with initial values of α = γ = a G = 1.0 in our computations).For the clan that corresponds to each jokbo (i.e., the grooms' side), we compute the number of administrative regions N admin in which it appears based on census data from 1985 and 2000.We use the census data to compute a radius of gyration r g (km) for both 1985 and 2000 and to estimate a diffusion constant D (km 2 /year) for diffusion of clans between those two years.As in Fig. 5   in the 1985 (respectively, 2000) census was 3 520 (respectively, 4 303).There are 3 481 clans in common in the two censuses: 39 clans disappeared and 822 appeared.
In Fig. 10, we indicate how many administrative regions the 822 "new" clans occupy.New clans might correspond to foreigners who obtained South Korean citizenship during the fifteen-year period 1985-2000, or these clans might simply have been missing from the 1985 census due to error.Figure 10 supports the hypothesis that these are genuinely new clans, because their members have not spread to a large number of administrative regions.This gives a total of 6 687 distinct clans after we also incorporate the 2 384 additional clans that are listed only in the jokbo.In Table I, we indicate the number of distinct clans in each of the ten jokbo.There are 162 clans that appear in all ten jokbo.For all calculations with the gravity and radiation models, we use the 4 303 clans listed in the 2000 census data.When we use the population-product model (for which γ = 0), we do not require geometrical information, so we also use the additional clans listed in each jokbo.In this case, we denote the number of clans by N γ=0 (see Table I).

Appendix C: Standardizing Administrative Regions in 1985 and 2000
For the administrative regions, we use municipal divisions that are composed of city (ᄉ ᅵ in Korean), county (ᄀ ᅮ ᆫ in Korean), and district (ᄀ ᅮ in Korean) [65].In South Korea, there were 232 (respectively, 246) such administrative regions in 1985 (respectively, 2000).The difference in the number of regions between the two years reflects a slight restructuring of the political units.
For our computations, we need to unify the two different partitionings to be able to systematically compare results from 1985 and 2000 and to compute diffusion constants.To do this, we manually extract 199 "standardized" regions that we use for all computations involving administrative regions.Our construction necessitates many instances of operations like the following: For each operation, the region on the right is the standardized one that we use in our computations.In a given line, each different region is represented by a different letter.Thus, in the first line, two distinct regions from the 1985 census have merged into one region (and correspond exactly to that region) from the 2000 census, and we use this last region as one of our 199 standardized regions.In other examples, such as in the third line above, the standardized region does not correspond exactly to a single region from either census.Finally, we remark that the above operations are examples of what we needed to do to reconcile the 1985 and 2000 administrative regions.This is not an exhaustive list (e.g., four regions in 1985 corresponding to six regions in 2000), and we treat these other cases similarly.
For each standardized region, we sum the associated areas and populations of the constituent regions to obtain the area and population values that we use in our computations.The list of standardized regions is available online at https://drive.google.com/file/d/0B6cEJA5TN6vfNFU3OWhjcmJkS3c.For each standardized region, this data includes the component region names (in Korean) in 1985 and 2000, the latitudes and longitudes (and UTM easting and northing coordinates; see Appendix E) of the component region administrative centers, the geographical areas of the component regions, and the populations of component regions in 1985 and 2000.The data are in a tab-delimited text file, for which we have used the UTF-16 (16-bit Unicode Transformation Format) encoding scheme [66] for the Korean characters.
The regional boundaries drawn in Fig. 1 are from the 2010 data downloaded from Ref. [38].There is a slight difference between the regional boundaries in 2000 and 2010, so we map the coordinates of administrative regions in 2000 to those in 2010 by checking which "polygon" in 2010 encloses the coordinates of administrative regions from 2000.FIG. 8.For each jokbo, we plot the number of distinct clans N c versus the total number of entries N e on a doubly logarithmic scale.We calculate the red line via a linear regression to Heaps' law [64] using the expression N c = 10 b J N a J e .This yields a slope of a J ≈ 0.55 (7) and an intercept of b J ≈ 0.6(3) (with 95% confidence intervals).

Appendix D: Obtaining Geographical Information from Google Maps
To obtain the coordinates of the clans' origins and the administrative regions, we wrote a Python script that returns the latitude and longitude given a clan origin location's name.We used a Python module for geocoding via Google Maps Application Programming Interface (API) [67][68][69].We were able to successfully retrieve 3 900 clan origin locations out of the total of 4 303 clans present in the 2000 census data (see Fig. 11).We excluded the remaining 403 clan origin locations as erroneous because in each case there is a distance of more than TABLE III.Gravity-model parameters α and γ in Eq. ( 1) calculated using the clan origin locations (instead of the population centroids from the census data) for all of the jokbo by minimizing the sum of squared differences using with the scipy.optimizepackage in Python [63] We confirmed by exhaustive checking that the modern administrative regions of South Korea are accurate.(The first author, who is South Korean, manually checked all of the locations.)However, as shown in Fig. 12, the clan origin locations are severely undersampled in the northern part of the Korean peninsula due to Google Maps' lack of information about North Korea.We hope to include more North Korean regions in future studies, and this might be possible because Google is adding details of North Korea to their mapping service [71].
In Table IV, we present our Python code using Google Maps API.It requires pygeocoder (we used version 1.1.4),which is available as of July 16, 2013 [73].The code returns coordinates in latitude and longitude, which can then be converted to metric units (see Appendix E).

Appendix E: Universal Transverse Mercator (UTM) Coordinates
All of the distance measures that we employ use the Universal Transverse Mercator (UTM) geographical coordinate system, which assigns a local two-dimensional Cartesian coordinate system to a given zone on the surface of Earth [74].We use the UTM Python module [75], which converts (ϕ, λ) coordinates for latitude (ϕ) and longitude (λ) to UTM coordinates (and vice versa), where the UTM standard revision used by this module is WGS84 [76].Conversion from (ϕ, λ) coordinates to UTM coordinates can be also done using the website [77].
A point (i E , i N ) defined by UTM coordinates has two components: easting i E and northing i N .For example, the mean UTM coordinates of our standardized regions are i E ≈ 381.3 and i N ≈ 4 017.7,where the numbers are in units of kilometers from a reference point.The UTM scheme splits Earth into 60 zones.Calculating distances between two points in different zones is complicated in general, but the Korean peninsula lies entirely in zone 52 [78], which simplifies the calculation considerably.For example, Seoul's (latitude, longitude) coordinates are (ϕ, λ) ≈ (37.58, 127.00), and its UTM coordinates are (i E , i N ) ≈ (323.4,4 161.5).For our computations, we use formulas from Ref. [79].In a given grid zone, these formulas are accurate to within less than a meter.

Appendix F: Czech Republic Surname Data
The census data for the Czech republic was derived from the 2009 Central Population Register (produced by the Czech Ministry of the Interior) by the authors of Ref. [25].The vast majority of the Czech Republic is within zone 33, although a small part of it is in zone 34 [74, 78].As with the Korean clan origin locations, we used Google's API to roughly geolocate each of 206 Czech administrative regions by searching for the name of the administrative region followed by the string ", Czech Republic".We then converted all of the resulting latitude and longitude coordinates to UTM, forcing zone 33 for consistency.
We use the surname concentration (i.e., surname population density) to define an "anomaly" that indicates the difference in value from what would be observed for an "ergodic" for surname, whose members are well-mixed with the population.First, we obtain the centroid coordinates as in Eq. 4. The surname density anomaly, similar to what we defined for the Korean clans, is where c i (k, t) is the population density of surname i in region k at time t, the quantity m i (t) is the total population of surname i in all regions at time t, the quantity N(t) is the total population of all surnames in all regions at time t, and ρ(k, t) is the population density of all surnames in region k.We use the same notational conventions as for Korean clans, so It is necessary to introduce the anomaly (F1) because the total population is not conserved.Otherwise, we could simply take J ∝ ∇c i as the flux of people with surname i.Instead, we take the flux to be

Appendix G: Estimating Diffusion Constants
To estimate a diffusion constant for each clan, we start with the known anomaly distribution based on 1985 census data.Using an initial guess for the diffusion constant D i (where i indexes the clan), we integrate forward in time until 2000.We then compare the numerical prediction to the known anomaly distribution based on 2000 census data using a single number: the relative error, which we define as the sum of square deviations divided by the sum of square anomalies from 2000.
After the relative error is known, we can adaptively change the "guess" for D i and repeat the above process until we find the optimal D i .In practice, we use Matlab's built-in numerical optimization routine fminbnd [80], which implements a Nelder-Mead downhill simplex search.

Some Subtleties
Because the census data is irregular, we first interpolate it to a regular grid before numerical integration of the diffusion equation.The grid that we use covers the UTM zone 52 rectangular region from 245 to 545 km easting and from 3800 to 4250 km northing with a uniform 2.5 km spacing between grid points.(We exclude Jeju Island, which is distant from mainland Korea and is located south of the mainland.)We employ a standard five-point stencil to approximate the Laplacian operator in space and integrate in time with a 4th/5th order adap-tive Runge-Kutta scheme.We impose Neumann conditions at the boundaries.
Because the numerical integration is unstable for negative values of D i , we restrict D i to be positive.We test for the possibility of "negative diffusion" by repeating the entire optimization procedure after interchanging the 1985 and 2000 data sets; that is, we start from 2000 and integrate backwards in time with a positive diffusion constant.

Testing against a Null Hypothesis
Because our hypothesis of clan diffusion is somewhat speculative, it is important to test it against a basic null hypothesis.We take the null hypothesis to be a model in which clans do indeed diffuse, but in which clan affiliation plays no role: members simply diffuse in accordance with the local population density ρ.Therefore, We accept the null model as preferable to the clan diffusion model whenever it yields a lower relative error using its bestfit D i .

Numerical Results
The results of our computational examination of diffusion are as follows.When we include all 3 481 clans for which data was available, we obtain a mean diffusion constant D = D i (where we average over the clans) of D ≈ 2.91 and a standard deviation of σ D ≈ 10.4.When we remove "small" clans (i.e., those with fewer than 2 000 members in the year 2000), the 707 remaining clans have a mean diffusion constant of D ≈ 0.79 and a standard deviation of σ D ≈ 4.8.When we also remove "ergodic clans" (by eliminating the 75% with the largest year-2000 radii of gyration), the 182 remaining clans have a mean diffusion constant of D ≈ 1.3 and a standard deviation of σ D ≈ 2.0.
When we use all 3 481 clans, our computations favor the clan-diffusion model over the null model in about 84% of the cases.Additionally, about 64% of all clans had both positive best-fit diffusion coefficients and have mobility patterns that are explained better by the clan-diffusion model.
When we exclude both small and ergodic clans, our computations favor the clan-diffusion model over the null model in about 81% of the cases (i.e., for 148 clans).Moreover, 78% of all clans had both positive best-fit diffusion constants and have mobility patterns that are explained better by the clandiffusion model than by the null model.
In Fig. 4, we show a histogram of diffusion constants for the subset of clans for which the clan-diffusion model appears to be valid.These 148 clans are non-ergodic, have a positive best-fit diffusion constant, and are fit better by the clandiffusion model than by the null model.They have a mean diffusion constant of D ≈ 1.6 and a standard deviation of ) As one can see in the inset (for which we use a doubly logarithmic scale), the maximum movement distance is larger than 400 km.However, as the main panel illustrates, most clans moved considerably shorter distances.

Appendix H: Construction of Population-Flow Network between Regions
Although it is impossible to track the movement of individual people from the census data (because it does not include such information), it is possible to construct a rough estimate of the population flow between a pair of regions by considering the movement of clans (i.e., of the smallest demographic unit that it is possible to resolve with our data) between 1985 and 2000.For each clan i, we define its population centroid [note the contrast with the clan anomaly centroid from Eqs. ( 5) and ( 6)] as where the normalization constant is Recall that k indexes the region in Korea, r(k) = [x(k), y(k)] gives the coordinate of that region's centroid, A k is its area, and c i (k, t) is the population density of clan i in region k at time t.Additionally, recall that N i (k, t) = c i (k, t)A k is the population of clan i in region k at time t (see Sec. III B) We considered approximating the movement of each clan i by but this would entail treating an entire clan population as a "point mass," so it neglects valuable information from the spatial variation (as illustrated by our calculations of clan ergodicity).In addition, as we show in Fig. 13, the amount of movement for the majority of clans is too small to proceed further if we took such an approach.(One can also infer the small scale of clan movements from Fig. 4.) As an alternative that avoids the undesirable point-mass approximation, we attempt to infer the flow of a clan between two regions from changes in population ratios.Let each of the 199 standardized administrative regions of South Korea (see Appendix C for details) be individual nodes of a populationflow network [81] between 1985 and 2000.To examine the central nature of Seoul in such a network, we merge the 17 regions that corresponding to different parts of Seoul into a single node that we call "Seoul."The resulting population-flow network has 183 nodes.The raw change in relative populations between regions k and k in our census data is where we exclude the regions with N i (k, t = 1985) = 0 or N i (k, t = 2000) = 0 to avoid singularities.We then define the normalized relative population change as The quantity U i,k→k is a proxy for a more finely-grained quantity, which we denote by W i,k→k , that describes the real population flow of clan i from region k to region k (where k → k is a directed edge) and would be desirable to construct from data for flows of individuals.The smallest units in the census data are clans, so we need to estimate population flows from them as our basic units.We thus instead calculate Wi,k→k and use the normalized relative population changes U i,k→k in Eq. (H4) as the adjacency-matrix elements of a weighted and directed population-flow network.
Our proxy is not guaranteed to be "correct," but it has several properties that are consistent with reasonable flow behavior: (1) If the ratio of clan populations N 1 /N 2 (in regions k = 1, k = 2) does not change from 1985 to 2000, then both the proxy flow U i,2→1 and the inferred flow W i,2→1 from region 2 to 1 are equal to 0; (2) if the ratio N 1 /N 2 increases from 1985 to 2000, then the proxy and inferred flows from region 2 to region 1 are both positive; and (3) if the ratio N 1 /N 2 decreases from 1985 to 2000, then the proxy and inferred flows from region 2 to region 1 are both negative.Naturally, both the proxy flow and the inferred flow are asymmetric (so W i,2→1 −W i,1→2 ).
As a downside, uniform population growth biases the proxy calculation slightly in favor of flow towards regions with smaller populations.Additionally, the proxy cannot capture circulating flows and is unlikely to do a good job when flow is strongly multipolar (i.e., if more than one area attracts a significant amount of flow).When flow is mostly between Seoul and other regions, we call it "unipolar." Using Eq. (H4), we define a population-flow network for each clan i.We include all clans by using the entire population density of region k in year t.That is, we calculate N tot (k, t) = In Fig. 14(a), we show box plots of the distribution of U tot,k→k using all pairs k, k .We also show the distributions of U i,k→k for the clans Kim from Gimhae and Lee from Hakseong.We show flows to Seoul separately from flows to other regions.Note that the values of U tot,k→k are distributed much more broadly for flows to Seoul than for flows to other regions, even though there are many more adjacency-matrix elements for the latter (182 for flow to Seoul and 33 124 for flow to other regions).One can also observe this feature in the shapes of distribution themselves [see Figs.14(d)-(f)].
One simple but intuitive way to check the centrality of Seoul is to extract a maximum relatedness subnetwork (MRS) [82] from each population-flow network.We construct an MRS as follows.For each node, we examine the weight of each of its edges and keep only the single directed edge with maximum weight.(When there are ties, we keep all edges that have the maximum weight.)We exclude out-edges from Seoul for the MRS in order to focus on the movement from other regions to Seoul.We will later compare the MRS to a nullmodel network which also disallow out-edges from the central node.In Figs.15(a,b), we show the MRS from the adjacency matrix with elements U tot,k→k .The central role of Seoul is apparent.As we indicate in Table V, Seoul's in-degree in this MRS is 109.This constitutes nearly 60% of the MRS edges and is consistent with the rapid growth of the Seoul area that we illustrate in Fig. 21.
We model the population flow using a simple rewirednetwork model inspired by the model in Ref. [83].We start with a "monocentric" network, with Seoul as the central node, in which all directed edges start in some region (aside from the center) and terminate at Seoul.We then rewire each edge with independent, uniform probability p.Each network that we construct in this way has one directed edge for each node aside from the central one, so we can use an ensemble of such networks as a null model for our MRSs.
As we indicate in Table V, the edges in the MRSs are distributed rather heterogeneously among the regions.For example, the region in Gyeonggi Province (which has the second largest in-degree) has about 39% of the edges for the MRS that we constructed using all clans.When constructing nullmodel networks, we use a rewiring probability of p = 0.4    [82] of the combined population-flow network for all clans, (b) a magnified portion of this MRS that includes all regions with nonzero in-degrees, and (c) an instance of a model (inspired by the one in Ref. [83]) of a rewired version of a monocentric network (with Seoul as the center) with a rewiring probability of p = 0.4.(See the discussion in the text.)We color the directed edges toward Seoul in blue, and we color the directed edges towards other regions in blue.
to ensure that about 60% of the directed edges terminate in Seoul on average (as suggested by the data when considering all clans).The null-model network ensemble generated from the rewiring process has a binomial (or approximately Poisson, as the MRS is rather sparse) in-degree distribution as a result of the given fraction p of edges that are redirected uni-formly at random except for the central node or Seoul [81].Therefore, the emergence of a second-largest hub comparable in size to the largest hub (Seoul) is extremely unlikely.We illustrate one instance of such a rewired network in Fig. 15(c), and the MRS for all clans that we constructed from empirical data differs significantly from the null-model network.(See 24.4% Fig. 16(b

Table V as well).
It is also instructive to examine the population-flow networks for individual clans.As with prior discussions, we will use Kim from Gimhae as an example of an ergodic clam and non-ergodic Lee from Hakseong as an example of a non-ergodic clan (see Fig. 1).
When we consider the population-flow network for the clan Kim from Gimhae [by using N i (k, t) with i corresponding to Kim from Gimhae in Eq. (H3)], we obtain a qualitatively similar result -namely, an abundance of edges terminating in Seoul -as what we obtained when using all clans.See Fig. 14(b), Fig. 16(a), and Table V.By contrast, we find that two different locations "attract" the population for Lee from Hakseong.Following the general trend in the population, one area is the Gyeonggi Province in the northwestern part of South Korea surrounding the Seoul area.(The name means "the area surrounding capital" in Korean and it is often construed to be essentially an "extended Seoul.")The other area is Ulsan/Busan in the southeastern part of South Korea (where the clan origin is located).See Fig. 14(c), Fig. 16(b), and Table V As we can see from Fig. 14(c), the Seoul region is not special for this clan.Therefore, we see that this young, non-ergodic clan has a different mobility pattern from the stabilized, ergodic clans that follow the general trend in population flow.

Appendix I: Other Results
In this section, we discuss several figures that illustrate additional results.In Fig. 17, we show the distribution of clan density anomalies for the clans of the two Korean authors of this publication.As we illustrated in Fig. 1, Kim from Gimhae appears to be ergodic, whereas Lee from Hakseong appears to be more localized.
In Fig. 18, we examine the correlation between the distance that a clan has moved and its current ergodicity.We consider two measures of ergodicity-radius of gyration and number of regions occupied-and we also show the correlation between these two diagnostics for all clans.Some clans do not exist in the 2000 census data, and other clans only exist in one administrative region in 2000.We thus obtain radii of gyration of r g = 0 for these clans.There are 3 120 clans with r g > 0, and our calculations involving radius of gyration only use these clans.In Fig. 19, we show the distribution of clan ergodicities-using both number of regions occupied and radius of gyration-for the Czech Republic.This is like Fig. 3, in which we showed this information for Korea.In Fig. 20, we use scatter plots to examine the possible correlation between the calculated diffusion constants and distance a clan has moved.We similarly illustrate the connection between the diffusion constants and the two measures of ergodicity.
In Fig. 21, we show two "cartograms" [62] of South Korea.In these images, we distort the administrative regions in proportion to the population of people who live there.The growth of the Seoul metropolitan area over the past 40 years is clearly visible.
To examine an alternative characterization of ergodicity as the fraction of ergodic clans (see Fig. 5 in the main text), we examine radii of gyration r g versus distance to Seoul and versus distance between clan origin location and present-day centroid.We show our results in Fig. 22, which we see are qualitatively similar to those in Fig. 5.As in Figs.18 and 20, we use the 3 120 clans instead of 3 900.Consequently, we repeat the computation from Fig. 5 using this smaller set of clans.As one can see in Fig. 23, we obtain the same qualitative result.clan origin location and present-day centroid.The Pearson correlation between the diagnostics is positive and statistically significant up to 170 km (r ≈ 0.86, p ≈ 0.01) and is negative and significant for larger distances (r ≈ −0.96, p ≈ 0.005).For all of the panels, we estimate r g separately in each of 11 equally-sized bins for the displayed range of distances.The gray regions give 95% confidence intervals.(We use the same bins as in the left panel.)The correlation between the variables is positive and significant up to 170 km (r ≈ 0.99, p ≈ 0.0001) and is negative and significantly for larger distances (r ≈ −0.92, p ≈ 0.01).For all of the panels, we estimate the fraction of ergodic clans in each of 11 equally-sized bins for the displayed range of distances.The gray regions give 95% confidence intervals.

FIG. 1 .
FIG. 1. Examples of (a) ergodic and (b) non-ergodic clans.We color the regions of South Korea based on the fraction of the population composed of members of the clan in the year 2000.We use arrows to indicate the origins of the two clans: Gimhae on the left and Ulsan ("Hakseong" is the old name of the city) on the right.In this map, we use the 2010 administrative boundaries [38].See the Appendices for discussions of data sets and data cleaning.

10 - 9
FIG.2.Flux predictions from the population-product model (i.e., the special case of the gravity model with γ = 0) with α = 1 and the radiation models for jokbo 1.(a) Scatter plot of the number of clan entries in jokbo 1 versus the corresponding centroid in 2000 using the population-product-model flux with α = 1.We compute the line using a linear regression to find the fitting parameter a G ≈ 6.55(4) × 10 −11 (with a 95% confidence interval) to satisfy the expression N i = a G G i j , where G i j is the population-product-model flux and N i is the total number of entries from clan i in the jokbo.(b) The same clan entries compared to the radiation model.We compute the line using a linear regression to find the fitting parameter a R ≈ 0.049(2) to satisfy the expression N i = a R R i j , where R i j is the radiation-model flux and N i is the total number of entries from clan i in the jokbo.In both panels, we color the points using the number of administrative regions that are occupied by the corresponding clans [see Figs.3(a) and (b)].The red markers (outliers) in both panels correspond to the clan of jokbo 1 (i.e., the case i = j).

FIG. 3 .
FIG. 3. Distribution of the number of different administrative regions occupied by clans.(a) Probability distribution of the number of different administrative regions occupied by a Korean clan in the year 2000.(b) Probability distribution of the number of different administrative regions occupied by the clan of a Korean individual selected uniformly at random in the year 2000.The difference between this panel and the previous one arises from the fact that clans with larger populations tend to occupy more administrative regions.Note that the rightmost bar has a height of 0.17 but has been truncated for visual presentation.(c) Probability distribution of radii of gyration (in km) for clans in 2000.(d) Probability distribution of radii of gyration (in km) for clans of a Korean individual selected uniformly at random in 2000.The difference between this panel and the previous one arises from the fact that clans with larger populations tend to occupy more administrative regions.Solid curves are kernel density estimates (from Matlab R2011a's ksdensity function with a Gaussian smoothing kernel of width 5).

FIG. 4 .
FIG.4.Distribution of estimated diffusion constants (in km 2 /year) computed using 1985 and 2000 census data and Eq.(3).The solid curve is a kernel density estimate (from Matlab R2011a's ksdensity function with default smoothing).See the Appendix for details of the calculation of diffusion constants.

7 FIG. 5 .
FIG. 5. Fraction of ergodic clans and distance scales of clans.(a)Fraction of ergodic clans versus distance to Seoul.The correlation between the variables is positive and statistically significant.(The Pearson correlation coefficient is r ≈ 0.83, and the p-value is p ≈ 0.0017.)For the purpose of this calculation, we call a clan "ergodic" if it is present in at least 150 administrative regions.We estimate this fraction separately in each of 11 equally-sized bins for the displayed range of distances.The gray regions give 95% confidence intervals.(b) Fraction of ergodic clans versus the distance between location of clan origin and the present-day centroid.We measure ergodicity as in the left panel, and we estimate the fraction separately for each range of binned distances.(We use the same bins as in the left panel.)The correlation between the variables is positive and significant up to 150 km (r ≈ 0.94, p ≈ 0.0098) and is negative and significant for larger distances (r ≈ −0.98, p ≈ 2.4 × 10 −4 ).
K. was supported by an NRF grant funded by Korean government (2011-0015731).D.M.A. acknowledges grant #220020230 from the James S. McDonnell Foundation.S.H.L. did the majority of his work at University of Oxford.

FIG. 7 .
FIG.7.Scatter plots of the number of clan entries in jokbo 2-10 versus the corresponding centroid in 2000 using the radiation-model flux.We show our results in numerical order of the jokbo in panels (a)-(i), so jokbo 2 is in panel (a), etc.In each panel, we calculate the line using a linear regression to determine the fitting parameter a R for N i = a R R i j , where R i j is the radiation-model flux and N i is the total number of entries from clan i in the jokbo.The parameter values are (a) a R ≈ 0.062(2) [jokbo 2], (b) a R ≈ 0.0098(7) [jokbo 3], (c) a R ≈ 0.040(3) [jokbo 4], (d) a R ≈ 0.075(7) [jokbo 5], (e) a R ≈ 0.23(2) [jokbo 6], (f) a R ≈ 0.0069(5) [jokbo 7], (g) a R ≈ 0.12(1) [jokbo 8], (h) a R ≈ 0.11(1) [jokbo 9], and (i) a R ≈ 0.11(1) [jokbo 10].The red markers in panels (a), (c), and (h) correspond to the clans of the depicted jokbo, and N i | i= j=own clan = 0 for all of the other jokbo.In each case, we use a 95% confidence interval and color the points according to the number of administrative regions occupied by the corresponding clans.

FIG. 9 .
FIG.9.Fraction of entries in each jokbo versus the birth year of brides using (a) linear and (b) semi-logarithmic scales.The sudden drop on the right of each panel simply reflects the fact that women who are too young are not married yet.

FIG. 10 .
FIG.10.Probability distribution for the number of different administrative regions occupied by the 822 "new" clans that are in the 2000 census data but are not in the 1985 census data.The solid curve is a kernel density estimate (from Matlab R2011a's ksdensity function with smoothing width 1.) FIG. 11.Probability distribution for how far clans have moved in terms of the distance from the historical clan origin location to the clan centroid from 2000.We geographically identified the origin and centroid for 3 900 clans among the 4 303 clans in the 2000 census data.The rightmost bar corresponds to ≥ 500 km, and the solid curve is a kernel density estimate (from Matlab R2011a's ksdensity function with default smoothing).

FIG. 13 .
FIG. 13.Distribution of clan-centroid movement between 1985 and 2000, which we compute using r pop i,C (t = 1985) − r pop i,C (t = 2000) from Eq. (H1).(The norm • is the Euclidean norm.)The curve marked by squares weighs each clan equally, and the curve marked by circles weighs each clan by its population.We indicate distance in units of populations in 2000.(The 1985 data gives a very similar distribution.)As one can see in the inset (for which we use a doubly logarithmic scale), the maximum movement distance is larger than 400 km.However, as the main panel illustrates, most clans moved considerably shorter distances.

FIG. 15 .
FIG. 15. (a)Maximum-relatedness subnetwork (MRS)[82] of the combined population-flow network for all clans, (b) a magnified portion of this MRS that includes all regions with nonzero in-degrees, and (c) an instance of a model (inspired by the one in Ref.[83]) of a rewired version of a monocentric network (with Seoul as the center) with a rewiring probability of p = 0.4.(See the discussion in the text.)We color the directed edges toward Seoul in blue, and we color the directed edges towards other regions in blue.

FIG. 16 .
FIG. 16.(a) MRS[82] of the population-flow network for (a) Kim from Gimhae and (b) Lee from Hakseong.In panel (a), we color directed edges towards Seoul in red and directed edges towards other regions in blue.In panel (b), we color directed edges towards Ulsan in purple and directed edges towards other regions in blue.
explore clan ergodicity in more detail, and Fig.21illustrates the "convective" effect of movement into the Seoul metropolitan area.

FIG. 17 .
FIG. 17. Distributions of clan density anomalies {φ i } in 2000 over the 199 standardized administrative regions for (a) Kim from Gimhae and (b) Lee from Hakseong.The leftmost and rightmost peaks correspond to ≤ −500 and ≥ 500, respectively.Solid curves are kernel density estimates (from Matlab R2011a's ksdensity function with a Gaussian smoothing kernel of width 20).
FIG. 19.(a) Probability distribution of the number of different administrative regions occupied a Czech family name in 2009.Note that the leftmost two bars have heights of 0.17 and 0.03 but have been truncated for visual presentation.(This data was initially analysed in Ref. [25].)(b) Probability distribution of the number of different administrative regions occupied by the clan of a Czech individual selected uniformly at random in 2009.The difference between this panel and the previous one arises from the fact that clans with larger populations tend to occupy more administrative regions (c) Probability distribution of radii of gyration (in km) of Czech family names in 2009.Note that the leftmost bar has a height of 0.11 but has been truncated for visual presentation.(d) Probability distribution of radii of gyration (in km) of Czech family names of a Czech individual selected uniformly at random in 2009.The difference between this panel and the previous one arises from the fact that clans with larger populations tend to occupy more administrative regions.Observe that the distributions in panels (a) and (b) are starkly different from the distributions in panels (a) and (b) from Fig. 3. Solid curves are kernel density estimates (from Matlab R2011a's ksdensity function with a Gaussian smoothing kernel of width 5).

FIG. 20 .
FIG. 20.Scatter plots of diffusion constants from the difference between the 1985 and 2000 census data versus (a) distance between origin and population centroid, (b) number of administrative regions, and (c) radius of gyration of clans in 2000.The corresponding Pearson correlation values are (a) r ≈ 0.13 (from 3 481 clans that include all of the required information; the p-value is p ≈ 1.9 × 10 −15 ), (b) r ≈ 0.088 (from 3 481 clans that include all of the required information; the p-value is p ≈ 1.7 × 10 −7 ), and (c) r ≈ 0.097 (from all 3 481 clans in the 1985 census data; the p-value is p ≈ 1.0 × 10 −8 ).

FIG. 21 .FIG. 22 .
FIG. 21.Density-equalizing population cartograms[62] for South Korea using population data from (a) 1970 and (b) 2010 censuses[37].The coordinates are longitude on the horizontal axis and latitude on the vertical axis.The growth of the Seoul metropolitan area over the past 40 years is clearly visible (compare this to a regular map of South Korea, such as the one in Fig.1).

8 FIG. 23 .
FIG.23.Fraction of ergodic clans and distance scales of clans using only the clans that we used for the calculations in Fig.22.We obtain the same qualitative result as in Fig.5in the main text.(a) Fraction of ergodic clans versus distance to Seoul.The correlation between the variables is positive and statistically significant.(The Pearson correlation coefficient is r ≈ 0.70, and the p-value is p ≈ 0.02.)For the purpose of this calculation, we call a clan "ergodic" if it is present in at least 150 administrative regions.(b) Fraction of ergodic clans versus the distance between location of clan origin and the presentday centroid.We measure ergodicity as in the left panel, and we estimate the fraction separately for each range of binned distances.(We use the same bins as in the left panel.)The correlation between the variables is positive and significant up to 170 km (r ≈ 0.99, p ≈ 0.0001) and is negative and significantly for larger distances (r ≈ −0.92, p ≈ 0.01).For all of the panels, we estimate the fraction of ergodic clans in each of 11 equally-sized bins for the displayed range of distances.The gray regions give 95% confidence intervals.
, the clans with N 2000 admin ≥ 150 is considered to be ergodic.Based on this definition, all ten clans in the jokbo data are ergodic.

TABLE IV .
[72]on code to obtain location coordinates in latitude and longitude from the Google Maps API.We set the delay of two seconds not to exceed Google Maps API's OVER QUERY LIMIT[72]in case of a large set of locations.

TABLE V .
List of regions (which we name based on the current administrative regions) with nonzero values of in-degree in our MRSs. )