Paper deep dive
Mean-Field Path-Integral Diffusion: From Samples to Interacting Agents
Michael Chertkov
Intelligence
Status: succeeded | Model: Gemma-4-26B-A4B | Prompt: intel-v1 | Confidence: 92%
Last extracted: 7/20/2026, 2:48:15 PM
Summary
The paper introduces Mean-Field Path-Integral Diffusion (MF-PID), a framework where diffusion samples act as interacting agents whose drift depends self-consistently on population density. This couples generative modeling with multi-agent control via McKean-Vlasov dynamics. Key theoretical results include an exact linear interpolation of global means for quadratic potentials and a reduction to Riccati ODEs for Linear-Quadratic-Gaussian (LQG) systems. Applied to demand-response control in energy systems, MF-PID achieves 19-24% energy savings compared to independent-agent baselines.
Entities (8)
Relation Signals (6)
MF-PID → achieves → Energy Savings
confidence 95% · MF-PID achieves 19--24% reductions in cumulative control energy over independent-agent baselines
MF-PID → appliesto → Demand Response
confidence 95% · Applied to demand-response control of energy systems... MF-PID achieves 19--24% reductions
MF-PID → governs → McKean-Vlasov Equation
confidence 92% · yields a closed McKean–Vlasov stochastic control problem
MF-PID → extends → Path Integral Diffusion
confidence 90% · We introduce Mean-Field Path-Integral Diffusion (MF-PID)... extend the Path Integral Diffusion (PID) framework
MF-PID → models → Stochastic Optimal Transport
confidence 90% · The coupling converts distribution matching into a McKean--Vlasov extension of the stochastic optimal transport problem
LQG → reducesto → Riccati ODEs
confidence 90% · LQG benchmark in which the infinite-dimensional mean-field system reduces to a finite set of Riccati and linear ODEs
Cypher Suggestions (0)
No Cypher suggestions yet.
Abstract
Abstract:Independent sample generation is the prevailing paradigm in modern diffusion-based generative models of AI. We ask a different question: can samples \emph{coordinate} through shared population statistics to transport probability mass more efficiently? We introduce Mean-Field Path-Integral Diffusion (MF-PID), a framework in which samples are promoted to interacting agents whose drift depends self-consistently on the evolving population density. The coupling converts distribution matching into a McKean--Vlasov extension of the stochastic optimal transport problem, unifying generative modeling and multi-agent control under the same Hamilton--Jacobi--Bellman/Kolmogorov--Fokker--Planck duality. We identify two analytically tractable regimes: a Linear--Quadratic--Gaussian (LQG) benchmark in which the infinite-dimensional mean-field system reduces to a finite set of Riccati and linear ODEs, and a Gaussian-mixture regime governed by a piecewise-constant protocol that preserves closed-form solvability. For a quadratic interaction potential with schedule $\beta_t$ and zero base drift we prove that the self-consistent MF guidance is the \emph{exact} linear interpolant between initial and target global means -- a result that holds for arbitrary initial and target densities and any $\beta_t$. Applied to demand-response control of energy systems, where agents aggregated into an ensemble are energy consumers (e.g.\ thermal zones within a building), MF-PID achieves 19--24\% reductions in cumulative control energy over independent-agent baselines while matching the prescribed terminal distribution exactly, and reveals how coordination redistributes actuation effort across heterogeneous sub-populations.
Tags
Links
- Source: https://arxiv.org/abs/2605.00007v1
- Canonical: https://arxiv.org/abs/2605.00007v1
Trouble viewing inline? Open PDF directly →
Full Text
95,413 characters extracted from source content.
Expand or collapse full text
Mean-Field Path-Integral Diffusion: From Samples to Interacting Agents Michael (Misha) Chertkov Graduate Interdisciplinary Program in Applied Mathematics & Department of Mathematics, University of Arizona, Tucson, AZ, USA Correspondence: chertkov@arizona.edu Abstract Independent sample generation is the prevailing paradigm in modern diffusion-based generative models of AI. We ask a different question: can samples coordinate through shared population statistics to transport probability mass more efficiently? We introduce Mean-Field Path-Integral Diffusion (MF-PID), a framework in which samples are promoted to interacting agents whose drift depends self-consistently on the evolving population density. The coupling converts distribution matching into a McKean–Vlasov extension of the stochastic optimal transport problem, unifying generative modeling and multi-agent control under the same Hamilton–Jacobi–Bellman/Kolmogorov–Fokker–Planck duality. We identify two analytically tractable regimes: a Linear–Quadratic–Gaussian (LQG) benchmark in which the infinite-dimensional mean-field system reduces to a finite set of Riccati and linear ODEs, and a Gaussian-mixture regime governed by a piecewise-constant protocol that preserves closed-form solvability. For a quadratic interaction potential with schedule βt _t and zero base drift we prove that the self-consistent MF guidance is the exact linear interpolant between initial and target global means — a result that holds for arbitrary initial and target densities and any βt _t. Applied to demand-response control of energy systems, where agents aggregated into an ensemble are energy consumers (e.g. thermal zones within a building), MF-PID achieves 19–24% reductions in cumulative control energy over independent-agent baselines while matching the prescribed terminal distribution exactly, and reveals how coordination redistributes actuation effort across heterogeneous sub-populations. The energy saving is independent of the number of zones per building (d=1d=1–3232 tested), confirming that the linear guidance formula broadcasts a single d-vector with (d)O(d) communication and grows mildly in compute (sub-cubically for d≤32d≤ 32, asymptotically (d3)O(d^3) for d≫1d 1). Mean-Field Path-Integral Diffusion (MF-PID) == Sample/Agent-Coordinated Optimal TransportIndependent Agentsdxt= score dt+dWtdx_t= [rgb]1,0,0 [named]pgfstrokecolorrgb1,0,0 score dt+dW_tp(in)p^(in)p(out)p^(out)Path Integral Diffusion score =∇xlogψt(x)= _x _t(x)"Schödinger Bridge" +Vt(x)+V_t(x)ℰIA(0)E^IA(0)Scen. A: 31.3 B: 17.2MFInteracting AgentsMF-PID = H-PIDwith νt=[xt] _t=E[x_t]p(in)p^(in)ρt _tp(out)p^(out)ut(∗)=∇logψt(x,ρt)u^(*)_t=∇\! _t(x, _t)Vt(eff)(x)=∫Vt(x−y)ρt(y)dyV^(eff)_t(x)\!\!=\!\! \!V_t(x\!-\!y) _t(y)dyℰMFE^MF↓ 11–23% energy savedExact LQGLinear f, Quadratic VVGaussian p(tar)p^(tar)⟶ Riccati ODEs Score is Explicitp(init)p^(init) & p(tar)p^(tar) are GMs+ βt _t is PWC⇒ ψt _t is GM (in x)& interval-analytic in ttApplyDemand ResponseTCL fleetℰE – energy savings ABIA(0)IA(0)IA(m¯)IA( m)MFMFExtends to Other Naturaland Engineered Systems Overview of the MF-PID framework. (a) Independent Agents/Samples: In the baseline, each sample evolves under a pre-computed drift; trajectories carry no information about one another, yielding no inter-agent coupling. (b) Interacting Agents: MF-PID promotes samples to interacting agents whose optimal drift ut(∗)u^(*)_t depends self-consistently on the evolving population density ρt _t through the effective potential Vt(eff)V^(eff)_t. (c) LQG: Two analytically tractable regimes— Linear–Quadratic–Gaussian (LQG), where the infinite-dimensional MF system reduces to Riccati ODEs with variance/mean decoupling; and Gaussian-Mixture (GM) representation for densities with Piece-Wise-Constant (PWC) in time protocol for βt _t results in an explicit – gradient-of-log-of-Gaussian mixture expression for the score function. (d) DR: Application to demand-response control of a building TCL fleet; MF coordination achieves 11.6 % (Scenario A) and 22.6 % (Scenario B) reductions in cumulative control energy ℰE relative to the unguided IA baseline. Introduction Generative AI has been transformed by diffusion models, which frame sample generation as a stochastic process steered from noise to data [1, 2, 3]. A key structural feature of these models — shared with other generative models, e.g. normalizing flows [4, 5] — is that samples are generated independently: the trajectory of one particle carries no information about any other. Similarly, stochastic optimal transport (SOT) and Schrödinger bridge formulations [6, 7, 8] cast distribution matching as an independent-particle path optimization, yielding tractable convolutions of Green functions but discarding inter-particle information; stochastic interpolants [9] construct flexible transport bridges between arbitrary densities via tunable continuous-time stochastic processes, recovering the Schrödinger bridge as a special limit — again in an independent-particle framework. A natural complementary question is whether coordinated generation can improve efficiency. In physical and engineered systems, collective behavior routinely outperforms individual action: flocks exploit aerodynamic coupling [10], synchronized HVAC fleets reduce peak demand [11, 12, 13, 14], and robotic swarms leverage formation geometry [15, 16]. In each case, coupling through a shared field — air pressure, a power grid, or relative position — allows the population to achieve its collective objective at lower individual cost. We make this intuition precise for generative modeling. Specifically, we extend the Path Integral Diffusion (PID) framework [17, 18, 19] — itself a strict generalization of the Schrödinger bridge [6, 7, 8, 20, 21], recovered as the special case of PID — to a Mean-Field (MF) setting in which an ensemble of agents, each performing a controlled diffusion, interacts through its evolving empirical distribution. In the limit of infinitely many agents this yields a closed McKean–Vlasov stochastic control problem: the optimal drift of each agent depends self-consistently on the population density, introducing a nonlinear coupling absent in Independent-Agent (IA) formulations. Four contributions. First, we derive the MF-PID equations in both terminal-cost and SOT (hard marginal constraint) formulations (SI §1), identifying the precise way in which mean-field coupling breaks global integrability of the IA case while preserving local analytical structure. Second, we show that when the base drift is linear, the interaction potential quadratic, and the target Gaussian, the infinite-dimensional mean-field system collapses to a finite set of coupled Riccati and linear ODEs (SI §2). This LQG benchmark is the first closed-form mean-field generalization of PID, providing explicit energetic comparisons between MF and IA strategies. The covariance dynamics are identical in the MF and IA cases; MF coordination operates entirely through the linear coefficient sts_t, giving a clean analytical separation of the energetic benefit. Third, we prove that for a quadratic interaction potential with zero base drift, the self-consistent MF guidance is the exact linear interpolant between initial and target global means (Theorem, SI §3.2). The result is independent of the β-schedule and the shape of both distributions. The proof rests on a structural cancellation in the Itô–HJB system: nonlinear score terms cancel exactly in the mean acceleration equation, leaving a zero-force condition that forces the mean to evolve linearly. This converts what appeared to be a nonlinear fixed-point problem into a one-shot explicit computation. Fourth, we extend analytic tractability to multi-modal Gaussian-mixture targets via a piecewise-constant (PWC) protocol (SI §3). With the guidance known analytically, the score function is assembled in a single pass through closed-form Green-function coefficients without iteration. The formulation naturally accommodates non-delta initial distributions, enabling closed-form transport bridges. Application. We demonstrate MF-PID on a physically motivated Demand Response (DR) scenario for large ensembles of multi-zone buildings. Post-curtailment recovery is cast as a Gaussian-mixture SOT bridge with K sub-population types and d zones per building. MF coordination reduces cumulative control energy by 19–24% relative to independent-agent baselines while exactly matching the terminal temperature distribution. Three scalability properties are verified numerically: (i) the per-zone energy saving is invariant across d=1d=1–3232 zones, confirming the theorem’s dimension-independence; (i) the saving grows consistently (19%–22%) as fleet heterogeneity increases from K=2K=2 to K=8K=8 sub-types; (i) adding AR(1) inter-zone thermal coupling (ρ=0ρ=0–0.80.8) leaves the saving unchanged at ≈21%≈ 21\%. The guidance is a single linearly-interpolated d-vector, broadcast once to the fleet with no iteration required. Relation to existing work. The value of inter-sample communication is not new to filtering. Ensemble Kalman filters (EnKF) [22] famously replace an intractable Gaussian update with a finite-ensemble approximation in which every particle is corrected by the empirical covariance of the full ensemble — a linear, observation-driven coupling that breaks particle independence at each assimilation step. Our setting differs in three respects: (i) there are no sequential observations; the coupling arises instead from a terminal-cost objective and an explicit interaction potential; (i) the interaction is nonlinear and self-consistent, governed by a McKean–Vlasov equation rather than a Kalman gain; and (i) optimality is measured by control energy rather than posterior approximation error. MF-PID can therefore be read as the continuous-time, generative-transport analogue of the EnKF idea: structured inter-sample communication in the service of a well-defined variational objective. Mean-field Schrödinger bridges have been formulated and their ergodic and propagation-of-chaos properties analyzed [23, 24]. Our work extends this line by placing mean-field coupling within the Path-Integral Diffusion (PID) framework introduced in [17], which strictly generalizes the classical Schrödinger bridge. The independent-agent bridge is recovered as the special case Vt=ft=At=0V_t=f_t=A_t=0 of PID. Harmonic PID (H-PID) further allows affine base drift and quadratic potentials while preserving Gaussian Green functions [17] (see also related analysis in [25]), and Guided PID introduces a time-dependent quadratic steering potential, Vt(x)=βt(x−ηt)2/2V_t(x)= _t(x- _t)^2/2, that shapes trajectories without coupling particles [19]. In all of these constructions, agents remain independent. MF-PID replaces the externally prescribed potential by a self-consistent effective interaction, Vt(eff)(x)=12∫βt(x−y)2pt(y)dy,V_t^(eff)(x)= 12 _t(x-y)^2\,p_t(y)\,dy, thereby promoting samples to interacting agents governed by a McKean–Vlasov stochastic control problem. Unlike existing mean-field Schrödinger formulations [23, 24], our setting incorporates both a nontrivial base drift and an explicit interaction potential. After the Cole–Hopf transform, the system reduces to quasi-linear HJB–KFP equations coupled through the evolving density, yielding a class of mean-field entropic stochastic optimal transport problems not previously analyzed. While linear–quadratic mean-field control and games are classical [26, 27], they have not been synthesized with entropic SOT or diffusion-based generative modeling. Stochastic interpolants [9] provide a broad independent-particle framework for flows and diffusions – MF-PID add structure (via PID construct) and then lifts this paradigm to the interacting regime. Finally, the piecewise-constant analytic machinery developed in [18, 19] is extended here to accommodate mean-field coupling without sacrificing closed-form tractability. Results The MF-PID Framework Setup. Consider N exchangeable agents, each evolving under the controlled Itô diffusion dxt(i)=(ft(xt(i))+ut(i))dt+dWt(i),x0(i)=0,dx_t^(i)= (f_t(x_t^(i))+u_t^(i) )dt+dW_t^(i), x_0^(i)=0, (1) where ftf_t is a base drift, ut(i)u_t^(i) is the control (score), and Wt(i)W_t^(i) are independent standard Brownian motions. As N→∞N→∞, the marginal density ptp_t is governed by the MF Kolmogorov–Fokker–Planck (KFP) equation, and each agent obeys the McKean–Vlasov SDE dxt=(ft(xt)+ut(xt,pt))dt+dWt,x0=0.dx_t= (f_t(x_t)+u_t(x_t,p_t) )dt+dW_t, x_0=0. (2) Optimal control. We minimize the mean-field cost Jt(x)=infu∫t1(12‖ut′‖2+∫Vt′(xt′−y)pt′(y)dy)dt′,J_t(x)= _u\,E\! _t^1\!\! ( 12\|u_t \|^2+\! \!V_t (x_t -y)\,p_t (y)\,dy )dt , (3) where Vt(⋅)V_t(·) is an interaction potential penalizing relative displacement. Under the Hopf–Cole substitution Jt=−logψtJ_t=- _t, the HJB equation for ψ becomes linear, and the optimal control is ut∗=∇xlogψtu_t^*= _x\! _t. In the SOT formulation, a hard constraint p1=p(tar)p_1=p^(tar) replaces the terminal penalty, and the optimal control is ut(∗)(x)=∇xlog∫p(tar)(y)Gt(−;(eff))(x;y)G1(+;(eff))(y;0)dy,u^(*)_t(x)= _x \!\!p^(tar)(y) G_t^(-;(eff))(x;y)G_1^(+;(eff))(y;0)\,dy, (4) where Gt(±;(eff))G_t^(±;(eff)) are Green functions of the HJB/KFP system evaluated under the effective potential Vt(eff)(x)=∫Vt(x−y)pt(y)dyV^(eff)_t(x)= V_t(x-y)\,p_t(y)\,dy. The MF coupling. In the IA case the Green functions can be computed independently of ptp_t; in the MF case, Gt(±;(eff))G_t^(±;(eff)) depend on ptp_t through Vt(eff)V^(eff)_t, and ptp_t in turn depends on the Green functions. This self-referential structure is the source of both the difficulty and the power of MF-PID. Full derivations are given in SI §1. LQG Benchmark: Analytic Closed-Form Solution Model. We specialize to a linear base drift ft(x)=Atx+Btf_t(x)=A_tx+B_t, quadratic interaction Vt(x)=12x⊤QtxV_t(x)= 12x Q_tx (Qt⪰0Q_t 0), and Gaussian target p(tar)=(m1(tar),Σ1(tar))p^(tar)=N(m_1^(tar), _1^(tar)). Gaussianity is preserved under the controlled dynamics, so pt=(mt,Σt)p_t=N(m_t, _t) for all t∈[0,1]t∈[0,1]. The effective potential is Vt(eff)(x)=12(x−mt)⊤Qt(x−mt)+12Tr(QtΣt)V^(eff)_t(x)= 12(x-m_t) Q_t(x-m_t)+ 12Tr(Q_t _t) and a quadratic ansatz for JtJ_t yields the closed system −S˙t - S_t =Qt+StAt+At⊤St−St2, =Q_t+S_tA_t+A_t S_t-S_t^2, (5) −s˙t - s_t =−Qtmt+At⊤st−Stst+StBt, =-Q_tm_t+A_t s_t-S_ts_t+S_tB_t, (6) m˙t m_t =(At−St)mt+Bt−st, =(A_t-S_t)m_t+B_t-s_t, (7) Σ˙t _t =(At−St)Σt+Σt(At−St)⊤+I, =(A_t-S_t) _t+ _t(A_t-S_t) +I, (8) with boundary conditions m0=0m_0=0, Σ0=0 _0=0, m1=m1(tar)m_1=m_1^(tar), Σ1=Σ1(tar) _1= _1^(tar). The optimal MF control is affine: ut∗(x)=−Stx−stu_t^*(x)=-S_tx-s_t. A key structural observation is that (5)–(8) decouple into two sequential sub-problems: (i) a variance block (St,Σt)(S_t, _t) that depends only on QtQ_t and AtA_t and is identical in the MF and IA cases; and (i) a mean block (mt,st)(m_t,s_t) whose source term −Qtmt-Q_tm_t couples to the current population mean rather than a fixed exogenous centre. MF coordination therefore operates entirely through the linear coefficient sts_t — a clean analytical separation of the energetic benefit. Scalar TCL example. We illustrate with d=1d=1, f(x)=−κxf(x)=-κ x (Ornstein–Uhlenbeck relaxation), and V(x)=q2x2V(x)= q2x^2. This models a population of Thermostatically Controlled Loads (TCLs) [11, 12, 28, 29, 30, 31], where xtx_t is a normalized temperature deviation. Let Δ=κ2+q = κ^2+q. A closed-form shooting procedure (SI §2) determines the unique ρ∈(−1,1)ρ∈(-1,1) such that the bridge constraint Σ1=(σ(tar))2 _1=(σ^(tar))^2 is satisfied: ρ=A−1Ar0−1,A=2Δ(σ(tar))21−r0,r0=e−2Δ,ρ= A-1Ar_0-1, A= 2 (σ^(tar))^21-r_0, r_0=e^-2 , (9) after which StS_t, Σt _t, and the mean mt=m(tar)sinh(κt)/sinh(κ)m_t=m^(tar) (κ t)/ (κ) are all explicit. We compare the MF bridge to the IA family parameterised by an exogenous centre m¯∈[0,m(tar)] m∈[0,m^(tar)]. Since StS_t and Σt _t are identical across schemes, the comparison reduces to the linear coefficient sts_t. Three performance metrics are used: the mean trajectory mtm_t, the instantaneous control power (t)=St2σt2+(Stmt+st)2P(t)=S_t^2 _t^2+(S_tm_t+s_t)^2, and the cumulative energy ℰ(t)=∫0t(u)duE(t)= _0^tP(u)du. Figure 1: MF vs. IA controls in the scalar TCL example. Top left: Trajectory ensembles for MF (blue), IA with m¯=0 m=0 (orange), and IA with m¯=m(tar) m=m^(tar) (green); thin curves: sample paths, thick: analytic means, shaded band: target variance. Top right: Instantaneous control power (t)P(t). Bottom left: Cumulative control energy ℰ(t)=∫0t(u)duE(t)= _0^tP(u)\,du. Bottom right: Mean trajectories mtm_t. All curves are obtained analytically. Fig. 1 shows that the MF mean follows a smooth hyperbolic sine arc, while IA means exhibit stronger curvature driven by Δ . The MF cumulative energy curve lies strictly below both IA curves throughout [0,1][0,1]: MF coordination achieves the same terminal accuracy with reduced total control effort. The advantage originates entirely from sts_t, which in the MF case adapts to the endogenous population mean rather than a fixed exogenous reference. Exact Linear MF Guidance For the isotropic quadratic potential Vt(x)=βt2‖x‖2V_t(x)= _t2\|x\|^2 with ft≐0f_t 0, the effective potential reduces to Vt(eff)(x)=βt2‖x−νt‖2+constV^(eff)_t(x)= _t2\|x- _t\|^2+const, where νt=mt=pt∗[x] _t=m_t=E_p_t^*[x]. MF-PID therefore becomes a guided H-PID whose guidance is the solution of the self-consistency condition νt(MF)=x∼pt∗(⋅|ν(MF))[x],t∈[0,1]. _t^(MF)=E_x p_t^*(\,·\,|\,ν^(MF))[x], t∈[0,1]. (10) Theorem (SI Thm. 3.1) Let ft≐0f_t 0, Vt(x)=βt2‖x‖2V_t(x)= _t2\|x\|^2 for any βt>0 _t>0, and let p(in)p^(in), p(tar)p^(tar) be probability measures on ℝdR^d with finite first moments m¯(in)=p(in)[x] m^(in)=E_p^(in)[x], m¯(tar)=p(tar)[x] m^(tar)=E_p^(tar)[x]. Then the self-consistent MF guidance satisfies νt(MF)=mt=(1−t)m¯(in)+tm¯(tar)∀t∈[0,1], _t^(MF)=m_t=(1-t)\, m^(in)+t\, m^(tar) ∀\,t∈[0,1], (11) exactly, independently of βt _t, the number and geometry of mixture components, and the shape of p(in)p^(in) and p(tar)p^(tar). Proof sketch. Applying Itô’s formula to ut∗(xt)u_t^*(x_t) and differentiating m˙t=[ut∗(xt)] m_t=E[u_t^*(x_t)] gives m¨t=[∂tut∗+(ut∗⋅∇)ut∗+12Δut∗] m_t=E[ _tu_t^*+(u_t^*·∇)u_t^*+ 12 u_t^*]. Taking the spatial gradient of the HJB equation yields ∂tut∗=∇xVt(eff)−(ut∗⋅∇)ut∗−12Δut∗ _tu_t^*= _xV^(eff)_t-(u_t^*·∇)u_t^*- 12 u_t^*, so the nonlinear score terms cancel exactly inside m¨t m_t, leaving m¨t=[∇xVt(eff)(xt)]=βt(mt−νt(MF)) m_t=E[ _xV^(eff)_t(x_t)]= _t(m_t- _t^(MF)). The self-consistency condition νt(MF)=mt _t^(MF)=m_t then forces m¨t=0 m_t=0, and the linear arc (11) follows from the boundary conditions m0=m¯(in)m_0= m^(in), m1=m¯(tar)m_1= m^(tar). Full proof in SI §3.2. Corollaries. When p(in)=δ(⋅)p^(in)=δ(·), equation (11) reduces to νt(MF)=tm¯(tar) _t^(MF)=t\, m^(tar) (SI Cor. 3.2). For ft(x)=−κxf_t(x)=-κ x the same Itô–HJB cancellation yields m¨t=κ2mt m_t=κ^2m_t, recovering the sinh -arc mt=m(tar)sinh(κt)/sinh(κ)m_t=m^(tar) (κ t)/ (κ) of the LQG benchmark (SI Cor. 3.3). Perspective (alternative viewpoint) Self-consistent guidance. Although we derive the linear guidance by formulating MF-PID and reducing it to guided H-PID, the logic can be reversed. Start from guided PID with prescribed νt _t and potential Vt(x)=βt2‖x−νt‖2V_t(x)= _t2\|x- _t\|^2. A natural “closure” is to choose νt _t so that it matches the mean it induces: νt=pt(ν)[x],t∈[0,1], _t=E_p_t^(ν)[x], t∈[0,1], where pt(ν)p_t^(ν) is the time-marginal under guidance ν. This turns guidance design into a fixed-point problem in trajectory space. Our main theorem shows that for ft≐0f_t 0 and quadratic interaction this fixed point is explicit: νt=(1−t)m¯(in)+tm¯(tar), _t=(1-t)\, m^(in)+t\, m^(tar), independent of βt _t and of the shapes of the endpoint distributions. In this view, MF-PID provides the variational/self-consistency principle that selects νt _t, while guided PID supplies the operational generative mechanism. Practical implication. Eq. (11) provides the MF guidance in closed form, without any iteration: one computes the two global means m¯(in) m^(in) and m¯(tar) m^(tar), sets νt(MF)=(1−t)m¯(in)+tm¯(tar) _t^(MF)=(1-t) m^(in)+t m^(tar), and proceeds directly to score function evaluation via the closed-form Green-function coefficients (SI §3.3–3.4). The self-consistency fixed-point iteration is therefore unnecessary for ft≐0f_t 0 and is retained in SI only for reference and for the ft≐̸0f_t 0 setting. Gaussian-Mixture Score and Demand Response Explicit score assembly. For a Gaussian-mixture target p(tar)=∑k=1Kπk(x;mk,Σk)p^(tar)= _k=1^K _kN(x;\,m_k, _k) with guidance (11) set analytically, the PWC Green-function coefficients within each interval [ti,ti+1)[t_i,t_i+1) are determined by the Riccati ODEs of SI §3.3. With the guidance pre-computed, all K sets of coefficients are evaluated in a single forward pass (no outer iteration). The score function then takes the closed form ut(∗)(x)=bt(−)(y^(t;x)−Υt(x)),u^(*)_t(x)=b_t^(-) ( y(t;x)- _t(x) ), (12) where y^(t;x) y(t;x) is the mixture-weighted posterior target-component mean and Υt(x) _t(x) an affine function of x; all quantities are explicit hyperbolic functions of the PWC coefficients (SI §3.3–3.4, eqs. S3.4–S3.9). No neural network or iterative solve is required. Demand-response setting. We consider a fleet of buildings, each with d thermal zones, partitioned into occupied (πocc=0.60 _occ=0.60, mocc(tar)=0.0m^(tar)_occ=0.0, σocc(tar)=0.20σ^(tar)_occ=0.20) and unoccupied (πunocc=0.40 _unocc=0.40, munocc(tar)=1.5m^(tar)_unocc=1.5, σunocc(tar)=0.30σ^(tar)_unocc=0.30) sub-populations (non-dimensional units; 0=^ 20∘C0\, =\,20\, C, unit =^ 3∘C =\,3\, C). A curtailment event displaces temperatures from setpoints; the recovery phase is cast as a GM-to-GM H-PID bridge. By Theorem (11) the guidance νt(MF)=(1−t)m¯(in)+tm¯(tar) _t^(MF)=(1-t) m^(in)+t m^(tar) is exact, requiring no iteration. We compare three strategies — IA(ν=0ν=0), IA(ν=m¯ν= m), and MF (11) — under two post-curtailment scenarios (Fig. 2; full diagnostics in SI §4): • Scenario A (wide): modes overlap heavily (σ(in)=3.0σ^(in)=3.0), aggressive curtailment. • Scenario B (narrow): modes well-separated (σ(in)=(0.5,0.7)σ^(in)=(0.5,0.7)), mild curtailment. Figure 2: Top: Initial and target distributions. Middle/Bottom: Sample trajectories (50 paths) for Scenarios A (middle) and B (bottom). Blue: occupied; orange: unoccupied. Black: ensemble mean; gray band: ±1σ± 1σ; shaded rectangles: target ±σk(tar)± _k^(tar) bands. Energy results. The ordering ℰMF<ℰIA(m¯)<ℰIA(0)E^MF<E^IA( m)<E^IA(0) holds in both scenarios. For Scenario A: ℰMF=27.67E^MF=27.67, ℰIA(m¯)=29.68E^IA( m)=29.68, ℰIA(0)=31.30E^IA(0)=31.30 (11.6% saving). For Scenario B: ℰMF=13.27E^MF=13.27, ℰIA(m¯)=15.47E^IA( m)=15.47, ℰIA(0)=17.15E^IA(0)=17.15 (22.6% saving). Table 1 decomposes these totals by sub-population: MF slightly increases cost for the easy (occupied) mode while substantially reducing it for the expensive (unoccupied) one — in Scenario B the unoccupied energy drops from 37.40 to 28.07 (−-25%), invisible at the level of the global mean. Table 1: Per-mode control energy. MF guidance (11) redistributes effort from the occupied to the unoccupied sub-population, yielding net savings. Scenario A Scenario B Method ℰoccE_occ ℰunoccE_unocc Total ℰoccE_occ ℰunoccE_unocc Total ν=0ν=0 13.89 56.92 31.30 3.38 37.40 17.15 ν=m¯ν= m 13.43 53.59 29.68 2.63 34.36 15.47 ν=νMFν=ν^MF 14.72 46.73 27.67 3.21 28.07 13.27 Mechanism. The MF advantage amplifies from Scenario A to B because fleet heterogeneity increases. In Scenario A the two modes overlap heavily, the ensemble behaves nearly as a single Gaussian, and a fixed constant guidance already captures most of the benefit. In Scenario B the modes are quasi-independent: the unoccupied cluster travels ≈4≈4 units while the occupied cluster moves only ≈1.5≈1.5. The self-consistent guidance (11) adapts to this asymmetry, whereas any constant ν cannot. The energy savings arise from the timing of guidance imposed by the βt _t schedule on the exact linear trajectory, not from a nonlinear displacement of νt _t itself. Mean-trajectory, power-spectrum, and guidance residual diagnostics for both scenarios are provided in SI §4. Multi-Zone Scalability We now lift the d=1d=1 restriction and test the method across three axes: zone count d, fleet heterogeneity K, and inter-zone thermal coupling. In all cases the MF guidance remains (11): a single linearly-interpolated d-vector requiring no iteration. Dimension sweep (d=1d=1–3232, K=2K=2). Each particle represents a building whose state x∈ℝdx ^d encodes temperature deviations of d zones (perimeter and interior, distinguished by a sinusoidal zone-type vector). We use Scenario-B parameters broadcast isotropically to d dimensions. Results are summarized in Table 2. Table 2: Dimension sweep (K=2K=2, Scenario B parameters). Per-zone energy ℰ(1)/dE(1)/d is invariant across d=1d=1–3232, confirming Theorem (11)’s dimension-independence. Wall-clock time (single CPU, B=4,000B=4,000 particles, nsteps=2,500n_steps=2,500) scales sub-cubically in this range; the asymptotic (d3)O(d^3) Cholesky cost dominates beyond d≫32d 32. d ℰ/dE/d (MF) ℰ/dE/d (IA0) ℰ/dE/d (IAm) Saving Time (s) 1 12.26 16.17 14.34 24.2% 88 2 13.07 16.96 15.13 23.0% 93 4 13.37 17.22 15.41 22.4% 105 8 13.57 17.40 15.59 22.0% 106 16 13.50 17.32 15.51 22.1% 152 32 13.57 17.19 15.42 21.1% 353 The key finding is that ℰ(1)/d≈13.4±0.2E(1)/d≈ 13.4± 0.2 (MF) and 17.2±0.317.2± 0.3 (IA0) across the entire range: the per-zone energy saving of ≈22%≈ 22\% is independent of d. Each additional zone costs the same to coordinate as the first. This directly validates Theorem (11): the guidance adds (d)O(d) compute (a single mean vector) regardless of how many zones share the building. The wall-clock time grows from 88 s (d=1d=1) to 353 s (d=32d=32) — a 4× increase for a 32× increase in dimension. For d≤8d≤ 8 the overhead from the Python time-stepping loop dominates the d×d× d Cholesky factorizations; the compute cost scales roughly as (d1.3)O(d^1.3) over d=1d=1–3232 and is expected to transition to the asymptotic (d3)O(d^3) regime for d≫32d 32. For real-time building control (d≤50d≤ 50, dispatch intervals ≥1≥\!1 min), this places MF-PID comfortably within operational time budgets. Fleet heterogeneity sweep (K=2K=2–88, d=4d=4). We fix d=4d=4 and vary the number of building sub-types K from 2 to 8. Component means are spaced uniformly and initial distributions are displaced by +4 units (aggressive curtailment), with weights decreasing geometrically so that efficient-to-move buildings are more numerous. The MF saving grows consistently: 19.3% (K=2K=2), 21.0% (K=3K=3), 21.6% (K=4K=4), 22.4% (K=8K=8). The trend is real but moderate because this design keeps the per-component displacement constant (|Δmk|=4| m_k|=4 for all k): only the weight asymmetry increases with K. The more dramatic contrast (11.6% vs 22.6%) seen in the d=1d=1 Scenarios A/B arises from geometric mode heterogeneity (different travel distances), which is the dominant driver of MF advantage. A structural observation from the K-sweep: by construction the target global mean m¯(tar)=0 m^(tar)=0 for all K, so the two IA baselines (ν≐0ν 0 and ν≐m¯(tar)ν m^(tar)) coincide exactly. This confirms a general property: for any two constant-guidance IA strategies that share the same time-invariant reference, the costs are identical regardless of K or d. The MF guidance, by contrast, uses a time-varying interpolant and achieves strictly lower cost in all cases. Inter-zone coupling (ρ=0ρ=0–0.80.8, d=8d=8, K=2K=2). Real buildings have spatially correlated zone temperatures through shared walls and HVAC ducts. We model this with a spatial AR(1) covariance: Σij(k)=σk2ρ|i−j|,i,j=0,…,d−1, ^(k)_ij= _k^2\,ρ^|i-j|, i,j=0,…,d-1, (13) where |i−j||i-j| is the integer distance between zone indices and ρ∈[0,1)ρ∈[0,1) controls how rapidly correlation decays with spatial separation. “AR(1)” refers to the first-order autoregressive structure in the spatial index: each zone’s temperature is most similar to its immediate neighbours, with influence decreasing geometrically with distance. Here ρ=0ρ=0 recovers independent zones, ρ=0.5ρ=0.5 gives moderate coupling (adjacent-zone correlation 0.5), and ρ=0.8ρ=0.8 represents strongly coupled perimeter-to-interior heat transfer (correlation 0.64 across two zones, 0.51 across three). (Note – as a remark towards future work – that the spatial AR(1) is distinct from accounting for a spatio-temporal correlation, which would additionally link the same zone across time steps. Within our framework, temporal dependence between successive states is naturally introduced by a nonzero base drift ft(x)f_t(x) — for example, the OU relaxation ft(x)=−κxf_t(x)=-κ x of the LQG section couples the current state to its past through the mean-reversion dynamics. Theorem (11) applies with spatial AR(1) covariances unchanged (it requires only finite first moments, not diagonal covariances).) Results: the MF saving is 22.0%22.0\% at ρ=0ρ=0, 21.8%21.8\% at ρ=0.5ρ=0.5, and 21.0%21.0\% at ρ=0.8ρ=0.8 — insensitive to zone coupling. Increasing ρ raises absolute energies (coupled zones collectively require more coordinated actuation) but the relative MF advantage is preserved, confirming the method needs no retuning for realistic building envelopes. Figure 3: Multi-zone scalability (K=2K=2). (a) Per-zone cumulative energy ℰ(1)/dE(1)/d vs. number of zones d. MF (blue diamonds) is flat at ≈13.5≈ 13.5; IA(ν=0ν=0) (orange) at ≈17.2≈ 17.2. (b) MF energy saving (%) vs. d: stable at 21–24% across the full range. (c) Wall-clock time on a single CPU core. The dashed line shows the asymptotic (d3)O(d^3) trend (Cholesky of d×d× d covariance matrices per time step); for d≤32d≤ 32 overhead dominates and scaling is milder. Figure 4: Spatial mechanism at d=10d=10, K=2K=2. Left: Zone×time heatmap of ensemble mean temperature x¯j(t) x_j(t) for MF (top) and IA(ν=0ν=0) (bottom). Under MF each zone mean follows a near-linear arc; under IA high-displacement zones (top/bottom rows) show stronger curvature. Right: Per-zone cumulative energy. MF redistributes effort: high-displacement zones (e.g. zone 5) see the largest reduction (≈14%≈ 14\%) while low-displacement zones change negligibly. Total MF saving: 22.3%. Discussion Samples as agents. MF-PID reveals a precise mathematical duality between generative modeling and coordinated stochastic control: what appears in generative AI as “sampling from a target distribution” is equivalent, in the MF limit, to “redistributing a population under minimal actuation.” The governing equations — HJB, KFP, and their Green-function kernels — do not distinguish between the two interpretations; only the context does. The duality between samples and interacting agents also has a distinguished antecedent in data assimilation: ensemble Kalman filters [22] achieve tractable Bayesian updates by having each particle absorb information from the empirical ensemble covariance. MF-PID extends this philosophy beyond the linear-Gaussian, observation-driven setting: the coupling is encoded in a potential, the interaction is self-consistent across continuous time, and the figure of merit is control energy rather than filtering accuracy. Energetic value of coordination. The LQG analysis provides a clean analytical proof that MF coordination is energetically superior to exogenously centred IA strategies. The advantage originates in the endogenous population feedback embedded in sts_t: the MF controller adapts its reference to the current ensemble mean rather than a pre-fixed value. In the Gaussian-mixture setting, the self-consistent field allocates actuation effort proportionally to per-mode transport difficulty. Exact linearity and its implications. The theorem (11) is the central new analytical result of this paper. What appeared empirically as “near-linearity” of the guidance is in fact exact linearity — a consequence of a structural cancellation in the Itô–HJB system that holds for any score function, any distribution, and any β-schedule. This has three direct practical payoffs. First, the guidance is known analytically before any simulation is run, eliminating online iteration and reducing the communication overhead to a single broadcast of a linearly-interpolated d-vector — ideal for real-time dispatch of large building fleets. Second, the exact linear guidance (11) provides a provably correct initializer for any hybrid neural–analytic extension of MF-PID, removing the need for warm-start iteration in high dimensions. Third, for the OU base drift ft(x)=−κxf_t(x)=-κ x, the same Itô–HJB argument yields the sinh -arc in closed form (Corollary, SI §3.2), so the LQG benchmark mean-dynamics result is a special case of a broader structural principle. Open-loop vs. closed-loop inference. The MF-PID construction is distributionally closed-loop at design time but open-loop at inference time: once νt(MF) _t^(MF) is fixed by (11), individual samples evolve independently under a pre-determined drift. A natural extension is to replace νt(MF) _t^(MF) in the drift with the empirical ensemble mean m^t(N)=N−1∑i=1Nxt(i) m_t^(N)=N^-1 _i=1^Nx_t^(i) during simulation, similarly to the DR implementation in [30]. Since Theorem (11) guarantees that the theoretical guidance equals the population mean, this closed-loop implementation introduces no additional design cost: agents broadcast their current mean and apply the linear schedule. In continuous time this corresponds to a McKean–Vlasov SDE; stability is suggested by the contraction properties established for Schrödinger bridge iterations in linear settings [21]. Neuralization and high dimensions. The present construction is deliberately free of neural networks: the analytic GM-to-GM setting serves as a mathematically controlled backbone. For high-dimensional targets beyond the Gaussian-mixture regime, the PWC inner-solve structure provides a natural scaffold for hybridization: learn score corrections beyond quadratic structure using neural networks while retaining analytic Green-function evolution within each PWC interval. Crucially, the exact linear guidance (11) provides a provably correct initialization for any such learning procedure. Neural SDE frameworks [32] provide the theoretical underpinning. Physical time and engineering implications. Unlike diffusion models operating in synthetic noise time, MF-PID operates in real physical time: recovery horizons, maneuvering windows, or reliability intervals. Control energy, peak actuation, and transient risk are first-class design quantities. The DR application demonstrates three properties of direct relevance to energy-system practitioners. Dimension invariance: the per-zone saving of ≈22%≈ 22\% is flat across d=1d=1–3232 zones per building. A dispatcher can scale from a scalar TCL model to a multi-zone office floor without retuning, and the guidance computation adds only (d)O(d) work (one matrix-vector product per dispatch interval). Fleet heterogeneity: the saving grows consistently as the number of distinct building sub-types K (number of components in Gaussian mixture) increases, because MF coordination naturally exploits the disparity in per-type transport distances. Zone coupling: AR(1) inter-zone thermal correlation (ρ=0ρ=0–0.80.8) leaves the saving unchanged, validating the approach for realistic building envelopes without any algorithmic modification. Taken together, these findings show that the (d)O(d) guidance broadcast is not merely a theoretical convenience but a practical enabler of real-time fleet control at building-relevant scales (d≤50d≤ 50, dispatch intervals ≥1≥\!1 min). Broader outlook. Three physics-deep research directions emerge. First, nonlinear MF coupling may induce collective phenomena: symmetry breaking, mode locking, or phase-transition-like dynamics in multi-modal targets, well-studied in statistical physics but unexplored in generative modeling. Second, MF-PID aligns naturally with the emerging framework of “sampling decisions” [33], which combines diffusion, transformer, and reinforcement learning under a unified stochastic control umbrella (see also last chapter of [34]). Third, the identification of samples with agents positions MF-PID as a bridge between passive data synthesis and active coordination of engineered systems under real-time constraints. Methods Theoretical framework. Full derivations of the MF-PID equations in terminal-cost and SOT formulations, the LQG closed-form reduction, and the proof of Theorem (11) are provided in SI §1–§3. The Hopf–Cole linearisation, Green-function machinery, and piecewise-constant (PWC) analytic formulas follow [17, 18, 19]. Simulations. Scalar DR scenarios (Scenarios A and B): B=8,000B=8,000 particles, nsteps=2,500n_steps=2,500 Euler–Maruyama steps on [0,1][0,1]. Multi-zone scalability experiments (d-sweep, K-sweep, AR coupling): B=4,000B=4,000–6,0006,000 particles, nsteps=2,500n_steps=2,500. In all cases the MF guidance is set analytically via (11). The β-schedule is geometrically decreasing: βj=β0γj−1 _j= _0γ^j-1 with β0=12.0 _0=12.0, γ=0.65γ=0.65, M=8M=8 intervals. Demand-response parameters. Scalar (d=1d=1) baseline. Target mixture: π=(0.60,0.40)π=(0.60,0.40), m(tar)=(0.0,1.5)m^(tar)=(0.0,1.5), σ(tar)=(0.20,0.30)σ^(tar)=(0.20,0.30). Scenario A: m(in)=(1.0,6.0)m^(in)=(1.0,6.0), σ(in)=(3.0,3.0)σ^(in)=(3.0,3.0). Scenario B: m(in)=(1.5,5.5)m^(in)=(1.5,5.5), σ(in)=(0.5,0.7)σ^(in)=(0.5,0.7). Units: 0 = = 20°C, one unit = = 3°C deviation. Global means: m¯A(in)=2.8 m^(in)_A=2.8, m¯B(in)=2.3 m^(in)_B=2.3, m¯(tar)=0.6 m^(tar)=0.6; MF guidance νt(MF)=2.8−2.2t _t^(MF)=2.8-2.2t (A) and 2.3−1.7t2.3-1.7t (B). d-sweep. Zone heterogeneity encoded via zj=sin(2πj/d)z_j= (2π j/d); target mode means m0(tar)=0.1⋅+0.15zm^(tar)_0=0.1·1+0.15z, m1(tar)=1.5⋅−0.15zm^(tar)_1=1.5·1-0.15z; initial means displaced by +1.5+1.5 (occupied) and +4.0+4.0 (unoccupied); diagonal covariances. d∈1,2,4,8,16,32d∈\1,2,4,8,16,32\; B=4,000B=4,000. K-sweep. d=4d=4; K target component means uniformly spaced on [−1,2][-1,2]; initial means == target means +4+4; weights ∝(K,K−1,…,1) (K,K-1,…,1). By construction m¯(tar)=0 m^(tar)=0 for all K, so the two IA baselines (ν≐0ν 0 and ν≐m¯(tar)ν m^(tar)) are identical; the comparison reduces to MF vs. a single constant-guidance IA. AR(1) coupling. d=8d=8, K=2K=2, Scenario B parameters; covariance Σij(k)=σk2ρ|i−j| ^(k)_ij= _k^2\,ρ^|i-j|, ρ∈0.0,0.5,0.8ρ∈\0.0,0.5,0.8\. Code is available at https://github.com/mchertkov/MeanFieldPID. Acknowledgements The author thanks the University of Arizona start-up programme for financial support. This work was initiated during sabbatical visits to the University of Michigan Institute for Computational Discovery and Engineering, the International Centre for Theoretical Physics (ICTP), the Technische Universität Ilmenau (Humboldt Fellowship), Lawrence Livermore National Laboratory (faculty mini-sabbatical program), and KAIST Graduate School of AI. Scientific engagement and encouragement from colleagues at all five institutions are gratefully acknowledged. Large language models (Claude, Anthropic; ChatGPT, OpenAI) assisted with text editing and code refactoring; all mathematical derivations, scientific claims, and code were independently verified by the author. Supplementary Information Overview. This Supplementary Information (SI) provides full mathematical derivations and implementation details supporting the main text. SI A derives the MF-PID governing equations in both terminal-cost and SOT formulations. SI B presents the complete LQG closed-form reduction, including the scalar TCL example and explicit performance metrics. SI C develops the theory for the quadratic interaction potential: it first reduces MF-PID to a self-consistent guided H-PID, then proves the central result that the MF guidance is the exact linear interpolant between initial and target means (Theorem C.2), and finally assembles the fully explicit Gaussian-mixture score function and marginal density. SI D provides additional experimental diagnostics for the demand-response application. SI E recalls the independent-agent PID foundation [17, 18, 19] used throughout. Appendix A MF-PID: Governing Equations A.1. Agent dynamics and mean-field limit Consider N exchangeable agents with dynamics dxt(i)=(ft(xt(i))+ut(i))dt+dWt(i),x0(i)=0,t∈[0,1],dx_t^(i)= (f_t(x_t^(i))+u_t^(i) )dt+dW_t^(i), x_0^(i)=0, t∈[0,1], (14) where xt(i)∈ℝdx_t^(i) ^d, ftf_t is a pre-specified base drift, ut(i)∈ℝdu_t^(i) ^d is the control, and Wt(i)W_t^(i) are independent standard Brownian motions. The per-agent density pt(i)p_t^(i) satisfies the Kolmogorov–Fokker–Planck (KFP) equation ∂tpt(i)+∇⋅(pt(i)(ft+ut(i)))=12Δpt(i),p0(i)=δ(x). _tp_t^(i)+∇· (p_t^(i)(f_t+u_t^(i)) )= 12 p_t^(i), p_0^(i)=δ(x). (15) As N→∞N→∞ pt(N)≐N−1∑ipt(i)p_t^(N) N^-1 _ip_t^(i) converges to the mean-field marginal ptp_t governed by ∂tpt+∇⋅(pt(ft+ut))=12Δpt,p0=δ(x), _tp_t+∇· (p_t(f_t+u_t) )= 12 p_t, p_0=δ(x), (16) where ut=ut(x,pt(⋅))u_t=u_t(x,p_t(·)). A representative agent then obeys the McKean–Vlasov SDE (main text, eq. (2)). A.2. Formulation A: Terminal-cost MF-PID The mean-field cost-to-go is Jt(x,pt)=infut→1[∫t1(12‖ut′‖2+∫Vt′(xt′−y)pt′(y)dy)dt′+φ(x1)|eqs. (16), (14),xt=x],J_t(x,p_t)= _u_t→ 1E\! [ _t^1\!\! ( 12\|u_t \|^2+ V_t (x_t -y)\,p_t (y)\,dy )dt + (x_1)\; |\;eqs.~ SI:eq:kfp-mf, SI:eq:sde-agent,\,x_t=x ], (17) where Vt(x−y)V_t(x-y) is an interaction potential and φ is a terminal cost. We work in the classical mean-field control setting where the representative agent solves a control problem with coefficients depending on the population law ptp_t; at equilibrium ptp_t is generated by the optimal control (self-consistency). Then the mean-field HJB equation is −∂tJt=∫Vt(x−y)pt(y)dy⏟Vt(eff)(x)+ft(x)⋅∇xJt+12ΔxJt−12‖∇xJt‖2,J1=φ.- _tJ_t= V_t(x-y)\,p_t(y)\,dy_V^(eff)_t(x)+f_t(x)· _xJ_t+ 12 _xJ_t- 12\| _xJ_t\|^2, J_1= . (18) Hopf–Cole linearization (exact - not an approximation). Setting Jt=−logψtJ_t=- _t in (18) yields the mean-field linear HJB (MF-lin-HJB) equation: −∂tψt+Vt(eff)(x)ψt=ft(x)⋅∇xψt+12Δxψt,ψ1=e−φ.- _t _t+V^(eff)_t(x)\, _t=f_t(x)· _x _t+ 12 _x _t, _1=e^- . (19) Remark A.1 (Linear conditional on ptp_t). (19) is linear in ψ for a given effective potential Vt(eff)V^(eff)_t, but the overall MF system remains nonlinear through the self-consistency Vt(eff)(x)=∫Vt(x−y)pt(y)dyV^(eff)_t(x)= V_t(x-y)\,p_t(y)\,dy. The optimal control and the optimal marginal density are ut∗(x) u_t^*(x) =∇xlogψt(x), = _x _t(x), (20) ∂tpt∗+∇⋅(pt∗(ft+ut∗)) _tp_t^*+∇· (p_t^*(f_t+u_t^*) ) =12Δpt∗,p0∗=δ(x). = 12 p_t^*, p_0^*=δ(x). (21) Green-function representation. Because (19) and (21) are (quasi-) linear in ψ and p respectively, their solutions can be expressed via effective Green functions Gt(−;(eff))(x;y)G_t^(-;(eff))(x;y) and Gt(+;(eff))(x;y)G_t^(+;(eff))(x;y) satisfying t∈[1→0]:−∂tGt(−;(eff))+Vt(eff)(x)Gt(−;(eff))=ft(x)⋅∇xGt(−;(eff))+12ΔxGt(−;(eff)), t∈[1→ 0]: - _tG_t^(-;(eff))+V^(eff)_t(x)\,G_t^(-;(eff))=f_t(x)· _xG_t^(-;(eff))+ 12 _xG_t^(-;(eff)), (22) G1(−;(eff))(x;y)=δ(x−y), 85.35826ptG_1^(-;(eff))(x;y)=δ(x-y), t∈[0→1]:∂tGt(+;(eff))+Vt(eff)(x)Gt(+;(eff))=−∇x⋅(ft(x)Gt(+;(eff)))+12ΔxGt(+;(eff)), t∈[0→ 1]: _tG_t^(+;(eff))+V^(eff)_t(x)\,G_t^(+;(eff))=- _x· (f_t(x)G_t^(+;(eff)) )+ 12 _xG_t^(+;(eff)), (23) G0(+;(eff))(x;y)=δ(x−y). 85.35826ptG_0^(+;(eff))(x;y)=δ(x-y). Then ψt(x)=∫e−φ(y)Gt(−;(eff))(x;y)dy,pt∗(x)=Gt(+;(eff))(x;0)ψt(x)ψ0(0). _t(x)= e^- (y)\,G_t^(-;(eff))(x;y)\,dy, p_t^*(x)=G_t^(+;(eff))(x;0)\, _t(x) _0(0). (24) The MF coupling enters through Vt(eff)(x)=∫Vt(x−y)pt(y)dyV^(eff)_t(x)= V_t(x-y)\,p_t(y)\,dy, which depends on ptp_t itself. Eqs. (21)–(23) therefore constitute a self-consistent system: in general and unlike the independent-agent case, the Green functions cannot be computed independently of the population density. Remark A.2 (Heads up: Quadratic Isotropic Potential). We will see below in Section C.2 that in the case of a quadratic potential Vt(x)=βtx2/2V_t(x)= _tx^2/2 and zero basic drift ft=0f_t=0, the MF reduces to guided H-PID with explicit guidance – an independent-agent case with the linear-in-time interpolant between initial mean and target mean taken as a guidance νt _t within H-PID with the potential Vt(x)=βt(x−νt)2/2V_t(x)= _t(x- _t)^2/2. A.3. Formulation B: Stochastic Optimal Transport (SOT) Replace the terminal cost φ with a hard marginal constraint p1=p(tar)p_1=p^(tar). The optimal control becomes ut∗(x)=∇xlog∫p(tar)(y)Gt(−;(eff))(x;y)G1(+;(eff))(y;0)dy,u_t^*(x)= _x \!\!p^(tar)(y)\, G_t^(-;(eff))(x;y)G_1^(+;(eff))(y;0)\,dy, (25) and the density evolves under (21) with p0∗=δ(⋅)p_0^*=δ(·) and p1∗=p(tar)p_1^*=p^(tar). The governing system is (22)–(23) with (25) substituted into (21). Equivalently, pt∗(x)∝ψt(x)Gt(+;(eff))(x;0),p_t^*(x) _t(x)G_t^(+;(eff))(x;0), where ψt _t solves the backward equation (19). Remark A.3 (Comparison with mean-field Schrödinger bridges). The mean-field Schrödinger problem studied in [35, 24] corresponds to ft=0f_t=0, Vt=0V_t=0 with endpoint constraints. Our formulation allows an arbitrary base drift ftf_t and an explicit running interaction potential Vt(x−y)V_t(x-y), and does not reduce to the previously studied mean-field Schrödinger systems. Appendix B LQG MF-PID: Closed-Form Reduction B.1. Model specification We specialize to: • Linear drift: ft(x)=Atx+Btf_t(x)=A_tx+B_t, At∈ℝd×dA_t ^d× d, Bt∈ℝdB_t ^d. • Quadratic interaction: Vt(x)=12x⊤QtxV_t(x)= 12x Q_tx, Qt⪰0Q_t 0. • Gaussian target: p(tar)(x)=(x;m1(tar),Σ1(tar))p^(tar)(x)=N(x;\,m_1^(tar), _1^(tar)), Σ1(tar)≻0 _1^(tar) 0. Assuming p0=δ(⋅)p_0=δ(·), the controlled SDE (14) with linear optimal control preserves Gaussianity: pt=(x;mt,Σt)p_t=N(x;\,m_t, _t) for all t. B.2. Effective potential and quadratic ansatz Substituting the Gaussian ptp_t into Vt(eff)V^(eff)_t: Vt(eff)(x)=12(x−mt)⊤Qt(x−mt)+12Tr(QtΣt).V^(eff)_t(x)= 12(x-m_t) Q_t(x-m_t)+ 12Tr(Q_t _t). (26) The trace term is x-independent and does not affect ∇xJt _xJ_t or the control. We seek a quadratic cost-to-go Jt(x)=12x⊤Stx+st⊤x+rtJ_t(x)= 12x S_tx+s_t x+r_t, giving the affine optimal control ut∗(x)=−Stx−st.u_t^*(x)=-S_tx-s_t. (27) B.3. The closed ODE system Substituting (26) and the quadratic ansatz into the MF-HJB (18) and matching polynomial terms in x yields the following closed system. Proposition B.1 (LQG MF-PID reduction). Under the LQG model specification of §B.1, the optimal mean-field cost-to-go is quadratic and the four coupled ODEs LQG MF-PID Equations −S˙t - S_t =Qt+StAt+At⊤St−St2, =Q_t+S_tA_t+A_t S_t-S_t^2, (S2.1) −s˙t - s_t =−Qtmt+At⊤st−Stst+StBt, =-Q_tm_t+A_t s_t-S_ts_t+S_tB_t, (S2.2) m˙t m_t =(At−St)mt+Bt−st, =(A_t-S_t)m_t+B_t-s_t, (S2.3) Σ˙t _t =(At−St)Σt+Σt(At−St)⊤+I, =(A_t-S_t) _t+ _t(A_t-S_t) +I, (S2.4) with boundary conditions m0=0m_0=0, Σ0=0 _0=0, m1=m1(tar)m_1=m_1^(tar), Σ1=Σ1(tar) _1= _1^(tar), characterise the unique optimal control via ut∗(x)=−Stx−stu_t^*(x)=-S_tx-s_t. Structural decoupling. Equations (S2.1)–(S2.4) separate into two sequential sub-problems: 1. Variance block (St,Σt)(S_t, _t): Equations (S2.1) and (S2.4) depend only on QtQ_t and AtA_t, and are identical in the MF and IA cases. They are solved first as a two-point boundary-value problem in Σ (shooting on S1S_1). 2. Mean block (mt,st)(m_t,s_t): Given StS_t, equations (S2.2)–(S2.3) are linear in (mt,st)(m_t,s_t) with boundary conditions m0=0m_0=0, m1=m1(tar)m_1=m_1^(tar). The MF and IA cases differ here: the MF source term is −Qtmt-Q_tm_t (coupling to the current population mean), whereas the IA source is −Qtm¯-Q_t m (a fixed exogenous centre). This decoupling is the key structural result: MF coordination operates entirely through the mean block. B.4. Scalar TCL example: complete closed-form solution Set d=1d=1, At=−κA_t=-κ (κ>0κ>0), Bt=0B_t=0, Qt=qQ_t=q (q≥0q≥ 0), p(tar)=(m(tar),(σ(tar))2)p^(tar)=N(m^(tar),(σ^(tar))^2). B.4.1 Step 1: Riccati equation With Δ≐κ2+q κ^2+q, the scalar Riccati equation S˙t=St2+2κSt−q S_t=S_t^2+2κ S_t-q is solved by St=−κ+Δ1+ρe−2Δ(1−t)1−ρe−2Δ(1−t),S_t=-κ+ \, 1+ρ\,e^-2 (1-t)1-ρ\,e^-2 (1-t), (28) where ρ=(S1+κ−Δ)/(S1+κ+Δ)ρ=(S_1+κ- )/(S_1+κ+ ). B.4.2 Step 2: Variance matching (bridge constraint) The covariance ODE Σ˙t=−2(κ+St)Σt+1 _t=-2(κ+S_t) _t+1, Σ0=0 _0=0 has the explicit solution Σ1=1−ρ1−ρr0⋅1−r02Δ,r0=e−2Δ. _1= 1-ρ1-ρ r_0· 1-r_02 , r_0=e^-2 . (29) Imposing Σ1=(σ(tar))2 _1=(σ^(tar))^2 yields the closed-form shooting parameter: ρ=A−1Ar0−1,A=2Δ(σ(tar))21−r0,ρ= A-1Ar_0-1, A= 2 (σ^(tar))^21-r_0, (30) and then S1=−κ+Δ(1+ρ)/(1−ρ)S_1=-κ+ (1+ρ)/(1-ρ). B.4.3 Step 3: Mean dynamics With StS_t determined, eliminating sts_t from (S2.2)–(S2.3) yields m¨t−κ2mt=0 m_t-κ^2m_t=0, whose solution satisfying m0=0m_0=0, m1=m(tar)m_1=m^(tar) is mt(MF)=m(tar)sinh(κt)sinh(κ).m_t^(MF)=m^(tar)\, (κ t) (κ). (31) The linear coefficient is then st=−m˙t−(κ+St)mt=−κm(tar)cosh(κt)sinh(κ)−(κ+St)m(tar)sinh(κt)sinh(κ).s_t=- m_t-(κ+S_t)\,m_t=- κ\,m^(tar) (κ t) (κ)-(κ+S_t) m^(tar) (κ t) (κ). (32) B.4.4 Step 4: IA baseline For the IA family with exogenous centre m¯ m, (S2.2) becomes s˙t(IA)=qm¯+(κ+St)st(IA) s_t^(IA)=q m+(κ+S_t)s_t^(IA), with solution st(IA)=g(t)(s0(IA)+qm¯J(t)),s_t^(IA)=g(t) (s_0^(IA)+q m\,J(t) ), (33) where g(t)=eΔt(1−ρe−2Δ)/(1−ρe2Δ(t−1))g(t)=e t(1-ρ e^-2 )/(1-ρ e^2 (t-1)) and J(t)=(1−e−2Δt−ρe−2Δ(e2Δt−1))/(2Δ(1−ρe−2Δ))J(t)=(1-e^-2 t-ρ e^-2 (e^2 t-1))/(2 (1-ρ e^-2 )). The initial condition s0(IA)s_0^(IA) is fixed by the bridge constraint m1(IA)=m(tar)m_1^(IA)=m^(tar): s0(IA)=−m(tar)+qm¯B~(ρ)A~(ρ),A~(ρ)=g(1)−1∫01g(u)2du,B~(ρ)=g(1)−1∫01g(u)2J(u)du.s_0^(IA)=- m^(tar)+q m\, B(ρ) A(ρ), A(ρ)=g(1)^-1\! _0^1\!\!g(u)^2\,du, B(ρ)=g(1)^-1\! _0^1\!\!g(u)^2J(u)\,du. (34) B.4.5 Performance metrics Three metrics differentiate MF from IA: μu(t) _u(t) =−(Stmt+st),σu2(t)=St2σt2 =-(S_tm_t+s_t), _u^2(t)=S_t^2 _t^2 (control mean and variance), (control mean and variance), (35) (t) (t) =St2σt2+(Stmt+st)2 =S_t^2 _t^2+(S_tm_t+s_t)^2 (instantaneous power), (instantaneous power), (36) ℰ(t) (t) =∫0t(u)du = _0^tP(u)\,du (cumulative energy). (cumulative energy). (37) Since StS_t and σt2 _t^2 are identical across schemes, the energy difference arises exclusively from the linear coefficient sts_t. The MF controller achieves strictly lower ℰ(1)E(1) than any exogenously centred IA strategy. Appendix C MF-PID with Quadratic Interaction Potential C.1. Reduction to self-consistent guided H-PID We now consider MF-PID in the case of an isotropic quadratic interaction potential Vt(x)=βt2‖x‖2,βt>0,V_t(x)= _t2\|x\|^2, _t>0, (38) with arbitrary base drift ftf_t and arbitrary initial and target distributions. Computing the effective potential (17) gives Vt(eff)(x)=∫Vt(x−y)pt(y)dy=βt2‖x−mt‖2+βt2TrΣt,V^(eff)_t(x)= V_t(x-y)\,p_t(y)\,dy= _t2\|x-m_t\|^2+ _t2Tr _t, (39) where mt=pt[x]m_t=E_p_t[x] is the population mean and Σt=Covpt[x] _t=Cov_p_t[x]. The trace term is x-independent and does not affect the control. Defining the guidance centre νt≐mt _t m_t, the effective potential reduces to Vt(eff)(x)=βt2‖x−νt‖2+constV^(eff)_t(x)= _t2\|x- _t\|^2+const, which is the quadratic guidance potential of the H-PID framework [17, 19]. Proposition C.1 (MF-PID as self-consistent guided H-PID). For the quadratic potential (38), MF-PID is equivalent to a guided H-PID whose guidance trajectory νt=νt(MF) _t= _t^(MF) is determined endogenously by the self-consistency condition νt(MF)=x∼pt∗(⋅|ν(MF))[x],t∈[0,1]. _t^(MF)=E_x p_t^*(\,·\,|\,ν^(MF))[x], t∈[0,1]. (40) In the SOT formulation, the corresponding optimal control is ut∗(x)=∇xlog∫p(tar)(y)Gt(−)(x;y|ν(MF))G1(+)(y;0|ν(MF))dy,u_t^*(x)= _x \!\!p^(tar)(y)\, G_t^(-)(x;y\,|\,ν^(MF))G_1^(+)(y;0\,|\,ν^(MF))\,dy, (41) where Gt(±)(⋅|ν(MF))G_t^(±)(\,·\,|\,ν^(MF)) are the Green functions of the H-PID system evaluated under the guidance ν(MF)ν^(MF). C.2. Linearity of the MF guidance The following theorem is the central analytical result for the zero-drift case. Theorem C.2 (Linear MF guidance). Let ft≐0f_t 0, let VtV_t be the quadratic potential (38) with an arbitrary schedule βt>0 _t>0, and let p(in)p^(in), p(tar)p^(tar) be any probability measures with finite first moments m¯(in)=p(in)[x] m^(in)=E_p^(in)[x] and m¯(tar)=p(tar)[x] m^(tar)=E_p^(tar)[x]. Also assume that the controlled process is initialized with x0∼p(in)x_0 p^(in) (i.e. p0=p(in)p_0=p^(in)), and that the MF fixed point exists and yields finite first moments mt=pt∗[x]m_t=E_p_t^*[x] for all t. Then νt(MF)=mt=(1−t)m¯(in)+tm¯(tar)for all t∈[0,1]. _t^(MF)=m_t=(1-t)\, m^(in)+t\, m^(tar) all t∈[0,1]. (42) That is, the self-consistent MF guidance is the exact linear interpolant between initial and target global means, independently of βt _t – for any measurable βt>0 _t>0 such that the MF bridge is well-posed and ‖xt‖<∞E\|x_t\|<∞ for all t – and of the shape of p(in)p^(in) and p(tar)p^(tar). Proof. We show that m¨t=0 m_t=0. Step 1 (Itô + mean acceleration). With ft≐0f_t 0, the McKean–Vlasov SDE is dxt=ut∗(xt)dt+dWtdx_t=u_t^*(x_t)\,dt+dW_t. Differentiating mt=[xt]m_t=E[x_t] gives m˙t=[ut∗(xt)] m_t=E[u_t^*(x_t)]. Applying Itô’s formula to ut∗(xt)u_t^*(x_t) and taking expectations: m¨t=[∂tut∗+(ut∗⋅∇x)ut∗+12Δut∗]. m_t=E\! [ _tu_t^*+(u_t^*· _x)u_t^*+ 12 u_t^* ]. (43) Step 2 (Differentiated HJB in space). With Jt=−logψtJ_t=- _t and ut∗=∇xlogψtu_t^*= _x _t, the HJB equation (18) with ft=0f_t=0 reads ∂tJt=−Vt(eff)+12|ut∗|2+12Δlogψt. _tJ_t=-V^(eff)_t+ 12|u_t^*|^2+ 12 _t. Taking the gradient in x and using the symmetry of ∇xut∗=Hess(logψt) _xu_t^*=Hess( _t): ∂tut∗=∇xVt(eff)−(∇xut∗)ut∗⏟=(ut∗⋅∇x)ut∗−12Δut∗. _tu_t^*= _xV^(eff)_t- ( _xu_t^*)\,u_t^*_=(u_t^*· _x)u_t^*- 12 u_t^*. (44) Step 3 (Cancellation). Substituting (44) into (43), the (ut∗⋅∇x)ut∗(u_t^*· _x)u_t^* and 12Δut∗ 12 u_t^* terms cancel exactly, leaving m¨t=[∇xVt(eff)(xt)]. m_t=E\! [ _xV^(eff)_t(x_t) ]. This identity holds for any interaction potential; the complex score structure of ut∗u_t^* does not appear. Step 4 (Quadratic potential + self-consistency). For the quadratic potential (39): ∇xVt(eff)(x)=βt(x−νt(MF)). _xV^(eff)_t(x)= _t(x- _t^(MF)). Therefore m¨t=βt[xt−νt(MF)]=βt(mt−νt(MF)). m_t= _t\,E[x_t- _t^(MF)]= _t\,(m_t- _t^(MF)). At the MF fixed point, the self-consistency condition (40) gives νt(MF)=mt _t^(MF)=m_t, so m¨t=βt(mt−mt)=0. m_t= _t\,(m_t-m_t)=0. Conclusion. With m¨t=0 m_t=0 and boundary conditions m0=m¯(in)m_0= m^(in), m1=m¯(tar)m_1= m^(tar), we obtain mt=(1−t)m¯(in)+tm¯(tar)m_t=(1-t)\, m^(in)+t\, m^(tar), and since νt(MF)=mt _t^(MF)=m_t, the claim (42) follows. ∎ Remark C.1 (Independence of βt _t, initial law, and target). The proof uses only three ingredients: the Itô–HJB cancellation (Step 3), which is a structural identity holding for any smooth ut∗u_t^*; the linearity of ∇xVt(eff) _xV^(eff)_t in x (Step 4), which is a consequence of the quadratic form of VtV_t; and the self-consistency relation νt(MF)=mt _t^(MF)=m_t. The βt _t schedule, the shape of p(in)p^(in) and p(tar)p^(tar) enter only through the score function ut∗u_t^*, which has already canceled out. Theorem C.2 therefore holds for arbitrary p(in)p^(in) and p(tar)p^(tar) and for any βt _t resulting in a well-defined densities. Corollary C.3 (Linearity for delta initial condition). When p(in)=δ(⋅)p^(in)=δ(·) (i.e. m¯(in)=0 m^(in)=0), Theorem C.2 gives νt(MF)=tm¯(tar) _t^(MF)=t\, m^(tar). Remark C.2 (Relation to the Brownian and Schrödinger bridge). Setting Vt≐0V_t 0 (i.e. βt=0 _t=0) in the H-PID framework removes the interaction potential entirely, recovering the standard Schrödinger bridge [7]: the problem of finding the most likely path of a Brownian motion that transports p(in)p^(in) to p(tar)p^(tar). When both marginals are delta distributions, this further specialises to the classical Brownian bridge (a Brownian motion pinned at both endpoints). In both cases the mean trajectory mt=[xt]m_t=E[x_t] is the linear interpolant (1−t)m¯(in)+tm¯(tar)(1-t) m^(in)+t m^(tar): for the Brownian bridge this is immediate from the explicit formula [xt]=(1−t)x0+tx1E[x_t]=(1-t)x_0+tx_1; for the general Schrödinger bridge it follows from the affine structure of the Doob h-transform (see [7], Remark 1.8). Theorem C.2 is a strictly stronger statement. It establishes the same linear interpolation for any βt>0 _t>0, however large or time-varying, and for arbitrary non-Gaussian, multi-modal marginals p(in)p^(in) and p(tar)p^(tar). In the limit βt→0 _t→ 0 our result trivially reproduces the Schrödinger bridge case (the mean acceleration equation m¨t=βt(mt−νt(MF))→0 m_t= _t(m_t- _t^(MF))→ 0 degenerately), but the proof mechanism is completely different: it does not rely on the h-transform or any special structure of the bridge kernel. Instead, it rests on the Itô–HJB cancellation (Step 3 of the proof), which is a structural identity for the interacting system, and holds precisely because the quadratic potential makes ∇xVt(eff) _xV^(eff)_t linear in x. The non-interacting Schrödinger bridge therefore corresponds to the special limit of our theorem in which the interaction is switched off, not the other way around. A related observation holds in the stochastic interpolant framework [9], where the base interpolant It=(1−t)x0+tx1+σzI_t=(1-t)x_0+tx_1+σ z has mean (1−t)[x0]+t[x1](1-t)E[x_0]+tE[x_1] by construction; in MF-PID with Vt>0V_t>0, the same linearity is not imposed but derived from the Itô–HJB cancellation, and holds for any βt _t and any initial and final densities. The following corollary recovers the LQG result of §B.4 as a special case and identifies the precise role of the base drift. Corollary C.4 (OU base drift). For ft(x)=Axf_t(x)=Ax with A=−κIA=-κ I and p0=δ(⋅)p_0=δ(·), the same Itô–HJB cancellation in Steps 1–3 applies. More generally, for a base drift ftf_t (under sufficient regularity), Steps 1–3 yield the identity m¨t=dt[ft(xt)]+[∇Vt(eff)(xt)]−[(∇ft(xt))⊤ut∗(xt)]. m_t= ddtE[f_t(x_t)]+E[∇ V^(eff)_t(x_t)]-E[(∇ f_t(x_t)) u_t^*(x_t)]. For linear ft(x)=Axf_t(x)=Ax (constant Jacobian ∇ft=A∇ f_t=A) and quadratic interaction, [∇Vt(eff)]=βt(mt−νt(MF))E[∇ V^(eff)_t]= _t(m_t- _t^(MF)) vanishes at the MF fixed point. Using m˙t=Amt+[ut∗(xt)] m_t=Am_t+E[u_t^*(x_t)] (hence [ut∗(xt)]=m˙t−AmtE[u_t^*(x_t)]= m_t-Am_t) gives m¨t=(A−A⊤)m˙t+A⊤Amt. m_t=(A-A ) m_t+A A\,m_t. Specializing to A=−κIA=-κ I yields m¨t=κ2mt m_t=κ^2m_t and hence, with m0=0m_0=0 and m1=m¯(tar)m_1= m^(tar), νt(MF)=mt(MF)=m¯(tar)sinh(κt)sinh(κ). _t^(MF)=m_t^(MF)= m^(tar)\, (κ t) (κ). (45) This departs from the linear interpolant by O(κ2)O(κ^2); Eq. (31) is recovered. Remark C.3 (Practical implication). Theorem C.2 provides the MF guidance in closed form, without iteration: for ft≐0f_t 0, one sets νt(MF)=(1−t)m¯(in)+tm¯(tar) _t^(MF)=(1-t)\, m^(in)+t\, m^(tar) and proceeds directly to the score function computation of §C.4. The self-consistency fixed-point iteration is therefore unnecessary in this case and is included only for reference and for the ft≐̸0f_t 0 setting. C.3. Green-function solution: PWC protocol Following [18, 19], we discretise [0,1][0,1] by 0=t0<t1<⋯<tM=10=t_0<t_1<·s<t_M=1 and represent the protocol by piecewise-constant (PWC) values (βi,νi)( _i, _i) on each interval [ti,ti+1)[t_i,t_i+1). By Theorem C.2, for ft≐0f_t 0: νi=(1−ti)m¯(in)+tim¯(tar),i=0,…,M. _i=(1-t_i)\, m^(in)+t_i\, m^(tar), i=0,…,M. (46) The guided Green functions take the Gaussian form Gt(−)(x|y) G_t^(-)(x|y) ∝exp(−at(−)2∥x−νt∥2+bt(−)(x−νt)⊤(y−νt) \! (- a_t^(-)2\|x- _t\|^2+b_t^(-)(x- _t) (y- _t) −ct(−)2∥y−νt∥2+(rt(−))⊤(x−νt)+(st(−))⊤(y−νt)), 56.9055pt- c_t^(-)2\|y- _t\|^2+(r_t^(-)) (x- _t)+(s_t^(-)) (y- _t) ), (47) Gt(+)(y|0) G_t^(+)(y|0) ∝exp(−at(+)2‖y−νt‖2+(st(+))⊤(y−νt)). \! (- a_t^(+)2\|y- _t\|^2+(s_t^(+)) (y- _t) ). (48) The scalar coefficients satisfy the Riccati equations ∓a˙t(±)+βt=(at(±))2,b˙t(−)=at(−)bt(−),c˙t(−)=(bt(−))2.∓ a_t^(±)+ _t=(a_t^(±))^2, b_t^(-)=a_t^(-)b_t^(-), c_t^(-)=(b_t^(-))^2. (49) Within each PWC interval these admit closed hyperbolic forms (§C.5). C.4. Score function in closed form Proposition C.5 (Gaussian-mixture score). For a Gaussian-mixture target p(tar)=∑k=1Kπk(x;mk,Σk),p^(tar)= _k=1^K _kN(x;\,m_k, _k), and guidance νt _t (given by (46) for ft≐0f_t 0), substituting (47)–(48) into the SOT formula (25) yields the score function ut∗(x)=bt(−)(y^(t;x)−Υt(x)),u_t^*(x)=b_t^(-) ( y(t;x)- _t(x) ), (50) where Υt(x) _t(x) =νt+at(−)(x−νt)−rt(−)bt(−), = _t+ a_t^(-)(x- _t)-r_t^(-)b_t^(-), (51) y^(t;x) y(t;x) =∑k=1Kπ¯k(t;x)m¯k(t;x), = _k=1^K π_k(t;x)\, m_k(t;x), (52) m¯k(t;x) m_k(t;x) =(Σk−1+KtI)−1(Σk−1mk+Ktμt(x)), = ( _k^-1+K_tI )^-1 ( _k^-1m_k+K_t\, _t(x) ), (53) π¯k(t;x) π_k(t;x) =πk(μt(x);mk,Σk+Kt−1I)∑ℓπℓ(μt(x);mℓ,Σℓ+Kt−1I), = _k\,N( _t(x);\,m_k,\, _k+K_t^-1I) _ _ \,N( _t(x);\,m_ ,\, _ +K_t^-1I), (54) with probe distribution parameters Kt=ct(−)−a1(+),μt(x)=bt(−)Kt(x−νt)+st(−)−s1(+)Kt+νt.K_t=c_t^(-)-a_1^(+), _t(x)= b_t^(-)K_t(x- _t)+ s_t^(-)-s_1^(+)K_t+ _t. (55) C.5. Closed-form PWC evolution within each interval Fix interval i with βt=βi _t= _i, νt=νi _t= _i, and set ωi=βi _i= _i, τ=t−tiτ=t-t_i. Scalar coefficients. Forward branch (a(+)a^(+), propagated forward from t=0+t=0^+ where a(+)∼1/ta^(+) 1/t): a(+)(t)=ωicoth(ωiτ+φi(+)),a^(+)(t)= _i ( _iτ+ _i^(+)), (56) with φ1(+)=0 _1^(+)=0 and subsequent phases set by continuity. Backward branch (a(−)a^(-), propagated backward from t=1−t=1^- where a(−)∼1/(1−t)a^(-) 1/(1-t)). On the terminal interval: a(−)(t)=c(−)(t)=ωM−1coth(ωM−1(1−t)),b(−)(t)=ωM−1csch(ωM−1(1−t)).a^(-)(t)=c^(-)(t)= _M-1 ( _M-1(1-t)), b^(-)(t)= _M-1\,csch( _M-1(1-t)). (57) On earlier interval i, with right-endpoint anchors (ai+1(−),bi+1(−),ci+1(−))(a^(-)_i+1,b^(-)_i+1,c^(-)_i+1) and τ=ti+1−tτ=t_i+1-t: a(−)(t) a^(-)(t) =ωiai+1(−)+ωitanh(ωiτ)ωi+ai+1(−)tanh(ωiτ), = _i\, a^(-)_i+1+ _i ( _iτ) _i+a^(-)_i+1 ( _iτ), (58) b(−)(t) b^(-)(t) =bi+1(−)βi−(a(−)(t))2βi−(ai+1(−))2, =b^(-)_i+1 _i-(a^(-)(t))^2 _i-(a^(-)_i+1)^2, (59) c(−)(t) c^(-)(t) =ci+1(−)+(bi+1(−))2βi−(ai+1(−))2(ai+1(−)−a(−)(t)). =c^(-)_i+1+ (b^(-)_i+1)^2 _i-(a^(-)_i+1)^2 (a^(-)_i+1-a^(-)(t) ). (60) Vector (linear) coefficients. Introduce the re-centred linear coefficients θt(+)=st(+)+at(+)νt,θx,t(−)=rt(−)+(at(−)−bt(−))νt,θy,t(−)=st(−)+(ct(−)−bt(−))νt. _t^(+)=s_t^(+)+a_t^(+) _t, _x,t^(-)=r_t^(-)+(a_t^(-)-b_t^(-)) _t, _y,t^(-)=s_t^(-)+(c_t^(-)-b_t^(-)) _t. (61) The ODEs for these quantities are linear with hyperbolic driving terms from (56)–(60). Their closed-form solutions on each interval are: θ(+)(t) θ^(+)(t) =sinhφi(+)sinh(ωiτ+φi(+))θ(+)(ti)+βiωicosh(ωiτ+φi(+))−coshφi(+)sinh(ωiτ+φi(+))νi, = _i^(+) ( _iτ+ _i^(+))\,θ^(+)(t_i)+ _i _i\, ( _iτ+ _i^(+))- _i^(+) ( _iτ+ _i^(+))\, _i, (62) θx(−)(t) _x^(-)(t) =b(−)(t)bi+1(−)θx(−)(ti+1−)+(a(−)(t)−b(−)(t)bi+1(−)ai+1(−))νi, = b^(-)(t)b^(-)_i+1\, _x^(-)(t_i+1^-)+ (a^(-)(t)- b^(-)(t)b^(-)_i+1a^(-)_i+1 ) _i, (63) θy(−)(t) _y^(-)(t) =θy(−)(ti+1−)+ci+1(−)−c(−)(t)bi+1(−)θx(−)(ti+1−)+((bi+1(−)−b(−)(t))−ai+1(−)ci+1(−)−c(−)(t)bi+1(−))νi. = _y^(-)(t_i+1^-)+ c^(-)_i+1-c^(-)(t)b^(-)_i+1\, _x^(-)(t_i+1^-)+ ((b^(-)_i+1-b^(-)(t))-a^(-)_i+1 c^(-)_i+1-c^(-)(t)b^(-)_i+1 ) _i. (64) C.6. Time-marginal density in closed form Proposition C.6 (Optimal marginal density). Define the derived time-continuous quantities Kt=ct(−)−a1(+),αt=bt(−)/Kt,d¯t=(θy,t(−)−θ1(+))/Kt,Sk(t)=Σk+Kt−1I,K_t=c_t^(-)-a_1^(+), _t=b_t^(-)/K_t, d_t=( _y,t^(-)- _1^(+))/K_t, S_k(t)= _k+K_t^-1I, (65) and Mk(t) M_k(t) =(at(+)+at(−)−(bt(−)2/Kt))I+αt2Sk(t)−1, = (a_t^(+)+a_t^(-)-(b_t^(-)^2/K_t) )I+ _t^2S_k(t)^-1, (66) hk(t) h_k(t) =θt(+)+θx,t(−)+bt(−)d¯t+αtSk(t)−1(mk−d¯t). = _t^(+)+ _x,t^(-)+b_t^(-) d_t+ _tS_k(t)^-1(m_k- d_t). (67) Then the optimal marginal density is the Gaussian mixture pt∗(x)∝∑k=1Kπk|Sk(t)|−1/2exp(−12x⊤Mk(t)x+hk(t)⊤x−12(mk−d¯t)⊤Sk(t)−1(mk−d¯t)).p_t^*(x) _k=1^K _k|S_k(t)|^-1/2 \! (- 12x M_k(t)x+h_k(t) x- 12(m_k- d_t) S_k(t)^-1(m_k- d_t) ). (68) C.7. Non-zero initial condition and mixture-to-mixture transport C.7.1 Deterministic start x0=zx_0=z A coordinate shift x~=x−z x=x-z reduces the problem to a zero-start problem with shifted guidance ν~t=νt−z ν_t= _t-z and shifted target means m~k=mk−z m_k=m_k-z. The scalar Riccati coefficients (at(±),bt(−),ct(−))(a_t^(±),b_t^(-),c_t^(-)) are unchanged; the linear coefficients acquire z-linear corrections via three scalar shift propagators λ(+)(t)λ^(+)(t), λx(−)(t) _x^(-)(t), λy(−)(t) _y^(-)(t) satisfying the same ODEs as (θ(+),θx(−),θy(−))(θ^(+), _x^(-), _y^(-)) but with the source βtνt _t _t replaced by βt _t: θ~t(+)=θt(+)−λ(+)(t)z,θ~x,t(−)=θx,t(−)+λx(−)(t)z,θ~y,t(−)=θy,t(−)+λy(−)(t)z. θ_t^(+)= _t^(+)-λ^(+)(t)\,z, θ_x,t^(-)= _x,t^(-)+ _x^(-)(t)\,z, θ_y,t^(-)= _y,t^(-)+ _y^(-)(t)\,z. (69) The shift propagators share the same PWC closed forms as (62)–(64) with νi _i replaced by 11. C.7.2 Stochastic initial condition x0∼p(in)x_0 p^(in) For a mixture initial condition p(in)(x)=∑j=1Jπj(in)(x;mj(in),Σj(in))p^(in)(x)= _j=1^J _j^(in)N(x;\,m_j^(in), _j^(in)), sample z(i)∼p(in)z^(i) p^(in) per trajectory and apply the shift transformation above. By Theorem C.2 with ft≐0f_t 0, the global mean evolves as mt=(1−t)m¯(in)+tm¯(tar),m_t=(1-t)\, m^(in)+t\, m^(tar), (70) so the guidance is fully explicit. The shifted scalar coefficients and propagators are computed once; per-trajectory cost is (d)O(d) additions. The optimal marginal density is the Gaussian mixture with J×KJ× K components pt∗(x)∝∑j=1J∑k=1Kπ~j,k(t)(x;m~j,k(t),Sk(t)),p_t^*(x) _j=1^J _k=1^K π_j,k(t)\,N(x;\, m_j,k(t),\,S_k(t)), (71) where the weights π~j,k π_j,k and means m~j,k m_j,k are given by the following Gaussian completing-the-square formulas. Component weights and means. Define the effective precision Λj(t)=(Σj(in))−1+(at(+)−a1(+))I _j(t)=( _j^(in))^-1+(a_t^(+)-a_1^(+))I and the linear vector ℓj(t)=(Σj(in))−1mj(in)+θt(+)−a1(+)mk−θ1(+) _j(t)=( _j^(in))^-1m_j^(in)+ _t^(+)-a_1^(+)m_k- _1^(+). Then: m~j,k(t) m_j,k(t) =(Mk(t)+Λj(t))−1(hk(t)+Λj(t)−1((Σj(in))−1mj(in)+θt(+)−a1(+)mk−θ1(+))), = (M_k(t)+ _j(t) )^-1 (h_k(t)+ _j(t)^-1 (( _j^(in))^-1m_j^(in)+ _t^(+)-a_1^(+)m_k- _1^(+) ) ), (72) π~j,k(t) π_j,k(t) ∝πj(in)πk|Λj(t)|−1/2|Σj(in)|−1/2|Sk(t)|−1/2|M~j,k(t)|−1/2exp(12h~j,k⊤M~j,k−1h~j,k−q~j,k(t)), _j^(in) _k\,| _j(t)|^-1/2| _j^(in)|^-1/2|S_k(t)|^-1/2| M_j,k(t)|^-1/2 \! ( 12 h_j,k M_j,k^-1 h_j,k- q_j,k(t) ), (73) where M~j,k=Mk+Λj M_j,k=M_k+ _j, h~j,k=hk+Λj−1((Σj(in))−1mj(in)+θt(+)−a1(+)mk−θ1(+)) h_j,k=h_k+ _j^-1(( _j^(in))^-1m_j^(in)+ _t^(+)-a_1^(+)m_k- _1^(+)), and q~j,k(t)=12(mk−d¯t)⊤Sk−1(mk−d¯t)+12(mj(in))⊤(Σj(in))−1mj(in) q_j,k(t)= 12(m_k- d_t) S_k^-1(m_k- d_t)+ 12(m_j^(in)) ( _j^(in))^-1m_j^(in). Appendix D Demand-Response Application: Additional Diagnostics D.1. Experimental parameters See Table 3. Table 3: Numerical parameters for the demand-response simulations. Parameter Value Number of particles B 8,000 Euler–Maruyama steps nstepsn_steps 2,500 PWC intervals M 8 β-schedule (β0,γ)( _0,γ) (12.0, 0.65)(12.0,\;0.65) MF tolerance ϵε 2×10−42× 10^-4 Target mixture (π,m(tar),σ(tar))(π,m^(tar),σ^(tar)) (0.60,0.40)(0.60,0.40); (0.0,1.5)(0.0,1.5); (0.20,0.30)(0.20,0.30) Scenario A (σ(in),m(in))(σ^(in),m^(in)) (3.0,3.0)(3.0,3.0); (1.0,6.0)(1.0,6.0) Scenario B (σ(in),m(in))(σ^(in),m^(in)) (0.5,0.7)(0.5,0.7); (1.5,5.5)(1.5,5.5) Both scenarios use ft≐0f_t 0. By Theorem C.2, the MF guidance is therefore the exact linear interpolant νt(MF)=(1−t)m¯(in)+tm¯(tar) _t^(MF)=(1-t)\, m^(in)+t\, m^(tar) in both cases, with no iteration required. The small deviations max|ν(MF)−νlin| |ν^(MF)- _lin| reported in the numerical diagnostics (0.078 for Scenario A, 0.030 for Scenario B; Fig. 13) are an artefact of the PWC temporal discretisation: within each of the M=8M=8 intervals the guidance is held constant, so the piecewise-constant approximation to a linear function is not exactly linear at the midpoints. These residuals shrink as M→∞M→∞ and are negligible for all energy comparisons. D.2. Scenario A: additional diagnostics Figures 5–7 show density snapshots, trajectory ensembles, and per-component energy decomposition for Scenario A. Figure 5: Scenario A: density snapshots at t∈0.10,0.30,0.50,0.70,1.0t∈\0.10,0.30,0.50,0.70,1.0\ for the three methods. The initially broad, overlapping distribution contracts and splits into the two target modes. Black solid: target ρ(tar)ρ^(tar); black dashed: initial ρ(in)ρ^(in). Figure 6: Scenario A: trajectory ensembles (50 sample paths, colored by initial component). Blue: occupied; orange: unoccupied. Black curve: ensemble mean; gray band: ±1σ± 1σ. Shaded rectangles at t=1t=1: target ±σk(tar)± _k^(tar) bands. Figure 7: Scenario A: per-component analysis. (a) Per-component mean trajectories: the unoccupied mode (orange) must travel ∼4.5 4.5 units while the occupied mode (blue) travels ∼1.0 1.0. (b) Cumulative per-component energy: the unoccupied mode absorbs most of the control cost; MF coordination reduces this asymmetry. D.3. Scenario B: additional diagnostics Figures 8–10 provide the same diagnostics for Scenario B. Figure 8: Scenario B (narrow initial). The MF energy saving increases to 22.6%. Figure 9: Scenario B: density snapshots. Unlike Scenario A, the initial law is already bimodal; the two peaks translate and sharpen simultaneously. Figure 10: Scenario B: trajectory ensembles (50 sample paths, colored by component). The two clusters remain well separated throughout the bridge. D.4. Score-field affinity The LQG benchmark (SI §2) predicts an optimal drift affine in x: u∗(t,x)=−S(t)x−s(t)u^*(t,x)=-S(t)\,x-s(t). For Gaussian-mixture targets the drift is no longer exactly affine; we quantify the departure by fitting u∗≈−S(t)x−s(t)u^*≈-S(t)\,x-s(t) via weighted least squares at each time slice and reporting the coefficient of determination R2(t)R^2(t). Figure 11 overlays the spring constant S(t)S(t), shift s(t)s(t), and affinity R2(t)R^2(t) for all three methods. Figure 11: Scenario A: score-field affine decomposition u∗(t,x)≈−S(t)x−s(t)u^*(t,x)≈-S(t)\,x-s(t). (a) Spring constant S(t)S(t): identical across methods, set by β(t)β(t) and ρ(tar)ρ^(tar). (b) Shift s(t)s(t): differs between methods; |sMF||s^MF| is largest at early times. (c) Affinity R2(t)R^2(t): drops from ∼1.0 1.0 to ∼0.03 0.03, marking the onset of genuinely non-affine Gaussian-mixture score structure. Vertical gray lines: β-interval boundaries. D.5. Cross-scenario comparison Figures 12–14 summarise the comparison between scenarios. Figure 12: Cross-scenario energy comparison. (Left) Total control energy ℰ(1)E(1) for all three methods in Scenarios A and B. (Right) Cumulative energy curves. The MF saving increases from 11.6% (A) to 22.6% (B). Energy. The absolute energy is roughly halved from Scenario A to B (the narrow initial law starts closer to the target in Wasserstein distance), but the relative MF advantage approximately doubles: 11.6% (A) vs. 22.6% (B). This amplification arises because, in Scenario B, the two modes are quasi-independent: the unoccupied cluster must travel ≈ 4≈\,4 units while the occupied cluster moves only ≈ 1.5≈\,1.5, and the self-consistent guidance of Theorem C.2 adapts to this asymmetry more effectively than any single constant ν. Figure 13: Self-consistent MF guidance ν(MF)(t)ν^(MF)(t) for Scenarios A (blue) and B (orange), with their respective linear interpolants (dotted). By Theorem C.2, the two curves coincide exactly in the continuous-time limit; the small residuals (max|ν(MF)−νlin|=0.078 |ν^(MF)- _lin|=0.078 for A, 0.0300.030 for B) are due to the PWC temporal discretisation (M=8M=8 intervals). MF guidance. Figure 13 confirms Theorem C.2 numerically: the self-consistent guidance ν(MF)(t)ν^(MF)(t) lies essentially on the linear interpolant between the initial and target global means for both scenarios. The residuals are smaller in Scenario B, where the modes overlap less and the PWC approximation to a linear profile is more accurate. Figure 14: Convergence of the fixed-point iteration to the linear interpolant. Left: Scenario A (12 iterations to tolerance 2×10−42× 10^-4). Right: Scenario B (10 iterations). (a) max|Δν| | ν| vs. iteration (dashed: tolerance). (b) Residual νi(MF)−νilin _i^(MF)- _i^lin at each PWC midpoint, confirming Theorem C.2 in the PWC limit. Appendix E Independent-Agent PID: Background This section summarises the independent-agent PID construction [17, 19] that forms the foundation for MF-PID. E.1. IA PID in the SOT formulation In the absence of inter-agent coupling, the cost-to-go satisfies the standard HJB equation, and after the Hopf–Cole substitution the optimal control is ut(IA)(x)=∇xlog∫p(tar)(y)Gt(−)(x;y)G1(+)(y;0)dy,u_t^(IA)(x)= _x \!\!p^(tar)(y)\, G_t^(-)(x;y)G_1^(+)(y;0)\,dy, (74) where Gt(±)G_t^(±) satisfy (22)–(23) with Vt(eff)V^(eff)_t replaced by a prescribed potential VtV_t (no ptp_t dependence). The Green functions can therefore be computed independently of the population density—this is the key structural simplification that MF breaks. E.2. H-PID with quadratic guidance potential For ft≐0f_t 0 and guided quadratic potential Vt(x)=βt2‖x−νt‖2V_t(x)= _t2\|x- _t\|^2, the Green functions take the Gaussian forms (47)–(48) with scalar Riccati coefficients satisfying (49). The score function for a Gaussian-mixture target is given by Proposition C.5 with a prescribed guidance νt _t. The PWC closed forms of §C.5 then yield a fully explicit, training-free generative model. MF-PID specialises this to the endogenously determined guidance of Theorem C.2. E.3. Cross-validation identity Proposition E.1 (Cross-validation identity). The optimal score satisfies ut∗(x)=∇xlogpt∗(x)−∇xlogGt(+)(x;0),u_t^*(x)= _x p_t^*(x)- _x G_t^(+)(x;0), (75) relating the optimal control, the optimal marginal density, and the forward Green function. When Vt=0V_t=0 and ftf_t is pre-trained to yield p(tar)p^(tar) as a stationary law, G1(+)(x;0)=p(tar)(x)G_1^(+)(x;0)=p^(tar)(x), and (75) implies ut∗≐0u_t^* 0—confirming that no corrective control is needed when the drift already realises the target. References [1] Sohl-Dickstein, J., Weiss, E., Maheswaranathan, N. & Ganguli, S. Deep Unsupervised Learning using Nonequilibrium Thermodynamics. In Bach, F. & Blei, D. (eds.) Proceedings of the 32nd International Conference on Machine Learning, vol. 37 of Proceedings of Machine Learning Research, 2256–2265 (PMLR, Lille, France, 2015). URL https://proceedings.mlr.press/v37/sohl-dickstein15.html. [2] Ho, J., Jain, A. & Abbeel, P. Denoising Diffusion Probabilistic Models (2020). URL http://arxiv.org/abs/2006.11239. ArXiv:2006.11239 [cs, stat]. [3] Song, Y. et al. Score-Based Generative Modeling through Stochastic Differential Equations (2021). URL http://arxiv.org/abs/2011.13456. ArXiv:2011.13456 [cs, stat]. [4] Rezende, D. & Mohamed, S. Variational Inference with Normalizing Flows. In Bach, F. & Blei, D. (eds.) Proceedings of the 32nd International Conference on Machine Learning, vol. 37 of Proceedings of Machine Learning Research, 1530–1538 (PMLR, Lille, France, 2015). URL https://proceedings.mlr.press/v37/rezende15.html. [5] Papamakarios, G., Nalisnick, E., Rezende, D. J., Mohamed, S. & Lakshminarayanan, B. Normalizing flows for probabilistic modeling and inference. J. Mach. Learn. Res. 22 (2021). [6] Pavon, M. & Wakolbinger, A. On Free Energy, Stochastic Control, and Schrödinger Processes. In Modeling, Estimation and Control of Systems with Uncertainty, 334–348 (Birkhäuser Boston, Boston, MA, 1991). URL http://link.springer.com/10.1007/978-1-4612-0443-5_22. [7] Léonard, C. A survey of the Schrödinger problem and some of its connections with optimal transport (2013). URL http://arxiv.org/abs/1308.0215. ArXiv:1308.0215 [math]. [8] Chen, Y., Georgiou, T. T. & Pavon, M. Optimal Transport Over a Linear Dynamical System. IEEE Transactions on Automatic Control 62, 2137–2152 (2017). URL http://ieeexplore.ieee.org/document/7549018/. [9] Albergo, M. S. & Vanden-Eijnden, E. Building Normalizing Flows with Stochastic Interpolants (2023). URL http://arxiv.org/abs/2209.15571. ArXiv:2209.15571 [cs, stat]. [10] Brambati, M., Celani, A., Gherardi, M. & Ginelli, F. Learning to flock in open space by avoiding collisions and staying together (2026). URL http://arxiv.org/abs/2506.15587. ArXiv:2506.15587 [cond-mat]. [11] Callaway, D. S. Tapping the energy storage potential in electric loads to deliver load following and regulation, with application to wind energy. Energy Conversion and Management 50, 1389–1400 (2009). URL http://dx.doi.org/10.1016/j.enconman.2008.12.012. [12] Callaway, D. S. & Hiskens, I. A. Achieving controllability of electric loads. Proceedings of the IEEE 99, 184–199 (2011). [13] Mathieu, J. L., Kamgarpour, M., Lygeros, J., Andersson, G. & Callaway, D. S. Arbitraging Intraday Wholesale Energy Market Prices With Aggregations of Thermostatic Loads. IEEE Transactions on Power Systems 30, 763–772 (2015). [14] Beil, I., Hiskens, I. A. & Backhaus, S. Frequency Regulation From Commercial Building HVAC Demand Response. Proceedings of the IEEE 104, 745–757 (2016). [15] Borra, F., Cencini, M. & Celani, A. Optimal collision avoidance in swarms of active Brownian particles. Journal of Statistical Mechanics: Theory and Experiment 2021, 083401 (2021). URL http://arxiv.org/abs/2105.10198. ArXiv:2105.10198 [cond-mat]. [16] Kachar, K. G. & Gorodetsky, A. A. Dynamic multi-agent assignment via discrete optimal transport (2019). URL http://arxiv.org/abs/1910.10748. ArXiv:1910.10748 [cs]. [17] Behjoo, H. & Chertkov, M. Harmonic Path Integral Diffusion. IEEE Access 13, 42196–42213 (2025). URL https://ieeexplore.ieee.org/document/10910146/. [18] Chertkov, M. & Behjoo, H. Adaptive Path Integral Diffusion: AdaPID (2025). URL http://arxiv.org/abs/2512.11858. ArXiv:2512.11858 [cs]. [19] Chertkov, M. Generative Stochastic Optimal Transport: Guided Harmonic Path-Integral Diffusion (2025). URL http://arxiv.org/abs/2512.11859. ArXiv:2512.11859 [cs]. [20] Caluya, K. F. & Halder, A. Reflected Schrödinger Bridge: Density Control with Path Constraints (2020). URL http://arxiv.org/abs/2003.13895. ArXiv:2003.13895 [math]. [21] Teter, A. M. H., Chen, Y. & Halder, A. On the Contraction Coefficient of the Schrödinger Bridge for Stochastic Linear Systems. IEEE Control Systems Letters 7, 3325–3330 (2023). URL https://ieeexplore.ieee.org/document/10293168/. [22] Evensen, G. Data Assimilation: The Ensemble Kalman Filter (Springer, 2009), 2nd edn. [23] Backhoff-Veraguas, J., Conforti, G., Gentil, I. & Léonard, C. The mean field Schrödinger problem: ergodic behavior, entropy estimates and functional inequalities (2019). URL http://arxiv.org/abs/1905.02393. ArXiv:1905.02393 [math]. [24] Hernández, C. & Tangpi, L. Propagation of chaos for mean field Schrödinger problems (2024). URL http://arxiv.org/abs/2304.09340. ArXiv:2304.09340 [math]. [25] Teter, A. M. H., Wang, W. & Halder, A. Weyl Calculus and Exactly Solvable Schr\"odinger Bridges with Quadratic State Cost (2024). URL http://arxiv.org/abs/2407.15245. ArXiv:2407.15245 [math-ph, stat]. [26] Huang, M., Caines, P. E. & Malhame, R. P. Large-Population Cost-Coupled LQG Problems With Nonuniform Agents: Individual-Mass Behavior and Decentralized epsilon-Nash Equilibria. IEEE Transactions on Automatic Control 52, 1560–1571 (2007). [27] Bensoussan, A., Frehse, J. & Yam, S. C. P. Mean Field Games and Mean Field Type Control Theory. SpringerBriefs in Mathematics (Springer, 2013). [28] Hao, H., Sanandaji, B. M., Poolla, K. & Vincent, T. L. Aggregate flexibility of thermostatically controlled loads. IEEE Transactions on Power Systems 30, 189–198 (2015). [29] Grammatico, S., Gentile, B., Parise, F. & Lygeros, J. A Mean Field control approach for demand side management of large populations of Thermostatically Controlled Loads. In 2015 European Control Conference, ECC 2015 (2015). [30] Métivier, D. & Chertkov, M. Mean-field control for efficient mixing of energy loads. Physical Review E 101, 022115 (2020). URL https://link.aps.org/doi/10.1103/PhysRevE.101.022115. [31] Valenzuela, L. F., Williams, L. & Chertkov, M. Statistical mechanics of thermostatically controlled multizone buildings. Physical Review E 107, 034140 (2023). URL https://link.aps.org/doi/10.1103/PhysRevE.107.034140. [32] Tzen, B. & Raginsky, M. Theoretical guarantees for sampling and inference in generative models with latent diffusions (2019). URL http://arxiv.org/abs/1903.01608. ArXiv:1903.01608 [cs, math, stat]. [33] Chertkov, M., Ahn, S. & Behjoo, H. Sampling Decisions (2025). URL http://arxiv.org/abs/2503.14549. ArXiv:2503.14549 [cs]. [34] Chertkov, M. Mathematics of Generative AI (2025). URL https://github.com/mchertkov/Mathematics-of-Generative-AI-Book. [35] Backhoff-Veraguas, J., Bartl, D., Beiglböck, M. & Eder, M. Adapted Wasserstein Distances and Stability in Mathematical Finance. Finance and Stochastics 24, 601–632 (2020). URL https://doi.org/10.1007/s00780-020-00437-1.