Time evolution and convergence of simple migration models

In this project, we consider two of the most fundamental mobility models, the Gravity and the Radiation models, and investigate their long-term trends. The analysis consists of determining the models’ steady states and investigating their temporal dynamics for different applications and scenarios. We find that a simple Gravity model results in two different long-term solutions, depending on its parametrization, which are independent of spatial population divisions and initial population distributions. The Radiation model on the other hand shows a strong dependency on spatial properties, due to its usage of intervening opportunities. We find that the dynamics differ significantly when it is applied to gridded population division or to population distribution divided into heterogeneous administrative units, like national counties or municipalities

In this project, we consider two of the most fundamental mobility models, the Gravity and the Radiation models, and investigate their long-term trends. The analysis consists of determining the models' steady states and investigating their temporal dynamics for different applications and scenarios. We find that a simple Gravity model results in two different long-term solutions, depending on its parametrization, which are independent of spatial population divisions and initial population distributions. The Radiation model on the other hand shows a strong dependency on spatial properties, due to its usage of intervening opportunities. We find that the dynamics differ significantly when it is applied to gridded population division or to population distribution divided into heterogeneous administrative units, like national counties or municipalities.

I. INTRODUCTION
With an increasing amount of data and computational power, researchers gain the ability to make increasingly accurate projections for future migration scenarios. Calibrated with past migration data, these mathematical models are able to make migration projections for several decades to come [1].
With the desire to run models for a growing number of time steps, the models' long-term dynamics grow more important. While models used for these kinds of tasks tend to be rather complex, they are often based on a more simple version of a mobility model, such as the Gravity or Radiation model [5,7].
The fundamental difference between the Gravity and the Radiation models is how the distance between two places impacts the migration flows magnitude. The Gravity model uses the distance as a direct variable with migration flows being inversely proportional to the distance [8]. The Radiation model on the other hand uses the number of intervening opportunities [9], correlated with the population living in the area between origin and destination.
In the first part of this paper, we investigate the convergence of both the Gravity and Radiation models. We try to determine the steady states both models will reach and investigate their stability. Furthermore, we investigate different geographical and geometrical scenarios, resembling internal and international mobility using either gridded population nets or intrinsic geometries like country or state borders. In the second part, we consider the time evolution of both models towards the steady states and interpret them with regard to real-world ap- * kluge@pik-potsdam.de plications. Last, we take a look at model setups that have been used in previous research and assign them to the different scenarios we found.

A. Gravity Model
We start our investigation with a simple Gravity model [15]. The model estimates bilateral migration flows between two places and is given by: with i and j being origin and destination respectively. The origin and destination population are given by m i and m j , whereas d ij denotes the distance between both places. A, α, β and γ indicate fitting parameters.
To investigate the population changes between two places we compare the bilateral flows in both directions, M ij and M ji , with each other: The flow-ratio, expressed in Eq.(2), is displayed in Fig. 2.
As we can see for the cases where α − β > 0 , a larger flow will originate from the higher populated areas to the lesser populated areas. Since this is the case for all times and each pair, this case will eventually lead to a homogeneous population distribution. https://doi.org/10.32388/46M58H The case where α − β < 0, results in the opposite behaviour, meaning that lesser populated areas will send more people towards higher populated areas than they receive in return. In this scenario, the system will eventually converge towards a population distribution with one area inhabiting all people.
FIG. 1: Displayed is the ratio between in-and outflows of two distinct locations, depending on the ratio of two populations (y-axis) and the difference between the two fitting parameters α and β (x-axis). The functional dependency is given by Eq.(2).
As we have seen the Gravity model can converge towards two different steady states depending on the choice of parameters α and β. The parameters A and d ij do not influence the final state of the system but only the convergence time.
It should be pointed out that both a homogeneous and a single-city distribution are steady states independently of the choice of parameters. If systems are initialized with one of these population distributions they will not change. For other mobility models, the geometrical layout of the investigated area, like population cell sizes and shape, can play an important role, these factors do not have any major impact on the Gravity model.
We would add, that hypothetically there is a third steady state for α = β. In that case, the initial population distribution is the steady state, since M ij = M ji ∀j, i.

B. Radiation Model
Next, we investigate the long-term results of the Radiation model introduced by Simini et al. [8]. Just as the Gravity model the Radiation model estimates bilateral migration flows and is given by: with i and j being origin and destination respectively. conditions. Red dots represent the origin and destination respectively. Blue dots represent population cells that will count into the intervening opportunities.
the radius being the distance between the origin and destination: with d ij being the distance between two points i and j. If we look at the same flow comparison that we used for the Gravity model, we can see that this expression can not be evaluated as easily as for the Gravity model.
As a result, we use classical stability analysis to investigate the long-term trend of the Radiation model. First, we consider the following expression which must be satisfied for the Radiation mode to be in a steady state: In the following steps, we try to find scenarios that fulfill Eq.(5) and investigate the resulting dynamics. In the Radiation model, the magnitude of the intervening opportunities has a major impact on the flows. We will consider two different versions of how these opportunities can be calculated.
First, we will consider the case with non-periodic boundaries. This case is displayed in Fig. 2(a). Areas close to the borders include a smaller number of populated areas, since the "s ij -circle" is partially outside of the investigated area. As a result, origins in the middle of the investigated area have statistically higher numbers of intervening opportunities, leading to smaller outflows.
Second, we consider a scenario with periodic boundary conditions, as displayed in Fig.2(b). Here, all "s ijcircles" include the same number of population cells (even though the actual population might vary). We consider this approach since it resembles the global migration case.
Both these setups will first be applied on a gridded population distribution with equally sized and spaced population cells. Afterwards, we repeat the simulations using heterogeneous setups, with population cells being randomly sized and distributed. Such a scenario resembles intrinsic country and state borders.

Discrete Boundary Conditions
In this subsection, we investigate the long-term dynamics of the Radiation model without periodic boundary conditions. To do so, we consider multiple, differentsized, rectangular population grids with random initial population distributions. We found that independently of size and shape, the system will always converge towards a population distribution with the highest population in the center and a decaying number of inhabitants when moving away from the center. As an example, we show in Fig. 3 multiple snapshots of a 10x10 population cell grid at different time steps. We also include graphics showing the population change.
To estimate the functional dependency between distance and population in their steady state we use: The distance d is always measured towards the center of the grid.   (6). The scenario column indicates the size of the population grid. Every population grid was randomly initialized. The right column denotes the number of time steps after which the Eq.(6) was fitted to the data.
As one can see, the parameters listed in Tab. I are of the same order and their difference decreases the larger the grids get. The reason the data points for the 10x10 and the 10x15 do not match the Gauss curve perfectly is caused by insufficient time steps. We observed that the longer the model runs, the more accurately the data fits a Gauss distribution. We assume that the small differences in the estimated parameters are the result of using discrete grids and assuming a continuous dependency between population and distance. We expect the coefficients to converge for C → ∞ with C being the number of grid cells.
We have numerically shown that a Gaussian population distribution, centered in the middle of the population grid, seems to be a steady state for the radiation model. In this section, we will have a closer look at the stability of this distribution. To investigate the stability of such a scenario, we use a linear stability analysis.
We consider a population grid with side length N × M , each grid-cell is assigned to a population m i (t) = m (x,y) (t), with x ∈ 1, . . . , N and y ∈ 1, . . . , M being the coordinates of the grid-cell. The time evolution of the population is given by: M ij represents the migration flows from origin i to destination j as defined in Eq. (3). Summarized the population at time t + 1 is given by the population at time t plus immigrants minus emigrants. The population change is therefore given by: We define the population change at place i as the function F i , which is dependent on the population of all other grid cells expressed with the vector ⃗ m, including m i . In order for a system to be in a steady state, ⃗ F = 0 must hold. As we can see in Fig.5, if a Gaussian population distribution is applied to the population grid, ⃗ F ≈ 0. The approximation is due to the discontinuity of our system. We expect that for N, M → ∞, ⃗ F = 0. We show examples of this behaviour in Fig. 5 After showing that the Gauss distribution is a steady state, we want to check its stability. We will do this in two ways, numerically and analytically. Numerically we will apply perturbation to the system and see if returns to the previously given Gauss distribution. Examples are shown in Fig. 5.
To analytically investigate the state's stability, we calculate the Lyapunov exponents for our steady state. To obtain these coefficients, we first investigate small perturbations from the Gauss distribution δ = ⃗ m(t) − ⃗ m Gauss . Afterwards consider perturbations time evolutionδ = J( ⃗ m Gauss )δ, with J( ⃗ m Gauss ) being the Jacobian matrix. If the real parts of the Jacobian eigenvalues are smaller than zero, the perturbation will decay and the system will return to the previous steady state. The Jacobian diagonals are given by with Our results show that if the Eigenvalues of the Jacobian are calculated after the Gaussian population distribution is applied are mostly positive. For larger grids, a few, comparably small, positive eigenvalues emerge. We assume that those positive eigenvalues are caused by the fact, that we use some approximation and we do not have the analytically exact steady state. Even though the state we initialize the system with is not the analytically exact steady state, we still claim that our approach can be assumed to be an approximate stable steady state.
We justify this claim by initializing the population with our Gauss function and letting the system run for a few time steps. The population change due to these additional time steps is marginal and does not change the shape of the overall state, the positive eigenvalues on the hand turn negative over time. Furthermore, we observe that by simulating more time steps, the positive eigenvalues consecutively turn negative.
Finally, we briefly investigate what happens if we apply the model to areas more complex than rectangles. The shape we use in Fig. 6 is supposed to resemble the rough shape of Mexico. One of the shapes we considered is a population distribution consisting of two square areas connected by a small corridor, see Fig. 6. As we can see, the distribution first shows two maximums, one in each of the squares but still close to the corridor. Considering longer times, the maximums will shift towards the middle of the corridor and start to resemble the Gaussian population distribution we observed for the previous cases. Other Shapes have yielded the same results.

Periodic Boundary Conditions
The reason we observe a Gauss distribution centered in the middle of the investigated area is due to the layout of the problem. Cells located at the outer parts of the investigated area have naturally smaller numbers of intervening cells compared to centred cells.
In this section, we investigate what happens if we apply periodic boundary conditions as displayed in Fig. 2.
The simulation we run with these new adjustments all converged towards a homogeneous population distribution. As one can easily confirm, a homogeneous population distribution is a steady state by itself. We assume that m i = m j , ∀i, j. As a result s ij = s ji and therefore We now show that the model converges towards a homogeneous population distribution for other initial con-ditions as well. Let us assume two populations m i and m j with m i < m j . Since there are no geographical limits for the problem anymore, s ij always contains the same number of population cells. To be able to do further calculations we assume that s ij ≈ s ji . Writing down the probability to migrate: We can assume that for at least larger distances because s ij includes a large amount of the same cells. Comparing the two probabilities we can see that P ij < P ji , meaning it's less likely for a person to migrate from a smaller city to a bigger city than vice versa. Now looking at actual migration numbers calculated from Eq. we obtain: As a result, we can see that even though the probability for an individual to migrate from a smaller city to a larger city is higher than the other way round, the actual migration numbers show that in total more people from a larger city to a smaller one than the other way round.

Heterogeneous population cell sizes
Previously, we used population areas that were divided into equally sized square cells. In reality, countries often consist of administrative areas that have neither equal shapes nor sizes.
To investigate, whether such a topological change impacts the dynamics of the Radiation model we use the following setup: We randomly create a fixed number of coordinates within a geographical range (compare Fig. 7). These random coordinates represent the center of some administrative area. Due to the randomization of the center points, such a setup could represent all kinds of differently shaped population cells. The results are displayed in Fig. 8.
As we can see, all heterogeneous population grids show different long-term solutions. Neither is there one distinct population maximum, nor is the majority of the population necessarily aggregated in the center of the grid. Furthermore, we can see even though the pattern of the t = 500 and t = 1000 plots show roughly the same population distribution, the net change in population (left plot) still show high variations even for later times (contrary to the steadily decreasing population change for the equally sized population grid). One should also point out, that even a homogeneous initial population distribution is not a steady state anymore but will instead converge to some final state depending on the given geometry of the population map.
Next, we look at whether the degree of heterogeneity in the cell location impacts the system's steady state. We create homogeneously distanced center points. Afterwards, we slightly move the center point to create a small and scaleable degree of heterogeneity. Even if only a few center points disturb the homogeneous grid, the steady state can already differ significantly from the unperturbed system yielding a Gauss distribution.
Last, we again use the randomly distributed center points, but instead of a spatially limited grid we investigate the dynamics under periodic boundary conditions. The setup for the periodic boundary conditions is the same as we used earlier. Our simulations show similar results to the scenarios without periodic boundary conditions. We are not able to find a general steady state independent from the grid's geometry, each setup has shown different steady states.

III. SUMMARY & DISCUSSION
In this paper, we investigated the long-term trends of a simple Gravity-type model and the Radiation model. We found that both models can reach different steady states depending on parameters and setting.
The Gravity model converges towards a homogeneous population distribution or towards a distribution where one cell holds the whole population. Whether the model converges towards one or the other state depends on the population parameters α and β. For the Gravity model, the geographical and geometrical setting does not impact the type of steady state the system will reach.
For the Radiation mode, on the other hand, the setting used to run the model does play a major role. We found that for equally sized population cells, the Radiation model can create a homogeneous population distribution as well. The second steady state is a twodimensional Gauss-shaped population distribution centered in the middle of the investigated area. Since the Radiation model does not have any tunable parameters, the decision of whether it creates a homogeneous or Gaussshaped population distribution is dependent on how one chooses the boundary conditions. If one chooses an isolated area to apply the model on with no population cells behind the national border or coastlines, the Radiation model converges towards the Gauss-shaped distribution. If periodic boundary conditions (resembling the geometry of the globe) are picked, the Radiation model results in a homogeneous population distribution. One could consider such a scenario when a gridded population net is available.
If one wants to use the natural border of countries or local administrative areas, one would usually not describe those as homogeneous and equally shaped pop- ulation cells. For different-shaped population cells, the Radiation model does not converge towards one independent steady state. The steady state is dependent on the geometry.
In the following, we discuss the significance of these findings on the application of both models for migration purposes. First, we point out that homogeneous population distributions are an unlikely case if we consider population distribution nowadays. The trend towards urbanization contradicts this finding as a realistic longterm result for population distributions.
Considering the Gravity model "metropol" and the Radiation model "Gauss" scenario, our biggest concern is that both of them are only capable of yielding a single population dense area/cell. It would be desirable if a model would be able to create multiple high-population areas. Both of the models have further downsides, the Gravity model projects that every person will eventually move towards a single population cell, and for areas with large numbers of population cells, e.g. the US on a county level, the result that the whole population is aggregated in a single cell seems unlikely to impossible. On the other hand, the location of this designated highpopulation cell is not predetermined and can therefore everywhere within the investigated area. The Radiation model on the hand shows a better-distributed population center, resembling a high-population city core with less populated suburbs. Contrary to the Gravity model, this population center is always located in the center of the investigated area, due to the way of calculating the intervening opportunities. As a result, the population of the border and coastal regions will always be depleted. This point already becomes important for shorter-time investigations. During the process of evolving towards a single high population area, the Radiation model is capable of having multiple "metropol" areas. These "metropol" areas will still be moving away from the border and coastal area, ignoring the fact, that these regions often tend to be highly populated.
When considering internal or international migration, neither a densely populated center nor a homogeneous population represents the reality, but for other fields, the application of these models may yield better results. Thinking about urban migration or commuting patterns within a single metropol region, the Radiation models Gauss distribution does not seem completely unreasonable. Even though there are cities with multiple city centers, the idea of obtaining a highly populated city center with less populated outskirts/suburbs seems more realistic than a county with a single center region.
This brings us to the last results we found, i.e., the Radiation model applied on a heterogeneously shaped topology. Here, we found that different results are possible. Densely populated areas can be located close to the borders of the grid as well as in the center. Furthermore, it is also possible that more than one densely populated area exists. The only similarity they show, is that large parts of the population grid are uninhabited except for smaller strongly urbanized regions.
Finally, we want to include the process of reaching a steady state in our discussion. The homogeneous Gravity model, as well as the Radiation models, show reasonable population changes and convergences. The single-city Gravity model, on the other hand, exhibits a few problems. First, it converges comparably quickly, therefore it might not be suitable for longer time simulations. Secondly, the dynamic of a high population change peak followed by an instant drop to no further population change IV. APPENDIX

A. Uniqueness of Radiation steady states in heterogeneously shaped areas
Previously, we only found that when the Radiation model is applied to a population grid with unequally shaped grid cells, the system will not end in the Gauss distribution we found earlier but will instead converge towards a steady state depending on the geometry of that population cells. Here, we briefly want to answer the question of whether the initial population distribution has an impact on the final state of the system.
We randomly choose a set of center coordinates and investigated several different randomly created population distributions. We also include a homogeneous initial population distribution. The scatter plots in Fig. 9 display the population distribution of each of those scenarios after 2000 time steps (homogeneous scenario = top left). The line plots display the difference between the random population scenarios and the homogeneous one.
As one can see, not all scenarios converge toward the same steady state. The line plot shows that some maintain a large, non-decreasing difference to the homogeneous case whereas others converge towards the same state as the homogeneous distribution does. The green circles in the two top left panels indicate some population cells where the two steady states differ. It should be pointed out that the two steady state do not differ on a large, they rather show small deviations compared to each other.
FIG. 9: Each scatter plot displays the same combination of center coordinates after 2000 time steps. The initial population was randomized but normalized to 1E6 so that the total population would remain constant for all scenarios. The top left grid was initialized with a homogeneous population distribution. In the bottom right panel, each line indicates the total difference of all population cells between the homogeneous case and a randomly initialized scenario over time.