Fluid-Spatiotemporal Stochastic Geometry for Information Flow in Non-Stationary Fields: Kinematics and Field Reconstruction

1

Fluid-Spatiotemporal Stochastic Geometry for Information Flow in Non-Stationary Fields: Kinematics and Field Reconstruction

arXiv:2610.23277v1 [cs.NI] 16 Oct 2026

Wen-Yu Dong, Weiwei Jiang, Senior Member, IEEE, Song Zhao, Qi Bi, Fellow, IEEE, Sheng Chen, Life Fellow, IEEE

Abstract—Classical stochastic geometry primarily addresses forward analysis by prescribing spatial point-process laws and mobility models. It does not directly address the inverse problem of recovering directional macroscopic transport from aggregate density evolution in non-stationary spatial networks. This paper develops the kinematic foundation of fluid-spatiotemporal stochastic geometry (F-STSG), representing an evolving mobile population by an intensity field under an assumed continuum limit while retaining the network infrastructure as a discrete spatial measure. We formulate macroscopic transport reconstruction as an inverse boundary value problem. Under a strictly positive bounded intensity and compatible source and boundary data, an instantaneous weighted minimum-energy criterion selects a unique velocity representative from the admissible class. The reconstruction reduces to a weighted elliptic Neumann problem and separates density-identifiable transport from weighted-solenoidal circulation that is not identifiable from density evolution alone. We further derive the information flux, interface transport, velocity-divergence relation, and centroid dynamics. Numerical experiments examine reconstruction consistency and sensitivity, while the public New Era of Urban Traffic Monitoring with Aerial Footage (pNEUMA) dataset provides held-out validation of internal transport. A separate unmanned aerial vehicle (UAV) study evaluates the centroid as a causal downstream control reference under observation delay, motion constraints, and execution perturbations. Index Terms—Fluid-spatiotemporal stochastic geometry, nonstationary wireless networks, inverse boundary value problems, information flux, macroscopic transport.

I. I NTRODUCTION TOCHASTIC geometry (SG) has provided a rigorous framework for relating random spatial configurations of wireless nodes to coverage, interference, and scaling behavior [1]–[4]. In its classical form, SG represents the network as a spatial snapshot governed by a prescribed point-process law. Mobility-aware extensions introduce temporal evolution by specifying displacement or mobility models between successive spatial configurations [5]–[7]. Both formulations primarily support forward analysis: once the spatial law and mobility

S

W.-Y. Dong, S. Zhao and Q. Bi are with Future Technology Research Center, China Telecom Research Institute, Beijing 102209, China (E-mails: dongwy@chinatelecom.cn; zhaosong1@chinatelecom.cn; qibi@chinatelecom.cn) W. Jiang is with the School of Information and Communication Engineering, Beijing University of Posts and Telecommunications, Beijing, 100876, China (Email: jww@bupt.edu.cn) S. Chen is with the School of Electronics and Computer Science, University of Southampton, Southampton SO17 1BJ, U.K., and also with Faculty of Information Science and Technology, Ocean University of China, Qingdao 266100, China (E-mail: sqc@ecs.soton.ac.uk).

mechanism are specified, the resulting wireless performance can be evaluated. In operational networks, however, individual trajectories and exact mobility laws are often unavailable. Network observations more commonly take the form of timeindexed occupancy, traffic, or load fields aggregated over spatial regions. In non-stationary spatial networks, aggregate demand may evolve coherently across space rather than through independently resampled configurations. We refer to this structured macroscopic evolution as the “Digital Tide,” represented mathematically by a time-varying macroscopic intensity field λ. It encompasses the propagation, accumulation, depletion, and boundary exchange of spatial demand, all of which can influence handover, interference, and infrastructure load. An instantaneous density field indicates where demand is concentrated, while a sequence of such fields reveals how the spatial distribution changes over time. What remains unobserved is the directional transport that carries demand between locations. Recovering this transport from aggregate density evolution is therefore necessary to move beyond a sequence of spatial snapshots toward a macroscopic description of network dynamics. We represent this macroscopic evolution by a time-varying intensity field λ(x, t) [8]. Its evolution obeys a conservation law relating the density to the information flux J = λv, together with an internal source term and boundary exchange. Given the density evolution and source term, the conservation law determines the divergence of J , while the prescribed boundary exchange constrains its normal component. These constraints do not, in general, uniquely determine the full vector field because divergence-free components with compatible boundary behavior may remain invisible to the observed density evolution. The central theoretical problem is therefore to select a physically meaningful macroscopic transport representative from these aggregate observations. The present paper develops the kinematic foundation of fluid-spatiotemporal stochastic geometry (F-STSG). A broader preliminary formulation of F-STSG is available in the arXiv preprint [9], whereas this paper focuses on its densityto-transport inverse problem. Specifically, we formulate a weighted inverse boundary value problem [10], define the admissible class of velocity fields consistent with the estimated density evolution, specified source term, and prescribed normal boundary flux, and select its instantaneous minimum-energy element. Under explicit positivity, compatibility, boundary, and weighted-space assumptions, this selection yields a unique macroscopic velocity representative and separates the density-

identifiable transport component from microscopic circulation that cannot be recovered from aggregate density observations. Once the transport field is reconstructed, the selected velocity specifies a balance-consistent macroscopic transport representative, its information flux quantifies the corresponding signed transport across spatial regions, and its divergence characterizes compression and expansion within that representation. Together, these quantities provide a kinematic description of how spatial demand moves, accumulates, and crosses network boundaries, thereby forming a basis for downstream tasks such as handover analysis, load migration, routing, and resource placement. This paper focuses on reconstructing and characterizing these macroscopic kinematic quantities from aggregate field observations, rather than developing application-specific decision policies.

tool for spatial wireless analysis [13]–[16], [18]. Mobilityaware extensions introduce prescribed displacement, mobility, or time-varying intensity models to characterize handover, interference correlation, queueing dynamics, and temporal wireless performance [19]–[22]. F-STSG considers the inverse direction: an estimated density evolution, together with source and boundary information, is used to reconstruct a densityidentifiable macroscopic transport representative. Demandaware unmanned aerial vehicle (UAV) deployment methods [23] are relevant to the subsequent control application, where they provide application-level references and baselines. 3) Continuum and mean-field models: Continuous-media models describe dense interacting systems through macroscopic states. Representative examples include vehicular-flow models [24], [25], pedestrian-flow models [26], [27], and mean-field game (MFG) formulations for wireless caching, computation offloading, and interference management [28]– [31]. These models typically prescribe a velocity–density relation, cost functional, or agent objective and then solve a forward evolution, equilibrium, or control problem. The present formulation instead treats the density evolution, source term, and boundary exchange as given inputs and performs a fixed-time inverse selection without identifying individual trajectories or microscopic utilities. 4) Spatiotemporal estimation and forecasting: Temporal point processes, complex-network models, and data-driven methods support event modeling, traffic estimation, forecasting, and resource decisions [32]–[36]. Representative spatiotemporal graph models, including the spatiotemporal graph convolutional network (STGCN) [37] and Graph WaveNet [38], learn temporal dependencies over a discrete graph. Physical distance or spatial adjacency can be incorporated through graph construction in specific formulations, but Euclidean geometry is not necessarily retained as an explicit state variable across the broader graph and data-driven class. Their estimated or forecast density fields can serve as inputs to the conservation-constrained reconstruction, while the source and boundary terms must be obtained from corresponding measurements or separately specified models.

A. Related Works The inverse reconstruction studied in this paper is closely related to dynamic optimal transport, mobility-aware stochastic geometry, continuum and mean-field modeling, and spatiotemporal field estimation. These lines of work address different aspects of density evolution, transport, mobility, and field estimation. 1) Dynamic optimal transport: Dynamic optimal transport characterizes density evolution through a continuity constraint and a kinetic-energy objective over a density–velocity path [11], [12]. Building on this variational structure, the present formulation conditions on an estimated wireless density evolution and reconstructs, at each observation time, the weighted minimum-energy element of the admissible velocity class under specified source and boundary data. The reconstructed continuum quantities are subsequently coupled to a discrete wireless infrastructure measure. The minimum-energy selection is performed independently at each observation time and does not require the observed density trajectory to be a finitehorizon Wasserstein geodesic. 2) Mobility-aware stochastic geometry: Following the foundational works in [2], [3], SG has become a standard

TABLE I S COPE - BASED COMPARISON OF F-STSG WITH CLOSELY RELATED TECHNICAL APPROACHES Technical Dimension 1. Euclidean spatial state Uses physical coordinates, distances, or spatial measures 2. Non-stationary evolution Represents spatial states that evolve continuously over time 3. Density-to-transport inversion Infers a transport field from observed density evolution 4. No prescribed microscopic law Does not require individual trajectories or agent utilities 5. Conservation-constrained kinematics Links density evolution to velocity, flux, and divergence 6. Field–infrastructure coupling Couples a continuous demand field to discrete wireless resources Representative References

Mobility-Aware SG

Continuum Flow (MFG/LWR)

Graph and DataDriven Models

Dynamic Optimal Transport

F-STSG (This Work)

✓

✓

△

✓

✓

△

✓

✓

✓

✓

×

△

△

✓

✓

×

△

✓

✓

✓

△

✓

△

✓

✓

△

△

△

×

✓

[2], [5]

[8], [24], [28]

[34], [37], [38]

[11], [12]

This Work

The entries summarize the typical formulations and primary objectives of the cited works rather than universal capabilities of every method in each broad class. ✓: explicitly represented as a characteristic component of the class; △: present only in some representative formulations, under additional assumptions, or as a partial component; ×: outside the primary formulation of the cited representative works. MFG denotes mean-field games, LWR denotes the Lighthill–Whitham–Richards model, and F-STSG denotes fluid-spatiotemporal stochastic geometry.

2

Table I summarizes how the representative formulations treat spatial state, temporal evolution, density-to-transport reconstruction, conservation-constrained kinematics, and wireless infrastructure coupling. F-STSG integrates these elements through a wireless density-to-transport inverse formulation, a weighted admissible-space analysis, and the resulting velocity, information-flux, and divergence quantities under the stated positivity, compatibility, and boundary assumptions. The reconstructed quantities provide a kinematic interface to subsequent wireless analysis and application-specific control.

causal control reference under observation delay, motion constraints, and execution-layer nonidealities. The remainder of this paper is organized as follows. Section II develops the continuum model from the underlying user constellation and formulates the conservation law. Section III establishes the inverse-reconstruction theory and its minimum-energy solution. Section IV derives the weakvariation Green’s-function and spectral representations. Section V develops the kinematic quantities used for network analysis. Section VI evaluates finite-population convergence, matched reconstruction, sensitivity to observation and model perturbations, real-trace transport reconstruction, the derived kinematic identities, and centroid-reference UAV control under planning and execution constraints. Finally, Section VII concludes the paper. Throughout the paper, bold lower-case letters denote vectors, ∥ · ∥ denotes the Euclidean norm, and ∇, ∇·, and ∇2 denote the gradient, divergence, and Laplacian, respectively. The brackets ⟨·, ·⟩Ω and ⟨·, ·⟩∂Ω denote the volume and boundary duality pairings, respectively, and 1 denotes the constant test function equal to one. When both arguments are integrable, these pairings reduce to the corresponding volume and boundary integrals. The superscript ∗ identifies the selected minimum-energy field, whereas b· denotes an observation-based or numerical estimate. The symbol δ denotes a first-order perturbation, except in the standard Dirac notation δx and δ(·), whereas ∆ applied to a field denotes an estimatedminus-reference discrepancy. In the variational formulation, v denotes a generic admissible field; reconstructed kinematic quantities use v ∗ after Proposition 1. The indicator function is denoted by I(·). Symbols used only within one subsection are defined at their first occurrence. Table II summarizes the principal notation used throughout the paper.

B. Main Contributions This paper establishes the kinematic and inversereconstruction foundation of F-STSG for non-stationary spatial networks. Starting from the conservation of an evolving intensity field, it develops a principled route from aggregate spatial observations to macroscopic transport, network-level kinematic quantities, and constrained control references. The main contributions are as follows: • Conservation-based kinematic foundation: We formulate non-stationary network evolution in terms of the intensity field λ(x, t), the information flux J (x, t) = λ(x, t)v(x, t), internal source terms, and boundary exchange. This formulation connects the microscopic user constellation to a continuum transport description and makes explicit the inverse question at the core of F-STSG: which part of the macroscopic motion is identifiable from aggregate density evolution? • Well-posed inverse reconstruction and identifiability: We cast density-to-transport reconstruction as a weighted inverse boundary value problem. In the corresponding weighted admissible space, an instantaneous kineticenergy criterion selects a unique velocity field under the stated positivity, compatibility, and boundary assumptions. The resulting weighted elliptic problem provides both a constructive reconstruction method and a precise separation between density-identifiable transport and circulation that is invisible to the density evolution. • Analytical characterization under spatial heterogeneity: For a positive background intensity with weak spatial variation, we derive a first-order perturbation model of the weighted potential equation, together with its effective forcing, Green’s-function representation, and spectral decomposition. The modal form reveals how spatial heterogeneity modifies the reconstructed transport through the eigenstructure of the background operator and extends the analysis beyond radially symmetric fields. • Network-level kinematic quantities and empirical evaluation: From the reconstructed transport field, we derive signed interface transport, fixed-cell population balance, velocity divergence, the material-derivative relation, and centroid dynamics. Controlled experiments examine finite-population convergence and reconstruction sensitivity, while the real New Era of Urban Traffic Monitoring with Aerial Footage (pNEUMA) mobility traces provide held-out validation of internal crossing rates. A separate UAV study evaluates the centroid as a

II. S YSTEM M ODEL We introduce the F-STSG framework. Adopting a fieldtheoretic perspective, this framework models the network as a coupled system comprising two distinct components: a static discrete random measure representing the network infrastructure, and a continuous intensity field representing the hydrodynamic limit of the mobile-user population. TABLE II P RINCIPAL NOTATION USED ACROSS SECTIONS Symbol

Definition

Ω, ∂Ω, n; ΦB

Observation domain, boundary, outward normal; base-station process Mobile-user process, empirical measure, and limiting density Macroscopic intensity, its estimate, and total population Source/sink, boundary flux; update time and interval Admissible velocity, selected velocity, and scalar potential Weighted admissible set and velocity space Information flux, velocity divergence, and centroid Oriented interface, unit normal, and signed transport rate Background fields and their first-order perturbations Background weighted operator and Green’s function Eigenfunctions and eigenvalues of L0 Mean spectral efficiency and area spectral efficiency

(n) ΦM (t), µn , p

b N (t) λ, λ, S, g; tq , ∆tupd v, v ∗ , ϕ Vt , L2λ (Ω)d J, D, C Γ, nΓ , JΓ λ0 , δλ, ϕ0 , δϕ L0 , G ψn , νn R, A

3

use a normalization for which µn (Rd , t) → 1 in probability at each fixed time. We assume that, as n → ∞, µn converges weakly in probability to a deterministic probability density p(x, t): λ(x, t) w dx, (4) µn (·, t) − → p(x, t) dx ≜ N (t) where x ∈ RdR denotes spatial location and t ∈ R+ denotes time, N (t) = Rd λ(x, t) dx is the total expected population, and λ(x, t) is the resulting macroscopic mobile-node intensity field. Under homogeneous per-user demand, it is proportional to the traffic-demand density; heterogeneous demand requires a correspondingly weighted intensity.

Fig. 1. Illustration of the field-measure coupling in F-STSG.

A. Network Skeleton: Discrete Random Measures

Definition 1. The macroscopic intensity field λ(x, t) : Rd × R+ → R+ is the nonnegative continuum field associated with the limit in (4), with normalized spatial density p(x, t) = λ(x, t)/N (t). We use the term “Digital Tide” to describe the structured spatiotemporal evolution of this field.

The network infrastructure is modeled as a discrete random measure ΦB on the Euclidean domain Rd [39]: X ΦB = δXi , (1) i

where {Xi } denotes the random locations of the base stations (BSs), and δx is the Dirac measure concentrated at x. In a Ktier heterogeneous network (HetNet), ΦB is the superposition of K P independent tier-specific point-process measures, i.e., K ΦB = k=1 Φk . While the numerical validation in Section VI assumes a Poisson point process (PPP) for tractability, the core kinematic theory developed in Section III applies to general simple, stationary ergodic point processes. A mobile user at location x ∈ Rd associates with the infrastructure based on the maximum biased received power. Let K(i) ∈ {1, . . . , K} denote the tier index of the BS at Xi . We adopt the power-law path-loss model r−αk , where αk > 2 is the path-loss exponent for tier k. The association cell Vi of a generic BS Xi is defined by the weighted Voronoi tessellation: n Vi = x ∈ Rd : PK(i) WK(i) ∥x − Xi ∥−αK(i) ≥ o PK(j) WK(j) ∥x − Xj ∥−αK(j) , ∀j ̸= i , (2)

According to the displacement theorem of SG [1], [2], when points of an initial Poisson process undergo admissible independent displacements driven by a prescribed mapping, the resulting process remains Poisson with a transformed intensity. This provides one tractable microscopic realization of the macroscopic intensity evolution. The inverse field reconstruction below requires an estimated intensity and a balance law, but it does not require every operational snapshot to follow an inhomogeneous Poisson point process. In the corresponding Poisson setting, the discrete realization at time t, denoted by ΦM (t), can be modeled as an inhomogeneous Poisson point process governed by the advected intensity. For any Borel set A ⊆ Rd , the expected number of users NA (t) in A is Z E[NA (t)] = λ(x, t) dx. (5) A

The macroscopic intensity field obeys a balance law driven by a macroscopic velocity field v(x, t), a local source or sink S(x, t), and possible boundary exchange. The sourcefree closed model used for the main analytical development is the special case in which local density changes are attributable to flux divergence alone.

where Pk and Wk denote the transmit power and association bias weight for tier k, respectively. F-STSG represents the mobile-user population by a continuum field under the assumed scaling limit, while retaining the BS infrastructure as a discrete random measure. Thus, ΦB represents the discrete infrastructure, whereas λ(x, t) and J (x, t) describe the evolving traffic field and its macroscopic transport. This coupling allows transport quantities such as user flux across cell boundaries to be interpreted relative to the discrete network infrastructure.

Remark 1. The continuum approximation in (4) is asymptotic and scale dependent. Under the Poisson sampling and smoothing assumptions used in Section VI-A, relative fluctuations decrease with the effective number of observations in a smoothing support. For a bin or compact kernel support Bh (x) of scale h, a useful local diagnostic is Z Nheff (x, t) ≜ λ(y, t) dy, (6)

B. Continuum Limit and Macroscopic Intensity Field To characterize the macroscopic transport of traffic demand, (n) let {ΦM (t)}n≥1 be a sequence of mobile-user point processes indexed by the scaling parameter n. The associated normalized empirical measure is X 1 I(Yj ∈ A), A ∈ B(Rd ), (3) µn (A, t) = n (n)

Bh (x)

which is the expected number of observations contributing locally under the Poisson model. A large total population N (t) does not prevent poor reconstruction where Nheff is small. Eq. (6) is therefore a diagnostic variable, not a universal validity threshold: sparse rural-IoT regimes may require coarser bins, stronger regularization, a parametric density model, or a hybrid discrete–continuous description.

Yj ∈ΦM (t)

where Yj is the location of the j-th user, B(Rd ) denotes the Borel σ-algebra, and n is the population-scaling parameter. We 4

izationRfactor determined by the population-normalization condition Rd λ(x, t) dx = N (t). Applying the non-separability criterion, the mixed derivative yields   ∂ β̇(t) x ∇x ln λ(x, t) = − . (9) ∂t d0 + ∥x∥ ∥x∥

C. Observation of the Macroscopic Intensity Field The inverse problem uses a regularized field estimate rather than assuming direct access to a continuous λ(x, t). Let tq denote the qth field-update time. Aggregate location samples, b occupancy counts, or load telemetry can be mapped to λ(x, tq ) by a documented spatial estimator, and neighboring regud larized field estimates provide ∂ t λ(x, tq ). External ingress– egress measurements across ∂Ω may inform the boundary-flux estimate gb, whereas handover counts across internal interfaces may inform or validate internal flux integrals and do not enter the outer-boundary condition. The estimated intensity must remain positive on the modeled support through a documented domain restriction or intensity floor. Source and boundary terms must be observed, modeled, or set to zero because they are not separately identifiable from density evolution alone. Empirical validation of the continuum approximation requires mobility or occupancy records aligned with the spatial grid and reporting interval. Validation of the reconstructed transport further requires trajectory or boundary-crossing measurements that are withheld from the reconstruction. The count-based observation models and their sensitivity are evaluated in Section VI-C, while their use in constrained centroid control is evaluated in Section VI-F.

For any x ̸= 0, this vector field is non-vanishing for β̇(t) ̸= 0. Thus, a time-varying decay exponent changes the spatial shape of the intensity field. In the closed source-free model, such non-separable evolution requires nontrivial macroscopic transport. E. Dynamic State Process: Fluid-Infrastructure Coupling To represent the interaction between aggregate demand and network resources, F-STSG couples the continuous intensity field P λ(x, t) with the discrete infrastructure ΦB . Let Ψ(t) = i Ai (t)δXi , where Ai (t) ∈ {0, 1} indicates whether BS i is active at time t. Thus, Ψ(t) is the active-infrastructure random measure obtained by thinning ΦB . In a K-tier HetNet, Ψ(t) = PK Ψ (t), where Ψk (t) denotes the active thinning of tier k k=1 k. The activation of on-demand infrastructure may depend on the local intensity λ(x, t), thereby coupling the continuous demand field to the discrete infrastructure measure. Given a realization of the active set Ψ(t), a user at x associates with the serving BS X ∗ ∈ Ψ(t) providing the maximum biased received power. Consistent with the association policy in Section II-A, the serving BS is determined by

D. Topological Dynamics: Non-Separable Intensity Fields The structural properties of the intensity field λ(x, t) characterize how its spatial profile evolves over time. To distinguish shape-preserving from shape-changing evolution, we classify the field according to the coupling between its spatial geometry and temporal variation.

X ∗ = arg

max

Xi ∈supp Ψ(t)

PK(i) WK(i) ∥x − Xi ∥−αK(i) .

(10)

This dynamic association rule retains tier prioritization while allowing the active infrastructure to vary with the macroscopic state. The inverse reconstruction in Section III operates on the continuum field, while the discrete infrastructure re-enters when the reconstructed transport is interpreted relative to cell boundaries and active resources.

Definition 2. A strictly positive spatiotemporal field is separable if it admits the multiplicative factorization λ(x, t) = a(x)b(t). This implies that the spatial gradient of its temporal logarithmic derivative vanishes identically over the entire domain:   ∂ ln λ(x, t) ≡ 0, ∀(x, t) ∈ Rd × R+ . (7) ∇x ∂t Conversely, a field is non-separable if this mixed derivative is not identically zero (i.e., ∇x ∂t ln λ ̸≡ 0).

III. F IELD R ECONSTRUCTION : T HE I NVERSE B OUNDARY VALUE P ROBLEM

A separable field preserves its spatial shape and permits only a global amplitude change. In the closed source-free setting with constant total mass, this amplitude is constant; hence, nontrivial spatial redistribution requires a non-separable field. As an example, consider a time-varying power-law topology modeling a dynamic traffic profile. At an illustrative level, this profile can represent a localized traffic-demand hotspot, such as transient concentration around a transport hub, event venue, or dense Internet of Things (IoT) cluster, whose spatial footprint evolves over time: a larger β(t) corresponds to stronger central concentration, whereas a smaller β(t) produces a broader spatial distribution. The intensity is a radially symmetric field with a time-varying decay exponent β(t):

Under the continuum-limit assumption introduced in Section II, the mobile-user distribution is represented at the selected resolution by the continuous intensity field λ(x, t). While λ(x, t) is the aggregate state, the underlying driver of network evolution includes macroscopic transport, local generation or depletion, and boundary exchange. We formulate the identifiable part of these latent dynamics as an inverse problem: reconstructing the minimum-energy macroscopic velocity field from the estimated density evolution together with specified source and boundary data.

λ(x, t) = Λ0 (t)(d0 + ∥x∥)−β(t) ,

The fundamental link between topological evolution and transport dynamics is mass conservation. Let S(x, t) denote the local rate of demand generation or depletion. For a conservative network with a vanishing source term (S(x, t) = 0), the

A. Conservation Laws and Underdetermined System

(8)

where d0 > 0 ensures finiteness at the origin, β(t) > d ensures integrability over Rd , and Λ0 (t) is a time-dependent normal5

local rate of change in intensity is balanced by the divergence of the flux vector:
∂λ(x, t) + ∇ · λ(x, t)v(x, t) = 0. (11) ∂t

B. Weighted Minimum-Energy Selection and Its OptimalTransport Connection We resolve this non-uniqueness by selecting, among the admissible fields, the one with minimum instantaneous weighted kinetic energy. The objective has the same kinetic-energy form as the integrand in dynamic optimal transport, whereas uniqueness is established directly in the weighted admissible space below. 2 d For a fixed time, let R Lλ (Ω)2 be the weighted Hilbert space 2 with norm ∥v∥L2 = Ω λ∥v∥ dx. For the closed source-free λ case, define the weak admissible set  Z 2 d Vt ≜ v ∈ Lλ (Ω) : λv · ∇w dx Ω  (15) = ⟨∂t λ, w⟩Ω , ∀w ∈ H 1 (Ω) .

More generally, on a bounded observation domain Ω the balance law is ∂λ + ∇ · (λv) = S in Ω, (12) ∂t with outward normal flux g(x, t) ≜ λ(x, t)v(x, t) · n

on ∂Ω.

(13)

Integration over Ω gives the compatibility identity Z Z Z d λ dx = S dx − g dS. (14) dt Ω Ω ∂Ω The conservative closed model in (11) is recovered by setting S = 0 and g = 0.

Here, H 1 (Ω) denotes the standard first-order Sobolev space, H −1 (Ω) its dual space, and H −1/2 (∂Ω) the corresponding dual trace space used below for boundary data. For sufficiently smooth fields, this weak constraint is equivalent to (11) with zero normal flux.

Assuming homogeneous per-user traffic demand, λ(x, t) serves as the macroscopic descriptor for the spatial load. In the closed source-free setting, continuous user trajectories lead to the conservative continuity equation (11). While total load dynamics are driven by both local generation (S ̸= 0) and spatial advection (−∇·(λv)), the closed conservative model is used first to isolate transport-induced topology reconfiguration. The extension to known S and nonzero boundary flux g is given below. If either term is unknown, its contribution is confounded with advection and the inverse problem is not identifiable from λ alone.

Proposition 1. Let Ω be a bounded connected Lipschitz domain, let 0 < λmin ≤ λ(x, t) ≤ λmax < ∞ almost everywhere, and let ∂t λ ∈ H −1 (Ω) have zero mean. Among all fields in Vt , the irrotational field v ∗ = −∇ϕ is the unique minimizer of the instantaneous weighted kinetic energy Z E(t) = λ(x, t)∥v(x, t)∥2 dx. (16) Ω

Here, ϕ is the zero-mean weak solution of the Neumann problem in Theorem 1. This is a pointwise-in-time minimumenergy selection. It uses the kinetic-energy integrand appearing in the Benamou–Brenier action, but it does not assert that the complete observed density path is a global Wasserstein geodesic.

Mathematically, (11) provides only a single scalar constraint for the d-dimensional vector field v(x, t). Consequently, the inverse problem of determining v from λ is non-unique up to a weighted-solenoidal component vsol satisfying ∇·(λvsol ) = 0.

Proof. Theorem 1 establishes a weak potential ϕ, unique after fixing the zero-mean gauge. With v ∗ = −∇ϕ, the weak form of (24) gives Z λv ∗ · ∇w dx = ⟨∂t λ, w⟩Ω , ∀w ∈ H 1 (Ω), (17)

Remark 2. It is important to differentiate F-STSG from classical vehicular traffic models, such as the Lighthill–Whitham– Richards (LWR) model. First, traffic theory typically solves a forward problem, prescribing a velocity-density relation to predict density evolution. In contrast, F-STSG solves the inverse problem, treating the observable intensity evolution λ(x, t) as the input to reconstruct the latent velocity field. Second, while traffic models focus solely on fluid hydrodynamics, F-STSG couples this continuous field with the discrete network infrastructure. This coupling places boundary flux and centroid drift in the same spatial coordinate system as the discrete infrastructure, so that these quantities may be used as inputs to downstream network-control models.

Ω

so v ∗ ∈ Vt . For any other v ∈ Vt , let h = v − v ∗ . Subtracting the two weak constraints yields Z λh · ∇w dx = 0, ∀w ∈ H 1 (Ω). (18) Ω R Taking w = ϕ gives Ω λv ∗ · h dx = 0. Hence, ∥v∥2L2 = ∥v ∗ ∥2L2 + ∥h∥2L2 . λ

λ

λ

(19)

Therefore v ∗ is the unique minimum-energy element of Vt . This weighted orthogonality argument gives the pointwise-intime minimizer of the kinetic-energy integrand used in the Benamou–Brenier formulation [11]; no finite-horizon geodesic claim is required.

The underdetermination has two levels. Even when S and g are known, one may add a weighted-solenoidal component vsol satisfying ∇ · (λvsol ) = 0 and zero normal flux without changing the observed density. If S or g is unknown, source activity and transport are additionally confounded. The minimum-energy principle below resolves the first ambiguity by selecting one representative; it cannot resolve the second without supplementary observations or a source model.

Remark 3. The reconstructed field v ∗ is a macroscopic transport representative rather than a reconstruction of individual user velocities. If a finer-scale field is written as vtrue = v ∗ + vsol , where ∇ · (λvsol ) = 0 and vsol has zero 6

normal flux, then vsol does not alter the observed density evolution and is not identifiable from λ alone. It may nevertheless affect microscopic mobility or channel performance, whose quantitative evaluation requires additional trajectory or channel information.

condition obtained by testing the balance with the constant function. Pointwise-in-time scope: Theorem 1 applies separately at each time for which the stated coefficient, forcing, compatibility, and boundary assumptions hold. It does not by itself establish temporal continuity or differentiability of t 7→ ϕ(·, t) or t 7→ v ∗ (·, t). Such temporal regularity requires additional assumptions on the time dependence of λ, S, and g.

Substituting the optimal solution v ∗ = −∇ϕ into (11), we derive the governing Poisson-like PDE for the flow potential: ∂λ(x, t) = ∇ · (λ(x, t)∇ϕ(x, t)) . (20) ∂t This equation forms the theoretical core of this paper, establishing a deterministic link between the observable density evolution and the latent scalar potential. For known source and boundary data, the same minimumenergy construction leads to ∇ · (λ∇ϕ) = ∂t λ − S

in Ω,

(21)

(λ∇ϕ) · n = −g

on ∂Ω.

(22)

Corollary 1. Under the same positivity and domain assumptions as Theorem 1, let f = ∂t λ − S ∈ H −1 (Ω) and g ∈ H −1/2 (∂Ω) satisfy ⟨f, 1⟩Ω = −⟨g, 1⟩∂Ω .

Then (21)–(22) admit a weak solution unique up to an additive constant. Thus, crossings of the outer observation boundary are represented by g ̸= 0, rather than by imposing a reflecting edge. An internal cell or tracking-area interface Γ ⊂ Ω is treated separately through the flux integral in Proposition 2; only crossings of ∂Ω enter g.

The boundary condition is understood in the weak conormal sense. For this case, let VtS,g denote all v ∈ L2λ (Ω)d satisfying Z λv·∇w dx = ⟨∂t λ−S, w⟩Ω +⟨g, w⟩∂Ω , ∀w ∈ H 1 (Ω). (23)

b fb, gb) satisfy the compatibility Remark 4. Let (λ, f, g) and (λ, b ≥ λ > 0 condition, where f = ∂t λ − S, and assume λ, λ almost everywhere. Let ϕ and ϕb be the corresponding zerob mean weak solutions, and define ∆λ = λ−λ ∈ L∞ (Ω), ∆f = 2 −1/2 fb− f ∈ L (Ω), ∆g = gb − g ∈ H (∂Ω), and ∆ϕ = ϕb − ϕ. Subtracting the two weak formulations and testing with ∆ϕ gives 1 ∥∇∆ϕ∥L2 (Ω) ≤ CP ∥∆f ∥L2 (Ω) + Ctr ∥∆g∥H −1/2 (∂Ω) λ
(28) + ∥∆λ∥L∞ (Ω) ∥∇ϕ∥L2 (Ω) .

Ω

Under the mass-compatibility condition (14), the same weighted-orthogonality argument applies, so v ∗ = −∇ϕ is the unique minimum-energy element. Omitting a nonzero S or g perturbs the forcing or boundary functional but not the operator. C. Existence and Uniqueness of Scalar Potential Equation (20) constitutes a linear elliptic PDE for the scalar potential ϕ(x, t). We now establish the mathematical wellposedness of this field reconstruction.

Here, CP and Ctr are constants associated with the Poincaré and trace inequalities on Ω, respectively. Temporaldifferentiation and source-estimation errors enter through ∆f , boundary-flux errors enter through ∆g, and density-estimation errors enter through the coefficient perturbation ∆λ. Source variations that are resolved and included in S do not affect elliptic well-posedness, whereas unresolved variations may produce aliasing in ∆f . If S is unknown, source activity and transport are not identifiable from λ alone.

Theorem 1. Let Ω ⊂ Rd be a bounded connected Lipschitz domain. At a fixed time, assume λ ∈ L∞ (Ω) satisfies 0 < λmin ≤ λ ≤ λmax < ∞ almost everywhere, and let ∂t λ ∈ H −1 (Ω) satisfy ⟨∂t λ, 1⟩Ω = 0. Consider ∇ · (λ∇ϕ) = ∂t λ

in Ω,

(24)

(λ∇ϕ) · n = 0

on ∂Ω,

(25)

(27)

where the homogeneous Neumann condition is interpreted in the same weak conormal sense as above. Then there exists a weak solution ϕ ∈ H 1 (Ω) unique up to an additive constant. Consequently, the selected field v ∗ = −∇ϕ is unique in L2λ (Ω)d .

Corollary 2. In the source-free, no-flux setting of Theorem 1, for a sufficiently regular radially symmetric intensity field λ(r, t) on a disk of radius rmax , the selected radial velocity vr (r, t) at radius r admits the closed-form integral representation: Z r 1 ∂λ(ρ, t) vr (r, t) = − ρ dρ. (29) rλ(r, t) 0 ∂t

1 1 RProof. Fix the gauge by working on H⋄ (Ω) = {u ∈ H (Ω) : u dx = 0}. The weak problem is Ω Z λ∇ϕ · ∇w dx = − ⟨∂t λ, w⟩Ω , ∀w ∈ H⋄1 (Ω). (26)

Proof. In polar coordinates, the conservative PDE (24) sim∂ plifies to 1r ∂r (rλvr ) = − ∂λ ∂t . Integrating both sides with respect to ρ from 0 to r and applying the non-singular boundary condition limr→0 rλvr = 0 directly yields (29). This completes the proof.

Ω

The bilinear form is bounded and, because λ ≥ λmin > 0, coercive on H⋄1 (Ω) by the Poincaré inequality. The right-hand side is a bounded functional on the same space. The Lax– Milgram theorem therefore gives a unique zero-mean solution. Adding an arbitrary spatial constant yields the usual Neumann equivalence class, while ∇ϕ and hence v ∗ remain unique. The zero-mean condition on ∂t λ is the Neumann compatibility

Theorem 1 establishes well-posedness for the stated bounded-domain, positive-density, compatibility, and boundary assumptions; it does not cover vanishing-density regions, unknown source terms, or unrestricted microscopic velocity 7

Algorithm 1 F-STSG field reconstruction at one update time

one density, source, and boundary-flux value, respectively. The post-aggregation payload and signaling rate are Bupd , (30) Bupd = Mg bλ +MS bS +M∂ bg +Bmeta , Rsig = ∆tupd where Bupd is the number of transmitted bits per update, Bmeta is the metadata payload in bits, and Rsig is the postaggregation signaling rate in bits per unit time. The metadata may contain timestamps, region identifiers, and report-format information. The source and boundary terms are omitted when they are known to be zero or are not separately reported. Equation (30) does not include the cost of collecting raw observations from terminals, which depends on the telemetry architecture, reporting protocol, and local preprocessing. In the following complexity statements, O(·) denotes asymptotic scaling with respect to the indicated problemsize variables. Let Nobs be the number of observations processed at one update. Under constant-time assignment of each observation to a grid cell, histogram construction requires O(Nobs ) operations. More general spatial estimators may require additional estimator-specific regularization. Let Lt be the number of regularized field snapshots retained by the temporal-derivative estimator. Evaluating the derivative over all grid cells requires O(Lt Mg ) operations and O(Lt Mg ) buffer memory. Let Nnz denote the number of nonzero entries in the discretized weighted elliptic operator and let Nit denote the number of iterative-solver iterations. When each iteration is dominated by sparse matrix operations, the inverse solve requires O(Nit Nnz ) operations. Storing the sparse operator, a constant number of solver vectors, and the temporal field buffer requires O(Nnz + (Lt + 1)Mg ) memory. For a fixed local discretization stencil, Nnz = O(Mg ). The iteration count Nit nevertheless depends on the grid resolution, operator conditioning, preconditioner, and stopping tolerance. Forming the reconstructed velocity and flux and evaluating grid-based metrics such as the centroid require an additional O(Mg ) operations. A reconstruction can be completed before the next update only if the combined aggregation, estimation, temporaldifferentiation, and inverse-solver latency is smaller than ∆tupd . Prediction and constrained control introduce separate application-dependent latency. As a representative wall-clock benchmark, for the 20 × 20 grid used in the pNEUMA experiment, the inverse-reconstruction stage requires 6.733 ms per update on average, with a 95th-percentile latency of 8.715 ms on a workstation equipped with an Intel(R) Core(TM) Ultra 7 155U CPU. The statistics are obtained from 2,095 timed executions over 419 real field updates after warm-up. The measured stage starts from the prepared aggregate field and includes compatibility correction, sparse-system construction and solution, and recovery of ϕ, v ∗ , and J ; it excludes telemetry transfer, raw-data preprocessing, field estimation, and downstream control. This workstation-level timing provides a concrete computational reference alongside the asymptotic scaling analysis, but it does not establish end-to-end realtime feasibility, which additionally depends on the telemetry architecture, preprocessing and estimation pipeline, hardware platform, and update interval.

Require: Aggregated observations at tq , domain Ω, source estimate b tq ), and boundary-flux estimate gb(·, tq ) S(·, b tq ) on the Mg -cell spatial grid 1: Estimate and regularize λ(·, d 2: Estimate ∂t λ(·, tq ) from the available regularized field estimates 3: Evaluate the discrete mass-compatibility residual associated with (14) and apply any prescribed compatibility correction b tq ) under the gauge 4: Solve the weighted elliptic problem for ϕ(·, R b tq ) dx = 0 ϕ(x, Ω b tq ) and J(·, b tq )b b tq ) = λ(·, b∗ (·, tq ) = −∇ϕ(·, 5: Compute v v ∗ (·, tq ) 6: Return the reconstructed fields and required application metrics to a separate controller

components. Corollary 2 further provides the analytical basis for radially symmetric systems. Intuitively, (29) states that the radial velocity is the negative of the enclosed accumulation rate, normalized by the local boundary density. Thus, net accumulation inside radius r corresponds to inward transport (vr < 0), whereas net depletion corresponds to outward transport (vr > 0). This links the radial transport to the conservation of traffic demand. D. Analytical Properties and Physical Interpretation The selected field v ∗ (x, t), established under the assumptions of Proposition 1 and Theorem 1, is used in the subsequent F-STSG constructions. It has three relevant properties. First, Corollary 2 provides a closed-form representation for the radial setting. Second, the field is consistent with the balance law through the relation between ∂t λ and ∇ · (λv ∗ ). Third, under the stated positivity and regularity assumptions, the weak solution depends continuously on the coefficient, forcing, and boundary data, as quantified in Remark 4. These properties establish a well-defined kinematic reconstruction but do not by themselves guarantee forecasting or control performance, which also depends on observation error, latency, prediction, and the downstream controller. This selected field v ∗ (x, t) serves as the basis for the perturbation analysis in Section IV and the definition of the information flux in Section V. E. Reconstruction Workflow and Computational Requirements Algorithm 1 summarizes the reconstruction performed at one update time, from aggregate-field estimation and compatibility checking to the velocity and flux fields supplied to a separate application-specific controller. The corresponding post-aggregation signaling and computational requirements are described below. Let Mg denote the number of spatial cells used by the reconstruction grid. Let tq denote the qth field-update time, with tq+1 − tq = ∆tupd . The continuous-time notation therefore does not imply continuous transmission of individual trajectories. At each update, the reconstruction module receives spatially aggregated observations rather than individual trajectories. Consider a fixed-length uncompressed report containing one density value for each of the Mg spatial cells, MS scalar source summaries, and M∂ scalar boundary-flux summaries. Let bλ , bS , and bg denote the numbers of bits used to represent 8

IV. P ERTURBATION A NALYSIS AND G REEN ’ S F UNCTION F ORMALISM

B. Green’s Function Integral Representation The linearity of (34) enables the formulation of a general solution via the Green’s function method. Let G(x, y; t) be the generalized Green’s function for the background operator L0 = ∇ · (λ0 ∇) at time t. Under Neumann boundary conditions, the standard definition L0 G(x, y; t) = δ(x − y) is ill-posed due to the global mass conservation constraint. We fix R the additive constant by imposing the zero-mean gauge G(x, y; t) dx = 0. The Neumann Green’s function is then Ω defined by 1 , (39) ∇x · (λ0 (x, t)∇x G(x, y; t)) = δ(x − y) − |Ω|

Theorem 1 establishes the well-posed inverse problem, but spatially varying coefficients generally preclude closed-form solutions. We therefore derive a first-order approximation for a positive background field in the weak-variation regime ϵ ≪ 1. The approximation is used only within this regime; outside it, the full weighted elliptic problem is solved directly. A. Linearized Field Equation We decompose the intensity field λ(x, t) into a specified background state λ0 (x, t) and aR mass-conserving first-order perturbation δλ(x, t) satisfying Ω δλ(x, t) dx = 0: λ(x, t) = λ0 (x, t) + ϵδλ(x, t),

(31)

ϕ(x, t) = ϕ0 (x, t) + ϵδϕ(x, t),

(32)

subject to ∂n G = 0 on ∂Ω, where δ(·) denotes the Dirac delta function. Because the effective source Seff satisfies R S (y, t) dy = 0, the constant offset −1/|Ω| vanishes eff Ω upon integration. The potential perturbation therefore has the Green’s-function integral representation Z δϕ(x, t) = G(x, y; t)Seff (y, t) dy. (40)

where ϵ ≪ 1 is the perturbation parameter. The background pair (λ0 , ϕ0 ) satisfies ∂λ0 = ∇ · (λ0 ∇ϕ0 ). (33) ∂t Theorem 2. Substituting (31) and (32) into (20) and retaining the O(ϵ) terms gives ∂δλ ∇ · (λ0 ∇δϕ) = − ∇ · (δλ∇ϕ0 ). (34) ∂t Proof. Substituting the asymptotic expansions λ(x, t) = λ0 + ϵδλ and ϕ(x, t) = ϕ0 + ϵδϕ into the governing PDE (20) yields: ∂(λ0 + ϵδλ) = ∇ · ((λ0 + ϵδλ)∇(ϕ0 + ϵδϕ)) . (35) ∂t Expanding the right-hand side (RHS) of (35) and grouping by powers of ϵ, we obtain:

Ω

Consequently, substituting Seff , the velocity perturbation δv = −∇(δϕ) becomes: Z δv(x, t) = − ∇x Ω

+ O(ϵ ).

h ∂δλ(y, t) ∂t

i + ∇y · δλ(y, t)v0 (y, t) dy.

(41)

Within the linearization regime, G maps the compatible firstorder effective source to the first-order velocity perturbation. Outside the small-perturbation and positive-background assumptions, the full weighted elliptic problem must be solved directly. C. Spectral Decomposition of the Weighted Operator

RHS =∇ · (λ0 ∇ϕ0 ) + ϵ (∇ · (λ0 ∇δϕ) + ∇ · (δλ∇ϕ0 )) 2

G(x, y; t)

The following spectral decomposition is performed at a fixed observation time, and the time arguments of L0 , G, ψn , and νn are suppressed in this subsection. The background operator L0 = ∇ · (λ0 ∇) acts as a generalized Laplacian and is self-adjoint on L2 (Ω) under the stated Neumann boundary condition. According to spectral theory for elliptic operators [40], there exists a countable set of orthonormal eigenfunctions {ψn (x)} and corresponding non-negative eigenvalues, 0 = ν0 < ν1 ≤ ν2 ≤ . . . , satisfying:

(36)

The left-hand side (LHS) expands linearly as: ∂(δλ) ∂λ0 +ϵ . (37) LHS = ∂t ∂t The background state identically satisfies the zeroth-order equation ∂t λ0 = ∇ · (λ0 ∇ϕ0 ). Equating the first-order O(ϵ) coefficients from both sides yields: ∂(δλ) = ∇ · (λ0 ∇δϕ) + ∇ · (δλ∇ϕ0 ). (38) ∂t Rearranging to isolate the weighted elliptic operator ∇ · (λ0 ∇δϕ) directly produces the linearized equation (34).

∇ · (λ0 (x)∇ψn (x)) = −νn ψn (x),

(42)

subject to the Neumann boundary condition ∂n ψn = 0 on ∂Ω. The zero eigenvalue ν0 = 0 corresponds to the constant eigenfunction ψ0 = |Ω|−1/2 , which reflects the global mass conservation. By the spectral theorem for self-adjoint elliptic operators, the Green’s function on the subspace orthogonal to ψ0 admits the expansion ∞ X ψn (x)ψn (y) . (43) G(x, y) = − νn n=1

Defining the background velocity as v0 (x, t) ≜ −∇ϕ0 (x, t), we write the effective first-order forcing as Seff ≜ ∂t δλ + ∇ · (δλv0 ). This term represents the firstorder mass imbalance induced by the density perturbation. The left-hand side gives the weighted elliptic response of the potential perturbation under the background field λ0 . The differential operator L0 = ∇ · (λ0 ∇) depends on the specified positive background state, permitting the Green’s-function representation under the same boundary and compatibility conditions.

Substituting (43) into (40), the potential perturbation becomes: 9

δϕ(x, t) = −

∞ X ψn (x) n=1

νn

seff n (t),

Remark 5. The reconstructed transport field can also be coupled to a separately specified wireless-performance model. Let SINR(x, t) denote the local signal-to-interference-plusnoise ratio, and let

(44)

R where seff n (t) ≜ Ω ψn (y)Seff (y, t) dy is the projection of the effective source onto the n-th eigenmode. The expansion shows how the inverse operator weights the spatial eigenmodes of the effective forcing. Each mode contributes to δϕ with a factor 1/νn , so modes with larger νn have a smaller contribution than modes with smaller νn . This is a property of the elliptic inverse operator and does not imply a temporal decay rate. This spectral property does not by itself establish robustness to observation noise; finite-count and estimator sensitivity are evaluated separately in Section VI-C.

R(x, t) = E[log2 (1 + SINR(x, t))] denote the local mean spectral efficiency under a specified association, activity, scheduling, channel, and interference model. When λ(x, t) represents the active-link density, the local area spectral efficiency is A(x, t) = λ(x, t)R(x, t). Combining this definition with (12) gives ∂R ∂A + ∇ · (RJ ) = RS + λ + J · ∇R. ∂t ∂t Thus, RJ provides a spectral-efficiency-weighted transport flux, while R itself is determined by the separate wirelessperformance model. If λ denotes demand rather than activelink density, an additional activity or scheduling map is required. For a compatible finer-scale field vtrue = v ∗ + vsol with ∇ · (λvsol ) = 0, the omitted component contributes

V. K INEMATIC A NALYTICAL T OOLBOX Based on the selected macroscopic velocity field v ∗ (x, t), we formulate a suite of kinematic metrics. These quantities describe the instantaneous transport implied by the balance law. Their use for prediction requires causal estimation and an explicit model of subsequent field evolution. For the classical forms of the differential identities and interface integrals below, we additionally assume that λ is continuously differentiable in space and time, v ∗ is continuously differentiable in space, and the interfaces are piecewise smooth. These additional assumptions are not required for the weak existence and uniqueness result in Theorem 1.

∇ · (Rλvsol ) − λvsol · ∇R = R∇ · (λvsol ) = 0 to the area-spectral-efficiency balance for a specified R. The same circulation may nevertheless alter R through microscopic channel or interference dynamics; such an effect is not identifiable from λ alone and requires additional mobility or channel information. Hence, this relation provides a kinematic interface to, rather than a replacement for, a wirelessperformance model.

A. Information Flux Vector Generalizing the concept of microscopic node kinetics, we define the information flux as a vector field quantifying the spatial advection of the intensity field.

Related fluid-kinetic formulations have also been used for loop-free routing in ultra-dense networks [41], physics-native digital-twin synchronization in kinetic swarms [42], and timevarying interference modeling in LEO mega-constellations [43].

Definition 3. The information flux, J (x, t), is defined as the product of the local intensity and macroscopic velocity field. For reconstructed metrics, the velocity is the selected field v ∗ : J (x, t) ≜ λ(x, t)v ∗ (x, t).

B. Topological Compression and Material Derivative

(45)

To characterize the formation of demand hotspots at the macroscopic level, we analyze the compressibility of the selected flow.

The information-flux field describes transport of the modeled intensity measure and defines signed flow across an oriented interface Γ, such as a cell edge or tracking-area boundary.

Definition 4. The congestion divergence, D(x, t), is the scalar field defined by the divergence of the macroscopic velocity:

Proposition 2. Let Γ ⊂ Ω be an oriented (d − 1)-dimensional interface with chosen unit normal nΓ . The signed transport rate across Γ in the direction of nΓ is Z JΓ (t) = J (x, t) · nΓ dS. (46)

D(x, t) ≜ ∇ · v ∗ (x, t).

(48)

Proposition 3. For a conservative system (S = 0), the congestion divergence is mathematically equivalent to the negative logarithmic rate of change of the intensity, as observed in the Lagrangian frame of reference: 1 Dλ(x, t) D(x, t) = − , (49) λ(x, t) Dt

Γ

Proof. The result follows by integrating the normal component of J over the oriented interface Γ. The flux integral JΓ (t) characterizes instantaneous transport across an internal interface. For a fixed association partition {Vi }, integrating the continuity equation over Vi ∩ Ω gives Z Z X d λ dx = S dx − Jij − Ji,ext , (47) dt Vi ∩Ω Vi ∩Ω

D ∂ ≜ ∂t + v ∗ · ∇ is the material (or Lagrangian) where Dt derivative, capturing the rate of change experienced by a moving fluid element.

Proof. Expanding the divergence term in the conservative continuity equation (11) using the scalar-vector product identity yields: ∂λ + v ∗ · ∇λ + λ(∇ · v ∗ ) = 0. (50) ∂t

j∈Ni

where Ni denotes the neighboring regions, Jij is the signed flux from Vi to Vj defined by Proposition 2, and Ji,ext is the outward flux through ∂Vi ∩ ∂Ω. 10

∂λ Applying the definition of the material derivative Dλ Dt ≜ ∂t + ∗ ∗ v ·∇λ and the congestion divergence D ≜ ∇·v , the equation simplifies to: Dλ + λD = 0. (51) Dt Assuming λ > 0 and solving for D directly yields (49).

10-1

10-2

For the general balance law (12), the same expansion gives S(x, t) 1 Dλ(x, t) D(x, t) = − . (52) λ(x, t) λ(x, t) Dt Thus, D < 0 identifies local compression of the selected macroscopic transport field, whereas density growth may also arise from S > 0. Forecasting or early-warning use additionally requires causal field estimation and an evolution or service model.

10-3

103

Definition 5. The network centroid, C(t), is the instantaneous expectation of the user location vector: Z 1 xλ(x, t) dx. (53) C(t) ≜ Ex∼p(·,t) [x] = N (t) Rd

107

∂Ω

Thus, in the source-free closed case (S = 0 and g = 0), ĊΩ reduces to the intensity-weighted mean velocity. More generally, source activity and boundary exchange modify the centroid dynamics through their first-moment contributions, although these contributions may vanish or cancel in special cases.

Proposition 4. In a source-free conservative system (S = 0) with sufficiently fast flux decay at infinity, the velocity of the centroid, VC (t) ≜ dC dt , is exactly the intensity-weighted spatial average of the macroscopic velocity field: Z VC (t) = Ex∼p(·,t) [v ∗ (x, t)] = v ∗ (x, t)p(x, t) dx. (54)

VI. N UMERICAL AND E MPIRICAL VALIDATION The experiments in this section serve different evidentiary purposes, summarized in Table III. The matched tests are used only for convergence, implementation, identifiability, or analytical-consistency checks; controlled perturbations evaluate sensitivity; the pNEUMA experiment provides held-out real-trace validation of macroscopic transport; and the UAV study evaluates a downstream control use case rather than reconstruction fidelity.

Rd

R

xλ(x, t) dx and Proof. Differentiating C(t) = N (t) Rd using dN (t)/dt = 0 gives Z ∂λ(x, t) N (t)VC (t) = x dx. (55) ∂t Rd Under the source-free balance law, substituting ∂t λ = −∇ · (λv ∗ ) gives Z N (t)VC (t) = − x∇ · (λv ∗ ) dx. (56)

A. Finite-Population Convergence

Rd

Rd

106

R NΩ (t)−1 Ω xλ(x, t) dx. Differentiating NΩ (t)CΩ (t) and using the mass-compatibility identity (14) gives Z Z
NΩ (t)ĊΩ (t) = λv ∗ dx + x − CΩ (t) S dx ΩZ Ω
− x − CΩ (t) g dS. (59)

To track the global trajectory of the demand, we define the first moment of the intensity field.

For the i-th component, integration by parts gives Z Z xi ∇ · (λv ∗ ) dx = − ∇xi · (λv ∗ ) dx,

105

Fig. 2. Finite-population convergence under the matched generator. Blue markers show the simulated mean relative density-estimation error E(Npop ), the shaded band indicates ±1 standard deviation across independent realiza−1/2 tions, and the red dashed line shows the theoretical O(Npop ) scaling.

C. Global Centroid Dynamics

−1

104

Section II-B introduces the weak-convergence limit of the normalized empirical measure µn as n → ∞. We examine its finite-population behavior using the Monte Carlo population size Npop and the relative error

(57)

Rd

where the boundary term vanishes when the flux decays sufficiently rapidly at infinity. Since ∇xi = ei , where ei is the ith Cartesian basis vector, Z N (t)VC (t) = λ(x, t)v ∗ (x, t) dx. (58)

bN (x, t) − λ(x, t) λ pop E(Npop ) ≜

∥λ(x, t)∥L2 (Ω)

L2 (Ω)

.

(60)

Here, λ(x, t) is the theoretical macroscopic intensity, and bN (x, t) is the empirical density field reconstructed from a λ pop Monte Carlo realization containing Npop users. We use a time-varying Gaussian hotspot to generate a nonseparable intensity field. To separate sampling error from gridintegration error, the theoretical intensity λ(x, t) is evaluated

Rd

Dividing by N (t) yields (54). For the general balance law R on a bounded observation domain, define NΩ (t) = Ω λ(x, t) dx and CΩ (t) = 11

TABLE III E VIDENTIARY ROLE AND SCOPE OF THE EXPERIMENTS IN S ECTION VI. Section VI-A

Evidence type Matched finite-population convergence

VI-B

Implementation and identifiability check

VI-C

Sensitivity and robustness analysis

VI-D

Real-trace empirical validation

VI-E

Matched analytical-consistency check

VI-F

Downstream control evaluation

Purpose Check convergence of the selected density estimator under the prescribed generator and smoothing rule. Verify minimum-energy selection and separation of the density-identifiable and weighted-solenoidal components. Isolate finite-count, temporal-estimation, missing-data, source, boundary, and operator perturbations. Compare reconstructed internal transport with held-out trajectory crossings in pNEUMA. Check the derived scalar kinematic identities and illustrate their temporal interpretation. Evaluate the centroid as a causal reference under observation delay and motion constraints.

by the cumulative-distribution-function difference method. For (N ) each realization m ∈ {1, . . . , M }, a constellation Φm pop is generated by drawing Npop coordinates from the instantaneous bN density profile. The empirical field λ pop is then estimated using the same smoothing rule for all population sizes. The results are averaged over M = 105 independent realizations for each Npop . Fig. 2 shows that the relative error follows an approximately −1/2 Npop trend over the tested range, while the snapshot variability decreases with population size. At the largest simulated population, Npop ≈ 107 , the error is below 10−3 for the selected estimator and Poisson generator. These finiterange results are consistent with the expected asymptotic convergence under the matched generator and smoothing rule, but they do not provide empirical validation for operational mobility traces. For a measured network, the adequacy of the continuum approximation depends on the spatial aggregation scale, reporting interval, local effective observation count, and estimator bias and must be assessed on the target dataset.

Evidentiary limit Does not establish reconstruction accuracy on operational mobility data. The 0.092% normalized root-mean-square error (NRMSE) is not an unseen-data accuracy result. Uses manufactured reference fields and does not measure empirical fidelity. Evaluates macroscopic transport, not microscopic vehicle velocities. Does not provide independent reconstruction-accuracy evidence. Does not validate the inverse reconstruction itself.

minimum-energy admissible transport component and excludes weighted-solenoidal modes that do not alter the observed density evolution. To verify this property and the resolution of the inverse problem’s non-uniqueness, we construct a reference velocity field vtrue as a superposition of a compressive component and a weighted-solenoidal component: vtrue = virr + vrot .

(61)

Here, virr is a radial gradient field that drives the prescribed density evolution. The injected rotational component vrot is constructed to satisfy ∇ · (λvrot ) = 0. In the present construction, this weighted-solenoidal condition follows from ∇ · vrot = 0 and vrot · ∇λ = 0. The rotational component contributes to the local flux and kinetic energy but does not change the density evolution. The inverse solver receives only the temporal derivative ∂t λ generated by vtrue . Fig. 3(a) shows the reference field, which contains both the radial transport component and a rotational component, while Fig. 3(b) shows that the reconstruction retains the density-identifiable radial component. Fig. 3(c) further shows agreement with this identifiable component, with an NRMSE of approximately 0.092%, where the NRMSE is ∥v ∗ − virr ∥L2 /∥virr ∥L2 . Because the reference field and inverse solver use matched numerical settings, this small error is interpreted only as an implementation and identifiability check rather than reconstruction accuracy on unseen data. Finally, Fig. 3(d) shows the remaining weighted-solenoidal component, which is invisible to the density evolution but may still affect microscopic mobility or user-level performance.

B. Matched Minimum-Energy Reconstruction Complementing the stochastic convergence analysis in Section VI-A, this subsection checks the inverse solver using deterministic manufactured continuum fields. Proposition 1 shows that the F-STSG reconstruction selects the

C. Sensitivity to Observation, Estimation, and Model Perturbations We next examine the robustness of the reconstruction pipeline under controlled perturbations using manufactured fields on a 96×96 grid and 30 random seeds. These controlled tests complement the real-trace evaluation in Section VI-D by providing known reference fields and allowing individual error sources to be isolated. The references provide λ, ∂t λ, vref , and Jref = λvref , while the reconstruction pipeline produces b ∂ d bv ∗ . b∗ , and Jb = λb λ, t λ, v Three aggregate-count observation models are considered. In O1, each grid-bin count is sampled independently from a Poisson distribution, so the aggregate count varies across snapshots. In O2, each independently sampled Poisson snap-

0.1

0.05

0

-0.05 NRMSE = 0.092% -0.1

Fig. 3. Matched inverse-problem check. (a) Reference velocity vtrue , formed by expansion and rotation. (b) Reconstructed minimum-energy field v ∗ . (c) Relative reconstruction-error map. (d) Residual vtrue − v ∗ , corresponding to the injected density-invisible circulation.

12

106

102

35 16

102 105

30

104

25

15 14 13 12 0.9

103

4800

5000

1

1.1

1.2 104

5200

10

1

102 101 101

103

104

105

106

100 103

104

105

106

100 103

104

105

106

Fig. 4. Sensitivity to aggregate-count observations. O1 uses independently sampled grid-bin Poisson counts, O2 rescales each O1 snapshot to the prescribed aggregate total, and O3 uses fixed-total multinomial counts. Solid and dashed curves show the median and 90th percentile, respectively, over 30 seeds.

e ∂t λ =

ev,λ =

∂t λ L2 (Ω) b∗ − vref L2 (Ω) v λ

vref L2 (Ω)

.

Local y (m)

(62a)

Jb − Jref L2 (Ω) Jref L2 (Ω)

.

800

600

600

400

400

200

0

t = 365 s

(b)

0.002

200

0

-200

-200

-400

-400

-600

0.001

-600 0

0

.

200

400

Local x (m)

(62b)

600

0

200

400

600

Local x (m)

Fig. 5. Real-trace reconstruction using the pNEUMA mobility dataset. (a) Observation domain and three fixed internal interfaces. (b) Aggregate vehicle density and reconstructed information flux at a representative time.

λ

eJ =

(a)

Estimated density (vehicles/m2 )

d ∂ t λ − ∂t λ L2 (Ω)

800

Local y (m)

shot is rescaled to the prescribed aggregate total. In O3, the grid-bin counts are jointly sampled from a multinomial distribution with a fixed total count. Model O1 is retained as an independent-snapshot stress test. When an O1 realization violates the closed-domain compatibility condition, the spatial mean of the estimated forcing is removed before the Neumann solve. The reconstruction errors are measured relative to the corresponding manufactured reference fields as

(62c) becomes too small, a behavior that is also observed in the real-trace experiment below.

Here, e∂t λ measures temporal-derivative estimation error. The velocity metric ev,λ uses the intensity-weighted norm of the minimum-energy formulation and gives greater weight to errors in high-intensity regions. The flux metric eJ includes b and v b∗ . All cases included in the reported errors in both λ relative-error summaries have nonzero reference norms. Medians and 90th percentiles are reported because the finite-count error distributions are skewed. With exact fields, exact temporal derivatives, and matched discrete operators, the reconstruction errors approach numerical precision, which serves only as an algebraic implementation check. Once finite-count observations are introduced, Fig. 4 shows that increasing the aggregate count reduces the reconstruction error, while temporal differentiation remains the dominant error source in the tested observation settings. The isolated source, boundary, and operator perturbations produce smaller errors than the finite-count observation cases over the tested ranges. In contrast, combining low observation counts, missing grid bins, source and boundary perturbations, and discretization mismatch can produce velocity and flux errors larger than their corresponding reference norms. These results identify a tradeoff between spatiotemporal resolution and observation noise. Finer updates do not necessarily improve reconstruction when the available count per update

D. Real-Trace Validation With pNEUMA Mobility Data We next test F-STSG on the public pNEUMA urbanmobility dataset [44]. Three sessions recorded on October 29, October 30, and November 1, 2018 are used, comprising 23,494 vehicle trajectories and more than 7.3 × 107 recorded positions. The trajectories are binned into spatial counts to esb timate λ(x, t). Thus, the reconstruction uses only the evolution of the aggregate density together with measured entries and exits at the outer boundary; individual vehicle velocities and internal crossing records are not used. Unless otherwise stated, a 20 × 20 grid and a 5-s update interval are used. A spatially uniform compatibility correction is applied to the estimated forcing before the Neumann solve so that the discrete massbalance condition is satisfied. Fig. 5(a) shows the observation region and the three predefined internal interfaces Γ1 , Γ2 , and Γ3 . Panel (b) gives a representative density snapshot together with the reconstructed information flux. The interfaces are fixed before any performance metric is calculated. To evaluate the reconstructed flux, the original high-rate trajectories are used only to count vehicles crossing these interfaces. For an 13

Rate (vehicles/s) Rate (vehicles/s) Rate (vehicles/s)

4

Measured crossing

F-STSG

2 0 -2 -4 100

200

300

400

500

600

100

200

300

400

500

600

100

200

300

400

500

600

4 2 0 -2 -4

4 2 0 -2

Fig. 7. Normalized center density, signed outward boundary flux, and center compression indicator in the selected periodic example. Peak compression and peak inward transport precede the center-density maximum, and both indicators cross zero at the density peak.

-4

Session time (s)

Fig. 6. Real-trace validation of reconstructed internal crossing rates on the pNEUMA dataset. Panels (a)–(c) compare the trajectory-counted and reconstructed signed crossing rates at Γ1 , Γ2 , and Γ3 , respectively.

deployment-level system performance.

TABLE IV AVERAGE REAL - TRACE RESULTS OVER THREE P NEUMA SESSIONS AND THREE INTERNAL INTERFACES . MAE AND RMSE ARE IN VEHICLES / S .

E. Validation and Temporal Interpretation of Kinematic Metrics

Method Pearson Spearman MAE RMSE F-STSG reconstruction 0.511 0.515 0.766 0.965 Centroid baseline -0.028 -0.040 0.877 1.076

We use a matched periodic example to verify the scalar kinematic identities and illustrate their temporal interpretation. The velocity divergence computed from ∇ · v ∗ agrees with the material-derivative relation in (49), while the centroid velocity obtained from the generated trajectory agrees with the intensity-weighted mean velocity in (54). Figure 7 illustrates the temporal relation among the kinematic quantities R for the selected periodic radial example. Let J∂Br0 (t) = ∂Br J · n dS denote the signed outward flux 0 across the boundary of a disk of radius r0 , so that J∂Br0 < 0 represents net inward transport and J∂Br0 > 0 represents outward transport. Radial symmetry gives v ∗ (0, t) = 0, and (49) therefore relates the center compression directly to the local density growth rate. Peak compression and peak inward transport precede the center-density maximum at t = 12 h, and both indicators cross zero at the density peak. The exact lead times are specific to this periodic example.

F-STSG: fluid-spatiotemporal stochastic geometry; MAE: mean absolute error; RMSE: root-mean-square error.

interface Γ, the signed crossing rate is JΓtrace (t) =

NΓ+ (t) − NΓ− (t) . ∆tupd

(63)

where NΓ+ (t) and NΓ− (t) are the numbers of crossings in the two opposite directions. The corresponding F-STSG crossing rate is obtained by integrating the reconstructed information flux across Γ according to Proposition 2. Fig. 6 shows that the reconstructed rates capture part of the temporal variation observed in the trajectory counts, while noticeable errors remain at some time instants. As a simple reference, we also consider a centroid-translation baseline in which the whole density field moves with its centroid velocity. This baseline retains only the global drift and ignores spatially varying transport within the observation region. Table IV reports the averages over the three sessions and three interfaces. The F-STSG reconstruction gives Pearson and Spearman correlations of 0.511 and 0.515, whereas the centroid baseline gives −0.028 and −0.040. The MAE decreases from 0.877 to 0.766 vehicles/s, and the RMSE decreases from 1.076 to 0.965 vehicles/s. Averaged over the three sessions and internal interfaces, F-STSG achieves higher temporal correlation and lower MAE and RMSE than the centroidtranslation baseline. Additional tests with different grid resolutions and update intervals show that very fine temporal sampling can become less reliable when too few vehicles are observed in each update. For the default 20 × 20 grid, the inverse-reconstruction stage takes 6.733 ms per update on average, with a 95th-percentile latency of 8.715 ms on an Intel(R) Core(TM) Ultra 7 155U CPU. These timings exclude telemetry transfer, raw-data preprocessing, field estimation, and downstream control. The experiment therefore evaluates the recovery of macroscopic transport from aggregate density evolution; it does not assess microscopic vehicle velocities or

F. Constrained Centroid-Reference UAV Control 1) Centroid Reference and Controller: Unlike our prior UAV-logistics application, which used information flux for workload prediction and proactive infrastructure provisioning [45], the demonstration here uses the intensity-weighted centroid as a causal control reference. Let p(x, t) = λ(x, t)/N (t) be the normalized demand density, and let C(t) denote its intensity-weighted centroid. Assuming that p(x, t) has a finite second moment, the instantaneous mean squared service distance from a mobile resource at P satisfies Z Q(P , t) ≜ ∥x − P ∥2 p(x, t) dx (64) Rd = Q(C(t), t) + ∥P − C(t)∥2 . R The second equality follows from Rd (x − C(t))p(x, t) dx = 0. Hence, C(t) uniquely minimizes this instantaneous unconstrained squared-distance objective and is used below as a control reference. 14

At control update k, let τ ≥ 0 denote the observation delay and assume that it is an integer multiple of the control interval ∆tctrl . Define dτ ≜ τ /∆tctrl ∈ N0 and tk = k∆tctrl . The bk+h|k−d denotes the centroid predicted for update quantity C τ k +h using only observations available through update k −dτ . Over a horizon of H steps, let Ph and uh denote the predicted UAV position and velocity command, respectively. Let Fk ⊂ Rd be the admissible position set used at update k, including the flight-region bounds and the local obstacle exclusion. The predicted positions and velocity commands are obtained from H X

min

H−1 {Ph }H h=0 ,{uh }h=0

2

bk+h|k−d Ph − C τ

+η

h=1

H−1 X

the grid cell with the largest count. We compare delayed density-peak control (C1), delayed reactive-centroid control (C2), causal predictive-centroid control (C3), and a clairvoyant constrained reference (C4). C1 repeats the latest available density peak over the prediction horizon, C2 repeats the latest available centroid, and C3 extrapolates the centroid by a least-squares linear fit using up to the latest five available estimates. C4 uses the future true centroid and is included only as a noncausal reference. All controllers use the same initial condition, motion constraints, objective weight, and solver. The main setting uses H = 5, dτ = 4, vmax = 0.0500, amax = 0.0180, and η = 0.0350 in normalized simulation units. Let Nctrl denote the number of control updates in one run. For run r and controller c, let pr (x, t) denote the true demand c density used for evaluation, and let Pr,k and ucr,k denote the realized UAV position and applied velocity command at update k, respectively. For controller c ∈ {C1, C2, C3} and run r, define the cumulative service-distance objective and control effort as N ctrl Z X 2 c Lcsvc,r = ∆tctrl x − Pr,k pr (x, tk ) dx,

∥uh ∥2

h=0

(65) s.t.

P0 = P (tk ), Ph+1 = Ph + ∆tctrl uh , h = 0, . . . , H −1, ∥uh ∥ ≤ vmax ,

h = 0, . . . , H − 1,

∥uh −uh−1∥ ≤ amax ∆tctrl , h = 0, . . . , H −1, Ph ∈ Fk ,

h = 1, . . . , H.

Here, u−1 is the velocity applied immediately before update k, η ≥ 0 balances the tracking and control-effort terms, and vmax and amax are the speed and acceleration limits. Only the first command u0 is applied, after which the problem is solved again at the next update. The numerical quadratic program uses conservative polygonal approximations of the speed and acceleration constraints and a local linear approximation of the obstacle boundary. The optimization in (65) is a planning-level kinematic reference generator: uh is a commanded velocity rather than a rotor-level control input. The constraints on vmax and amax limit the planned motion but do not model command-response dynamics, execution-layer speed or acceleration saturation, or aerodynamic disturbances. These effects are examined separately below. 2) Evaluation Setup and Results: The closed-loop evaluation uses three demand-evolution scenarios representing hotspot translation, turning, and abrupt migration. At each update, a fixed total of Nobs = 50000 counts is sampled from the normalized demand field on a 25 × 25 spatial grid. The empirical centroid is the count-weighted mean of the grid coordinates, and the empirical density peak is the center of

Lcctrl,r = ∆tctrl

d k=1 R Nctrl X−1

(66) 2 ucr,k .

k=0

The first quantity is evaluated by the corresponding spatial-grid sum in the simulation. The true demand field pr is used only to evaluate all controllers against the same noise-free service objective and is not supplied to C1–C3 during control. The second quantity measures the accumulated squared velocity commands. For Fig. 8, all service-distance values are normalized by the mean C2 value at zero delay. The percentages reported in the text are instead obtained from the paired run-level gain. For Npair scenario-and-seed-matched runs, the gain of C3 relative to C2 is Npair C3 1 X LC2 svc,r − Lsvc,r G 3,2 = . (67) Npair r=1 LC2 svc,r The 95% confidence intervals use the normal approximation across the corresponding run-level quantities. For (67), the

(b) Control effort

(a) Service performance 1.8

(c) Tracking lag

1.6 1.4 1.2 1

10 0.04

Tracking lag (steps)

Cumulative control effort

0.045

0.035 0.03 0.025

8 6 4

0.02 2 0

1

2

4

6

8 C1 delayed peak

0

1

2

4

C2 reactive centroid

6

8

0

1

2

4

6

8

C3 predictive centroid

Fig. 8. Sensitivity to observation delay: (a) normalized cumulative service distance, (b) cumulative control effort, and (c) tracking lag. Markers indicate the mean, and shaded bands represent pointwise normal-approximation 95% confidence intervals across 30 scenario-and-seed-matched runs at each observation delay. Panel (a) shows normalized group means, whereas the gains reported in the text are paired run-level averages. C4 is omitted because it uses future information and serves only as a clairvoyant reference.

15

interval is computed from the paired gains rather than from the difference between group means. For controller c in run r, the tracking lag is the nonnegative temporal shift ℓ ∈ {0, . . . , Lmax } that minimizes the RMSE between the realized UAV positions and the delayed truecentroid trajectory, where Lmax = min(12, Nctrl − 3). Ties are resolved using the smallest shift. The main comparison uses 30 observation-count seeds in each of the three scenarios, giving 90 paired runs. C3 achieves a 12.80% mean paired service-distance gain over C2 (95% confidence interval: 11.52%–14.08%) and a 16.47% gain over C1. The C3-versus-C2 win/tie/loss counts are 90/0/0. The mean tracking lag is 3.33 steps for C3 and 7.00 steps for C2. This improvement is accompanied by a 9.59% increase in pooled control effort. Under nominal direct execution, no solver failure or speed, acceleration, obstacle, or flight-region violation is observed. Execution nonidealities are evaluated separately below. The delay sweep uses 10 seeds in each scenario, giving 30 paired runs at every delay. As shown in Fig. 8, the mean paired C3-over-C2 gain increases from 1.29% at zero delay to 21.94% at an eight-step delay, together with higher control effort. In the abrupt-migration scenario, the gain decreases to 5.81%, indicating that linear centroid extrapolation is less effective when the demand trajectory changes abruptly. The study covers aggregate-count centroid estimation, causal centroid prediction, and constrained UAV tracking, but not the complete PDE-reconstruction-to-control chain. 3) Execution-Layer Stress Test: We further examine the planning-to-execution gap by passing the commanded velocity uk through a surrogate UAV execution layer. Let wk denote the executed velocity. The command-response, acceleration, and position updates are modeled as uk − wk , acmd = k τa
cmd aexe + dk , k = sataphys ak (68)

TABLE V E XECUTION - LAYER STRESS TEST FOR C3. G AIN IS THE PAIRED SERVICE - DISTANCE GAIN OVER C2. S PEED SATURATION , ACCELERATION SATURATION , AND TRACKING FAILURE ARE REPORTED AS PERCENTAGES OF CONTROL UPDATES ; OBSTACLE VIOLATIONS ARE REPORTED AS UPDATE COUNTS . Case RMSE Gain Speed sat. Accel. sat. Track fail. Obstacle Nominal 0.0834 12.80% 0 0 6.11% 0 Lag 0.0936 12.96% 0 0 7.19% 0 Saturation 0.0897 12.46% 8.48% 20.61% 6.67% 0 Wind 0.0842 12.81% 2.83% 7.37% 6.11% 0 Combined 0.0958 12.68% 8.52% 5.07% 7.74% 24

Table V shows that execution nonidealities produce measurable tracking degradation. A response time of τa = 2∆tctrl increases the C3 tracking RMSE by 12.18% relative to nominal execution, while no saturation or constraint violation occurs in the lag-only case. Under tightened physical limits, the speed and acceleration bounds are activated in 8.48% and 20.61% of the C3 updates, respectively. Wind alone increases the C3 tracking RMSE by 0.89% under the tested disturbance level. Under the combined condition, the paired C3-over-C2 service gain remains 12.68%, compared with 12.80% nominally, but 24 obstacle-violation updates occur in 13 of the 90 C3 runs. In the abrupt-migration scenario, the tracking-failure update rate is 14.56%, and the C3-over-C2 gain is 4.44%. No flightregion violation or C3 solver failure is observed. These results show that execution mismatch can produce effective speed or acceleration saturation, tracking degradation, and local feasibility violations under the tested surrogate model; they do not constitute six-degree-of-freedom flight-dynamics or hardware validation. VII. C ONCLUSION This paper established the kinematic foundation of FSTSG for non-stationary wireless networks. Based on population conservation, macroscopic transport reconstruction was formulated as a weighted inverse boundary value problem. Under the stated positivity and compatibility conditions, the minimum-energy principle selects a unique admissible velocity representative, from which information flux, divergence, and intensity-weighted centroid are derived. The formulation also establishes the identifiability boundary of density-invisible circulation. Numerical experiments examined reconstruction consistency and sensitivity, while real pNEUMA mobility traces provided empirical validation of internal transport against heldout crossing measurements. A separate UAV study evaluated centroid-based control under observation delays and execution perturbations, revealing both performance gains and practical limitations. Future work will extend the kinematic framework to wireless interference and performance analysis, investigate causal prediction of evolving demand fields, and evaluate end-to-end implementations in operational cellular and UAV systems.

wk+1 = satvphys (wk + ∆tctrl aexe k ), Pk+1 = Pk + ∆tctrl wk+1 , where τa represents the finite command-response time, vphys and aphys impose effective speed and acceleration saturation, and dk is a bounded temporally correlated wind-acceleration disturbance. For a vector z, let satc (z) ≜ z min{1, c/∥z∥} for z ̸= 0, with satc (0) = 0, so that the Euclidean norm is limited to c. This surrogate model is intended to test execution mismatch rather than reproduce full six-degree-of-freedom quadrotor dynamics. The same three demand scenarios, 30 seeds, controller parameters, and initial conditions are retained. Five execution conditions are considered: nominal direct execution; a response lag with τa = 2∆tctrl = 2 s; tightened physical limits vphys = 0.8vmax and aphys = 0.8amax ; bounded correlated wind with root-mean-square magnitude 0.1amax , maximum magnitude 0.2amax , and correlation time 3∆tctrl ; and a combined condition containing all three effects. Paired controllers use the same wind realization. A tracking-failure update is recorded when the centroid-tracking distance exceeds the threshold 0.16.

R EFERENCES [1] S. N. Chiu, D. Stoyan, W. S. Kendall, and J. Mecke, Eds., Stochastic Geometry and Its Applications, 3rd ed. Hoboken, NJ, USA: Wiley, 2013.

16

[2] F. Baccelli and B. Błaszczyszyn, “Stochastic geometry and wireless networks: Volume I theory,” Foundations and Trends in Networking, vol. 3, nos. 3–4, pp. 249–449, 2009. [3] M. Haenggi, Stochastic Geometry for Wireless Networks. Cambridge, U.K.: Cambridge Univ. Press, 2012. [4] P. Gupta and P. R. Kumar, “The capacity of wireless networks,” IEEE Trans. Inf. Theory, vol. 46, no. 2, pp. 388–404, Mar. 2000. [5] M. Haenggi, “The local delay in Poisson networks,” IEEE Trans. Inf. Theory, vol. 59, no. 3, pp. 1788–1802, Mar. 2013. [6] S. Krishnan and H. S. Dhillon, “Spatio-temporal interference correlation and joint coverage in cellular networks,” IEEE Trans. Wireless Commun., vol. 16, no. 9, pp. 5659–5672, Sep. 2017. [7] H. Tabassum, M. Salehi, and E. Hossain, “Fundamentals of mobilityaware performance characterization of cellular networks: A tutorial,” IEEE Commun. Surveys Tuts., vol. 21, no. 3, pp. 2288–2308, 3rd Quart., 2019. [8] J. M. Lasry and P. L. Lions, “Mean field games,” Jpn. J. Math., vol. 2, pp. 229–260, Jan. 2007. [9] W.-Y. Dong, W. Jiang, S. Zhao, Q. Bi, and S. Chen, “Fluidspatiotemporal stochastic geometry: Information flow in non-stationary fields,” arXiv preprint arXiv:2607.00616, Jul. 2026. [10] V. Isakov, Inverse Problems for Partial Differential Equations, 2nd ed. New York, NY, USA: Springer, 2006. [11] J.-D. Benamou and Y. Brenier, “A computational fluid mechanics solution to the Monge-Kantorovich mass transfer problem,” Numerische Math., vol. 84, no. 3, pp. 375–393, Jan. 2000. [12] C. Villani, Topics in Optimal Transportation. Providence, RI, USA: Amer. Math. Soc., 2003. [13] W.-Y. Dong, S. Yang, and S. Chen, “Uplink performance analysis of heterogeneous non-terrestrial networks in harsh environments: A novel stochastic geometry model,” IEEE Trans. Commun., vol. 73, no. 8, pp. 6734–6747, Aug. 2025. [14] W.-Y. Dong, S. Yang, P. Zhang, and S. Chen, “Stochastic geometry based modeling and analysis of uplink cooperative satellite-aerial-terrestrial networks for nomadic communications with weak satellite coverage,” IEEE J. Sel. Areas Commun., vol. 42, no. 12, pp. 3428–3444, Dec. 2024. [15] W.-Y. Dong, S. Yang, P. Zhang, and S. Chen, “Modeling and performance analysis of IoT-over-LEO satellite systems under realistic operational constraints: A stochastic geometry approach,” IEEE Internet Things J., vol. 12, no. 15, pp. 30576–30593, Aug. 2025. [16] W.-Y. Dong et al., “Proper sea surface roughness enhances the performance of near-shore maritime networks,” IEEE Trans. Commun., vol. 74, pp. 13514–13530, Sep. 2026. [17] W.-Y. Dong et al., “Outage probability analysis of uplink heterogeneous non-terrestrial networks: A novel stochastic geometry model,” in Proc. IEEE GLOBECOM, Cape Town, South Africa, Dec. 8–12, 2024, pp. 2588–2593. [18] W.-Y. Dong et al., “Coverage probability density: A spatially-resolved performance analysis for non-homogeneous maritime networks,” in Proc. IEEE ICC, Glasgow, U.K., May 24–28, 2026, pp. 1–6. [19] X. Lu, M. Salehi, M. Haenggi, E. Hossain, and H. Jiang, “Stochastic geometry analysis of spatial-temporal performance in wireless networks: A tutorial,” IEEE Commun. Surveys Tuts., vol. 23, no. 4, pp. 2753–2801, 4th Quart., 2021. [20] M. Z. Win, P. C. Pinto, and L. A. Shepp, “A mathematical theory of network interference and its applications,” Proc. IEEE, vol. 97, no. 2, pp. 205–230, Feb. 2009. [21] H. H. Yang and T. Q. S. Quek, “Spatio-temporal analysis for SINR coverage in small cell networks,” IEEE Trans. Commun., vol. 67, no. 8, pp. 5520–5531, Aug. 2019. [22] Y. Zhong, X. Ge, T. Han, Q. Li, and J. Zhang, “Tradeoff between delay and physical layer security in wireless networks,” IEEE J. Sel. Areas Commun., vol. 36, no. 7, pp. 1635–1647, Jul. 2018. [23] S. Zhang, Y. Zhu, and J. Liu, “Multi-UAV enabled aerial-ground integrated networks: A stochastic geometry analysis,” IEEE Trans. Commun., vol. 70, no. 10, pp. 7040–7054, Oct. 2022. [24] M. J. Lighthill and G. B. Whitham, “On kinematic waves II. A theory of traffic flow on long crowded roads,” Proc. Roy. Soc. London, Ser. A, vol. 229, no. 1175, pp. 317–345, May 1955. [25] B. D. Greenshields, “A study of traffic capacity,” Highway Res. Board Proc., vol. 14, no. 1, pp. 448–477, 1935. [26] D. Helbing and P. Molnár, “Social force model for pedestrian dynamics,” Phys. Rev. E, vol. 51, no. 5, pp. 4282–4286, 1995. [27] R. L. Hughes, “A continuum theory for the flow of pedestrians,” Transp. Res. B: Methodol., vol. 36, no. 6, pp. 507–535, 2002.

[28] H. Kim, J. Park, M. Bennis, S.-L. Kim, and M. Debbah, “Mean-field game theoretic edge caching in ultra-dense networks,” IEEE Trans. Veh. Technol., vol. 69, no. 1, pp. 935–947, Jan. 2020. [29] R. Zheng et al., “Dynamic computation offloading in ultra-dense networks based on mean field games,” IEEE Trans. Wireless Commun., vol. 20, no. 10, pp. 6551–6565, Oct. 2021. [30] Y. Zhang, C. Yang, J. Li, and Z. Han, “Distributed interference-aware traffic offloading and power control in ultra-dense networks: Mean field game with dominating player,” IEEE Trans. Veh. Technol., vol. 68, no. 9, pp. 8814–8826, Sep. 2019. [31] Y. Wang, F. R. Yu, H. Tang, and M. Huang, “A mean field game theoretic approach for security enhancements in mobile ad hoc networks,” IEEE Trans. Wireless Commun., vol. 13, no. 3, pp. 1616–1627, Mar. 2014. [32] J. A. González, F. J. Rodrı́guez-Cortés, O. Cronie, and J. Mateu, “Spatiotemporal point process statistics: A review,” Spatial Statist., vol. 18, pp. 505–544, Nov. 2016. [33] S. Boccaletti et al., “Complex networks: Structure and dynamics,” Phys. Rep., vol. 424, nos. 4–5, pp. 175–308, 2006. [34] Y. Wang, M. Yan, G. Feng, S. Qin, and F. Wei, “Autonomous on-demand deployment for UAV assisted wireless networks,” IEEE Trans. Wireless Commun., vol. 22, no. 12, pp. 9488–9501, Dec. 2023. [35] R. Ding et al., “A novel PODMAI framework enhanced by user demand prediction for resource allocation in spectrum sharing UAV networks,” IEEE Trans. Commun., vol. 73, no. 9, pp. 8243–8259, Sep. 2025. [36] C. Li, Y. Gan, Y. Zhang, and Y. Luo, “A cooperative computation offloading strategy with on-demand deployment of multi-UAVs in UAVaided mobile edge computing,” IEEE Trans. Netw. Service Manag., vol. 21, no. 2, pp. 2095–2110, Apr. 2024. [37] B. Yu, H. Yin, and Z. Zhu, “Spatio-temporal graph convolutional networks: A deep learning framework for traffic forecasting,” in Proc. 27th Int. Joint Conf. Artif. Intell. (IJCAI), Stockholm, Sweden, Jul. 2018, pp. 3634–3640. [38] Z. Wu, S. Pan, G. Long, J. Jiang, and C. Zhang, “Graph WaveNet for deep spatial-temporal graph modeling,” in Proc. 28th Int. Joint Conf. Artif. Intell. (IJCAI), Macao, China, Aug. 2019, pp. 1907–1913. [39] D. J. Daley and D. Vere-Jones, An Introduction to the Theory of Point Processes, Volume I: Elementary Theory and Methods, 2nd ed. New York, NY, USA: Springer, 2003. [40] R. Courant and D. Hilbert, Methods of Mathematical Physics, Volume II: Partial Differential Equations. Hoboken, NJ, USA: Wiley, 2008. [41] W.-Y. Dong et al., “Vorticity dissipation based routing: A fluid-kinetic framework for loop-free transport in ultra-dense networks,” IEEE Trans. Mobile Comput., vol. 25, no. 11, pp. 22370–22387, Nov. 2026. [42] W.-Y. Dong et al., “Thermodynamic digital twins: Physics-native semantic synchronization for 6G kinetic swarms,” arXiv preprint arXiv:2610.05050, Oct. 2026. [43] W.-Y. Dong et al., “Fluid-dynamic interference modeling for LEO megaconstellations: A spatiotemporal kinetic field approach,” arXiv preprint arXiv:2608.20651, Aug. 2026. [44] E. Barmpounakis and N. Geroliminis, “On the new era of urban traffic monitoring with massive drone data: The pNEUMA large-scale field experiment,” Transp. Res. Part C: Emerg. Technol., vol. 111, pp. 50–71, Feb. 2020. [45] W.-Y. Dong, S. Zhao, R.-S. Han, Q. Bi, and S. Chen, “Digital tides: A fluid-dynamic framework for flux-aware infrastructure provisioning in UAV logistics networks,” IEEE Trans. Mobile Comput., vol. 25, no. 10, pp. 16705–16720, Oct. 2026.

Wen-Yu Dong received the B.S. degree in electronic and information engineering from Sichuan University, Chengdu, China, in 2019, and the Ph.D. degree in information and communication engineering from the School of Information and Communication Engineering, Beijing University of Posts and Telecommunications, Beijing, China, in 2025. He is currently a Researcher with the Future Technology Research Center, China Telecom Research Institute, Beijing, China. His current research interests include spaceair-ground integrated networks, fluid-spatiotemporal stochastic geometry (F-STSG), and wireless communications.

17

Weiwei Jiang (Senior Member, IEEE) received the B.Sc. Degree of Electronic Engineering and Ph.D. Degree of Information and Communication Engineering from the Department of Electronic Engineering, Tsinghua University, Beijing, China, in 2013 and 2018, respectively. He is currently an associate professor with the School of Information and Communication Engineering, Beijing University of Posts and Telecommunications, and Key Laboratory of Universal Wireless Communications, Ministry of Education. His current research interests include artificial intelligence, machine learning, big data, wireless communication and edge computing. He has published more than 100 academic papers in IEEE Trans and other journals, with more than 5700 citations in Google Scholar. He was included in Stanford’s World’s Top 2% Scientists list in 2022–2026.

Song Zhao is a senior engineer with the Future Technology Center, China Telecom Research Institute, Beijing, China. He received his Ph.D. degree in Telecommunications and Information Systems from Beijing University of Posts and Telecommunications, Beijing, China. He serves as a delegate of China Telecom in the 3GPP SA2 and SA5 working groups, and as rapporteur for multiple 3GPP work items related to network automation and intelligence. His research interests include wireless channel modeling, radio resource management, proximity services, and network AI.

Qi Bi (Fellow, IEEE) is Chief Scientist of China Telecom and CTO of its Research Institute, specializing in 5G and 6G. He received his M.S. from Shanghai Jiao Tong University and Ph.D. from Pennsylvania State University. He was awarded Bell Labs Fellow in 2002, the Bell Labs President’s Gold Awards in 2000 and 2002, the Asian American Engineer of the Year in 2005, the Beijing Outstanding Contribution Award for Innovation and Entrepreneurship for Overseas Scholars in 2019, and three Patent Silver Prizes from the China National Intellectual Property Administration from 2022 to 2024.

SHENG CHEN (IEEE Life Fellow) is an Emeritus Professor, University of Southampton, Southampton, U.K. He received the B.Eng. degree in control engineering from the East China Petroleum Institute, Dongying, China, in January 1982, the Ph.D. degree in control engineering from City University, London, in September 1986, and the higher doctoral (D.Sc.) degree from the University of Southampton, Southampton, U.K., in August 2005. From 1986 to 1999, he held research and academic appointments with the University of Sheffield, the University of Edinburgh, and the University of Portsmouth, all in U.K. In September 1999, he joined the School of Electronics and Computer Science, University of Southampton, and retired in June 2026. Since then, he is a University Visitor associated with the School of Electronics and Computer Science, University of Southampton. In 2023, Dr. Chen was appointed a Distinguished Professor at Ocean University of China. Professor Chen has established an over 40-year distinguished academic career in the fields of adaptive signal processing, wireless communications, neural networks and machine learning. He has authored/coauthored over 750 research papers. He has 23,000+ Web of Science citations with h-index 67, and 44,000+ Google Scholar citations with h-index 88. Professor Chen is IEEE ComSoc Signal Processing and Computing for Communications (SPCC) 2025 Technical Recognition Award recipient. Professor Chen was elected as a Fellow of the United Kingdom Royal Academy of Engineering in 2014. He is a fellow of Asia-Pacific Artificial Intelligence Association, a fellow of Industry Academy of the International Artificial Intelligence Industry Alliance, a fellow of IET, and one of the original 200 ISI Highly Cited Researchers in engineering (March 2004).

18

Related documents

• Fluid-Spatiotemporal Stochastic Geometry: Information Flow in Non-Stationary Fields#329045
• Transport of spherical microparticles in a three-dimensional vortex flow#1079351
• AnyFlow: Any-Step Video Diffusion Model with On-Policy Flow Map Distillation#182755
• AnyFlow: Any-Step Video Diffusion Model with On-Policy Flow Map Distillation#180694
• Velocities for moving random surfaces#1204156