# Decomposing geographical and universal aspects of human mobility

Louis Boucherie  DTU Compute, Technical University of Denmark  Center for
Social Data Science, University of Copenhagen  Benjamin F. Maier  DTU Compute,
Technical University of Denmark  Center for Social Data Science, University of
Copenhagen  Sune Lehmann Email<EMAIL_ADDRESS>DTU Compute, Technical University
of Denmark  Center for Social Data Science, University of Copenhagen

###### Abstract

Driven by access to large volumes of detailed movement data, the study of
human mobility has grown rapidly over the past decade [1]. This body of work
has argued that human mobility is scale-free [2, 3, 4], it has proposed models
to generate scale-free moving distance distributions [5, 6], and explained how
the scale-free distribution arises from aggregating displacements across
scales [7]. The field of human mobility, however, has not explicitly addressed
how mobility is structured by the constraints set by geography – that mobility
is shaped by the outlines of landmasses, lakes, and rivers; by the placement
of buildings, roadways, and cities. Based on unique datasets capturing
millions of moves between precise locations, we show how separating the effect
of geography from mobility choices, reveals a universal power law spanning
five orders of magnitude ($10\,\mathrm{m}$ to $1,000,000\,\mathrm{m}$). To do
so we incorporate geography via the ‘pair distribution function’, a
fundamental quantity from condensed matter physics that encapsulates the
structure of locations on which mobility occurs. We show that this
distribution captures the constraints that geography places on human mobility
across length scales. Our description conclusively addresses debates between
distance-based and opportunity-based perspectives on human mobility: By
showing how the spatial distribution of human settlements shapes human
mobility, we provide a novel perspective that bridges the gap between these
previously opposing ideas.

## Introduction

Human mobility is a fundamental component of societies, with constant flows of
individuals and goods circulating through cities and landscapes. This flow of
individuals moving from place to place underpins our social, economic, and
cultural lives [8]. Therefore, understanding mobility patterns enables us to
mitigate the spread of infectious diseases [9, 10], reduce pollution [11], or
organize urban planning [12]. Despite the fact that human movement is both
shaped by and shaping our geography [13, 14], these two processes have not
heretofore been considered together [15].

Geography trivially constrains human mobility. Consider two simple examples:
Two houses cannot occupy the same physical space, thereby imposing a minimum
distance for any movement between them. On a larger scale, residents of an
isolated island of a $10\,\mathrm{km}$ radius cannot travel distances greater
than $20\,\mathrm{km}$: they would end up in the sea.

Previous work has shown that human mobility is highly structured, with power
law distributions of distance [2, 3] and trip frequency [16], inflows between
administrative boundaries [5, 6], conservation laws [17], and scales
corresponding to structures of the built environment [7]. Simultaneously, the
structure of the built environment is well understood, from the hierarchical
organization of human settlements [18] at the scales of countries [19, 20] to
street-level networks [21, 22], as well as the fractal geometries of cities
[23, 24], and how they developed [13].

A key concern, however, is that the work cited above arises from two
literatures that have evolved separately. The recent surge of data-driven work
on human mobility has not explicitly addressed how the observed mobility data
is structured by the constraints set by geography [15]. In that sense, there
is currently a disconnect between the large-scale data-driven human mobility
research [1] and the community of geographers [13, 20].

The absence of a comprehensive geographical framework in mobility studies has
led to the emergence of two main theories for human movement [1]: The first
theory considers physical distance to deter mobility [25, 26], a concept often
modeled through gravity models [5]. The alternative viewpoint instead suggests
that human movement does not depend on distance, but rather on the quantity of
opportunities nearby [27, 28].

There have been scattered efforts to incorporate a notion of geography in
human mobility, but this work has been limited by incomplete data, e.g.
relying on coarse-grained data such as data in grids [29], administrative
units [30], or artificial geographical space [31], obscuring smaller scales
and reducing geography to discrete density maps [32]. The use of neural
embeddings also achieves a similar effect [33], and similarly leads to the
emergence of gravity laws [34].

Therefore, the question of how geography shapes human mobility still lacks a
comprehensive and satisfying answer. Below, we address this gap and present a
new framework to understand geography and the distribution of locations
(addresses, points of interest, stop-locations). We argue here that all one
needs to know about geography – the fine structure, rivers and lakes, the
country outlines – is captured in the locations where people stop. This set of
locations represents a kind of effective geography.

Our primary analysis is based on a previously unused dataset comprised of 36
years of registry data on 39M migration between addresses in Denmark,
pinpointed with uncertainty of only $2\,\mathrm{m}$. This resolution far
exceeds the state-of-the-art $50\,\mathrm{m}$ accuracy typical of GPS data.
Below, we extend the analysis to day-to-day mobility and point-of-interests.

Analyzing these datasets using tools from statistical physics that are
completely novel in the context of human mobility (e.g. the pair distribution
function), and by considering locations as particles, our approach unveils
mobility patterns at scales that were previously concealed.

##  1 Mobility takes place on the physical structure of the built environment

The concept of the pair distribution function. To characterize geography we
study the pair distribution function (‘pair distribution’ below) between
locations. The pair distribution is a powerful tool from statistical mechanics
developed to understand the structural properties of materials [35]. Here, we
argue that the pair distribution is able to capture the structural properties
of locations in a given geography. To get a sense of this distribution,
consider a hypothetical spatial arrangement of concentric circles of $n$
locations (stop-locations, point of interest, addresses) on a disk of radius
$R$ (Fig. 1a). For an individual positioned at the center of these circles (a
white $\times$), the potential moves are constrained by the position of other
locations. For instance, in the absence of a location at distance $R/2$, a
movement of that exact distance is impossible. To understand the potential
movements from the disk’s center, we compute the distribution of distances of
every location to the center (Fig. 1b) which can be interpreted as what we
would observe if an individual was to move many times from the center to one
of the locations on the disk, choosing their new location independently of
distance. We can place each of the $n$ locations in the focus of such an
analysis, finding a unique distance distribution for each of them, and
calculate an average distance distribution. The result quantifies the expected
number of locations found at distance $r$ from an _average_ location and is
mathematically equivalent to the pair distribution, which simply encodes the
number of location pairs found at distance $r$ (Fig. 1c).

Figure 1: Illustration of the influence of geography on human mobility, (a)
spatial arrangement of concentric circles of locations, (b) distance
distribution of points to center of the disk (white cross). (c) Pair
distribution of the points in the concentric circles, (d) simulated moving
distance distribution for different cost of distance. Examples of more shapes
and distributions of locations are studied in SI. 6.1. (e,f,g) Three
geographies: (e) Real Denmark with 3.3M addresses, (f) Disk Denmark, where the
macro structure has cities distributed uniformly at random on a disk, but the
microstructure of the cities is preserved, and (g) Uniform Denmark with the
macro structure of Denmark (shape of landmasses and city positions) but no
local microstructure, i.e. uniformly distributed locations within cities. (h)
Pair distribution of the three geographies, and (i) the simulated observed
moving distance for an intrinsic distance cost of $1/r$. In (h) and (i), I,
II, III respectively represent the intra-city scale, the city scale, and the
inter-city scale.

A formal definition of the pair distribution function. Formalizing this
notion, we consider points $x_{i}\in\Omega$ distributed in a $d$-dimensional
space $\Omega\subseteq\mathbb{R}^{d}$ according to a density $\varrho^{d}(x)$,
then the pair distribution $p(r)$ is uni-dimensional and determined by both
the shape of the space as well as $\varrho^{d}(x)$. Thus, the pair
distribution is given by,

$\displaystyle p(r)=\int_{\Omega}dx\ \int_{\Omega}dy\
\varrho^{d}(x)\varrho^{d}(y)\delta(r-\lVert x-y\rVert).$ (1)

with $\lVert x-y\rVert$ the distance between two points $x$ and $y$. When we
count flows between pairs of points of $\Omega$, we observe a distance-
dependent distribution, which we refer to as $f(r)$ or the observed movement
distance. The observed distribution $f(r)$ can be represented as composite of
the pair distribution $p(r)$ and what we call the intrinsic distance cost
$\pi(r)$ (Methods: Pair distribution) [36, 37].

$\displaystyle f(r)\propto\pi(r)p(r).$ (2)

Observed data is a product of two separate components. This simple observation
has a profound consequence: Any study of the observed movement distances
$f(r)$ – the entire study of human mobility – is focused on a quantity that
includes the geometry of space $\Omega$ and the density $\varrho^{d}(x)$, both
reflected in the pair distribution $p(r)$. We argue that, any ‘universal’
behavior of the system, i.e. a geometry-independent law that encapsulates the
behavior of the system, must be captured entirely in $\pi(r)$, with the
remainder depending on the specific configuration of locations. This intrinsic
distance cost depends only on the distance between two locations. Thus, to
uncover any geometry-independent behavior of our system, we must adjust our
observation $f(r)$ by the geometry encoded in $p(r)$, i.e.

$\displaystyle\pi(r)\propto\frac{f(r)}{p(r)}.$ (3)

To illustrate the impact of the intrinsic distance cost, $\pi$, we simulate
movements on the concentric circles example (Fig. 1a,b,c) and demonstrate the
effect of three different potential distance cost functions:
$\pi(r)=\mathrm{const}$, $\pi(r)=\sqrt{r}$ and $\pi(r)=1/\sqrt{r}$. Note that
for $\pi(r)=\mathrm{const}$, i.e. when the choice of moving is independent of
distance, the distribution of moving distance is equal to the pair
distribution (Fig. 1d).

Illustrating scale via three Denmarks. To illustrate how geography shapes the
pair distribution and the observed moving distances, we consider the mobility
trace simulation for three different geographies (Fig.1e,f,g). The first
geography, ’Real Denmark,’ uses the actual geography of Denmark, including
3.3M precise locations (addresses) (Fig. 1e). The second, ’Disk Denmark,’
maintains the microstructure of cities but alters city positions and landmass
shapes, distributing real city centers uniformly on a disk (Fig. 1f). The
third, ’Uniform Denmark,’ keeps the macro structure of landmasses and city
positions from ’Real Denmark’ but features uniformly dense cities (Fig. 1g).

The distribution of pair distribution effectively encapsulates the nuances of
each geographical layout (Fig. 1h). On the local scales (less than
$1\mathrm{km}$), the pair distribution of the ‘Real Denmark‘ geography aligns
with the one of ‘Disk Denmark‘, reflecting city layouts. However, at broader
scales (more than $10\mathrm{km}$), the pair distribution of ’Real Denmark’
matches ’Uniform Denmark’, indicating the influence of the shape of land
masses and city positions.

To generate mobility traces we impose the same intrinsic distance cost on each
geography, and make the ansatz that it follows

$\pi(r)=1/r$ (4)

Remarkably, despite having the same intrinsic distance cost, the resulting
observed mobility traces are distinct across the three Denmarks (Fig.,1i).
This outcome is significant as it shows that even under a similar movement
law, the variations in geography alone can lead to distinct mobility patterns.

Next, we focus on the geographic regularities manifested through the power
laws and the fluctuations in the pair distribution of real Denmark (Fig.1 h).

##  2 Multi-scale urban structure shapes the distribution of pair distribution

The condensed matter physics of locations. Having illustrated the fundamental
way geography is encoded in a country’s location pair distribution and how it
affects the dyadic process of human mobility, we now focus on a deeper
analysis of its shape and origin. Specifically, we show how its properties
naturally emerge from statistical physics arguments, and how understanding the
pair distribution allows us to identify the key attributes of geography in
terms of shaping human mobility. Thus, in this section, we go beyond simply
appropriating the concept of the pair distribution function for analyses of
human mobility, but use the tools from condensed matter physics [35] to create
simple models for the micro-, meso-, and macro-structure of the geography of
an entire country.

Regimes of the pair distribution. The pair distribution between residential
locations (buildings and addresses) in Denmark shows the behavior displayed in
Fig. 2a. In this view buildings are simply addresses that stack on top of each
other (e.g. apartment buildings). On the micro-scale (I) within distances of
$r\lesssim 25\,\mathrm{m}$, the ‘neighborhood’, we observe a linear onset of
neighbor density, oscillatorily modulated. For larger distances at the meso-
scale (II), this growth assumes a scaling of approximately $p(r)\propto
r^{0.67}$ between $25\,\mathrm{m}$ and $200\,\mathrm{m}$ (until
$4000\,\mathrm{m}$ for addresses). Afterward, on the macro-scale (III), this
growth curbs, decays super-exponentially, and eventually approaches zero as we
reach the limit of the finite system.

Regime I: The density of cities as an ideal gas in a potential. Our basis for
modeling cities is a simple particle model in two dimensions – where buildings
are modeled as particles distributed uniformly at random, corresponding to an
ideal gas kept in place by an external potential. The key constraint to
identify the properties of this 2D confining potential (Fig. 2b) is the
observation from the literature that population density per unit area decays
exponentially with distance to its center (Fig. 2c) [38, 39, 40]. A
generalized Gamma distribution of shape
$p(r)=(r/R^{2})\exp[-(r/R)^{m}]/\Gamma(2/m)$ accurately describes the pair
distribution of this ensemble, where variable $R$ allows us to associate every
city with a radius that can be estimated from the data via a maximum-
likelihood fit, without having to rely on a definition of city center (Fig. 2c
and SI: City Sizes).

Regime I: Forces between locations. To explain the oscillatory behavior
modulating the onset of the pair distribution, we take the condensed matter
approach a step further. We argue that there is both a repulsive and
attractive effective force between locations. The repulsive force originates
from the simple fact that two buildings cannot occupy the same space. The
effective attractive force results from the lower cost of placing buildings
near one another as it reduces infrastructural cost for shared amenities [41].

In implementing this idea, we recognize that – unlike gases – cities are not
constructed according to strict laws, so there is not a single, ‘true’
description of the system. Thus, we illustrate the influence of such effective
forces on the pair distribution by modeling the system as two distinct
ensembles of repulsive-attractive forces: (i) an ensemble of Lennard-Jones
(LJ) disks, a canonical model for sphere-like particles [42], and (ii)
attractive hard disks of heterogeneous size, both in the presence of a linear
external potential (Methods: Locations as interacting particles). Both models
reproduce the oscillatory features observed at the micro scale in our high-
resolution dataset (Fig. 2d).

We can study the oscillatory structure in terms of the pair-correlation
function $g(r)=p(r)/i(r)$, where $i(r)$ is the pair distribution of an ideal
gas. This function measures the over- and under-representation, respectively,
of neighbors at distance $r$ in relation to the expectation if no interaction
forces were present [43]. In both models and the data, there is a considerable
lack of neighboring buildings for distances $r\lesssim 4\,\mathrm{m}$,
suggesting presence of repulsive forces, later reaching a peak $g(r)>1$ at
$r\approx 6\,\mathrm{m}$ which is typically observed for systems with
attractive forces. This peak becomes wider when we consider heterogeneous
building radii. The discrepancy between the data and model(s) below
5$\,\mathrm{m}$ (distance between black and teal/pink lines in Fig. 2d) is due
to the definition of building location in the data as the location of its
front door, rather than the building center (Methods: Data).

Figure 2: (a) Pair distribution of residential locations in Denmark. (I)
Neighborhood-scale with a linear onset and an oscillating modulation for the
immediate neighborhood, (II) city-scale showing a power-law growth with
exponent $\alpha=0.67$, and (III) country-scale with a slower growth and
eventual fast decay. (b) With external potential $V(|\bm{x}|)$ that grows
linearly with distance to the city center $|\bm{x}|$, an ideal gas of
buildings will form a city-like structure with radial building density
following an Erlang distribution of scale $R_{0}$ (squares) (Methods: Model of
locations). The corresponding pair distribution approximately assumes a
generalized Gamma distribution with scale $R$ providing a definition for the
radius of a city (circles). (c) Maximum-likelihood fits of the generalized
pair distribution against pair distribution of buildings in three Danish towns
(upper) and all Danish towns, rescaled (lower). (d) Pair distributions for
four systems of building positions on the micro-scale, (i) buildings in
Denmark, (ii) random positions (ideal gas), (iii) Lennard-Jones disks, and
(iv) attractive hard disks of heterogeneous size. All systems display linear
onset with oscillatory modulations. The pair-correlation function $g(r)$
quantifies the under- and over-presence, respectively, of neighbors at
distance $r$ when comparing to the ideal gas. Repulsive and attractive forces
lead to valleys and peaks in $g(r)$ (see Methods: Model of locations, Hard
disks and Molecular Dynamics simulation) (e) Pair distribution for all towns
in Denmark for (i) random positions, (ii) non-attractive hard disks, (iii)
real data. Systems (ii) and (iii) display similar behavior with slow onset,
quickly approaching the ideal gas. Weighting with the numbers of building
pairs per pair of towns (pink and green lines) replicates the macro regime
behavior (III) of panel (a) but introduces slight deviations to the model
system (ii). (f) (i) Fractal model for building positions in cities, as a
hierarchically nested structure of self-similar patches. On each hierarchy
layer, a patch of buildings consists of $b$ sub-patches distributed non-
overlapping with constant density and packing fraction $\theta$. At the lowest
level, a sub-patch is a building. (ii) A sample configuration of this model
with $b=4$, $z_{1}=3\,\mathrm{m}$, $\theta=0.8$, and $L=8$. (iii-iv)
Independent-patch model. Equating the concept of “town” with “patch of radius
$R$ containing buildings”, the number of buildings per patch grows as $R^{2}$
and the number of pairs of buildings per patch grows as $R^{4}$ (with $R$ from
the analysis in panel c). At the same time, the tail of the patch radius
distribution falls as $R^{-(4-\alpha)}$ with $\alpha=0.67$. (v) For both the
fractal model and the independent-patch model, we recover the pair
distribution scaling observed in the data.

Regime III: The spatial distribution of cities. We now use the condensed
matter tool-set to understand the macro-scale (III) of the pair distribution.
For large distances, its shape is dominated by the locations of cities. We
compare an ideal gas of cities (i.e. random positions uniformly distributed in
the landmass of Denmark) to the empirical pair distribution of all Danish
towns (Fig. 2e). While the ideal-gas distribution has a linear onset, there is
a considerable lack of neighboring towns for small distances $r\lesssim
1\,\mathrm{km}$, but after that, the observed town pair distribution rapidly
approaches and matches the ideal-gas distribution remarkably well,
demonstrated by the fast approach of $g(r)\rightarrow 1$.

The fact that there is not a clear peak in $g(r)$ suggests the absence of a
strong attractive force. When we compare the observed pair distribution to
that of a non-attractive hard disks-model (Methods: City pair distribution) we
observe a similar behavior, indicating that the positions of towns are almost
statistically indistinguishable from random locations except for the condition
that towns may not overlap in territory.

To obtain the shape of the pair distribution for buildings, next we weigh
every contributing pair distance for a pair of towns by the product of their
respective numbers of buildings (Fig. 2e) This shifts the ideal gas pair
distribution to larger distances and introduces a peak in $g(r)$ at around
$r=10\,\mathrm{km}$. While the shifted onset is replicated in a building pair-
weighted hard disk model, we do not recover a peak. This behavior, however,
can be explained by means of central place theory [44], which states that
towns that are located near large cities tend to be substantially smaller,
hence introducing an effective attractive force between types of cities
(Methods: Non-overlapping disks)

Regime II: Fractal structures explain the meso-scale. As described above, the
initial linear growth of the pair distribution quickly approaches a sub-linear
scaling law of $\alpha\sim 0.67$ (Fig. 2a). We attribute the emergence of this
scaling law to the non-regular shapes of cities composed of smaller ‘patches’
of buildings that form larger clusters in a space where occupation is limited
by local geography (e.g. bodies of water, hills) or for developmental reasons
(e.g. industrial areas, parks, agricultural use).

Specifically, we adapt a model that was recently used to generate surrogate
city positions [32, 45]. Assume that a city is constructed as a hierarchically
nested structure of patches containing $b$ buildings with packing fraction
$\theta$ (Fig. 2f.i, Methods: Fractal model). For a range of realistic
parameters, this model yields $\alpha=0.67$ for both single realizations (Fig.
2f.ii) and analytically (SI. 7.7). Our model replicates the behavior observed
for the whole country and in agreement with $\alpha=0.69\pm 0.04$ for
individual cities. Furthermore, this nested description is consistent with
earlier observations that settlement size distributions are scale-free [23,
46].

By modeling the built environment as a multiverse of independent patches, we
can infer the patches’ area distribution from the distributions of radii using
the per-city pair distributions (Fig. 2f.iii and Methods: Patch Model), we
obtain $\alpha=0.71$ as well as a patch area distribution of $f(A)\propto
1/A^{2.15}$, a result that matches the exponents reported in [46]. The meso-
scale scaling law is robust as all these similar values of $\alpha$ arising
from different models are consistent with our empirical finding (Fig.3a-Fig.5)
and the literature [46].

##  3 Intrinsic Distance Cost function

A power law that spans five orders of magnitude. Having established a
fundamental understanding of the physics of the pair distribution, we now
demonstrate the impact of taking into account the information encoded within
it, by analyzing an exhaustive dataset—39M residential moves between 3.3M
Danish addresses over 36 years. The dataset presents reduced bias as it covers
all residential moves (Methods: Data). To uncover the intrinsic distance cost
of residential mobility, we first eliminate the geographical component by
normalizing the observed moving distance distribution by the pair distribution
(Fig. 3a). This normalization reveals that the intrinsic distance cost
function (Eq. (3)) follows a power law remarkably well:

$\displaystyle\pi(r)=\frac{f(r)}{p(r)}=\frac{1}{r^{s}}.$ (5)

This power law is consistent across scales ranging from 10m to 500km, covering
five orders of magnitude (Fig. 3b). A maximum likelihood fit provides an
exponent of $s=0.98$, which validates the ansatz of Eq. (4); the intrinsic
distance cost follows $\pi(r)=1/r$. We consider this power law to be an
extension of the gravity model to the continuous domain.

Figure 3:  (a) Empirical moving distance distribution and pair distribution in
Denmark. Dividing the mobility by the pair distribution yields an intrinsic
distance cost (b) that follows a power law from 10m to 500km (maximum
likelihood estimation). (c) This framework can be represented as a continuous
gravity law, if we take a toy model of two distant cities, (d) the pair
distribution at the distance $d$ between the two cities collapses into a Dirac
distribution, with the height equal to the product of cities population,
$N_{1}N_{2}$. (e) The toy intrinsic distance cost of the simulation model. For
a comparison with a continuous radiation model, see SI. 8.1.2. (f-g)
Residential mobility analysis in a small (f) and a large (g) Danish town
reveals a consistent piece-wise power law across all cities. For more cities,
see SI Fig.28-31, for details on the piece-wise power-law definition and
maximum likelihood estimation see SI. 8.4. (h) Distribution of the piece-wise
power laws’ exponents, for the first and second parts, respectively center
around $0.6\pm 0.20$ and $2.0\pm 0.21$ (mean $\pm$ standard deviation). (i)
The transition point in the piece-wise power law aligns with a city size
distribution. We fit different distributions (log-normal, power-Pareto) via
maximum-likelihood (see Methods: Power-law estimation). (j) We recover the
intrinsic distance cost from simulating the piece-wise power law (f,g) with
the mean values of $0.6$ and $2$ (h), over the city-scale (i). We assume
random position of cities as in Fig. 2e. For alternative derivation of the
intrinsic distance cost and details see SI. 8.5.

A ‘geography-free’ gravity law. In the gravity law for human migration, the
migration probability between cities is proportional to the product of their
population and inversely proportional to their distance. This model is
inherently discrete as it requires administrative boundaries to define the
population size.

However, the concept of pair distribution closely aligns with the idea behind
the gravity law. For example, consider two cities with different populations
and radii (Fig. 3c). When the distance between the cities far exceeds their
radii, the pair distribution tends to a Dirac distribution of height equal to
the product of the number of addresses in each city (Fig. 3d-e). Assuming that
the number of addresses is proportional to the population, the ’mass-product’
term of the gravity model emerges and Eq. (2) coincides with the gravity model
(Methods: A ‘geography-free’ gravity law). Our framework is continuous as it
operates at the granularity of individual addresses, the finest possible
scale.

Piece-wise local power laws. Considering the intrinsic distance cost starting
from individual Danish cities (as opposed to globally), we observe a piece-
wise process (SI Fig. 28-32) [47]. Specifically, for each city, we consider
the distribution of moves originating in the city of interest to the rest of
the country. The pair distribution is limited to pairs that include at least
one address in the city of interest. At the city scale, we find that the
intrinsic distance cost is a piece-wise power law (Fig. 3f-g) with the first
exponent centered around $0.60\pm 0.20$(SD) and the second exponent around
$2.0\pm 0.21$(SD) (Fig. 3h). The process is universal in the sense that the
exponents of the piece-wise power laws are consistent across the 1400 cities
of Denmark. The transition point between the two power laws defines a mobility
city radius (Fig. 3i) with a distribution that lies between a log-normal and a
power-Pareto distribution (Methods: Power-law estimation).

Putting the pieces together. So how does the global power-law (Fig. 3a) relate
to the local piece-wise picture (Fig. 3f-g) starting in each city? Here, our
statistical physics based understanding of the components of the pair
distribution can be used to complete the picture of human mobility, explaining
the intrinsic distance cost starting from the city-level description. We
simulate the piece-wise power law over a toy geography that reproduces the key
characteristics of geography: scale of cities (Fig. 3.i), local pair
distribution that follows a generalized gamma distribution (Eq. (47), Fig.
2c), and random position of cites (Fig. 2e.III), we recover the empirical
intrinsic distance cost, $\pi(r)=1/r$, validating our geographical model and
the piece-wise power law.

Generalizing to other geographies and types of mobility [48]. Finally. we show
that our results are not particular to either the geography of Denmark or
residential mobility. Figure 4a illustrates that the normalization by the pair
distribution unveils a power law for residential mobility in France. We also
find the same patterns for to day-to-day mobility across the diverse
geographies of Houston, Singapore, and San Francisco (Fig. 4b-d). The fact
that our results generalize highlights the general nature of the power law
distribution of intrinsic mobility cost.

Figure 4: (a) Observed moving distance distribution and pair distribution for
residential mobility in France. (b) Observed moving distance distribution and
pair distribution for day-to-day mobility between points of interest in San
Francisco, Houston (c), and Singapore (d).

##  4 Discussion

Geography trivially constrains human mobility. While the existing human
mobility literature has extensively documented the structure of human
mobility, it often remains disconnected from studies focusing on the
organization of the built environment. This led to the emergence of two
distinct viewpoints on human mobility [1]: one based on distance, similar to
the gravity model [5], and another focused on the availability of
opportunities [27]. Here, we have proposed the pair distribution of locations
as a key normalization variable of human mobility. The normalization reveals
an intrinsic power law consistent across various geographies and mobility
types. Additionally, we proposed a model inspired by statistical physics to
account for the shape of the pair distribution between locations. Our model
replicates the characteristics of the pair distribution at the scales of
dwelling and country. While our analysis primarily focuses on addresses,
further work could study the pair distribution between point-of-interests or
between pairs of residential and workplace locations. Such an extension could
help us understand commuting patterns, a significant component of day-to-day
mobility [48, 49, 50]. Finally, we explained our normalization as a continuous
gravity model that reconciles the two paradigms of human mobility. We also
show that this description does not hold at the city scale and emerges from
the aggregation of scales [7]. At the city level, other demographic factors
could explain relocation decisions, such as housing prices, employment
opportunities, and family evolution [47]. Our framework enhances the
understanding of human mobility and provides a bridge connecting two rich
bodies of literature.

## Methods

###  4.1 Data

####  4.1.1 Uncertainty Quantification

Estimates reported in the main text are reported as mean $\pm$ standard
deviation.

####  4.1.2 Residential Mobility in Denmark

The data on residential mobility within Denmark is derived from the Befolkning
database of Danmarks Statistik [51]. This database contains 39,297,646
residential moves between different addresses in Denmark from 1986 to 2020.
The location data for the 2,345,453 buildings that make up the 3,251,464
addresses are also provided by Danmarks Statistik, with an accuracy better
than $\pm 2$ meters [52]. The location of an address is defined as the
location of the building’s entrance door, we investigate the effect of this
definition of pair distribution function in SI. 7.4.

####  4.1.3 Bias in residential mobility in Denmark

In contrast to studies using digital tools such as mobile phones, which may
introduce biases [53], the Danish residential mobility dataset we study is
inherently representative of the entire population. In fact, the dataset
includes every residential move, a requirement for all residents to report,
thus eliminating potential biases associated with selective data collection
methods [54]. Demographic information (gender and year of birth) are reported
in Fig. 6. Cities are defined according to the definition of Danmarks
Statistik; to ensure the robustness of the results, we compare their
definition with a hierarchical clustering of the addresses (HDBSCAN [55] in
SI. 9)

####  4.1.4 Residential mobility in France

The location data for France is compiled from the Base Adresse Nationale,
accessible via https://adresse.data.gouv.fr/donnees-nationales, which contains
precise location data for 21,567,447 buildings. To match the housing unit
data, we utilized information from https://www.data.gouv.fr (FiLoSoFi), which
segments the country into a grid of 200m squares, each annotated with the
count of housing units and collective housing entities. We augment the data by
uniformly distributing the housing unit count amongst each collective housing
entity, resulting in 34,041,910 individual addresses. The residential mobility
data was obtained from https://www.insee.fr/fr/statistiques. It contains
40,465,288 inter-city migrations from 2016 to 2020. The distances are
calculated from one city center to another. The city centers are computed as
the centroid of contained addresses. As we only have access to inter-city
moves, migration distances that are less than the largest city’s dimensions in
France were disregarded.

####  4.1.5 Day-to-Day Mobility

The data on daily mobility was obtained from user check-ins on the location-
centric social network Foursquare over six-months (May 27, 2010 to November 3,
2010) as documented [28]. The data comprises a total of 239,788 movements
across 43,395 distinct locations. The regional breakdown of the data is,
47,996 movements occurred across 11,808 locations in Houston, 79,624 movements
across 15,617 locations in Singapore, and 112,168 movements across 15,970
locations in San Francisco.

####  4.1.6 GeoBoundaries

We obtain the country, state, and municipality shapes from the geoboundaries
project [56]. Throughout the paper, we work with the EPSG:23032 projection for
all points that lie within Denmark and the EPSG:27561 projection for points
that lie within France. We then perform analyses on the Euclidean geometry of
the projection.

###  4.2 Pair distribution function

This paragraph provides further details on how the pair distribution function
emerges from studying the dyadic process of human mobility between locations.
The expected number of movements from location $x$ to location $y$ is given
by,

$\displaystyle f(x,y)\mathrm{d}^{2}x\mathrm{d}^{2}y$ (6)

with $f(x,y)$ being a continuous number density. This number is influenced by
two functions. First, the number of observations is proportional to the number
of locations at $x$ as well as the number of locations at $y$, i.e., depending
on the product of location density
$\varrho(x)\varrho(y)\mathrm{d}x\mathrm{d}y$. Second, the number of movements
is proportional to the movement propensity, i.e., the tendency for a move that
begins at a location at $x$ to end at a location at $y$, which we denote as
$\pi(x,y)$,

$\displaystyle
f(x,y)\mathrm{d}^{2}x\mathrm{d}^{2}y=\pi(x,y)\varrho^{d}(x)\varrho^{d}(y)\mathrm{d}^{2}x\mathrm{d}^{2}y.$
(7)

The number of moves spanning distance $r$ within the domain $\Omega$ is

$\displaystyle
f(r)=\int_{\Omega}\pi(x,y)\varrho^{d}(x)\varrho^{d}(y)\delta(r-\lVert
x-y\rVert)dxdy.$ (8)

If we make the anstaz that the propensity to travel from $x$ to $y$ solely
depends on the distance between the two, i.e,

$\displaystyle\pi(x,y)=\pi(\lVert x-y\rVert))=\pi(r)$ (9)

as $\pi(r)$ is now independent of location $x$ and $y$, it can be factored out
of the integral, we obtain

$\displaystyle f(r)\propto\pi(r)p(r)$ (10)

with $p(r)$ the pair distribution function as in Eq. (1), and we recover the
main text’s Eq. (2) (see SI. 6.2 for more details). The ansatz of Eq. (9) is
the key point to reconcile the distance-based and opportunity-based
perspective on human mobility. The distance term is $\pi(\lVert x-y\rVert)$,
the opportunity term is the expected number of pairs of locations between
location $x$ and location $y$,
$\varrho^{d}(x)\varrho^{d}(y)\mathrm{d}^{2}x\mathrm{d}^{2}y$ (SI. 6.2, SI.
8.3) [57, 58].

Moreover, the pair distribution function characterizes the geometry of the
space as it represents the joint probability of finding two locations at
particular positions in the system, (SI. 6.1 and SI. 6.3). The computation of
the pair distribution function scales as $\mathcal{O}(n^{2})$, we use a k-d
tree to compute it for large number of locations (see SI. 7.5)

###  4.3 Locations as interacting particles

Our basis for modeling cities is a simple particle model in two
dimensions—where buildings are modeled as particles uniformly distributed at
random, corresponding to an ideal gas kept in place by an external potential.
Consider the city center, e.g., the central business district, as the primary
point of interest due to its proximity to amenities. This attractiveness
implies a higher cost of being located further away from the center,
represented by the distance $|\bm{x}|$.

At the same time, two simple mechanisms will prevent buildings from
accumulating at the exact center of the city. First, buildings have a certain
average radius $z$, so they cannot reside too close (at distance $<2z$). There
is an advantage to buildings not being too far apart, as they can share local
amenities. Rather than explicitly modeling this phenomenon, we assume that
there is an inherent temperature, $T$, in the system. This temperature
determines the distribution of house locations, which is influenced by
interactions and external potential.

From a statistical physics perspective, we describe the system with a simple
external potential that linearly increases with distance.

$\displaystyle V^{\mathrm{ext}}(\bm{x})=\gamma|\bm{x}|$ (11)

where we assume that the origin of the coordinate system is located at the
city center. In order to model repulsion and attraction between houses, we
further assume a Lennard-Jones interaction potential

$\displaystyle
V^{\mathrm{int}}(\bm{x}_{i},\bm{x}_{j})=\varepsilon\left[\left(\frac{2z}{|\bm{x}_{j}-\bm{x}_{i}|}\right)^{12}-2\left(\frac{2z}{|\bm{x}_{j}-\bm{x}_{i}|}\right)^{6}\right],$
(12)

which is commonly used to model intermolecular attraction and repulsion [42,
35].

In total, a system with these properties evolves according to the Hamiltonian

$\displaystyle\mathcal{H}(\\{\bm{x},\bm{p}\\})=\frac{1}{2}\sum_{i=1}^{N}\bm{p}_{i}^{2}+\gamma\sum_{i=0}^{N}|\bm{x}_{i}|+\frac{1}{2}\sum_{i=1}^{N}\sum_{j\neq
i}^{N}V^{\mathrm{int}}(\bm{x}_{i},\bm{x}_{j})$ (13)

with two-dimensional momenta $\bm{p}_{i}$ and locations $\bm{x}_{i}$. In a
canonical-ensemble formulation of the system, i.e. at constant inverse
temperature $\beta$, the probability of finding a configuration
$\\{\bm{x},\bm{p}\\}$ in volume-element $d^{2}\bm{x}d^{2}\bm{p}$, is given by

$\displaystyle\varrho[\\{\bm{x},\bm{p}\\}]d^{2}\bm{x}d^{2}\bm{p}$
$\displaystyle=\exp[-\beta\mathcal{H}(\\{\bm{x},\bm{p}\\})]d^{2}\bm{x}d^{2}\bm{p}.$
(14)

We will initially consider an ideal gas first where $V^{\mathrm{int}}=0$,
which will enable us to say something about the density of particles around
the center of city, i.e. we want to find the particle dwelling probability
$p(r)dr$ with $r$ being the particle’s distance to the center. For the sake of
simplicity, we set $N=1$ without loss of generality. This is justified by the
ergodicity of the system, which implies that the trajectory of a single
particle will eventually follow the density of the entire distribution. This
can be expressed as taking the $N$th root of the $N$-particle density.
Integrating over the momenta yields

$\displaystyle\varrho[\bm{x}]d^{2}\bm{x}$
$\displaystyle=\exp\left(-\beta\gamma|\bm{x}|\right)d^{2}\bm{x}$ (15)

With a change of variables to polar coordinates, we find

$\displaystyle\varrho[r,\phi]rdrd\phi=r\exp\left[-\beta\gamma
r\right]drd\phi,$ (16)

i.e.

$\displaystyle p(r)dr=r\exp\left[-\beta\gamma r\right]dr.$ (17)

This implies that in the absence of interactions, the distribution of houses
around the city center should follow a $k=2$ Erlang distribution with a scale
parameter of $\lambda=\beta\gamma=1/R_{0}$ where $R_{0}$ is half the city
radius. The role of the inverse temperature term, $\beta$, is explained in SI.
7.1.4.

To gain insight into the radial particle density in the context of strong
repulsion, we return to the Lennard-Jones perspective. In the limit of
$\varepsilon/T\gg 1$, particles will have a strong tendency to be found in
their respective potential minimum, i.e. at distance $2z$ from each other.
Effectively, we can think of them as hard disks of radius $z$ that have a
tendency to form clusters. If the temperature is low, the effective radius of
the city will be small, we expect a crystal to form in the center (see SI.
7.1.4).

####  4.3.1 Molecular Dynamics simulation

We are interested in finding configurations $\mathcal{C}$ that accurately
depict the canonical ensemble with number of particles $N$, an average-
constant temperature $T$, and a constant but irrelevant volume. We assume that
the particles are constrained to a radially symmetric external potential, with
the total volume of the system containing the particles being irrelevant.
Finally, the total potential must be

$\displaystyle
V=\sum_{i=1}^{N}\left[V^{\mathrm{ext}}(\bm{x}_{i})+\frac{1}{2}\sum_{j\neq
i}V^{\mathrm{int}}(\bm{x}_{i},\bm{x}_{j})\right].$ (18)

To this end, we integrate the equations determined by the system’s Hamiltonian
numerically using the velocity-Verlet algorithm [59]. Furthermore, we rescale
particle velocities according to the stochastic Berendsen thermostat with
relaxation time $\tau$ [60]. The details and parameters of the molecular
dynamics simulation are available in SI. 7.1.5.

####  4.3.2 Hard disks

We generate a configuration of hard disks with a non-specified attractive
force (i.e. we do not explicitly integrate the equations of motion for a hard-
sphere interaction potential with an additional attractive force). To do so,
we first generate $N=10^{4}$ random positions according to a radial Erlang
distribution of scale parameter $R=675$ (i.e. an ideal gas). Subsequently, we
draw a random radius for each position from a heterogeneous distribution with
power law tail. We first draw values $\hat{z}_{i}$,

$\displaystyle
p(\hat{z})=\frac{a^{2}-1}{2a}\times\begin{cases}\hat{z}^{a},&\hat{z}\leq 1,\\\
\hat{1}/z^{a},&\hat{z}>1,\end{cases}$ (19)

and then assign disk radius $z_{i}=3\hat{z}_{i}/\left<\hat{z}\right>$. We
choose $a=4$ to obtain a heterogeneous distribution with non-finite variance
in disk area. Afterwards, we run the collision algorithm outlined in SI. 7.2.

The modelling of interaction between the house in the external potential is
further developed in the SI. 7.3.

###  4.4 Non-overlapping disks

To emulate the position of cities, we randomly distribute $N$ hard disks of
radii $R_{i}$ in shape $\Omega$. We start with largest disk of radius
$R_{\mathrm{max}}$ and iterate over all disks, ordered decreasingly in size.
For every disk $j$, we generate a random position $\bm{x}\in\Omega$ until the
condition $|\bm{x}-\bm{x}_{i}|>R_{j}+R_{i},\,\forall j<i$ is met, i.e. drawing
new random positions until there are no overlaps with other, already placed
disks. Then, assign $\bm{x}_{j}\leftarrow\bm{x}$ and continue with the next
disk.

###  4.5 City pair distribution and radius

The pair distribution of each city is best represented by a generalized Gamma
distribution with linear onset

$\displaystyle
p(r,R,m)=\frac{r}{R^{2}\,\Gamma(2/m)}\exp\left[-\left(\frac{r}{R}\right)^{m}\right],$
(20)

as depicted in Fig. 2c. Here, $r$ is the distance between two locations, $R$
is the scale parameter corresponding to city-radius and $m>0$ is a shape
parameter controlling the decay of the tail. Assuming that population density
decays exponentially with distance to the city center and that consequently
the pair distribution of buildings assumes a generalized pair distribution of
the form of Eq. (20), we perform maximum-likelihood fits to infer parameters
$R$ and $m$ for each city in Denmark with more than 30 buildings, we generally
find strong correspondence to the empirical distribution, as illustrated by
the correlation of $r/R$ and the inverse of Eq. (20),

$\displaystyle(r/R)^{*}=-(\log[p(r)R^{2}\Gamma(2/m)/r])^{1/m},$ (21)

albeit with small deviations at smaller distances (cf. Fig. 2c) that can be
associated with the sub-linear scaling behavior observed on the meso-scale of
Fig. 2a. Remarkably, measuring the ‘radius’ $R$ of a city in terms of its pair
distribution’s scaling parameter comes with the advantage of not having to
rely on any definition of ‘city center’.

Moreover, no assumptions about urban growth are necessary; instead, the
observed behavior emerges from the simple ansatz that being placed at distance
$|\bm{x}|$ from the center of a town comes at a cost $\propto|\bm{x}|$, which
can be overcome for a multitude of reasons encoded in temperature $T$.

####  4.5.1 Fractal Model

The fractal model we adopt was originally devised to explain the positions of
galaxies and was recently used to generate surrogate city positions [32, 45].
We assume that a city is constructed as a hierarchically nested structure of
patches containing buildings (cf. Fig. 2f.i and SI. 7.7). At layer $k=1$, a
patch of radius $R_{1}$ contains $b$ buildings of radius $z_{1}$ placed in a
non-overlapping manner with packing fraction $\theta$. At layer $k=2$, a patch
of radius $R_{2}$ contains $b$ building-patches of radius $z_{2}=R_{1}$, again
placed without overlap and packing fraction
$\theta=\theta_{k}=bz_{k}^{2}/R_{k}^{2}$. Continue in that fashion until
reaching layer $k=L$. Due to the constant packing fraction $\theta<1$, larger
areas of absent buildings will form at each hierarchy layer, mirroring similar
observations in real urban structures. Note that demanding $\theta$ to be
constant across hierarchy layers also implies a predetermined patch radius of
$R_{k}=(\theta/b)^{k/2}z_{1}$. A sample configuration with $b=4$,
$z_{1}=4\,\mathrm{m}$ and $\theta=0.8$ can be seen in Fig. 2f.ii.

To evaluate the overall pair distribution of such a fractal structure, we
assume that for a patch of hierarchy layer $k$, the $b^{k-1}$ buildings that
are contained in each of its $b$ sub-patches are all located at their
respective centers. Then, the pair distance number density (pdnd) of an
average such patch will approximately assume Eq. (20), weighted with a total
of $b^{2k-2}b(b-1)/2$ building pairs (between subpatches). Moreover, at every
hierarchy layer $k$ there will be $b^{L-k}$ such patches. In total, each
hierarchy layer therefore contributes a pdnd of

$n_{k}(r)\propto b^{k}p(r,R_{k},m),$ (22)

implying that the observed pdnd of the whole structure is proportional to

$n(r)=\sum_{k=1}^{L}n_{k}(r)\propto\sum_{k=1}^{L}\frac{b^{2k}r}{\theta^{k}z_{1}^{2}}\exp\left[-\left(\frac{r}{z_{1}}\right)^{m}\left(\frac{b}{\theta}\right)^{mk/2}\right].$
(23)

In the limit of $L\rightarrow\infty$, we use Laplace’s method to find

$\log
n(r)\propto\underbrace{\left({1-\frac{2\log\theta}{\log(\theta/b)}}\right)}_{=\alpha}\log
r,$ (24)

i.e. a sub-linear scaling law emerges from considering the pair-weighted sum
of single-scale pair distributions with linear onset. The pair distance number
density $n(r)\propto r^{\alpha}$ with $\alpha=1-2\log\theta/\log(\theta/b)$.
Notably, the result is independent of the parameter $m$ that controls the tail
of the single-scale pair distribution, suggesting that the exact shape of the
composite pair distributions are irrelevant, which we further demonstrate in
the SI. 7.1 by discussing another functional form of the pair distribution.

####  4.5.2 Patch Model

We consider a _multiverse_ consisting of an infinite amount of universes
indexed $i$, each of which is inhabited by a patch of size $R_{i}$ with
$N_{i}$ buildings inside, contributing to the multiverse with a pair
distribution number density of $n(r,R_{i})\approx N_{i}^{2}p(r,R_{i},m)$ (each
universe contributes a number of building pairs amounting to $N_{i}^{2}$). We
also assume that there’s a constant population density $\rho_{0}\propto
N/R^{2}$, demanding that $N\propto R^{2}$.

For each of these universes (or rather, each of these patches), we calculate
the number of building pairs at distance $r$ and assume that patch sizes $R$
are distributed according to some—at this time still unknown—distribution
$f(R)$. Then, the joint distribution of house pairs at distance $r$ is given
by

$\displaystyle n(r)$ $\displaystyle\propto\int\limits_{0}^{\infty}dR\
n(r,R)f(R)$ (25) $\displaystyle\propto\int\limits_{0}^{\infty}dR\
f(R)R^{4}\frac{r}{R^{2}}\exp\left(-\left(\frac{r}{R}\right)^{m}\right).$ (26)

To simplify the integral, we change variables to $\beta=1/R$, such that
$\mathrm{d}R=-\beta^{2}\mathrm{d}\beta$ and

$\displaystyle n(r)$ $\displaystyle\propto\int\limits_{0}^{\infty}d\beta\
f(\beta^{-1})\frac{r}{\beta^{4}}\exp\left(-(\beta r)^{m}\right).$ (27)

This integral yields a solution if $f=\beta^{4-\alpha}$ with $\alpha<1$. Let’s
assume that this is the case, which means we have to solve the integral

$\displaystyle
n(r)=C\int_{0}^{\infty}\beta^{-\alpha}r\exp\left(-(r\beta)^{m}\right)d\beta$
(28)

where $C$ is a normalization constant. We recognize the integral of the gamma
function (see SI. 7.6) and therefore find,

$\displaystyle n(r)\propto
r^{\alpha}C\,\Gamma\left(\frac{1-\alpha}{m}\right).$ (29)

That is, if the radius of patches in the multiverse would be distributed as
$f(R)\propto R^{-(4-\alpha)}$, the initial growth of the joint pair
distribution would follow a sub-linear power law. Notably, the scaling
exponent does not depend on the tail parameter $m$.

Let’s see what this would mean in terms of the area of the patches. The area
would scale as $A\propto R^{2}$ (or $R\propto\sqrt{A}$), i.e. the distribution
of the area would follow

$\displaystyle\tilde{f}(A)$
$\displaystyle=\frac{\mathrm{d}R}{\mathrm{d}A}f(R=\sqrt{A})$ (30)
$\displaystyle\propto\frac{1}{A^{(5-\alpha)/2}}.$ (31)

In the data, we observe $\alpha=0.67$. That would mean that the area of the
patches would have to be distributed according to a power law with exponent
$\mu=-(5-\alpha)/2\approx-2.17$, which is well within what has been found
empirically for cities [46]. Furthermore, we find $f(R)\propto 1/R^{3.29}$
from our city radius inference analysis (cf. Fig. 21) which leads to
$\mu=-2.15$, showing that the results of these two separate analyses are
consistent.

###  4.6 A ‘geography-free’ gravity law

In its original form, the gravity law for human migration states that the
probability of moving between two cities is proportional to the product of
their population and inversely proportional to the distance separating them
[5],

$\displaystyle P_{1\to 2}\propto\frac{N_{1}N_{2}}{d_{1\to 2}}.$ (32)

This relationship requires administrative units to define cities and their
respective populations, thereby rendering it inherently discrete and subject
to arbitrary choices regarding administrative boundaries. Yet, constructing a
product of origin$\times$destination populations closely aligns with the idea
of normalizing by the pair distribution. To illustrate, consider two cities of
populations $N_{1},\,N_{2}$, with radii $r_{1}$ and $r_{2}$ separated by a
distance $d$. If we assume a uniform location density within these cities (as
illustrated in Fig.3c), the density of locations is given by

$\displaystyle p^{d}(x)=\frac{N_{1}}{\pi
r_{1}^{2}}\delta_{\Omega_{1}}(x)+\frac{N_{2}}{\pi
r_{2}^{2}}\delta_{\Omega_{2}}(x),$ (33)

with $\Omega_{j}$ the domain of city $j$. This leads to the pair distribution
(using Eq. 1) :

$\displaystyle p(r)=N_{1}N_{2}\int_{\Omega}\int_{\Omega}dxdy\frac{1}{\pi
r_{1}^{2}}\frac{1}{\pi
r_{2}^{2}}\delta_{\Omega_{i}}(x)\delta_{\Omega_{j}}(y)\delta(r-R(x,y)).$ (34)

In the limit case where $d$ significantly exceeds both $r_{1}$ and $r_{2}$,
when the city radii are negligible compared to the distance between cities,
the pair distribution around $d$ tends towards a Dirac distribution of height
equal to the product of the number of addresses in each city,
$p(r)\xrightarrow[d>>r_{1},r_{2}]{}N_{1}N_{2}\delta_{d}(r)$, therefore
(Fig.3d,e),

$\displaystyle
f(r)=p(r)\pi(r)\xrightarrow{}\frac{N_{1}N_{2}\delta_{d}(r)}{r}.$ (35)

Equation 35 illustrates how in the context of the the gravity model for
migrations, the ‘mass-product’ $N_{i}N_{j}$ emerges. Given that the number of
addresses in a city is proportional to its population, we recover the discrete
gravity law. We refer to our finding as ‘continuous’ since it operates at the
granularity of individual addresses, the finest possible scale (see SI.
8.1.1). A comparison between the gravity model, the radiation model, and the
pair distribution framework is available on figure 7.

We can also interpret the intrinsic distance cost $pi(r)=1/r$ as an emergent
property of random utility theory, where individuals appreciate distance with
a logarithmic scale. This interpretation is consistent with the logarithmic
utility function described in SI. 8.2.

###  4.7 Power-Law Estimation

The parameters for power-law distributions are determined using maximum-
likelihood estimation, as outlined in [61]. For continuous distributions, the
probability density function (pdf) for a power-law distribution is given by

$\displaystyle p(x)=Cx^{-\alpha}$ (36)

i where $x>x_{\text{min}}>0$ and $C$ is a normalization constant. The
parameter $\alpha$ is typically greater than 1, as the probability density
function described in equation 36 does not integrate over
$x\in[x_{min},+\infty)$ for $\alpha\leq 1$. However, in our analysis, the
intrinsic distance cost has an exponent $\alpha\sim 1$. Moreover, the maximum-
likelihood estimation method from [61] is known to be biased when
$\alpha<1.5$, a limitation that has been confirmed by [62]. Consequently, for
estimating power laws with $\alpha<1.5$, we employ the maximun-likelihood
estimator developed in [62, 63], with a maximum boundary for the power law, so
that the pdf in equation 36 is integrabled over
$x\in[x_{\text{min}},x_{\text{max}}]$ for any $\alpha>0$. It is reasonable to
assume that power-law distributions are bounded within the context of this
study, given that geographical limitations inherently bound the movement
distances in our data set. We define $x_{\text{max}}$ as the diameter of each
domain $\Omega$, which corresponds to the largest distance between two points
(Table.,1). We assess the power-law parameter estimates through goodness-of-
fit tests employing the Kolmogorov–Smirnov (KS) statistic and likelihood
ratios. The Kolmogorov–Smirnov test, presented with p-values in Table. 1.
Likelihood ratios provide evidence on whether the power-law distribution is a
superior fit compared to other distributions. The methodology for estimating
piece-wise power laws for each Danish city (Fig. 3f-g) is based on a maximum
likelihood estimator derived in [64], the details of the fit are available SI.
8.4 and SI. 10. In the analysis of the Mobility City Radius the power Pareto
distribution [64] is identified as the optimal model with a slightly better
fit, measured by an Akaike Information Criterion (AIC) score of 24046.03. This
is compared to the log-normal distribution, which has an AIC of 24161.17. A
comparison with additional distributions is available in Table 2.

The mesoscale power law of the pair distribution function (Fig. 2a) and the
simulation of the fractal and patch model (Fig. 2fv) are also estimated with a
maximum likelihood estimator.

##  5 Figures

Figure 5: Pair distribution of buildings and addresses in France shows three
scales: (I) neighborhood-scale with linear onset and oscillating modulation,
(II) city-scale with power-law growth ($\alpha=0.67$, for addresses), and
(III) country-scale with slower growth and rapid decay (see SI. 6.4)

Figure 6: Demographic characteristics of the Danish Residential Mobility
dataset. The gender imbalance is the first panel due to the $1.07$ male/female
ratio rate at birth [65]. Indeed we study the total population, not the
currently alive population, that is balanced for gender due higher death rate
for men. Therefore although there is an imbalance in gender, the dataset is
not biased. The second panel shows the birth year of every individual in the
dataset, it closely follows the birth rate of Denmark, and the differences are
due to migrations [66].

Figure 7: Comparison between our framework, the gravity model [5] and the
radiation model [6] for cities.

Table 1: Power law estimation Type | Dataset | $n$ | $x_{\text{min}}$ | $x_{\text{max}}$ | Exponent | KS | p
---|---|---|---|---|---|---|---
Residential Mobility | Denmark | 39M | 10 | 475032 | 0.98 | 0.01 | 0.99
| France | 40M | 3000 | 1134667 | 1.07 | 0.01 | 0.98
Day-To-Day Mobility | San Francisco | 112168 | 10 | 55639 | 0.94 | 0.07 | 0.96
| Houston | 47996 | 43 | 71209 | 0.95 | 0.03 | 0.98
| Singapore | 15167 | 10 | 79674 | 0.89 | 0.12 | 0.67
Geography | City Radius | 1367 | 2878 | N/A | $3.05$ | $0.58$ | 0.11

The table presents information about the estimation of the exponent of the
power law. The estimation of the power law is based on a maximum-likelihood
estimator as in [62]. The values of $x_{\text{min}}$ and $x_{\text{max}}$ are
reported in meters. The table reports also the p-values associated with the
Kolmogorov-Smirnov test for the power law exponent estimators [61]. According
to [61], the suitability of a power law model is considered statistically
implausible if the p-value is less than or equal to 0.1. This threshold
indicates a less than or equal to $10\%$ probability that the observed
deviation of the data from the model predictions is due to random variation
alone. Therefore the power-law model is statistically significant if $p>0.1$
(in bold).

Table 2: Log-likelihood Ratio between power law and other heavy-tailed distributions Type | Dataset | Log-Normal | Exponential | Weibull
---|---|---|---|---
Residential Mobility | Denmark | $128697^{***}$ | $19209^{***}$ | $89676^{***}$
| France | $5042^{***}$ | $411^{***}$ | $18535^{***}$
Day-to-day Mobility | San Francisco | $1091^{***}$ | $91^{***}$ | $2001^{***}$
| Houston | $446^{***}$ | $35^{**}$ | $2946^{**}$
| Singapore | $729^{***}$ | $104^{***}$ | $3651^{***}$
Geography | Mobility City Radius | $-0.96$ | $1.91^{*}$ | $1.18$

The table presents the log-likelihood ratio $\mathcal{R}$ (ref. [61])
comparing the power law (cf. Table 1) to other heavy-tailed distributions (one
per column). When $\mathcal{R}$ is positive, the power law distribution has a
higher likelihood compared with the alternative. When $\mathcal{R}$ is
negative, the other distribution has a higher likelihood compared with the
power law. The table reports also the p-values associated with $\mathcal{R}$
$(^{*}P\leq 0.05,\,^{**}P\leq 0.01,\,^{***}P\leq 0.001)$.

#### Data availability

The data that support the findings of this study are available on the Zenodo
repository: However, restrictions apply to the availability of these data on
in some cases only the anonymized and aggregated data is publicly available.
Source data are provided with this paper, see code repository.

#### Code availability

Code is available at https://github.com/LCB0B/role-of-geo/. The repository
contains the geography model, the statistical estimation codes, the source
data, the code for the figures, and additional figures.

#### Acknowledgement

The authors thank J. Dzubiella and YY. Ahn for helpful comments in the early
development of this study, as well as L. Alessandretti and S. De Sojo Caso for
providing insightful comments on the manuscript. The work was supported in
part by the Villum Foundation and the Danish Council for Independent Research.

#### Author contributions

L.B., B.M. and S.L. designed the study and the model. L.B. and B.M. performed
the analyses and implemented the model. L.B., B.M. and S.L. analyzed the
results and wrote the paper.

## Supplemental Information

###### Contents

  1. 1 Mobility takes place on the physical structure of the built environment
  2. 2 Multi-scale urban structure shapes the distribution of pair distribution
  3. 3 Intrinsic Distance Cost function
  4. 4 Discussion
  5. 4.1 Data
    1. 4.1.1 Uncertainty Quantification
    2. 4.1.2 Residential Mobility in Denmark
    3. 4.1.3 Bias in residential mobility in Denmark
    4. 4.1.4 Residential mobility in France
    5. 4.1.5 Day-to-Day Mobility
    6. 4.1.6 GeoBoundaries
  6. 4.2 Pair distribution function
  7. 4.3 Locations as interacting particles
    1. 4.3.1 Molecular Dynamics simulation
    2. 4.3.2 Hard disks
  8. 4.4 Non-overlapping disks
  9. 4.5 City pair distribution and radius
    1. 4.5.1 Fractal Model
    2. 4.5.2 Patch Model
  10. 4.6 A ‘geography-free’ gravity law
  11. 4.7 Power-Law Estimation
  12. 5 Figures
  13. 6 The concept of pair distribution function distribution
    1. 6.1 Pair distribution examples
    2. 6.2 The pair distribution function emerges as a normalization for human mobility
    3. 6.3 Conditions for the pair distribution function to identify a set of points uniquely
    4. 6.4 Pair distribution function for France
  14. 7 Models for pair distribution function
    1. 7.1 Pair distribution functions
      1. 7.1.1 Generalized pair distribution function
      2. 7.1.2 Circle
      3. 7.1.3 Parabola
      4. 7.1.4 Building locations by external and interaction potentials
      5. 7.1.5 Molecular Dynamics simulation
      6. 7.1.6 Emulating attractive hard disks
    2. 7.2 Collision-resolving algorithm
    3. 7.3 Random non-overlap positioning algorithm of disks
    4. 7.4 Influence of building location definition on pair distribution
    5. 7.5 Measuring the pair distribution function for large dataset
      1. 7.5.1 Small distances
      2. 7.5.2 Larger distances
      3. 7.5.3 Between cities
    6. 7.6 Pair distribution from independent patches of heterogeneous size
    7. 7.7 Pair distribution of a self-similar modular hierarchical model of building locations
  15. 8 Studying mobility with the pair distance function
    1. 8.1 Comparison with other Mobility models
      1. 8.1.1 Gravity
      2. 8.1.2 Radiation
    2. 8.2 The intrinsic distance cost as the inverse of distance
    3. 8.3 Reconciling Opportunity and Distance-Based Mobility Paradigms
    4. 8.4 Local piece-wise power law
      1. 8.4.1 The remarkable case of islands
    5. 8.5 Simulation of piece-wise intrinsic distance cost
  16. 9 The definition of city
  17. 10 Maximum likelihood estimation of power laws

###### List of Figures

  1. 1 Illustration of the influence of geography on human mobility, (a) spatial arrangement of concentric circles of locations, (b) distance distribution of points to center of the disk (white cross). (c) Pair distribution of the points in the concentric circles, (d) simulated moving distance distribution for different cost of distance. Examples of more shapes and distributions of locations are studied in SI. 6.1. (e,f,g) Three geographies: (e) Real Denmark with 3.3M addresses, (f) Disk Denmark, where the macro structure has cities distributed uniformly at random on a disk, but the microstructure of the cities is preserved, and (g) Uniform Denmark with the macro structure of Denmark (shape of landmasses and city positions) but no local microstructure, i.e. uniformly distributed locations within cities. (h) Pair distribution of the three geographies, and (i) the simulated observed moving distance for an intrinsic distance cost of $1/r$. In (h) and (i), I, II, III respectively represent the intra-city scale, the city scale, and the inter-city scale.
  2. 2 (a) Pair distribution of residential locations in Denmark. (I) Neighborhood-scale with a linear onset and an oscillating modulation for the immediate neighborhood, (II) city-scale showing a power-law growth with exponent $\alpha=0.67$, and (III) country-scale with a slower growth and eventual fast decay. (b) With external potential $V(|\bm{x}|)$ that grows linearly with distance to the city center $|\bm{x}|$, an ideal gas of buildings will form a city-like structure with radial building density following an Erlang distribution of scale $R_{0}$ (squares) (Methods: Model of locations). The corresponding pair distribution approximately assumes a generalized Gamma distribution with scale $R$ providing a definition for the radius of a city (circles). (c) Maximum-likelihood fits of the generalized pair distribution against pair distribution of buildings in three Danish towns (upper) and all Danish towns, rescaled (lower). (d) Pair distributions for four systems of building positions on the micro-scale, (i) buildings in Denmark, (ii) random positions (ideal gas), (iii) Lennard-Jones disks, and (iv) attractive hard disks of heterogeneous size. All systems display linear onset with oscillatory modulations. The pair-correlation function $g(r)$ quantifies the under- and over-presence, respectively, of neighbors at distance $r$ when comparing to the ideal gas. Repulsive and attractive forces lead to valleys and peaks in $g(r)$ (see Methods: Model of locations, Hard disks and Molecular Dynamics simulation) (e) Pair distribution for all towns in Denmark for (i) random positions, (ii) non-attractive hard disks, (iii) real data. Systems (ii) and (iii) display similar behavior with slow onset, quickly approaching the ideal gas. Weighting with the numbers of building pairs per pair of towns (pink and green lines) replicates the macro regime behavior (III) of panel (a) but introduces slight deviations to the model system (ii). (f) (i) Fractal model for building positions in cities, as a hierarchically nested structure of self-similar patches. On each hierarchy layer, a patch of buildings consists of $b$ sub-patches distributed non-overlapping with constant density and packing fraction $\theta$. At the lowest level, a sub-patch is a building. (ii) A sample configuration of this model with $b=4$, $z_{1}=3\,\mathrm{m}$, $\theta=0.8$, and $L=8$. (iii-iv) Independent-patch model. Equating the concept of “town” with “patch of radius $R$ containing buildings”, the number of buildings per patch grows as $R^{2}$ and the number of pairs of buildings per patch grows as $R^{4}$ (with $R$ from the analysis in panel c). At the same time, the tail of the patch radius distribution falls as $R^{-(4-\alpha)}$ with $\alpha=0.67$. (v) For both the fractal model and the independent-patch model, we recover the pair distribution scaling observed in the data.
  3. 3  (a) Empirical moving distance distribution and pair distribution in Denmark. Dividing the mobility by the pair distribution yields an intrinsic distance cost (b) that follows a power law from 10m to 500km (maximum likelihood estimation). (c) This framework can be represented as a continuous gravity law, if we take a toy model of two distant cities, (d) the pair distribution at the distance $d$ between the two cities collapses into a Dirac distribution, with the height equal to the product of cities population, $N_{1}N_{2}$. (e) The toy intrinsic distance cost of the simulation model. For a comparison with a continuous radiation model, see SI. 8.1.2. (f-g) Residential mobility analysis in a small (f) and a large (g) Danish town reveals a consistent piece-wise power law across all cities. For more cities, see SI Fig.28-31, for details on the piece-wise power-law definition and maximum likelihood estimation see SI. 8.4. (h) Distribution of the piece-wise power laws’ exponents, for the first and second parts, respectively center around $0.6\pm 0.20$ and $2.0\pm 0.21$ (mean $\pm$ standard deviation). (i) The transition point in the piece-wise power law aligns with a city size distribution. We fit different distributions (log-normal, power-Pareto) via maximum-likelihood (see Methods: Power-law estimation). (j) We recover the intrinsic distance cost from simulating the piece-wise power law (f,g) with the mean values of $0.6$ and $2$ (h), over the city-scale (i). We assume random position of cities as in Fig. 2e. For alternative derivation of the intrinsic distance cost and details see SI. 8.5.
  4. 4 (a) Observed moving distance distribution and pair distribution for residential mobility in France. (b) Observed moving distance distribution and pair distribution for day-to-day mobility between points of interest in San Francisco, Houston (c), and Singapore (d).
  5. 5 Pair distribution of buildings and addresses in France shows three scales: (I) neighborhood-scale with linear onset and oscillating modulation, (II) city-scale with power-law growth ($\alpha=0.67$, for addresses), and (III) country-scale with slower growth and rapid decay (see SI. 6.4)
  6. 6 Demographic characteristics of the Danish Residential Mobility dataset. The gender imbalance is the first panel due to the $1.07$ male/female ratio rate at birth [65]. Indeed we study the total population, not the currently alive population, that is balanced for gender due higher death rate for men. Therefore although there is an imbalance in gender, the dataset is not biased. The second panel shows the birth year of every individual in the dataset, it closely follows the birth rate of Denmark, and the differences are due to migrations [66].
  7. 7 Comparison between our framework, the gravity model [5] and the radiation model [6] for cities.
  8. 8 Example of pair distribution function for different toy geographies and density of points. (a) Four different shapes of landmasses and their respective pair distribution functions (same maker shape and greyscale). (b) For the circle shape, four different toy layouts of the built environment and their respective pair distribution function (same greyscale).
  9. 9 Given a ring with uniformly distributed points, the pair distance distribution of points will be constant.
  10. 10 Map of 34,041,910 individual addresses in France (a), and zoom on the Paris area (b).
  11. 11 pair distribution of buildings and addresses in France. Neighborhood-scale with a linear onset and an oscillating modulation for the immediate neighborhood, then city-scale showing a power-law growth with exponent $\alpha=0.67$, and finally country-scale with a slower growth and eventual fast decay.
  12. 12 Energy evolution of a single MD simulation of a canonical ensemble of LJ particles in a linear external potential, thermalized by a stochastic Berendsen thermostat. Energies are displayed per particle in units of the LJ potential depth $\varepsilon$. Time is units of the stochastic Berendsen relaxation time $\tau$.
  13. 13 Three models for positioning buildings within an external potential of shape $V(\bm{x})=\gamma|\bm{x}|$. (Left) An ideal gas, forming with radial density $p(|\bm{x}|)=(|\bm{x}|/R_{0}^{2})\exp(-|\bm{x}|/R_{0})$. (Middle) A configuration of Lennard-Jones disks. (Right) A configuration of heterogeneously sized hard disks with abstract attractive force.
  14. 14 Three models for city positioning within Denmark. (Left) An ideal gas, i.e. positions chosen uniform at random within the administrative boundaries. (Middle) Hard disks are placed in a non-overlapping, random manner (see Sec. 7.3). Radiuses are equal to the inferred $R$ values from the pair distribution MLE of Eq. 47 to Denmark’s towns with more than 30 buildings. (Right) Cities placed at the centroid of their respective administrative boundaries with radius equal to $R$ inferred from the respective city’s empirical pair distribution .
  15. 15 How the LJ pair distribution and pair-correlation function changes when building location is redefined to (i) randomly within the disk of radius $z$ and (ii) randomly on the rim of disk of radius $z$.
  16. 16 Inferred values of the sub-linear scaling exponent $\alpha$ in $p(r)\propto r^{\alpha}$ for $r\in(25\,\mathrm{m},200\,\mathrm{m})$ for the 30 largest Danish cities. (Left) buildings. (Right) addresses.
  17. 17 Building pair distribution functions and pair-correlation functions $g(r)$ with $r<200\,\mathrm{m}$ for (top) every building of Denmark, (second from top) Hovedstadens (capital region), (second from bottom) Aarhus, and (bottom) every building not located with the administrative boundaries of any city/town.
  18. 18 Address pair distribution functions and pair-correlation functions $g(r)$ with $r<200\,\mathrm{m}$ for (top) every building of Denmark, (second from top) Hovedstadens (capital region), (second from bottom) Aarhus, and (bottom) every address not located with the administrative boundaries of any city/town.
  19. 19 Empirical building pair distribution functions for the largest 100 cities in Denmark and fits of $p(r,R,m)=(r/R^{2})\exp(-(r/R)^{m})/\Gamma(2/m)$. We mark the inferred radius $R$.
  20. 20 Empirical address pair distribution functions for the largest 100 cities in Denmark and fits of $p(r,R,m)=(r/R^{2})\exp(-(r/R)^{m})/\Gamma(2/m)$. We mark the inferred radius $R$.
  21. 21 (Upper left) empirical pair distributions for all Danish towns with more than 30 buildings, rescaled by $R$ and $m$ obtained from (upper right) MLE fits of $p(r,R,m)=(r/R^{2})\exp(-(r/R)^{m})/\Gamma(2/m)$ to the empirical data. (Lower left) Distribution of inferred $R$ values and power-law fit to its tail, giving $p(R)\propto R^{-\alpha}$ with $\alpha=3.29$ and $x_{\min}=824\,\mathrm{m}$. (Lower right) Distribution of inferred $m$ values with $\left<m\right>=1.61$ and $\mathrm{Std}[m]=0.44$.
  22. 22 building-pair-weighted pair distribution function of cities, both for real city positions as well as a hard disk model without attractive force. Overlaid are the building pair distribution and the ideal gas pair distribution.
  23. 23 Integrating over the multiverse of patches of size $R$, summing up the contributions to the pairwise-distance distribution of all patches, for three assumptions of (i) generalized pair distribution, (ii) pair distribution of a uniform, circular patch of buildings, and (iii) a parabola pair distribution, with varying $\alpha$. Here, we used a minimum patch size of $R=10\mathrm{m}$, which fixes a scale and leads to linear growth of the pairwise-distance distribution for small distances.
  24. 24 Upper two panels: Theoretical pairwise-distance distributions of the hierarchical building placement model with $N=4$ and $\theta=0.8$, in color the single contributions Eq. (109) of scale $k$, in black the total joint distribution function over all scales Eq. (110), including a power-law fit showing the resulting scaling exponent $\alpha\approx 2/3$. Lower two panels: A sample from the hierachical placement model with $L=8$ hierarchy layers and the corresponding resulting pairwise-distance distribution of the whole structure, showing the desired scaling behavior on the mesoscale.
  25. 25 Gravity model as in [5] for migration between (left) municipalities in Denmark, (right) parishes in Denmark. The normalized flow represented on the $y$-axis is the number of moves from $i$ to $j$ divided by the product of population between $i$ and $j$, $T_{i,j}/(N_{i}N_{j})$.
  26. 26 Piecewise intrinsic distance for the city of Rønne and Nexø both situated on an island.
  27. 27 For each city in Denmark, the distribution of observed moving distance (orange) and the relative pair distribution function. Pairs are restricted to those containing at least one address from the city of interest. The second plot shows the intrinsic distance cost (black), as the ratio between the observed moving distance and the relative pair distribution function. The pink dashed line is the inferred mobility city radius, the light teal corresponds to intra-city moves, and the dark teal to inter-city moves.
  28. 28 For each city in Denmark, the distribution of observed moving distance (orange) and the relative pair distribution function. Pairs are restricted to those containing at least one address from the city of interest. The second plot shows the intrinsic distance cost (black), as the ratio between the observed moving distance and the relative pair distribution function. The pink dashed line is the inferred mobility city radius, the light teal corresponds to intra-city moves, and the dark teal to inter-city moves.
  29. 29 For each city in Denmark, the distribution of observed moving distance (orange) and the relative pair distribution function. Pairs are restricted to those containing at least one address from the city of interest. The second plot shows the intrinsic distance cost (black), as the ratio between the observed moving distance and the relative pair distribution function. The pink dashed line is the inferred mobility city radius, the light teal corresponds to intra-city moves, and the dark teal to inter-city moves.
  30. 30 For each city in Denmark, the distribution of observed moving distance (orange) and the relative pair distribution function. Pairs are restricted to those containing at least one address from the city of interest. The second plot shows the intrinsic distance cost (black), as the ratio between the observed moving distance and the relative pair distribution function. The pink dashed line is the inferred mobility city radius, the light teal corresponds to intra-city moves, and the dark teal to inter-city moves.
  31. 31 For each city in Denmark, the distribution of observed moving distance (orange) and the relative pair distribution function. Pairs are restricted to those containing at least one address from the city of interest. The second plot shows the intrinsic distance cost (black), as the ratio between the observed moving distance and the relative pair distribution function. The pink dashed line is the inferred mobility city radius, the light teal corresponds to intra-city moves, and the dark teal to inter-city moves.
  32. 32 For each city in Denmark, distribution of observed moving distance (orange) and the relative pair distribution function. Pairs are restricted to those containing at least one address from the city of interest. The second plot shows the intrinsic distance cost (black), as the ratio between the observed moving distance and the relative pair distribution function. The pink dashed line is the inferred mobility city radius, the light teal corresponds to intra-city moves, the dark teal to inter-city moves.
  33. 33 For each city of the piecewise fit, Figure.28-32, we compare the log-likelihood ratio between the piecewise power law distribution of equation 152 and a Pareto distribution. The left plot shows the distribution of p values with a pink line at $p=0.05$. The right plot shows the distribution of the log-likelihood ratio, with a pink line at the decision boundary $R=0$. A positive log-likelihood ratio indicates that the piecewise models describe the data better, the p-value give the significance of the result.
  34. 34 For each city of the piecewise fit, Figure.28-32, we compare the log-likelihood ratio between the piecewise power law distribution of equation 152 and a lognormal distribution. The left plot shows the distribution of p values with a pink line at $p=0.05$. The right plot shows the distribution of the log-likelihood ratio, with a pink line at the decision boundary $R=0$. A positive log-likelihood ratio indicates that the piecewise models describe the data better, the p-value give the significance of the result.
  35. 35 For each city of the piecewise fit, Figure.28-32, we compare the log-likelihood ratio between the piecewise power law distribution of equation 152 and a piecewise exponential-Pareto distribution. The left plot shows the distribution of p values with a pink line at $p=0.05$. The right plot shows the distribution of the log-likelihood ratio, with a pink line at the decision boundary $R=0$. A positive log-likelihood ratio indicates that the piecewise models describe the data better, the p-value give the significance of the result.
  36. 36 Comparison between the Danmarks Statistik definition of cities (on the left) and the cluster obtained from HDBSCAN (on the right). The colors are ordered according to the size of the cluster/city.
  37. 37 Heatmap of the NMI between the clustering of addresses locations and their official city and the values for different combinations of HDSCAN parameters.
  38. 38 Kolmogorov-Smirnov (KS) test for Houston data, the blue dots represent the distribution of the test statistics over the sampled data. The mean is represented by a dashed black line. The 5th and 95th percentiles are represented by dashed green lines, indicating the range within which 90% of the observed values fall.
  39. 39 Kolmogorov-Smirnov test for Singapore data, the blue dots represent the distribution of the test statistics over the sampled data. The mean is represented by a dashed black line. The 5th and 95th percentiles are represented by dashed green lines, indicating the range within which 90% of the observed values fall.
  40. 40 Kolmogorov-Smirnov test for San Francisco data, the blue dots represent the distribution of the test statistics over the sampled data. The mean is represented by a dashed black line. The 5th and 95th percentiles are represented by dashed green lines, indicating the range within which 90% of the observed values fall.
  41. 41 Kolmogorov-Smirnov test for France data, the blue dots represent the distribution of the test statistics over the sampled data. The mean is represented by a dashed black line. The 5th and 95th percentiles are represented by dashed green lines, indicating the range within which 90% of the observed values fall.

##  6 The concept of pair distribution function distribution

###  6.1 Pair distribution examples

This section presents additional examples of toy geographies that vary both
the shape of the landmasses and the density distribution of points. These
examples intend to demonstrate that the pair distribution distribution can
encode the geography of the built environment. Figure 8a, illustrates how the
pair distribution function captures the shape of uniformly distributed points.
Figure 8b illustrates the impact of spatial point distribution on a fixed
shape.

h

Figure 8: Example of pair distribution function for different toy geographies
and density of points. (a) Four different shapes of landmasses and their
respective pair distribution functions (same maker shape and greyscale). (b)
For the circle shape, four different toy layouts of the built environment and
their respective pair distribution function (same greyscale).

###  6.2 The pair distribution function emerges as a normalization for human
mobility

Whenever points $x_{i}\in\Omega$ are distributed in a $d$-dimensional subspace
$\Omega\subseteq\mathbb{R}^{d}$ according to some density $p^{d}(x)$, their
pair distribution function $p(r)$ is one-dimensional and determined by both
the geometry of the space as well as $p^{d}(x)$, given that there exists a
(pseudo-) metric $r\equiv R(x,y)$ that determines a distance between two
points $x$ and $y$. This pair distribution is given by

$\displaystyle p(r)=\int_{\Omega}dx\ \int_{\Omega}dy\
p^{d}(x)p^{d}(y)\delta(r-R(x,y)).$ (37)

When observations are made between two distinct points (i.e. when a quantity
is counted) in the domain of interest, a distance-dependent distribution of
the observed entities is observed. This distribution is hereafter referred to
as the function $f(r)$. These observations can be quantified as the number of
packages sent from location A to location B at distance r, the number of
people moving from A to B at distance r, or the number of power line
connections that connect substations at distance r in a power grid. The
observed distribution function $f(r)$ is the result of two factors: the number
of pairs of points that exist in the domain of interest, $\Omega$, with
associated point density $p^{d}(x)$, and the one-dimensional probability
$\pi(r)$ with which pairs of distance $r$ are manifested in the real world to
be observed. In short, the number $f(r)$ of observations of distance $r$, is
determined by the number of possible pairs $p(r)$ at distance $r$ and the
probability $\pi(r)$ that they would exist at this distance, i.e.

$\displaystyle f(r)\propto\pi(r)p(r),$ (38)

bar a normalization constant.

When we measure $f(r)$, we always measure with it the geometry of the subspace
$\Omega$ as well as the density $p^{d}(x)$, manifested in the pair
distribution function $p(r)$. The “universal”, i.e. geometry-independent law
that encapsulates the behavior of the system we want to study, is encoded in
$\pi(r)$. To properly deduce this geometry-independent behavior of our system,
we need to adjust our observation $f(r)$ by the geometry encoded in $p(r)$,
i.e.

$\displaystyle\pi(r)\propto\frac{f(r)}{p(r)}.$ (39)

A pertinent question is which reference topology (geometry and point
distribution) would result in an observation $f$ that is directly proportional
to the behavior $\pi$. This question holds practical significance, as
considering behavior benefits greatly from a conceptual framework for the
distribution of points within a space. Indeed, our goal is to understand the
mechanisms that link two points within this space. Without defining these
points or the space itself, this endeavor becomes somewhat pointless.

In our formalism, this implies that we are seeking a topology where,

$\displaystyle f(r)\propto\pi(r),$ (40)

which implies

$\displaystyle p(r)=\mathrm{const.}$ (41)

Looking at Eq. (37), we have to find $\Omega$ and $p^{d}(x)$ such that
$p(r)=\mathrm{const}.$ While there might be a multitude of solutions, the
simplest one is a one-dimensional ring, which can be conceptualized as a box
of length $L$ and domain $x\in[0,L)$ with periodic boundary conditions, i.e.
an associated distance

$\displaystyle R(x,y)=\begin{cases}|x-y|,&0\leq|x-y|\leq L/2,\\\
L-|x-y|&L/2\leq|x-y|<L\end{cases}$ (42)

and uniformly distributed points, i.e.

$\displaystyle p^{d=1}(x)=\begin{cases}L^{-1},&\ 0\leq x<L,\\\ 0,&\
\mathrm{otherwise}.\end{cases}$ (43)

then the pair-wise distance distribution evaluates to,

$\displaystyle p(r)$
$\displaystyle=L^{-2}\int_{0}^{L}dx\left(\int_{x-L/2}^{x}dy\delta(r-(x-y))+\int_{x}^{x+L/2}dy\delta(r-(y-x))\right)$
(44) $\displaystyle=\begin{cases}L^{-2}\int_{0}^{L}dx\left(1+1\right),&\ 0\leq
r\leq L/2,\\\ 0,&\ \mathrm{otherwise}\end{cases}$ (45)
$\displaystyle=\begin{cases}(L/2)^{-1}&\ 0\leq r\leq L/2,\\\ 0,&\
\mathrm{otherwise},\end{cases}$ (46)

for an illustration see Fig. 9. Hence, the observed distance distribution
$f(r)$ occurring between pairs of locations at distance $r$ on this topology
will be proportional to the geometry-independent behavioral part $\pi(r)$.
This technique has been used in [36, 58, 67, 37].

Figure 9: Given a ring with uniformly distributed points, the pair distance
distribution of points will be constant.

###  6.3 Conditions for the pair distribution function to identify a set of
points uniquely

According to [68], in general, point configurations can be uniquely determined
by their pair distribution distributions, up to a rigid transformation.
However. there are counterexamples, for example when two distances are equal.

In the context of geographic analysis, for a set of 2D points representing
addresses, some distances between points will inevitably be repeated. If we
consider $1\times 10^{6}$ points located in a $1000\,\mathrm{km}\times
1000\,\mathrm{km}$ square and whose coordinates are known with a
$1\,\mathrm{m}$ precision, by a simple combinatorial argument, the number of
pair distances is $1\times 10^{12}.$ Although the maximum distance in the
square is $\sqrt{2}\times 1000\,\mathrm{km}$, or $1.414\times
10^{6}\,\mathrm{m}$, which means that the set of possible distances contains
$1.414\times 10^{6}$ different values (due to the limited precision), but
there are $1\times 10^{12}$ instances, so some distances must be equal, and
the pair distribution function does not uniquely identify a set of 2d points.

However, in the case of geography, due to the regularities and scaling law of
the pair distribution, one can adopt a coarse-grained view of the problem.
Instead of considering each point individually, we can first examine the pair
distribution function between urban centers or other areas and iteratively
reconstruct the geography (set of 2d points) in iterative. First, each urban
center position should be a unique configuration, then independently on the
local geography around each city. This should lead to a unique configuration
(up to isometries that are of the second order for the pair distribution
function). The coarse-grained construction would be similar to the quad-tree
partitioning or the HDSCAN clustering (see section 9).

###  6.4 Pair distribution function for France

Figure 10: Map of 34,041,910 individual addresses in France (a), and zoom on
the Paris area (b).  Figure 11: pair distribution of buildings and addresses
in France. Neighborhood-scale with a linear onset and an oscillating
modulation for the immediate neighborhood, then city-scale showing a power-law
growth with exponent $\alpha=0.67$, and finally country-scale with a slower
growth and eventual fast decay.

The pair distribution function between residential buildings in France shows
similar patterns as the one of Denmark, as shown in Fig.11. On the micro-scale
(i.e. within distances of $r\lesssim 25\,\mathrm{m}$, the “neighborhood”). We
observe a linear onset of neighborhood density, oscillatory modulated. For
larger distances (mesoscale), this growth assumes a scaling of approximately
$p(r)\propto r^{0.67}$ between $25\,\mathrm{m}$ and $200\,\mathrm{m}$, which
is even more pronounced in the distribution of address distances (up to
$10000\,\mathrm{m})$. The mesoscale has a larger amplitude for France than for
Denmark ($10\,\mathrm{km}$ vs $4\,\mathrm{km}$)due to larger urban areas
(Paris $2853\,\mathrm{km}^{2}$ vs $526\,\mathrm{km}^{2}$ for Copenhagen).
After that, at the macro scale, this growth slows down, decays rapidly, and
finally approaches zero as we reach the limit of the finite system.

##  7 Models for pair distribution function

###  7.1 Pair distribution functions

####  7.1.1 Generalized pair distribution function

In the main text, we defined the generalized pair distribution function model
as

$\displaystyle
p_{m}(r,R)=C\frac{r}{R^{2}}\exp\left(-\left(\frac{r}{R}\right)^{m}\right).$
(47)

It has a linear onset and a tail that falls as a stretched exponential and is
a special case of the generalized Gamma distribution. We find $C$ as follows.
Begin with the Gamma function

$\displaystyle\Gamma(z)=\int\limits_{0}^{\infty}t^{z-1}\exp\left(-t\right)\mathrm{d}t$
(48)

and substitute $t=(r/R)^{m}$ such that
$\mathrm{d}t=(m/R)(r/R)^{m-1}\mathrm{d}r$ and

$\displaystyle\Gamma(z)=\frac{m}{R}\int\limits_{0}^{\infty}(r/R)^{mz-m+m-1}\exp\left[-(r/R)^{m}\right]\mathrm{d}r.$
(49)

Now we demand $mz-1=1$, i.e. $z=2/m$, to find

$\displaystyle\Gamma(2/m)=m{C^{-1}}\int\limits_{0}^{\infty}C\frac{r}{R^{2}}\exp\left[-(r/R)^{m}\right]\mathrm{d}r.$
(50)

Due to the normalization condition, we have

$\displaystyle C=\frac{m}{\Gamma(2/m)}.$ (51)

To obtain the first moment, we demand $mz-1=2$ such that $z=3/m$ so we find

$\displaystyle\Gamma(3/m)$
$\displaystyle=\frac{mC^{-1}}{R}\int\limits_{0}^{\infty}C\frac{r^{2}}{R^{2}}\exp\left[-(r/R)^{m}\right]\mathrm{d}r$
(52) $\displaystyle=\frac{\Gamma(2/m)}{mR}\left<r\right>$ (53)
$\displaystyle\left<r\right>$
$\displaystyle=\frac{mR\Gamma(3/m)}{\Gamma(2/m)}.$ (54)

We can fit this distribution to data by using the per-sample log-likelihood

$\displaystyle{\mathcal{L}}=\frac{1}{n}\ln L$
$\displaystyle=\ln{m}-\ln\Gamma(2/m)+2\ln w+\left<\ln
r\right>-w^{m}\left<r^{m}\right>.$ (55)

where $w=1/R$ and we have an observational set $\\{r_{i}\\}$ of empirical
pairwise distances with sample size $n$. From $\partial\mathcal{L}/\partial
w=0$ we find that the inverse city scale that maximizes the likelihood is
given by

$\displaystyle\hat{w}(m)=\left(\frac{2}{m\left<r^{m}\right>}\right)^{1/m}.$
(56)

With $\hat{w}(m)$, we can find the zero of

$\displaystyle\frac{\partial\mathcal{L}}{\partial
m}=\frac{1}{m}+\frac{2}{m^{2}}\psi(2/m)-\Big{<}[\hat{w}(m)r]^{m}\ln[\hat{w}(m)r]\Big{>}$
(57)

numerically, which gives $\hat{m}$. Here,
$\psi(z)=\Gamma^{\prime}(z)/\Gamma(z)$ is the digamma function.

####  7.1.2 Circle

The pair distribution of a uniform distribution of random points within a disk
of radius $R$ is given by

$\displaystyle p(r,R)=\begin{cases}\frac{4r}{\pi
R^{2}}\arccos\left(\frac{r}{2R}\right)-\frac{2r^{2}}{\pi
R^{3}}\sqrt{1-\frac{r^{2}}{4R^{2}}},&\qquad r\leq 2R,\\\
0&\qquad\mathrm{otherwise,}\end{cases}$ (58)

see [69].

####  7.1.3 Parabola

Looking at Eq. (58), we see that this distribution looks somewhat close to a
parabola with zeros at $r=0$ and $r=2R$. A parabolic pair distribution with
such properties is given as

$\displaystyle
p(r,R)=\begin{cases}-\frac{3}{4R}\left(\frac{r}{R}\right)^{2}+\frac{3r}{2R^{2}}&\qquad
r\leq 2R,\\\ 0&\qquad\mathrm{otherwise.}\end{cases}$ (59)

Note that this is the Beta distribution with $\alpha=2$ and $\beta=2$ for
random variable $x=r/R$.

####  7.1.4 Building locations by external and interaction potentials

Consider the location of a city as the literal center of interest, for
example, the “central business district”. Since it might be attractive for
individuals to reach this center as quickly as possible (because amenities
will be close to the center), we assume that there is an increased cost of
living at a distance $|\bm{x}|$ from the center. For example, consider having
to commute to a job within the central business district, which has a cost
that increases with $|\bm{x}|$. At the same time, two simple mechanisms will
prevent buildings from accumulating in the exact center of the city. First,
buildings have a certain average radius $z$, so they cannot be too close
together (at a distance $<2z$). There is an advantage to buildings not being
too far apart, as they can share local amenities. Instead of modeling this
explicitly, we simply assume that there is an inherent temperature $T$ in the
system, according to which house locations are distributed following an
interaction and an external potential.

From a statistical physics point of view, we describe the system with the
simplest external potential, which increases linearly with distance

$\displaystyle V^{\mathrm{ext}}(\bm{x})=\gamma|\bm{x}|$ (60)

where we assume that the origin of the coordinate system is in the center of
the city. To model repulsion and attraction between houses, we also assume a
Lennard-Jones interaction potential

$\displaystyle
V^{\mathrm{int}}(\bm{x}_{i},\bm{x}_{j})=\varepsilon\left[\left(\frac{2z}{|\bm{x}_{j}-\bm{x}_{i}|}\right)^{12}-2\left(\frac{2z}{|\bm{x}_{j}-\bm{x}_{i}|}\right)^{6}\right],$
(61)

which is commonly used to model simultaneous attraction and repulsion between
molecules in chemical solutions [42, 35].

In total, a system with these properties evolves according to the Hamiltonian

$\displaystyle\mathcal{H}(\\{\bm{x},\bm{p}\\})=\frac{1}{2}\sum_{i=1}^{N}\bm{p}_{i}^{2}+\gamma\sum_{i=0}^{N}|\bm{x}_{i}|+\frac{1}{2}\sum_{i=1}^{N}\sum_{j\neq
i}^{N}V^{\mathrm{int}}(\bm{x}_{i},\bm{x}_{j})$ (62)

with two-dimensional momenta $\bm{p}_{i}$ and locations $\bm{x}_{i}$. In a
canonical-ensemble formulation of the system, i.e. at constant inverse
temperature $\beta$, the probability of finding a configuration
$\\{\bm{x},\bm{p}\\}$ in volume-element $d^{2}\bm{x}d^{2}\bm{p}$, is given by

$\displaystyle\varrho[\\{\bm{x},\bm{p}\\}]d^{2}\bm{x}d^{2}\bm{p}$
$\displaystyle=\exp[-\beta\mathcal{H}(\\{\bm{x},\bm{p}\\})]d^{2}\bm{x}d^{2}\bm{p}.$
(63)

For now, we restrict ourselves to an ideal gas with $V^{\mathrm{int}}=0$,
which will allow us to say something about the density of particles around the
center of the city, i.e., we want to find the probability of a particle being
present $p(r)dr$. Without loss of generality, we set $N=1$, because due to
ergodicity the trajectory of a particle will eventually follow the density of
the whole distribution (think of it as taking the $N$th root of the $N$
particle density). Integrating over the momenta yields

$\displaystyle\varrho[\bm{x}]d^{2}\bm{x}$
$\displaystyle=\exp\left(-\beta\gamma|\bm{x}|\right)d^{2}\bm{x}$ (64)

Changing the variables to polar coordinates, where $r$ is the distance of the
particle from the center, we find

$\displaystyle\varrho[r,\phi]rdrd\phi=r\exp\left[-\beta\gamma
r\right]drd\phi,$ (65)

i.e.

$\displaystyle p(r)dr=r\exp\left[-\beta\gamma r\right]dr.$ (66)

This means that if there are no interactions, the distribution of houses
around the city center should follow an Erlang distribution with scale
parameter $\lambda=\beta\gamma=1/R_{0}$, where $R_{0}$ is half the city
radius. Note that we have postulated that all particles have the same mass
$m=1$, so $\gamma r$ has the dimension of energy. The definition of the
inverse temperature $\beta=1/T$ implies that the temperature also has an
energy dimension.

Relying on the arguments of the kinetic theory of gases, we can relate the
temperature to the momenta of a particle with the identity

$\displaystyle K=N_{f}N\frac{T}{2}$ (67)

where $K=(1/2)\sum_{i}p_{i}^{2}$ is the kinetic energy and $N_{f}$ is the
degree of freedom of each particle, i.e. for single-atom particles in two
dimensions, $N_{f}=2$ (two translational degrees of freedom, no rotations, no
oscillations). Note that this relates to the root-mean-square velocity of a
single particle as

$\displaystyle v_{0}\equiv\sqrt{\left<v^{2}\right>}=\sqrt{T/2}.$ (68)

In this sense, the instantaneous temperature plays the role of a particle’s
ability to overcome the potential energy. If $V^{\mathrm{ext}}(r)=\gamma r$
represents a _cost_ of being at distance $r$ from the center, the temperature
gives a measure of how well particles in the system can overcome that cost. If
the temperature is low, particles cannot overcome this cost and the density in
the city center will be high. If the temperature is high, particles can
overcome the cost easily because of the larger amount of kinetic energy
available in the system.

Increasing complexity by going back to the Lennard-Jones perspective, we want
to obtain an intuition about how the radial particle density changes when they
strongly repel each other. In the limit of $\varepsilon/T\gg 1$, particles
will have a strong tendency to be found in their respective potential minimum,
i.e. at distance $2z$ from each other. Effectively, we can think of them as
hard disks of radius $z$ with a tendency to form clusters. If the temperature
is low, the effective radius of the city will be small (remember that the
Erlang shape parameter is $\lambda=\gamma/T=1/R_{0}$). In this case, it may
happen that the number of particles that we would expect to lie within radius
$r$ from the center (according to the ideal gas) will be greater than the
number of hard disks that can fit within a circle of radius $r$. When this
happens, we expect a crystal to form at the center.

The maximum number of disks that can fit within a circle of size $r$ can be
approximated by

$\displaystyle N_{\mathrm{disks}}(r)=\theta
A_{\mathrm{circle}}(r)/A_{\mathrm{disk}}.$ (69)

Here, $A_{\mathrm{circle}}(r)$ is the area of a circle of radius $r$,
$A_{\mathrm{disk}}$ is the area of a circle of radius $z$, and $\theta$ is an
optimal packing fraction, where we can approximately assume $\theta\approx
0.9$ for optimal hexagonal packing. Then,

$\displaystyle N_{\mathrm{disks}}(r)=\theta r^{2}/z^{2}.$ (70)

From the ideal gas distribution, we expect to find

$\displaystyle N_{\mathrm{id}}(r)$ $\displaystyle=NP(r)$ (71)
$\displaystyle=N\left(1-\mathrm{e}^{-\gamma r/T}\left(1+\frac{\gamma
r}{T}\right)\right)$ (72)

within radius $r$ (where $P(r)$ is the Erlang cumulative density function).

Following the aforementioned argumentation, we expect a nucleation effect when
$N_{\mathrm{id}}(r)>N_{\mathrm{disks}}(r)$ for $gr/T\ll 1$. Linearizing the
exponential factor, that happens when

$\displaystyle N_{\mathrm{id}}(r)$ $\displaystyle>N_{\mathrm{disks}}$ (73)
$\displaystyle N\left(1-\left(1-\frac{\gamma r}{T}\right)\left(1+\frac{\gamma
r}{T}\right)\right)$ $\displaystyle>\theta\frac{r^{2}}{z^{2}}$ (74)
$\displaystyle N\left(1-\left(1-\frac{(\gamma r)^{2}}{T^{2}})\right)\right)$
$\displaystyle>\theta\frac{r^{2}}{z^{2}}$ (75) $\displaystyle\frac{N(\gamma
r)^{2}}{T^{2}}$ $\displaystyle>\theta\frac{r^{2}}{z^{2}}$ (76)
$\displaystyle\frac{N\gamma^{2}z^{2}}{\theta T^{2}}$ $\displaystyle>1,$ (77)

or in terms of the reference velocity $v_{0}$ and the Lennard-Jones distance
$d$,

$\displaystyle\frac{N\gamma^{2}d^{2}}{8\theta v_{0}^{2}}$ $\displaystyle>1.$
(78)

The radius of this cluster will approximately be given by the solution to the
equation $N_{\mathrm{id}}(r)=N_{\mathrm{disks}}(r)$, which can be obtained
numerically.

####  7.1.5 Molecular Dynamics simulation

We are interested in finding configurations $\mathcal{C}$ that accurately
represent the canonical ensemble with number of particles $N$, constant but
irrelevant volume 111because we force the particles to be confined within a
radially symmetric external potential, the total volume containing the
particles does not matter., and average-constant temperature $T$ with total
potential energy of

$\displaystyle
V=\sum_{i=1}^{N}\left[V^{\mathrm{ext}}(\bm{x}_{i})+\frac{1}{2}\sum_{j\neq
i}V^{\mathrm{int}}(\bm{x}_{i},\bm{x}_{j})\right].$ (79)

To this end, we integrate the equations determined by the system’s Hamiltonian
numerically using the velocity-Verlet algorithm [59]. We also rescale the
particle velocities according to the stochastic Berendsen thermostat [60] with
relaxation time $\tau$. To initiate the system in a state of sufficiently low
potential energy, we assign initial particle positions according to the
corresponding ideal gas ensemble. Then we run the collision algorithm
described in Sec. with collision strength $u=1$. Moreover, we set the
temperature by defining the root-mean-square initial velocity $v_{0}$ per
particle and assigning a random velocity vector drawn from a two-dimensional
Gaussian distribution with standard deviation $v_{0}$.

As parameters, we choose $v_{0}=6$, $\gamma=0.08$, $N=10^{4}$, $\Delta
t=0.01$, $\varepsilon=20$, $z=3$, $\tau=100\Delta t$. To speed up the
numerical integration, we only consider pairs of particles that lie within
distance $r\leq 6z$, found by constructing and querying a k-d-tree for each
time step. When calculating the interaction energies, we therefore shift the
potential $V_{\mathrm{LJ}}$ so that $V_{\mathrm{LJ}}(r=6z)=0$. We integrate
the equations of motion until $t=10^{4}\tau$. The energy time series for a
single run can be seen in Fig. 12. The final configuration of this run is
shown in Fig. 13.

Figure 12: Energy evolution of a single MD simulation of a canonical ensemble
of LJ particles in a linear external potential, thermalized by a stochastic
Berendsen thermostat. Energies are displayed per particle in units of the LJ
potential depth $\varepsilon$. Time is units of the stochastic Berendsen
relaxation time $\tau$. Figure 13:  Three models for positioning buildings
within an external potential of shape $V(\bm{x})=\gamma|\bm{x}|$. (Left) An
ideal gas, forming with radial density
$p(|\bm{x}|)=(|\bm{x}|/R_{0}^{2})\exp(-|\bm{x}|/R_{0})$. (Middle) A
configuration of Lennard-Jones disks. (Right) A configuration of
heterogeneously sized hard disks with abstract attractive force.

####  7.1.6 Emulating attractive hard disks

We generate a configuration of hard disks with an unspecified attractive force
(i.e. we do not explicitly integrate the equations of motion for a hard-sphere
interaction potential with an additional attractive force). To do this, we
first generate $N=10^{4}$ random positions according to a radial Erlang
distribution with scale parameter $R=675$ (i.e. an ideal gas). Then, we draw a
random radius for each position from a heterogeneous distribution with power-
law tail. We first draw values $\hat{z}_{i}$

$\displaystyle
p(\hat{z})=\frac{a^{2}-1}{2a}\times\begin{cases}\hat{z}^{a},&\hat{z}\leq 1\\\
\hat{1}/z^{a},&\hat{z}>1\end{cases}$ (80)

and then assign disk radius $z_{i}=3\hat{z}_{i}/\left<\hat{z}\right>$. We
choose $a=4$ to obtain a heterogeneous distribution with non-finite variance
in disk area. Afterward, we run the collision algorithm outlined in Sec.7.2
with collision strength $u=1$. This leads the initially overlapping disks to
take positions where their boundaries touch, i.e. an unlikely configuration to
be found in the absence of an attractive force that would cause the disks to
lie right at each other’s boundaries.

An example configuration of this method is displayed in Fig. 13.

Figure 14: Three models for city positioning within Denmark. (Left) An ideal
gas, i.e. positions chosen uniform at random within the administrative
boundaries. (Middle) Hard disks are placed in a non-overlapping, random manner
(see Sec. 7.3). Radiuses are equal to the inferred $R$ values from the pair
distribution MLE of Eq. 47 to Denmark’s towns with more than 30 buildings.
(Right) Cities placed at the centroid of their respective administrative
boundaries with radius equal to $R$ inferred from the respective city’s
empirical pair distribution .

###  7.2 Collision-resolving algorithm

We have $N$ disks with initial positions $\bm{x}_{i}$, radii $R_{i}$, and
diameters $D_{i}=2R_{i}$, respectively. With the distance vector
$\bm{r}_{ij}=\bm{x}_{i}-\bm{x}_{j}$ and distance $r_{ij}=|\bm{r}_{ij}|$, let

$\displaystyle\mathcal{E}=\left\\{(i,j):(i<j)\wedge(r_{ij}<\max(D_{i},D_{j}))\wedge(r_{ij}<R_{i}+R_{j})\right\\}$
(81)

be the set of pairs of disks that overlap. This set can be found by iterating
over all disks $i$, finding all neighbors within distance $D_{i}$, for
instance using a k-d-tree. With this definition, let

$\displaystyle\mathcal{J}_{i}=\\{j:(i,j)\in\mathcal{E}\vee(j,i)\in\mathcal{E}\\}$
(82)

be the set of disks that overlap with disk $i$. Then,

$\displaystyle\tilde{\Delta}_{i}=u\sum_{j\in\mathcal{J}_{i}}\frac{\bm{r}_{ij}}{r_{ij}}\big{(}R_{i}+R_{j}-r_{ij}\big{)}\frac{m_{j}}{m_{i}+m_{j}}.$
(83)

is the demanded initial displacement, with the displacement rate $0<u\leq 1$
(imagine two overlapping disks—with $u=1$, the collision would be resolved
after one update). The masses $m_{j}$ control the strength of the
displacement. Suppose that disk $i$ has s small mass and a colliding disk $j$
has a great mass. The colliding disk $j$ should move less. Hence, the
influence on disk $i$ from disk $j$ should be proportional to $m_{j}$.
Assuming homogeneous density of all disks, we can set

$\displaystyle m_{i}=R_{i}^{2}.$ (84)

We also want to avoid that disks move too far per one single update to avoid
large jumps. Therefore, we move each disk by the final displacement vector

$\displaystyle\Delta_{i}={\tilde{\Delta}_{i}}\times\min(1,R_{i}/|\tilde{\Delta}_{i}|).$
(85)

i.e. the disk shouldn’t move more than its radius per update.

We update the whole ensemble of disks step by step until either

$\displaystyle|\mathcal{E}|=0$ (86)

or

$\displaystyle\max\limits_{i}|\Delta_{i}|<\epsilon.$ (87)

With a default value of $\epsilon=10^{-10}\,\mathrm{m}$.

If all disks are of equal radius $R$, we instead choose the second stop
condition as

$\displaystyle\min\limits_{(i,j)\in\mathcal{E}}r_{ij}>(1-\epsilon)2R$ (88)

and $\epsilon=10^{-3}$.

###  7.3 Random non-overlap positioning algorithm of disks

We want to randomly distribute $N$ hard disks of radii $R_{i}$ in shape
$\Omega$. We start with the largest disk of radius $R_{\mathrm{max}}$ and
iterate over all disks in decreasing order of size. For each disk $j$, we
generate a random position $\bm{x}\in\Omega$ until the condition
$|\bm{x}-\bm{x}_{i}|>R_{j}+R_{i}\ \forall j<i$ is satisfied, i.e. drawing new
random positions until there are no overlaps with other, already placed disks.
Then, assign $\bm{x}_{j}\leftarrow\bm{x}$ and continue with the next disk.

###  7.4 Influence of building location definition on pair distribution

In the dataset, the location of a building is defined as the location of the
building’s entrance door. We want to check how the Lennard-Jones ensemble’s
pair distribution and the pair-correlation function $g(r)$ change when the
location of a building is not associated with its center. To this end, we take
the configuration of Lennard-Jones disks as shown in Fig. 2d in the main text
and redefine a disk’s location to be (i) randomly within the disk and (ii)
randomly on the rim of the disk. The resulting pair distribution and $g(r)$
are shown in Fig. 15. We see that the sharp peak almost disappears for both.
At the same time, the onset of $g(r)$ becomes less abrupt, approaching a shape
similar to that observed in the data.

Figure 15: How the LJ pair distribution and pair-correlation function changes
when building location is redefined to (i) randomly within the disk of radius
$z$ and (ii) randomly on the rim of disk of radius $z$.

###  7.5 Measuring the pair distribution function for large dataset

Calculating the pair distribution function for large datasets requires a
substantial amount of memory. The number of different pairs among $N$ points
is $N(N-1)/2$. Consequently, for the dataset of 34 million address coordinates
within France (Fig.10), we need to store $10^{15}$ distance values in memory,
which is equivalent to 4,000 terabytes of data using 32-bit floats. Such a
large amount of data storage is practically unfeasible. To overcome this
computational hurdle, we proceed in two steps using a k-d-tree.

####  7.5.1 Small distances

Let $\mathcal{D}$ be the two-dimensional, non-contiguous shape containing
every building (or address, respectively) in Denmark and $\delta$ the set of
buildings (addresses) $i$ with locations $\bm{x}_{i}\in\mathcal{D}$. For each
predefined subset of this form $\Omega\subseteq\mathcal{D}$ (for example, the
official boundaries of a city), we iterate over all building (address)
positions $\omega=\\{i\in D:x_{i}\in\Omega\\}$ located within $\Omega$ and
compute the distances to all buildings (addresses) $j\in\delta\wedge j\neq i$
with $0<|x_{i}-x_{j}|<r_{\max}=200\,\mathrm{m}$. That is, for each predefined
subset of $\mathcal{D}$ (e.g. city), we find all distances of every one of its
buildings with respect to _all_ buildings (addresses) that lie within radius
$r_{\max}$, not just to those that also lie within its boundaries. We compute
this histogram with a resolution (i.e. bin width) of $1\,\mathrm{m}$. For this
task, we use a k-d-tree on all locations of $\delta$.

For each city $i$ in Denmark, let $\Omega_{i}$ be the shape that is defined by
its administrative boundaries. Note that none of the $\Omega_{i}$ overlap. We
define as

$\displaystyle\overline{\Omega}=\mathcal{D}\Big{\backslash}\left(\bigcup_{i}\Omega_{i}\right)$
(89)

the shape that includes land that is not associated with a city. Thus,
iteration over all $\Omega_{i}$ and $\overline{\Omega}$ allows us to analyze
the small-scale structure of every city, every building (address) that is not
located in a city, and - by combining all these histograms - of every building
(address) in Denmark.

We inferred the mesoscale scaling parameter $\alpha$ by fitting
$p(r)=Cr^{\alpha}$ against the respective empirical pair distribution with 1m
resolution in the range $r\in(25\,\mathrm{m},200\,\mathrm{m})$, using least
squares. We find $\left<\alpha\right>=0.69$ and $\mathrm{Std}[\alpha]=0.04$
for the building pair distribution s of the 30 largest cities as well as
$\left<\alpha\right>=0.67$, and $\mathrm{Std}[\alpha]=0.05$ for addresses, cf.
Fig. 16. Here, “largest” refers to the number of registered residential
buildings with locations within the administrative boundaries of the city. For
all buildings in Denmark, we find $\alpha=0.67$ and for all addresses, we have
$\alpha=0.69$. Example analyses can be seen in Fig. 17 for buildings and in
Fig. 18 for addresses. Note that the local environment around buildings and
addresses that are not within a city/town boundary grows much slower with
$\alpha=0.30$ for both.

Figure 16: Inferred values of the sub-linear scaling exponent $\alpha$ in
$p(r)\propto r^{\alpha}$ for $r\in(25\,\mathrm{m},200\,\mathrm{m})$ for the 30
largest Danish cities. (Left) buildings. (Right) addresses.

Figure 17: Building pair distribution functions and pair-correlation functions
$g(r)$ with $r<200\,\mathrm{m}$ for (top) every building of Denmark, (second
from top) Hovedstadens (capital region), (second from bottom) Aarhus, and
(bottom) every building not located with the administrative boundaries of any
city/town.

Figure 18: Address pair distribution functions and pair-correlation functions
$g(r)$ with $r<200\,\mathrm{m}$ for (top) every building of Denmark, (second
from top) Hovedstadens (capital region), (second from bottom) Aarhus, and
(bottom) every address not located with the administrative boundaries of any
city/town.

####  7.5.2 Larger distances

To compute the pair distribution for larger distances, we proceed as follows.
For each shape $\Omega_{i}$ (and $\mathcal{D}$, respectively), we find the set
of its building (address) locations $\omega_{i}$ (and $\delta$, respectively).
From this set we sample $n=\min(|\omega_{i}|,3\times 10^{4})$ unique locations
(without replacement). Then, we find the pairwise distances of all pairs of
these sampled locations and bin them to find histograms.

Note that to obtain the country-wide pair distribution displayed in Fig. 2a in
the main text, we combine the respective pair distributions from the small-
distance and large-distance analyses by requiring that they take the same
value at $r=184.5\,\mathrm{m}$.

We show the pair distributions for buildings and addresses for the 100 largest
Danish cities in Fig. 19 and 20. In Fig. 21 we show the empirical and fit pair
distributions of all Danish cities with more than 30 buildings as well as the
distributions of the inferred values of $R$ and $m$, respectively. The tail of
the radius distribution $R$ scales as $p(R)\propto 1/R^{3.29}$, with
$x_{\min}=824\,\mathrm{m}$, inferred by the MLE technique in refs. [70, 71].
The tail decay $(r/R)^{m}$ has values of $m$ with $\left<m\right>=1.61$ and
$\mathrm{Std}[m]=0.44$.

Figure 19: Empirical building pair distribution functions for the largest 100
cities in Denmark and fits of
$p(r,R,m)=(r/R^{2})\exp(-(r/R)^{m})/\Gamma(2/m)$. We mark the inferred radius
$R$. Figure 20: Empirical address pair distribution functions for the largest
100 cities in Denmark and fits of
$p(r,R,m)=(r/R^{2})\exp(-(r/R)^{m})/\Gamma(2/m)$. We mark the inferred radius
$R$. Figure 21: (Upper left) empirical pair distributions for all Danish towns
with more than 30 buildings, rescaled by $R$ and $m$ obtained from (upper
right) MLE fits of $p(r,R,m)=(r/R^{2})\exp(-(r/R)^{m})/\Gamma(2/m)$ to the
empirical data. (Lower left) Distribution of inferred $R$ values and power-law
fit to its tail, giving $p(R)\propto R^{-\alpha}$ with $\alpha=3.29$ and
$x_{\min}=824\,\mathrm{m}$. (Lower right) Distribution of inferred $m$ values
with $\left<m\right>=1.61$ and $\mathrm{Std}[m]=0.44$.

####  7.5.3 Between cities

To compute the pair distribution and $g(r)$ between cities, we define the
‘city center‘ as the centroid of a city’s (multi-) polygon, i.e. its geometric
center: the center of mass of its shape.

Weighting each inter-city distance $r_{ij}$ by the number of pairs of
buildings $m_{i}m_{j}$ it contains, we find a pair distribution that
approximates the empirical building pair distribution (see Fig. 22). This pair
distribution has a clear peak in $g(r)$, suggesting an abstract attractive
force between cities of a certain size. However, the onset of $g(r)$ is
consistent with the weighted pair distribution of the hard-disk model
configuration.

Figure 22: building-pair-weighted pair distribution function of cities, both
for real city positions as well as a hard disk model without attractive force.
Overlaid are the building pair distribution and the ideal gas pair
distribution.

###  7.6 Pair distribution from independent patches of heterogeneous size

We consider a _multiverse_ consisting of an infinite amount of universes
indexed by $i$, each of which is inhabited by a patch of size $R_{i}$ with
$N_{i}$ buildings inside, contributing to the multiverse with a pair
distribution number density of $n(r,R_{i})\approx N_{i}^{2}p(r,R_{i},m)$ (each
universe contributes a number of building pairs amounting to $N_{i}^{2}$).
Furthermore, we assume that there is a constant population density, which is
proportional to the number of individuals per unit area, and that this density
is inversely proportional to the square of the radius, $\rho_{0}\propto
N/R^{2}$. This implies that the number of individuals in a given area is
proportional to the area itself, $N\propto R^{2}$.

For each of these universes (or rather, each of these patches), we calculate
the number of building pairs at distance $r$ and assume that the patch sizes
$R$ are distributed according to some distribution $f(R)$, which is currently
unspecified. The joint distribution of house pairs at distance $r$ is then
given by

$\displaystyle n(r)$ $\displaystyle\propto\int\limits_{0}^{\infty}dR\
n(r,R)f(R)$ (90) $\displaystyle\propto\int\limits_{0}^{\infty}dR\
f(R)R^{4}\frac{r}{R^{2}}\exp\left(-\left(\frac{r}{R}\right)^{m}\right).$ (91)

To simplify the integral, we change variables to $\beta=1/R$, such that
$\mathrm{d}R=-\beta^{2}\mathrm{d}\beta$ and

$\displaystyle n(r)$ $\displaystyle\propto\int\limits_{0}^{\infty}d\beta\
f(\beta^{-1})\frac{r}{\beta^{4}}\exp\left(-(\beta r)^{m}\right).$ (92)

The integral in question has a solution in the case where $f=\beta^{4-\alpha}$
with $\alpha<1$. We assume that this is the case, the integral now becomes,

$\displaystyle
n(r)=C\int_{0}^{\infty}\beta^{-\alpha}r\exp\left(-(r\beta)^{m}\right)d\beta$
(93)

where $C$ is a normalization constant.

We begin with the Gamma function

$\displaystyle\Gamma(z)=\int_{0}^{\infty}t^{z-1}\exp(-t)dt.$ (94)

This integral converges for $z>0$ if $z\in\mathbf{R}$. Substituting
$t=(r\beta)^{m}$ yields $\mathrm{d}t=mr^{m}\beta^{m-1}\mathrm{d}\beta$ and
therefore

$\displaystyle\Gamma(z)$
$\displaystyle=\int_{0}^{\infty}(r\beta)^{mz-m}\exp\left(-(r\beta)^{m}\right)m(r\beta)^{m}\beta^{-1}d\beta$
(95)
$\displaystyle=mr^{mz-1}\int_{0}^{\infty}\beta^{mz-1}r\exp\left(-(r\beta)^{m}\right)d\beta$
(96)
$\displaystyle=mr^{-\alpha}C^{-1}C\int_{0}^{\infty}\beta^{-\alpha}r\exp\left(-(r\beta)^{m}\right)d\beta.$
(97)

Here, we introduce $z=(1-\alpha)/m$, which, due to the lower bound of $z>0$,
leads to an upper bound of $\alpha<1$. We recognize the integral $n(r)$ on the
right and therefore find

$\displaystyle n(r)\propto
r^{\alpha}C\,\Gamma\left(\frac{1-\alpha}{m}\right).$ (98)

If the radius of patches in the multiverse were distributed according to the
function $f(R)\propto R^{-(4-\alpha)}$, the initial growth of the joint
pairwise distribution would follow a sublinear power law. Notably, the scaling
exponent would not depend on the tail parameter $m$.

Figure 23: Integrating over the multiverse of patches of size $R$, summing up
the contributions to the pairwise-distance distribution of all patches, for
three assumptions of (i) generalized pair distribution, (ii) pair distribution
of a uniform, circular patch of buildings, and (iii) a parabola pair
distribution, with varying $\alpha$. Here, we used a minimum patch size of
$R=10\mathrm{m}$, which fixes a scale and leads to linear growth of the
pairwise-distance distribution for small distances.

Let us consider the implications of this for the area of the patches. The area
would scale as $A\propto R^{2}$ (or $R\propto\sqrt{A}$), indicating that the
distribution of the area would follow

$\displaystyle\tilde{f}(A)$
$\displaystyle=\frac{\mathrm{d}R}{\mathrm{d}A}f(R=\sqrt{A})$ (99)
$\displaystyle\propto\frac{1}{A^{(5-\alpha)/2}}.$ (100)

In the data, we observe $\alpha=0.67$. That means that the area of the patches
would have to be distributed according to a power law with exponent
$\mu=-(5-\alpha)/2\approx-2.17$, which is well within what has been found
empirically for cities [72]. Furthermore, our city radius inference analysis
(cf. Fig. 21) indicates that $f(R)\propto 1/R^{3.29}$, which leads to
$\mu=-2.15$. This demonstrates that the results of these two separate analyses
are consistent. As previously demonstrated, the sublinear scaling exponent
$\alpha$ is independent of the tail parameter $m$. We extend our analysis by
using the parabola pair distribution model Eq. 59. We have

$\displaystyle n(r)$ $\displaystyle\propto\int\limits_{0}^{\infty}dR\
n(r,R)f(R)$ (101) $\displaystyle=\int\limits_{r/2}^{\infty}dR\
f(R)R^{4}\left(-\frac{3}{4R}\left(\frac{r}{R}\right)^{2}+\frac{3r}{2R^{2}}\right).$
(102)

Here, the lower bound in the integral comes from the condition that $r\leq
2R$. Now, as above, we demand $f(R)\propto 1/R^{4-\alpha}$ to find

$\displaystyle n(r)$ $\displaystyle\propto\int\limits_{r/2}^{\infty}dR\
R^{\alpha}\left(-\frac{3}{4R}\left(\frac{r}{R}\right)^{2}+\frac{3r}{2R^{2}}\right)$
(103) $\displaystyle=\int\limits_{r/2}^{\infty}dR\
\left(-\frac{3r^{2}}{4R^{3-\alpha}}+\frac{3r}{2R^{2-\alpha}}\right)$ (104)
$\displaystyle\propto r^{\alpha}.$ (105)

We solve the respective integrals $n(r)$ for the circle, generalized pair
distribution, and parabola pair distribution models numerically and find that
the above derivation holds for all three (see Fig. 23).

###  7.7 Pair distribution of a self-similar modular hierarchical model of
building locations

One limitation of the multiverse approach is that it is only applicable if the
patches are truly independent or sufficiently separated so that the scale of
the pair distribution number density $n(r)$ does not affect the outcome.
However, it is plausible that patches may be in close proximity to each other.
This is supported by findings in [72]. Additionally, they discovered that a
collection of these patches forms a fractal, or self-similar structure. While
it is possible to demonstrate that a fractal dimension of building location
does not necessarily result in sub-linear growth of the pair distribution
function, it is certainly possible to investigate the consequences of such
self-similar patch location on it.

Figure 24: Upper two panels: Theoretical pairwise-distance distributions of
the hierarchical building placement model with $N=4$ and $\theta=0.8$, in
color the single contributions Eq. (109) of scale $k$, in black the total
joint distribution function over all scales Eq. (110), including a power-law
fit showing the resulting scaling exponent $\alpha\approx 2/3$. Lower two
panels: A sample from the hierachical placement model with $L=8$ hierarchy
layers and the corresponding resulting pairwise-distance distribution of the
whole structure, showing the desired scaling behavior on the mesoscale.

We assume a self-similar modular hierarchical structure comprising patches of
buildings. We start with a single unit of $b$ buildings, each of which has a
radius $z_{1}$. These are located within a patch of size $R_{1}$. Our
objective is to regulate the building density in such a way that the number of
buildings in a patch is given by $b=\theta(R_{1}^{2}/z_{1}^{2})$. Here,
$\theta$ represents the packing fraction. This implies that the patch radius
is computed as follows: $R_{1}=z_{1}\sqrt{b/\theta}$.

Now, consider that there are $b$ of these patches of radius $R_{1}$, located
in a larger patch of higher order, which is (self-)similar to the basic patch.
We posit that each of the lower-order patches has a radius of $z_{2}=R_{1}$,
while the higher-order patch has a radius of $R_{2}=z_{2}\sqrt{b/\theta}$.

Subsequently, we add $b-1$ similarly constructed patches to form an even
larger patch. This entails constructing a self-similar structure of patches
where the size of each patch of hierarchical order $k$ is given by

$\displaystyle R_{k}=z\left(\frac{b}{\theta}\right)^{k/2}$ (106)

and the size of each sub-patch of a patch is given by

$\displaystyle z_{k}=z_{1}\left(\frac{b}{\theta}\right)^{(k-1)/2}.$ (107)

Given a maximum number of $L$ orders (layers), the total number of houses is
eventually $b^{L}$. In each layer $k$, there are $b^{L-k}$ patches of order
$k$.

Now assume that for each patch $i$ of order $k$ and location $\bm{s}_{i,k}$,
we distribute the location of its $b$ sub-patches randomly within this patch
such that the pair distribution function of sub-patches leads to a pair
distribution $p(r,R_{k}$ of scale $R_{k}$. To prevent the patches from
overlapping excessively, a collision algorithm is employed for the $b$ sub-
patches (preventing collisions between disks with radius $z_{k}$). We repeat
this process recursively for each patch until the depth of the $N$-ary
hierarchy tree reaches $L$. The leaves of the tree represent buildings. As
these buildings will overlap, another collision algorithm is run until they no
longer overlap.

We now estimate the joint pairwise-distance distribution of the entire
structure. Consider a container patch that contains $b$ sub-patches. If we
assume that the total the $B$ buildings within a sub-patch of scale $k-1$ are
sufficiently concentrated within its center, the contribution of this patch of
radius $R_{k}$ will be proportional to the pair distribution $p(r,R_{k})$,
weighted with the total number of pairs of buildings within this container-
patch, except the pairs of buildings within the same sub-patch. Consequently,
the number of pairs therefore scales as $\propto B^{2}b(b-1)$, or, neglecting
the linear contribution, as $\propto(bB)^{2}$. For each patch in layer $k$,
there are going to be $B=b^{k-1}$ buildings in a sub-patch. Therefore, the
number of pairs of buildings in a patch of scale $k$ grows as $\propto
b^{2k}$.

Each of these patches contributes to approximately $b^{2k}$ pairs to the joint
pairwise-distance distribution. There are $b^{L-k}$ patches of order $k$ with
radius $R_{k}$. Therefore, the total contribution to the joint pairwise-
distance distribution of this scale is

$\displaystyle n_{k}(r,R_{k})$ $\displaystyle\propto b^{k}p(r,R_{k})$ (108)
$\displaystyle=b^{k}p(r,z(b/\theta)^{k/2})$ (109)