Structural and Dynamic of Global Population Migration Network

People are the most important factors of production and the primary carriers of social culture. In addition to the impacts on the origin and destination, the volume of cross-border migrants can reflect the relationship between two countries. In fact, the migration relationship of countries is complex and multilateral, and network theories could provide a better description and more clearly show the structure and statistical characteristics. Based on the estimated bilateral migration data and disparity filter algorithm, we extract the global migration networks describing the multilateral migration relationships between 200 countries over fifty years. The results show that the global migration network during 1960-2015 exhibits obvious clustering and a disassortative property and that it is experiencing globalized and multipolarized changes during these years. Moreover, the global migration network has a typical"core-periphery"structure, and we certify its hyperbolic geometry and embed the network into a Poincar\'e disk, which could present the hierarchical structure of the global migration network and describe the status, role and evolution trends of different countries. Finally, we analyze the correlation and evolution of communities. The results show the stability of most communities, but there are still structural changes, such as the influence of the collapse of the Soviet Union on the Eurasian pattern, the intimacy between France and Africa, Canada's relative affinity between America and the British Commonwealth, and the combination of the countries around the Indian Ocean.


Introduction
In this globalization wave, the scale and diversity of international migration are substantially increasing [4,41]. In 2015, approximately 244 million people, or 3.3% of the world's population, lived in a country other than their birthplace [58,3], and this value is forecasted to double by 2050 [46]. Population migration could bring about important effects on both importing and exporting countries [22,10,23], and some scholars have used quantitative models to analyze the influencing factors, evolution patterns and trends of global population migration [33,41].
However, most of the early methods focus on the bilateral migration flow and relations between two countries, such as the conventional gravity model [51,62,48], the RUM model [30,45], and the self selection method [52,13,14]. In reality, potential migrants usually face multiple destinations, and the conventional bilateral model cannot effectively describe the actual migration selection behavior. Subsequent scholars put forward the definition of a multilateral migration barrier and introduced the structural gravity model [27,15,6] and multilateral probability model [41,40,39]. Such improvements from bilateral to multilateral analysis are meaningful yet still insufficient because the migration relationships between countries are complex and exhibit coupling. Therefore, it is necessary to take global population migration as a whole and comprehensively describe the relationship and status of countries, where the influences of other countries should not be ignored when discussing the migration flow between any two countries [41,39].
In fact, complex networks can better describe the multilateral relationships, and they show the overall structure and statistical characteristics of the system more clearly. As a new method of systematic research, complex networks have been successfully applied in many fields, including biological systems [47], cortical circuits [28], geographic maps [7] and so on. In recent years, some researchers have used the complex network method to study international migration and have achieved some preliminary results [46]. Fagiolo and Mastrorillo were the first scholars to study migration from a complex-network perspective, and they indicated that the global migration network was organized with a small-world binary pattern displaying the characteristics of disassortativity and high clustering [25]. This finding was later certificated by other scholars [20,49]. Tranos analyzed the pull and push factors behind international migration flows between OECD countries with the network method and gravity model [56]. Porat and Benguigui analyzed the degree distribution and connectivity of global migration networks and classified 145 destination countries into three classes [50]. Some scholars decomposed the world migration network into communities and analyzed its structure evolution [46], along with the glocalization, polarization and globalization of the network [19]. These research results have encouraged us to study population migration more deeply and comprehensively with the recent complex network method.
In the existing literature, some scholars used the immigration stock data to represent past flow quantities [46,20,25], while others subtracted the stock data of different eras to obtain the simplified flow data [49,50]. As an improvement, this paper extracts the networks based on the estimated bilateral flow data by a Pseudo-Bayes (PB) method for every ten years during 1960-2015, which is more recognized in the field of population research [3]. The results should thus be more reliable and scientific. This paper studies the statistical characteristics and topological structure of global migration networks composed of 200 countries/districts, with an evolution time greater than 50 years.
Geometric features have also become the focus of recent network theories. With the discovery of hyperbolic features of some real-world networks [12,29], here, we take the initiative in studying the geometric features of population migration networks. In addition to proposing the hyperbolic characteristics of the population migration network, the geometric configuration is helpful for intuitively analyzing the regional and global structures of the whole system. This paper is organized as follows: Section 2 introduces the data source and method to extract the backbone network of global bilateral migration, which is called GMN in following sections. Section 3 analyzes the skeletal construction and community dynamics of GMNs, including the changes in network statistical characteristics during 1960 to 2015, and the structure evolution. The results confirm that the GMN is a disassortative and highly clustering network, exhibiting globalized and multipolarized changes during 1960-2015. Additionally, the network represents the hyperbolic and hierarchical characteristics of international migration by embedding the countries on a 2-D Poincar disk. The positions on the disk show the status of each country/district in the global migration network, and the hyperbolic distance can indicate the migration relations between countries. Section 4 provides the conclusions and discussion.

Data source
Here, the analysis requires the data of bilateral migration flow between countries. However, from the perspective of statistics, authoritative institutions generally only provide data on the composition of immigrants (immigration stock data), such as the 'UN Global Migration database' [57] and 'World Bank Global Bilateral Migration database' [59], which cover most of the countries in the world. Some existing global migration networks are directly based on the immigrant stock data [46,20,25].
In addition, there are three common methods to estimate the bilateral migration flows based on the immigrant stock data published by the World Bank or United Nations [3]: 1. use the differences in successive bilateral stocks to estimate the corresponding migration flows [49,50,41]; 2. approximate the migration flow rates, which are then multiplied by additional data to obtain the estimated global migration flows [21]; and 3. frame the changes in migrant stocks as the residuals in a global demographic account [4,2]. Among the literature, the third method, called "demographic accounting, could estimate migration flows to match increases or decreases in the reported bilateral stocks with births and deaths during the period. Some scholars consider "demographic accounting with a Pseudo-Bayesian method as the most effective estimating method [3,8].
This paper uses the estimated bilateral migration flows with a Pseudo-Bayesian method provided by Abel between 200 countries or districts during 1960-2015 [1]. To reduce the impact of contingency, we separate the data into 6 periods: 1960-1969, 1970-1979, 1980-1989, 1990-1999, 2000-2009 and 2010-2015. Since there are only five complete years of data in 2010-2015, we have doubled the estimated flow data to match other periods.

Global Migration Network (GMN)
The topological information of networks helps us to understand the relationship structure and evolution of population migration between countries. We first construct an undirected complex network based on estimated bilateral migration flows for each period. Here, the nodes are the countries/districts, and the connection between nodes depends on whether there are migrant flows between them. The weight of the edge indicates the volume of migrants, which is the sum of emigrant and immigrant flows between two nodes.
Global migration is a complex system with complicated microstructure and evolutionary characteristics [25,50,49]. The backbone of the network offers a perspective that allows us to uncover stable large-scale patterns, such as the small-world property, heterogeneous distributions of the number of migration partners (degree), and high levels of transitive relationships (clustering) [54]. From a practical point of view, we extract the backbone network with a disparity filter algorithm (in Appendix A). This backbone network, which is called the global migration network (GMN) in the following sections, could help us better analyze the structural characteristics and evolution trends of the world's population migration system.
To illustrate the feasibility of the extracting method, we first assess the effect of inhomogeneity at the local level; for each country i with k migration routes, we calculate the Herfindahl-Hirschman index kY i (k) = k j ( wi,j si ) 2 [54,29] (in Appendix B). w i,j is the weight of the edge connecting i and j; s i = j w i,j . The local heterogeneity in the distribution of migration reveals that not all migration channels are equally significant (in Fig. 1(a), most blue triangles are below the red line of y = x, i.e., kY < k), and thus, the disparity filter can be applied in order to select only migration channels that are significant to at least one of the countries at the end of the channel.
To keep more countries, more weights and fewer links in the backbone network, with a higher fraction of remaining nodes N BB /N and weights W BB /W but a lower fraction of remaining links L BB /L, we choose a position (in Fig.  1(b)) where the number of nodes begins to be lower than the initial value as our specific indicator α s for extracting the backbone networks. N BB , L BB and W BB are the numbers of nodes, links and the weights in the backbone network, respectively, while N , L and W are those in the original flow network. The extracted backbone network is shown in Fig. 2. The color of the node indicates the community of the country/district, which is consistent with Sections 3.3 and 3.4. The color of the edge is mixed with the colors of the origin and destination nodes. The results show that the GMN became denser in the 2010s.

Hyperbolic geometry and embedding methods
In contemporary network science, hyperbolic spaces have started to receive attention because they are well suited to model hierarchical data. Some scholars have proposed that many real-world networks exhibit natural hyperbolic geometry [38,29], ranging from biology to economics, finance and trade [16,5,38,9,12,29,60]. Some machine learning methods, for instance, embeddings of graphs such as latent space embeddings, Node2vec, and Deepwalk, have found important applications for community detection and link prediction in social networks [44]. Maximilian Nickel and Douwe Kiela present an efficient algorithm to learn the embeddings based on Riemannian optimization [44]. This method is based on the Poincar ball model, as it is well suited for gradient-based optimization. Here, we embedded the GMN in all time periods , using a 2-dimensional Poincar disk, with a learning rate of 0.1 and a negative sample size of 30. This setup produced hyperbolic embeddings in which each node i-a country in the migration embedding-has radius r i and angle θ i . Nodes of small radius hold central positions in the circularly arrayed hierarchy. The hyperbolic distance (depending on angle and radius) between two nodes quantifies their migration relation. Please find the comparative evaluation of hyperbolic embedding and Euclidean embedding in Section 3.3.

Clustering and disassortative network
Analysis from the perspective of complex network attributes can help us understand the basic characteristics of global migration networks (GMNs). Fig. 3 (a) shows the evolution of the strength and numbers of edges for the GMN from 1960 to 2015. Obviously, in general, the number of edges and the sum of weights exhibit a growth trend over time, which can also be observed in Fig. 2.  shows the degree distribution of the backbone networks for all periods in gray and blue. This figure reveals that only a few nodes in the network have a large degree, while most other nodes have a small degree, which appears to be the power law distribution. Thus, we use the backbone network of 2010-2015 as an example to perform power rate fitting for the degree distribution by the nonequidistant bin method, shown in Fig. 3 (b) with red dots and line. It was found that the network essentially conformed to this feature. In other words, most countries have a single direction of population migration, while a small number of countries have population exchanges with many countries; this means that population migration exhibits local concentration.
The clustering coefficient measures how connected a nodes neighbors are to one another; the average shortest path reflects the difficulty of one node connecting to another node in the network. The clustering coefficients and the average shortest path (in the maximum connected subnet) of GMNs are shown in Fig. 3 (c). The clustering coefficient exhibits an upward trend, while the shortest path follows a downward trend, exhibiting the enhancement of small world attributes. The mobility barriers to international migration appear to be decreasing in recent years, while the global migration relations have become closer, which also indicates the increase in network clustering during 1960-2015.
The degree correlation in complex networks reflects the connection preference of nodes in the network, as defined below: where k nn indicates the average degree of the first neighbors of nodes with degree k. If µ > 0, then this is an assortative network; similarly, µ = 0 indicates a neutral network, while µ < 0 indicates a disassortative network. We calculate the degree correlation of each network by degree correlation function, and their µ values are all negative during 1960-2015. The degree correlation of the 2010-2015 GMN is shown in Fig. 3 (d), with µ = −0.43 as the example. The GMNs in other years exhibit similar characteristics. This method shows that all of the networks have the characteristics of negative matching, indicating that nodes with lower degrees are more likely to be connected with nodes with higher degrees. The reason may be that the migration of most small countries shows preferential movement to several large countries that have survival advantages or economic advantages. In fact, the disassortative characteristic of the international migration network has been verified in some existing studies [49,25]. In addition, the international oil trading network shows a disassortative feature in which countries with fewer trading partners tend to develop oil trading relations with countries with more trading partners [61]. Table 1 shows the statistical characteristics and evolution trends of global migration networks(Note that: APL: Average path length; CC: Clustering coefficient; ND: Node degree; NS: Node strength). We can observe an increase in the edges and weights, along with increased network density and connectivity, which means that many countries have become closer to each other in the GMNs.

Globalized and multipolarized networks
The cumulative distribution of degrees and weights shows that, on the whole, they both tend to be relatively flat. This means that many top countries are  1960-1969 1970-1979 1980-1989 1990-1999 2000-2009   reducing their proportions of edges and flows, while other countries with low flows are experiencing relatively rapid development ( Fig. 4 (a) and (b)). The Gini coefficients of degree and weight both exhibit an overall downward trend, and the entire network becomes more balanced over time (Fig. 4 (c)).
In recent years, some scholars have proposed the "globalization of migration" hypothesis and emphasized both the progressively increasing number of countries involved in global migration and the diversification of origins and destinations [42,19]. Some other scholars offered an alternative understanding, suggesting that in recent years, globalization did not contribute to an overall increase in mobility possibilities but instead widened the gap between rich and poor countries [53,34,31], leading to polarization of the global migration networks.
From the network perspective, we could analyze the "globalization" or "polarization" trends based on comparing the global and local connectivity in migration communities. There are many ways to explore the communities of a complex network, such as the GN algorithm [43] based on network topology and Potts model [26] based on network dynamics. In this paper, the Louvain algorithm [11] based on modularity, which is rapid and exhibits an obvious clustering effect, is adopted. The algorithm divides each round of calculation into two steps: in the first step, the algorithm scans all nodes, traverses all neighbor nodes of the node, and measures the modularity benefit of adding the node to the community of its neighbor node; it then selects the corresponding neighbor node with the highest modularity gain and joins its community. This process is repeated until the results are stable. During 1960-2015, the modularity value Q is within 0.66-0.76, which proves the validity of clustering (green line in Fig. 4(d)). Obviously, the result of clustering is relatively stable over the most recent fifty years. According to the composition of members, we define ten typical communities. 1. America: including the countries in North America, Central America, South America and the Caribbean; 2. French related: including France and other neighboring European countries, as well as French territories and former colonies around the world; 3. British Commonwealth: including some original and current Commonwealth countries, including Australia, New Zealand, Canada, etc.; 4. Indian Ocean: centered on India, including South Asia, North Africa, Southeast Asia and other countries close to the Indian Ocean; 5-7. most sub-Saharan African countries are divided into three communities: East-Middle, Western and Southern Africa; 8. Europe: some countries in western, central and eastern Europe; 9. former Soviet Union: including Russia and some former Soviet countries; and 10. East Asia: mainly East Asian countries, including some Southeast Asian countries in early times, and Russia and some former Soviet countries in recent years. Here we use the external-internal index (E-I index) to measure the comparison of local and global cohesion, which is widely used in group embeddedness [36,32,19]. We define the E-I index of GMNs as The "internal" edge connects the two nodes in the same community, and the "external" edge connects the nodes from the different communities. EK and EW are the sums of external degrees and weights for all nodes, respectively; IK and IW are the sums of internal degrees and weights for all nodes, respectively. The E-I index ranges from 0 to 1. Smaller E-I index values indicate stronger connectivity between communities; larger E-I index values indicate stronger connectivity within the community and show that the community is more independent. Fig. 4 (d) shows a downward trend of the E-I index with respect to both degrees and weights. This figure indicates the continuous growth trend of crosscommunity connection and to some extent proves the significant trend of globalization in GMNs over the past fifty years.  Fig. 6. E-I index (degree) for each community from the 1960s to 2010s. The countries/districts having the largest degrees in the communities are listed in the table.
In addition, Fig. 6 shows E-I index values for ten typical communities; the countries/districts possessing the largest degree in the communities are listed, which could be regarded as the center of communities. For most communities, the central country/district is relatively stable and unchanged for these 50 years or is only adjusted between neighboring countries. It is worth noting that there are mergers and splits of communities 8-10, Europe, the former Soviet Union and East Asia, which will be described in detail in Section 3.4. Here, dark green indicates the smaller E-I index values for the community in the corresponding time period.
The results present that over the past 50 years, the migration relation between different communities has become closer, which is reflected in the overall decline in E-I index values presented in Fig. 6. In particular, the two communities of the former Soviet Union and Indian Ocean are typically introverted, with most of the migration flow coming from the "internal edges" connecting the community members; the communities of America, French related and Europe (since the 2000s) are extroverted, where the cross-community migration relation is greater than that within the communities. On the whole, the number of extroverted communities is increasing over time, and this also shows that for the potential immigrants, the possible moving routes among communities become more abundant. In contemporary times, the number of communities with E-I index values below 50% (dark green grids) has increased from ONE to FOUR, which indicates that the GMNs became more globalized and multipolarized from the 1960s to 2010s.

Hyperbolic and hierarchical network
To intuitively analyze the topological structure of the network, we try to separately embed GMNs into Euclidean and hyperbolic planes. We consider the least squares error function used in [18]. After unified measurement, the cumulative errors of two spaces are revealed in Fig. 7 with dotted lines (details in Appendix C). In addition, we compute the embedding score in both Euclidean space and hyperbolic space for comparison (considering the error function performance, we only embed the network in Euclidean space by nonclassical MDS). The result is shown in Fig. 7 with solid lines. This figure shows that all errors in hyperbolic embedding are lower than those in Euclidean embedding; according to the Score, hyperbolic embedding also offers a more professional performance in the expression of data size relations. This result indicates that the GMNs exhibit a significant hyperbolic characteristic during 1960-2015. Here, we embed the GMNs into hyperbolic space and finally acquire the locations of 200 countries/districts on a 2-dimensional Poincar disk.
Taking the GMNs of 1960-1969 and 2010-2015 as an example (Fig. 8), the Score of over 90 represents the effectiveness of the embedding (Fig. 7). In Fig.  8, each node represents a country/district. The size of the node expresses its degree in the GMN, and the locations of countries originate from the mobility distance matrix. Briefly, the distances represent their relations in the global migration network, and the edge and its weight are the same as those in the backbone network. For clarity, only the names of nonperipheral countries are shown, including the nodes with distance less than 0.98 from the origin of the coordinates. The color of the node indicates the geographical location of the country/district. First, the figure indicates that the GMN presents an obvious hierarchical structure. In the 1960s, the network is sparse, while in the 2010s, 200 countries have closer and more complex migration relations, which also conforms to the common law of global integration. The migration of population in the 1960s primarily occurred within the regions or communities, and there are fundamentally fewer typical countries in the center of the Poincar disk ( Fig. 8 (a)). The D.R. of the Congo (COD) is the second largest country in Africa. Once a Belgian colony, it became independent in 1960 and became a link in the global migration network between France (FRA) in Europe and African countries such as Sudan (SUD) and Madagascar (MDG).
Furthermore, for the period of 2010-2015, the network structure is more complex, and the hierarchy is more obvious. Here, the United States and Canada (with large degrees of k = 53 and k = 30, respectively), which have closer population migration relations with other countries in various regions of the world, are more centrally located in the disk. France and the United Kingdom, which mainly connect local communities such as European and some African countries, have been slightly more marginal since the 2010s. Although the degrees of Portugal and Yemen are not large (k = 8 and k = 2, respectively), their locations indicate contributions to connecting the migrations of several continents.
Additionally, it is shown that the hyperbolic distance between countries/districts is not entirely determined by the geographical location. It represents the correlation of embedding distances and geographic distances, which is significant during 1960-2015 (with sig.≈ 0). This means that hyperbolic distance d h is positively correlated with geographical distance d g but encodes more than purely geographical information.

Structure evolution of global migration network
To assess the network structure more clearly, we analyze the correlation of the network communities with time. Fig. 9 shows the matrix of Jaccard similarity coefficients of ten typical communities covering 90-97% of the countries/districts. In general, the composition of the members of each community is stable, and the characteristics related to geographical location are shown (Fig. 10). Over more than 50 years, the central countries/districts of most communities have not changed. The green color indicates greater correlation, along with the higher coincidence of the members for the cluster between two eras. In contrast, the yellow color indicates that the structure of the communities changed greatly during this time.
Focusing on the yellow grids in Fig. 9, combined with the specific composition of each community in Figs. 5 and 10, we found some structural evolution of global migration networks during the past 50 years.
The community of former Soviet Union: the ninth community, centered on Russia, was an independent cluster in the 1960s-1970s; in the 1980s-1990s, it merged into the community of Europe centered on Germany. Such structural changes may be related to the collapse of the Soviet Union in 1991, when it pursued the policy of deporting the nonnative population, together with the boom of immigrants from the East into Western Europe [37,35]. After 2000, the former Soviet Union cluster left the Germany group and merged into the community of East Asia. In fact, since then, Russia gradually replaced Hong Kong SAR as the new center of the community. The map also shows that Russia and these former Soviet Union countries have had a closer relationship with East Asia in GMNs since 2000 ( Figure 10).  The country of Canada: once belonging to the British Commonwealth, Canada had a close relationship with Hong Kong SAR, and they were in the same community in the 1960s-1970s. Canada changed during 2000-2010 to join the community of America, which contained most of the American countries (Fig. 10). In fact, some scholars have certified the relationships of the countries in Latin America and North America, including the United States [24,17]. However, beginning in 2010, Canada left the United States community and became the new center of the British Commonwealth community; it is also the third closest country to the center on the hyperbolic plane, after the United States and French Guiana (Fig. 8).
The community centered on France: the structure of the second community centered on France also substantially changed. In the 1990s, 68% of its members belonged to European countries or their territories, but in the 2010s, the European members only accounted for 58%. In contrast, in the 1960s, African countries accounted for only 15% of this community, but after 2010, the proportion of African countries increased to 42%, which also shows that France, as the representative and center of the community, became increasingly close to African countries in the global migration network. This trend should be closely related to the influence of language and historical colonies [55,41].
The countries including Malaysia, Singapore and Indonesia: in the 1960s-1970s, Malaysia, Singapore and Indonesia were in the East Asia community, with Hong Kong SAR as the center, but since the 1980s, these three countries have been transferred to the community centered on India, which has greatly impacted the East Asia community and greatly reduced the number of its members. After 2000, the former Soviet Union community was merged and the center of the cluster was adjusted to Russia, which greatly changed the structure of the East Asia community once again.

Conclusion
Global population migration is a typical complex system. At the micro level, each potential migrant makes a rational decision on "whether" and "where" to migrate according to the diversity utility function. Although the individuals are heterogeneous, specific migration patterns and evolution rules are continually emerging on the macro level.
The migration relationship between countries is complex and multilateral, and network theories could provide better description and more clearly exhibit its structure and statistical characteristics. This paper constructs an undirected global migration network (GMN) based on estimated bilateral migration flows during 1960-2015. The GMNs display the characteristics of disassortativity and high clustering with a typical power law in degree distribution. In the most recent fifty years, the network density and clustering have been increasing; the Gini coefficient of the degree and weight both exhibit an overall downward trend; and the entire network becomes more balanced and exhibits greater connectivity with time.
From the network perspective, we analyze the evolution trend of international migration by comparing the global and local migration connectivity in communities. On the whole, the number of extroverted communities is increasing over time. This observation indicates the continuous growth trend of cross-community connection and, to some extent, proves the significant trends of "globalization" and "multipolarization" in the global migration network since the 1960s.
The existing literature does not discuss the geometric characteristics of the population migration network. This paper indicates that the GMNs exhibited a significant hyperbolic characteristic and hierarchical structure during 1960-2015, which is becoming more obvious these years. We embed the GMNs into hyperbolic space and finally obtain the locations of 200 countries/districts on a 2-dimensional Poincar disk.
Finally, we analyze the correlation between network communities and the structural evolution of GMNs with time. In general, from 1960 to 2015, the composition of the members of each community remained stable, and the central countries/districts of most communities did not change. In addition, we still find some changes: the former Soviet Union community merged into Germany during the 1980s-1990s, which could be related to the collapse of the Soviet Union, and it left the German group and replaced Hong Kong SAR as the center of the East Asian community after 2000; the community centered in France reduced the proportion of members from Europe, and more African countries, especially those with colonial relations and labor contracts with France, gradually joined the group beginning in the 1960s; some southeastern Asian countries such as Singapore, Malaysia and Indonesia were in the East Asian community but transferred to the Indian Ocean community centered in India since the 1980s.
This paper provides a creative way to analyze the structural, statistical, and geometric characteristics and hierarchical structure of the population migration network. With respect to complex human migration behavior, it is far from sufficient to analyze only the migration flow data. In future research, we will consider the economic, social and policy factors that affect the decision-making of potential migrants and research the features and evolution of population migration patterns more comprehensively and scientifically.

Conflicts of Interest
The authors declare no conflict of interest including any financial, personal or other relationships with other people or organizations.

Funding Statement
This work was supported by the Chinese National Natural Science Foundation (71701018, 61673070), the National Social Sciences Fund, China (14BSH024), and the Beijing Normal University Cross-Discipline Project.

A Disparity filter algorithm
The disparity filter proceeds as follows. We first normalize the weights of edges linking node i with its neighbor j as p i,j = ω i,j /s i , with s i = j ω i,j being the strength of node i and ω i,j the weight of the edge connecting i and j. For each migration channel of a given country i, we compute the probability α i,j that the link takes the observed value p i,j according to the purely random null model. By imposing a significance level α, we can determine the statistical significance of a given migration channel by comparing α i,j to α. Therefore, if α i,j > α, the flow through that migration channel can be considered compatible with a random distribution (with the chosen significance level α) and is thus discarded. The statistically relevant channels are those that satisfy for at least one of the two countries i and j. k represents the degree of node i. By applying this selection rule to all of the links in the network, we find the backbone, a new graph containing, in general, fewer links and nodes, as the GMN in this paper. However, the number of links and nodes removed depends on the value of the significance level α. To find the appropriate value of α, it is convenient to plot the fraction of remaining nodes N BB /N and the fraction of remaining weights W BB /W in the backbone vs. the fraction of remaining links L BB /L for different values of α. As the filter becomes more restrictive, the number of links decreases while keeping almost all nodes until a certain critical point, after which the number of nodes begins a steep decay. To retain more countries, more weights and fewer links, we choose a point where the number of nodes begins to be lower than the initial value as our specific indicator α s for extracting the backbone networks [54].

B Inhomogeneities of the networks
To calculate inhomogeneities at the local level, for each country i with k migration routes, we calculate the Herfindahl-Hirschman index (HHI) kY i (k), which is extensively used as an economic standard indicator of market concentration, and it is also denoted as the disparity measure in the complex networks literature where ω ij is the total flow between countries i and j and s i = j ω i,j is the strength (aggregated migration) of country i. If country i distributes its migration homogeneously between its migration partners, then kY i (k) = 1; in the opposite case, if all its migration is concentrated on a single link, then kY i (k) = k.

C Effectiveness of hyperbolic embedding method
The method proposed by Maximilian Nickel and Douwe Kiela is based on the Poincar ball model, as it is well suited for gradient-based optimization [44]. In particular, let B d = x ∈ R d | x < 1 be the open d -dimensional unit ball, where · denotes the Euclidean norm. The Poincar ball model of hyperbolic space then corresponds to the Riemannian manifold (B d , g x ), i.e., the open unit ball equipped with the Riemannian metric tensor where x ∈ B d and g E denotes the Euclidean metric tensor. Furthermore, the distance between points u, v ∈ B d is given as Note that Equation 6 is symmetric and that the hierarchical organization of the space is solely determined by the distance of nodes to the origin. Due to this self-organizing property, Equation 6 is applicable in an unsupervised setting where the hierarchical order of objects such as text and networks is not specified in advance. Remarkably, Equation 6 therefore allows us to learn embeddings that simultaneously capture the hierarchy of objects (through their norms) as well as their similarity.
We embedded our GMN in all time periods , using a 2-dimensional Poincar disk, with a learning rate of 0.1 and a negative sample size of 30. This setup produced hyperbolic embeddings in which each node i-a country in the migration embedding-has radius r i and angle θ i . Nodes of small radius hold central positions in the circularly arrayed hierarchy. The hyperbolic distance (depending on angle and radius) between two nodes quantifies their migration difference. To further explain the embedding effect, the general form of the function is as below: where δ j,k is the dissimilarity between nodes j and k and d j,k represents the embedded distances. To transfer the migration matrix to the dissimilarity matrix, for every weight ω i,j , we use 1 − (ω i,j /ω max ) (ω max denotes the maximum weight in the matrix) to replace the original data. In Euclidean space, we use two kinds of regular MDS (multidimensional scaling) methods, namely, the nonmetric MDS (hereinafter referred to as NMM) and nonclassical MDS (hereinafter referred to as NCM), to embed the data for comparison.
The error function described in the previous section calculates the cumulative difference between the embedded distance and the actual data. However, it also concerns another issue: whether the two countries with closer relations are actually closer to each other than other countries after embedding. Here, we propose a scoring scheme to assess this possibility. For any two edges l i,j and l m,n that exist in the network, where i = m and j = n, and with corresponding embedding distances d i,j and d m,n , we calculate that We repeat this random selection n times and obtain the following scores Score = 100 · S N . Since there are almost no consistent original migration data or embedding distances between different countries, (l i,j −l m,n )(d i,j −d m,n ) = 0 is almost nonexistent. Thus, Score indicates the probability of meeting the required rules for any two links. The closer the Score is to 100, the better the embedding distance can interpret the size relationship in the original data.

D Results of community division
We use the Louvain algorithm to divide the community of all networks. For example, the complete community results from 2010 to 2015 are presented in Table 2.

Group
Name Members