Demography and the emergence of universal patterns in urban systems

The study addresses why and when city size distributions follow Zipf’s law and why empirical data often deviate from it, especially for very small and very large cities. Prior models often invoke proportional growth (Gibrat’s law) or preferential attachment at a mesoscopic level, without deriving distributions from fundamental demographic processes. The authors propose an evolutionary population-dynamics framework grounded in births, deaths, and migration between cities. They hypothesize that Zipf’s distribution emerges as a neutral law when selection is absent and key symmetries in demographic rates (vital and migration rates) are restored through long-run stochastic fluctuations. Deviations from Zipf’s law then encode information about preferences and selection in demographic choices. The purpose is to derive city size distributions from first principles of demography, reconcile observed deviations, and clarify the conditions under which universal patterns appear in urban systems.

Power-law type frequency distributions appear widely in complex systems and have been linked to multiplicative processes and network structures. Urban Zipf’s law has been observed since Auerbach (1913) and popularized by Zipf (1930s–1940s). Empirical studies show substantial and systematic deviations from Zipf’s law at distribution tails and variation of the Pareto exponent across countries and time (e.g., China, US, Europe). Critiques emphasize statistical ambiguity and non-universality. Prior generative models rely on Gibrat’s law, preferential attachment, or other mesoscopic mechanisms. Cristelli et al. (2012) stress systemic properties like coherence and screening beyond mere power-law form. The literature also includes gravity-type migration models and urban growth histories indicating temporally shifting regional growth advantages. This work integrates these strands by deriving size distributions from demographic rates and migration flows, identifying neutrality and symmetry as core explanatory principles.

- Demographic framework: City i’s population evolves by births, deaths, and migration: N_i(t+1) = N_i(t) + ν_i N_i(t) + Σ_{j≠i} J_{ij}, with ν_i = b_i − d_i and J_{ij} migrants moving i→j per period. Set Δt = 1. Migration currents are expressed as rates J_{ij} = m_{ij} N_i(t). The system is written in matrix form N(t) = A(t) N(t−1), with A_ii = 1 + ν_i − Σ_{k≠i} m_{ik} and A_{ij} = m_{ji} for i≠j. The population structure vector x has components x_i = N_i / N_T. - Ergodic results: For strongly connected, aperiodic environments, the strong ergodic theorem implies convergence of structure to the leading eigenvector of A when A is static. In stochastic environments, distributions of x converge to stationary distributions independent of initial conditions (weak/strong stochastic ergodic theorems). - Static environment solution: In constant A, N(t) ∝ λ_0^t e_0, where e_0 is the positive leading eigenvector (eigenvector centrality). Convergence rate is set by τ = 1/ln(λ_0/λ_1). Using US intercity migration data, τ is several centuries using census-based flows and shorter using IRS tax-return flows. - Migration decomposition and diagonalization: Decompose flows into symmetric and antisymmetric parts by writing J_{ij} = [(s_{ij} + δ_{ij} N_j(t))/2] N_i(t), with s_{ij} = s_{ji} and δ_{ij} = −δ_{ji}. This connects to the (symmetric) gravity law when δ_{ij}=0 and s_{ij}=G N_i f(d_{ij}). Only the antisymmetric component affects net growth: N_i(t+1) = [1 + v_i − δ_i] N_i(t), with δ_i = Σ_j δ_{ij}. - Self-consistency and rotational symmetry: There exists a structure x* such that v_i = v* = Σ_j δ_{ij} x*_j, with Σ_i δ_{ij} x*_i = 0. The antisymmetric matrix δ acts as an infinitesimal generator of rotations on x. Specific differences in relative growth rates break rotational symmetry and yield non-Zipfian steady hierarchies in deterministic settings. - Stochastic dynamics and symmetry restoration: Let ε_i(t) = v_i(t) − δ_i(t) be zero-mean relative growth fluctuations with variance σ^2 over long horizons; temporal averaging can restore rotational symmetry (neutrality) so that ⟨ε_i⟩→0 for each city. Under these conditions, the structure components follow geometric Brownian motion without drift: d x_i / x_i = σ dW. The associated Fokker–Planck equation for P(x_i,t) admits stationary solutions; with zero probability current (J=0) and size-independent σ^2 (Gibrat’s law), P(x) ∝ 1/x^2 (Zipf). With constant nonzero current, P has exponential form over x. If σ^2 depends on N (e.g., σ^2 ∼ N^{−γ}), the Pareto exponent shifts to 2−γ. - Boundary conditions: Stationarity and zero-current require reflecting/absorbing-like conditions at lower and upper bounds x_l, x_u. The lower bound x_m is crucial; without conserved current, solutions drift to a lognormal concentrated at small x. Implementing a ‘ghost’ of small cities at the lower boundary preserves probability current and Zipfian scaling for observable cities. - Numerical experiments: Simulate urban systems (e.g., 100 cities) under (i) static, (ii) nonlinear deterministic, and (iii) stochastic environments, comparing trajectories and rank-size distributions. Convergence to e_0 in static A; stronger fixed hierarchy under nonlinear deterministic dynamics; and long-run fluctuations around Zipf’s law under stochastic symmetry restoration. KL divergence D_KL(P||P_Z) quantifies deviations. - Data construction (US): Build A matrices using US Census Bureau (USCB) migration/mobility tables (decennial) and IRS annual migration proxies (tax filings), aggregated from counties to MSAs via NBER crosswalk. Net vital rates inferred from residual population change not explained by domestic migration. Empirically, annual inter-MSA migration ≈1.8%; migration probabilities are lognormally distributed. Lower-bound handling resets x_i slightly above x_m when violated and renormalizes Σ_i x_i=1.

- Static environments lead to non-Zipfian steady states: City sizes converge to the leading eigenvector of the environment matrix (eigenvector centrality), independent of initial conditions, with convergence time τ = 1/ln(λ_0/λ_1) estimated to be several centuries using US flows. - Deterministic nonlinear environments (with antisymmetric migration components) also yield stationary hierarchies distinct from Zipf’s law and unfold more slowly, producing stronger urban hierarchies. - Stochasticity can restore rotational symmetry (neutral dynamics): When city-specific growth advantages reverse over long times so that ⟨ε_i⟩=0, the structure follows geometric Brownian motion without drift. The zero-current stationary solution is Zipf’s law (P ∝ 1/x^2), emerging as a long-run time average; instantaneous snapshots fluctuate around it. - Pareto exponent depends on fluctuation scaling: If growth-rate variance σ^2 depends on city size (σ^2 ∼ N^{−γ}), the exponent changes (P ∼ 1/N^{2−γ}), potentially destroying scale invariance if not power-law. - Two attractor classes: (i) zero probability current (Zipf), and (ii) constant current (exponential in x), both reached at long times t ≫ 1/σ^2. - Boundary conditions are decisive: The lower cutoff x_m controls tail behavior; maintaining zero current effectively requires a Zipfian ‘ghost’ of many small cities, explaining empirical deficits of very small towns and deviations in high ranks. - Empirical support (US): Decennial data for the largest 100 US cities (1790–1990) show that cumulative temporal averaging approaches Zipf’s law, while single-year distributions deviate; in recent decades, persistent deviations grow, largely due to faster growth of mid-sized and large metros relative to the very largest. - Timescales: Neutrality-induced symmetry and Zipf’s emergence operate on very long timescales (centuries to millennia), consistent with slow reversal of regional advantages.

The findings explain Zipf’s law as a neutral outcome of demographic dynamics when selection (persistent growth advantages) averages out and rotational symmetry in the structure vector is restored by stochastic fluctuations. In integrated systems with strong intercity migration, symmetry need not be perfect; asymmetric flows can maintain coherent but non-Zipfian hierarchies. Thus, Zipf’s law is not a unique signal of integration but one among several possible integrated states. Deviations from Zipf’s law capture information about selection—preferences for certain places due to socio-economic opportunities, shocks, or demographic changes—and can be quantified by KL divergence. The framework clarifies systemic properties such as coherence and screening: zero-current conditions make subsamples non-representative and protect the rank of the largest city. Boundary conditions, especially for small cities, critically shape observed distributions and reveal hidden departures from true scale invariance. Empirically, US city sizes averaged over long periods approximate Zipf’s law, while recent decades exhibit increasing, informative deviations consistent with selective growth of specific regions and size classes.

Starting from general demographic equations for births, deaths, and migration, the paper derives urban size distributions and explains how universal patterns, including Zipf’s law and gravity-like flows, emerge as long-run averages when demographic symmetries are restored by stochasticity. In static or selectively biased environments, city size hierarchies converge to eigenvector-defined structures rather than Zipf’s law. Zipf’s law acquires the status of a neutral, maximum-entropy distribution under zero probability current, with deviations measuring the information content of preferences and selection. The approach generalizes to other complex systems governed by multiplicative growth and type transitions. Future research should: (i) quantify selection signals via information measures across countries and epochs; (ii) refine boundary-condition modeling for small settlements; (iii) incorporate spatial frictions and network evolution endogenously; and (iv) test timescale predictions and variance-scaling effects with high-frequency, long-span migration and vital-rate data.

- Dependence on long-run averaging: Emergence of neutrality and Zipf’s law requires very long timescales (centuries or more), which may not hold in rapidly transforming systems. - Boundary-condition sensitivity: Zero-current Zipfian solutions depend critically on lower-bound handling and the implied ‘ghost’ of small cities, which is unrealistic; mis-specification affects tail behavior. - Assumptions on stochasticity: Gaussian growth-rate fluctuations and stationarity over long periods are approximations; size-dependent variances can alter exponents. - Model scope: The framework abstracts from economic mechanisms, policy, and endogenous network evolution; asymmetric, persistent shocks or structural changes can sustain non-neutral dynamics. - Data aggregation and coverage: US analyses aggregate counties to MSAs and rely partly on IRS tax data (≈85% coverage), potentially biasing migration estimates; international migration is absorbed into vital rates implicitly.