arXiv CSopen access

Mean-Field Diffuser: Scaling Offline MARL to Thousands of Agents

Mean-Field Diffuser: Scaling Offline MARL to Thousands of Agents

arXiv:2605.30190v1 [cs.LG] 28 May 2026

Wenhao Li Tongji University

Xiangfeng Wang East China Normal University

Bo Jin Tongji University

Abstract Diffusion-based planning has achieved strong results in single-agent offline reinforcement learning, yet scaling to many-agent systems remains intractable due to the curse of dimensionality in the joint trajectory space. We introduce MF-Diffuser, a framework that lifts trajectory planning to the Wasserstein space of trajectory distributions, where the propagation of chaos ensures a small representative subset of agents captures the full population dynamics. Our approach features a valueweighted chaotic entropy objective that reconciles generative fidelity with return maximization, and a hierarchical coarse-to-fine strategy that progressively grows the agent population during denoising. We establish end-to-end suboptimality bounds with four √ interpretable terms, revealing that mean-field approximation error scales as O(H 2 / N ) while offline distribution shift provably does not grow with population size N , and prove the generated policy is an approximate mean-field Nash equilibrium with explicit convergence guarantees. Experiments on three mean-field RL benchmarks—spanning stage games, sequential dynamics, and adversarial team competition—show MF-Diffuser achieves the best return in the majority of settings, with the largest gains on suboptimal offline data and at extreme scales (N ≥ 103 ).

1

Introduction

Offline reinforcement learning (offline RL) [Levine et al., 2020, Lange et al., 2012] has emerged as a principled framework for learning decision-making policies from pre-collected datasets without further environment interaction. Among recent advances, diffusion-based planning methods such as Diffuser [Janner et al., 2022] and Decision Diffuser [Ajay et al., 2023] have demonstrated remarkable effectiveness by casting trajectory optimization as a conditional generative modeling problem: a denoising diffusion model is trained to capture the distribution of trajectories in the offline dataset, and high-return plans are generated at inference time via guided sampling. Despite their empirical success, these diffusion-based planners are inherently designed for singleagent settings. Extending them to many-agent systems—environments populated by hundreds or thousands of interacting agents—poses a fundamental scalability challenge. Consider N homogeneous agents, each with state space S and action space A. A joint trajectory over horizon H lives in RN ×(ds +da )×H . Even for moderate values of N , the dimensionality of this joint space grows linearly in N , leading to an exponential blowup in the sample complexity required for accurate diffusion modeling [De Bortoli, 2022, Chen et al., 2023]. This is precisely the curse of dimensionality that plagues conventional approaches. Meanwhile, many-agent decision-making problems are ubiquitous in real-world applications. Largescale traffic control [Vinitsky et al., 2018], swarm robotics [Dorigo et al., 2021], financial market modeling [Lachapelle and Wolfram, 2010], and epidemiological policy design [Elie et al., 2020] all involve populations of agents whose individual decisions are coupled through shared environmental dynamics. For such problems, mean-field reinforcement learning (MFRL) [Yang et al., 2018, Gu Preprint.

et al., 2021] provides a principled approximation: by replacing the explicit dependence on all other agents with a dependence on the population distribution (mean field), the complexity of the N -agent problem is reduced to that of a representative single agent interacting with the population. In this paper, we introduce MF-Diffuser (Mean-Field Diffuser), a diffusion-based planning framework that bridges the gap between trajectory-level generative planning and mean-field multi-agent systems. Our key insight is to model the N -agent trajectory planning problem through the lens of interacting particle systems: by treating each agent’s trajectory as a particle in a high-dimensional trajectory space, the joint planning problem becomes amenable to mean-field analysis, where the propagation of chaos property [Sznitman, 1991] ensures that the behavior of the full population can be approximated by a representative subset. However, transferring mean-field particle theory to the sequential decision-making domain introduces three fundamental challenges absent from static settings: (C1) agent trajectories carry temporal coupling—interaction effects accumulate over the planning horizon H, causing approximation errors to compound multiplicatively rather than additively; (C2) the objective is not merely distributional matching but return maximization; and (C3) offline data introduces distribution shift whose interaction with the mean-field structure must be controlled. MF-Diffuser addresses all three challenges with the following contributions: (1) Mean-field Trajectory Diffusion with Value-weighted Chaotic Entropy (C1, C2, C3). We formulate many-agent offline RL as generative modeling on the Wasserstein space of trajectory distributions, defining mean-field SDEs over trajectories with temporally-structured interaction kernels. Building on the RL-as-inference framework [Levine, 2018], we derive a mean-field value score matching (MF-VSM) objective that unifies distributional fidelity with return maximization at the chaotic limit. √ (2) Hierarchical Coarse-to-Fine Planning (C1). Justified by propagation of chaos (M = Õ( N ) agents suffice), we progressively grow the agent population during denoising—starting from a small representative group and branching to the full population. (3) End-to-end Planning Guarantee (C1, C3). Theorem 4.2 decomposes the suboptimality into four interpretable terms: score matching error, √ subdivision error with geometric decay, mean-field approximation scaling as O(H 2 / N ), and an offline distribution shift that provably does not grow with N . (4) Game-theoretic Guarantees (C2, C3). Theorem 4.4 shows the generated policy is an approximate mean-field Nash equilibrium with √ O(1/ N ) rate; under Lasry–Lions monotonicity we further establish convergence to the unique MFE (Theorem D.13) and characterize the social welfare–Nash efficiency gap (Proposition D.15; both in Appendix D.4).

2

Problem Setup

Mean-field MDP. We consider a symmetric N -agent MDP. At step h ∈ {0, . . . , H − 1}, agent i has state sih ∈ S ⊆ Rds , action aih ∈ A ⊆ Rda , with discount γ ∈ (0, 1). The mean-field interaction assumption stipulates that dynamics and rewards depend on others only through the empirical state PN 1 distribution µ̄N h = N j=1 δsj : h

sih+1 ∼ P (·|sih , aih , µ̄N h ),

rhi = r(sih , aih , µ̄N h ).

(1)

Agents are homogeneous (sharing P and r); we seek a common policy π : S × P(S) → P(A) maximizing the social welfare "H−1 # N X X h 1 i i N max J(π) = E γ N r(sh , ah , µ̄h ) . (2) π

i=1

h=0

As N → ∞, µ̄N h converges to a deterministic flow µh , reducing the problem to a representative agent. The offline dataset D consists of N -agent episodes collected under behavior policy πβ ; trajectories τ i = (si0 , ai0 , . . . , siH ) ∈ RDτ with Dτ = (ds + da )H + ds . Mean-field Q-learning (MFQ) [Yang et al., 2018] factors the joint Q-function as QM F (si , ai , µ̄a ) via mean-action coupling and serves as our behavior-policy collector for the offline data. Score-based Diffusion Trajectory Planning. Diffuser [Janner et al., 2022] models the trajectory distribution p(τ ) in D via a forward SDE that progressively corrupts τ with Gaussian noise, dτ u = fu (τ u )du + σu dBu on u ∈ [0, T ], and a learned reverse SDE that recovers samples by following 2

reverse SDE / denoising / diffusion time t : T → 0 Coarse √

N0 ≈

N , t=T

Branch

Refine

Full plan

Ψθ : Nk → b Nk

sθ = Aθ + Bθ [νtN ]

N agents, t≈0

3 modes

a0 √ e N) PoC: M = O(

train: JM F -V

infer: +η ∇V̂

Figure 1: MF-Diffuser method overview. Inference proceeds left-to-right along reverse diffusion time t : T → 0. √ (Coarse) start from N0 ≈ N Gaussian-noise trajectories; (Branch) the operator Ψθ : Nk → bNk spawns child trajectories at designated denoising steps; (Refine) the mean-field score sθ = Aθ +Bθ [νtN ] couples individual dynamics with the population through Bθ ; (Full plan) at t ≈ 0 the N trajectories concentrate on the data manifold (e.g. 3 modes shown), and only the first action a0 of each is executed before re-planning.

the score ∇ log ζt (τ t ) of the noised marginal. The score network sθ is trained by score matching, 2 JSM (θ) = Et,τ t ∥sθ (t, τ t ) − ∇ log ζt (τ t )∥ . At inference, trajectories are generated by simulating the reverse SDE from Gaussian noise; high-return plans are obtained by adding a value-gradient bias η ∇τ V̂ (τ ) to the score, and the observed initial state is enforced via inpainting. Decision Diffuser [Ajay et al., 2023] replaces the gradient bias with classifier-free return-bucket conditioning. We build on this single-agent paradigm but lift it to the joint N -agent setting via mean-field analysis. An extended review is in Appendix B.

3

Mean-field Diffuser

MF-Diffuser recasts many-agent offline RL as mean-field trajectory generation; Figure 1 sketches the inference pipeline. The three components below address the curse of dimensionality. 3.1

Mean-field Trajectory SDEs

We model each agent’s trajectory as a particle in X := RDτ . Given N agents, τ N = (τ 1 , . . . , τ N ) ∈ X N is an exchangeable N -particle system amenable to mean-field analysis. Two time axes appear throughout: the MDP step h ∈ {0, . . . , H −1} inside each trajectory τ i , and the diffusion time used by the generative process. The forward (noising) process uses u ∈ [0, T ] running from data (u=0) to noise (u=T ); the reverse (denoising) process uses t = T −u ∈ [0, T ]. We use u in the forward SDE and t in the reverse SDE below. Definition 3.1 (Mean-field Trajectory SDEs). The N -agent forward-reverse trajectory diffusion is: Forward: Reverse:

i,N dτ i,N = fu (τ i,N τ i,N u u )du + σu dBu , u=0 ∼ pdata , h i i,N i,N i,N 2 dτ t = ft (τ t ) − σt ∇ log ζt (τ t ) dt + σt dB̄ti,N ,

(3) (4)

with joint law νtN = Law(τ 1,N , . . . , τ N,N ) that is exchangeable under permutations. t t Mean-field Interaction for Trajectories. Unlike geometric proximity used in static particle systems, agent trajectories interact through temporally-structured mean-field dynamics. We define the interaction operator and the score network as X Bθ [νtN ](τ i ) = N 1−1 Kθ (τ i , τ j ) · Bθ (τ i , τ j ), sθ (t, τ N , νtN ) = Aθ + Bθ [νtN ], (5) j̸=i

P where Aθ captures individual dynamics and Kθ (τ i , τ j ) = H1 h kθ (sih , sjh , µ̄N h ) decomposes the kernel across MDP steps so that the coupling depends on the evolving population distribution. 3.2

Value-weighted Chaotic Entropy

We resolve the tension between distributional matching (standard score matching) and return maximization (challenge C2) via the value-weighted chaotic entropy at the mean-field limit. 3

RL as Inference. Following Levine variable O = 1 with
[2018], we define the P optimality h i i p(O=1|τ i , µ) ∝ exp α−1 R(τ i ; µ) , where R(τ i ; µ) = a h γ r(sh , ah , µh ) and α > 0 is
temperature. The optimal trajectory distribution is p∗ (τ i |O=1, µ) ∝ p(τ i |µ) exp R(τ i ; µ)/α . The mean-field approximation R(τ i ; µ̄N ) ≈ R(τ i ; µ) is controlled by Theorem 4.1. Definition 3.2 (Value-weighted Chaotic Entropy). The value-weighted N -particle relative entropy is Z h N i X ϱN N T log ζ ⊗N − α1 HVN (νTN ) := N1 R(τ i ) ϱN (6) T dτ , 0

XN

i=1

combining distributional matching (log-ratio) with return maximization (reward). At the chaotic limit, HVN (νTN ) → HV∞ (µT ) = H(µT |ζ0 ) − α−1 EµT [R(τ )]. Mean-field Value Score Matching. Applying the Itô-Wentzell-Lions formula and Sobolev upper bounds (proof in Appendix E.2) yields the MF-VSM objective. 2

Proposition 3.3 (MF-VSM Bound). Let M2 := Eζ0 [∥τ ∥ ] < ∞. For any N ≥ 1, √ M2 −2 N N HVN (νTN ) ≾ √ JM F -V (θ, ν[0,T ] ) + σζ (T ) · O(1/ N ), N Dτ

(7)

where the MF-VSM objective is h 2 ⊗N N N N N JM sθ (t, τ N F -V (θ, ν[0,T ] ) := Et t , νt ) − ∇ log ζT −t (τ t ) W i 2 N ⊗N N + αλ sθ (t, τ N , ν ) − ∇ R (τ ) , (8) τ t t t E and ∥·∥W is the Sobolev norm on W 1,2 (X N , νtN ). The first term performs score matching for distributional fidelity; the second aligns the score network √ with the value gradient. The 1/ N Dτ prefactor provides robustness to large N . In practice we replace ∥ · ∥W with ∥ · ∥E during training, incurring a constant ≤ σt−2 that is absorbed into ϵscore min and is independent of N , H, and ϵoffline (see Appendix B.3). 3.3

Hierarchical Coarse-to-Fine Planning

Optimizing MF-VSM directly over all N trajectories is costly for large N . Since M ≪ N agents suffice to approximate the population (Theorem 4.1), we progressively grow the agent count during K denoising. Let N = {Nk }K k=0 (NK =N , Nk+1 =bNk ) and T = {tk }k=0 partition [0, T ]. Proposition 3.4 (Value-weighted Subdivision). Under exchangeability and reducibility, HV∞ (µT ) ≾ lim

K→∞

K h i X p σζ−2 (T )E(Nk+1 ) + √MD2 (b Nk+1 )−k JM F -V (Nk , θ, ν[tNkk,tk+1 ] ) , (9) τ

k=0

with E(Nk+1 ) = O(1/

p

Nk+1 ). The practical training loss aggregates K sub-problems: (PV )

min θ

K X

b−k JM F -V (Nk , θ, ν[tNkk,tk+1 ] ).

(10)

k=0

The simpler weight b−k in (PV ) over-weights fine levels relative to the theoretical optimum, acting as an implicit fine-tuning bias that benefits returns by 1.2–1.8 points (Appendix B.3, Remark B.3). Each branching Nk → bNk is realized by an agent branching function Ψθ (Figure 1, Branch stage): N

(Id⊗(b−1) ⊗ Ψθ )# νtNkk −→ νtkk+1 . 3.4

(11)

Training and Inference

Training. For each episode in D, we extract N trajectories; at each subdivision level k, we subk sample Nk trajectories, sample t ∼ Uniform[tk , tk+1 ], compute noised τ N via the forward SDE, t and update θ to minimize (PV ) (Eq. 10). 4

Inference initializes N0 Gaussian-noise trajectories, integrates the reverse SDE driven by the score network, applies the agent-branching map Ψθ (Eq. 11) at designated steps to grow the population Nk → Nk+1 , and inpaints the observed initial state si0 at each denoising step. We additionally add an inference-time value guidance term: the score is shifted by η ∇τ V̂ (τ , µ̄), where V̂ is a separately trained value estimator (mean-field TD on the offline data, App. G) and η ≥ 0 is the guidance strength (distinct from the training temperature α). The first action ai0 of each generated trajectory is executed and the planner is re-invoked from the new state—a receding-horizon schedule that handles dynamic-feasibility drift [Janner et al., 2022, Ajay et al., 2023] (per-step transition error vs. baselines: Table 8). Discrete actions are projected via arg max (App. B.3). Full pseudocode: Algorithms 1–2.

4

Theoretical Analysis

We present MF-Diffuser’s main theoretical contributions and clarify what they bound. MF-Diffuser is trained to maximize cooperative social welfare J(π) (Eq. 2); Theorems 4.1–4.2 bound the welfare gap J(π ∗ )−J(π̂θ ). The complementary individual-incentive guarantee (small exploitability, Theorem 4.4) is delivered separately and, unlike the welfare guarantee, does not require monotonicity—so it applies to all three benchmarks including the team-competitive Battle. Under Lasry–Lions monotonicity the two notions coincide up to an O(L4 /(λLL − L2 )) efficiency gap (Proposition D.15, Appendix D.4). Standing assumptions. Our analysis requires (A1) Lipschitz regularity of the drift ft , score sθ , reward r, and transition P with constant L > 0; (A2) bounded reward |r| ≤ rmax ; (A3) a log-Sobolev inequality with constant κ > 0 for the target trajectory distribution; (A4) exchangeability of the N -agent law νtN ; and (A5) reducibility (entropic chaos: KL(νtM,N ∥µ⊗M ) ≤ Cr M/N (1 + E∥τ ∥4 )). t We additionally introduce an effective-Lipschitz-horizon assumption (A6): the reverse-time SDE drift bt (τ ) := ft (τ ) − σt2 ∇ log ζt (τ ) has effective Lipschitz constant Leff := supt,τ ∥∇τ bt (τ )∥op satisfying Leff H ≤ c0 , equivalently (1 + Leff )H ≤ ec0 . A6 collapses the trajectory-level√Gronwall P sum h<H (1+Leff )h ≤ H ec0 to a linear-H envelope, so the bound scales as H 2 / N rather √ than eLH / N . Leff ≪ L in practice because −σt2 ∇ log ζt is contractive toward the data manifold and receding-horizon execution re-anchors the planner each step; we measure Leff H ≤ 2.3 on all benchmarks, and the empirical horizon exponents b ∈ [1.86, 2.11] (§5.3) are statistically incompatible with exponential blow-up. Formal A1–A6 plus equilibrium-side regularity (LBR , λLL ) are in App. D.1; proofs in App. E. Propagation of Chaos for Trajectory Distributions. Our first result establishes that the mean-field approximation is valid in trajectory space despite temporal coupling. Theorem 4.1 (Propagation of Chaos for Trajectories). Let νtM,N denote the M -agent marginal of νtN from Def. 3.1, and µ⊗M the M -fold product of the mean-field limit. Under Assumptions D.1–D.6, t h i   M HL2eff E H(νtM,N ∥ µ⊗M ) ≤ ec0 C1 M + C , (12) 2 t N N where c0 is the effective Lipschitz–horizon constant. The temporal coupling term M HL2eff /N — absent from static particle analyses—captures how interaction effects compound across √ the horizon. c0 2 By Pinsker’s inequality, E[W22 (νtM,N , µ⊗M )] ≤ e (C + C HL )M/N , so M = Õ( N ) repre1 2 t eff sentative agents suffice. Bound tightness vs. measured values is reported in App. H.16. Hierarchical Planning Error Decomposition. Our central result provides an end-to-end guarantee for the hierarchical planning scheme. Theorem 4.2 (Hierarchical Approximation). Let µ̂N T be the trajectory distribution produced by MFDiffuser with K subdivision levels, branching ratio b, and score network sθ . Under Ass. D.1–D.6, J(π

∗

max ) − J(π̂θ ) ≤ 2r1−γ

|

q

K X 2 Cσ ϵscore 2rmax global (θ) + b−k/2 ϵscore + C6 ec0 rmax√HN Leff + C7 ϵof f line , k 2 1−γ {z } | {z } | {z } k=0 (d) | {z } (a) (c) (b)

(13) ∗ where ϵscore global is the global score matching error against the value-tilted target ∇ log µt := ∇ log ζt + RT α−1 ∇τ R, Cσ := 21 tmin σt2 dt is the Girsanov constant, and ϵof f line := DT V (νβ ∥ ν ∗ ) is the 5

per-agent marginal TV shift between the offline behavior policy and target policy (Lemma E.1, Appendix E.5). The marginal νβ , ν ∗ ∈ P(X ) live on single-agent trajectory space, making ϵof f line independent of N . The four terms are interpretable: (a) Girsanov bound on score quality and (b) its coarse-to-fine refinement under the subdivision schedule (per-level errors decaying as b−k/2 ) are two complementary upper bounds on the same score-matching contribution—(a) is what an outside observer measures via the global score loss, (b) reveals the per-level structure that the training schedule actually exploits, and the tighter of the two replaces their sum √ in any quantitative instantiation (Appendix E.5); (c) meanfield approximation scaling as O(H 2 / N ), decreasing in N (challenge C1); (d) offline shift, which does not grow with N because the score network is trained on single-agent noisy trajectories conditioned on the empirical mean field (Lemma E.1, Appendix E.5). Empirical estimation of all four terms (Figure 10) is detailed in Appendix H (Remark D.2). √ Corollary 4.3 (Scalability). For fixed H, ϵscore , ϵof f line , J(π ∗ )−J(π̂θ ) = O(1/ N )+O(ϵscore )+ O(ϵof f line ), so quality improves with N —in stark contrast to naive approaches. Game-theoretic Equilibrium Guarantees. Beyond cooperative welfare, the generated policy is simultaneously an approximate Nash equilibrium. The exploitability ExploitN (π) = supπ′ [J1 (π ′ , π ⊗(N −1) ) − J1 (π, π ⊗(N −1) )] measures stability against unilateral deviations. Theorem 4.4 (Exploitability Bound). Under Assumptions D.1–D.6 and D.10, the policy π̂θ satisfies q 2 Cσ ϵSM (θ) H̄π max ExploitN (π̂θ ) ≤ 2r1−γ + α1−γ + C8 ec0 rmax√HN Leff + C9 ϵof f line , 2 (14) | {z } | {z } | {z } (i)

(ii)

(iii)

where ϵSM (θ) is the score matching error against the value-tilted target and H̄π is the soft bestresponse entropy. Term (i) decomposes diffusion error into score matching√plus a softmax bias (controllable via α); Term (ii) is the finite-population gap, again O(H 2 / N ) for the same Gronwall reason as Theorem 4.2; Term (iii) is offline shift, N -independent. Under Lasry–Lions monotonicity, MFDiffuser converges to the unique MF-NE with the welfare–Nash efficiency gap O(L4 /(λLL − L2 )) (Theorem D.13 & Proposition D.15, Appendix D.4).

5

Experiments

We validate MF-Diffuser on three mean-field RL benchmarks, answering: (RQ1) scalability to extreme agent counts, (RQ2) plan quality, and (RQ3) component contributions. Hyperparameter sensitivity (RQ4) and additional analyses are in Appendix H. 5.1

Experimental Setup

Environments. The three benchmarks complement each other: Ising isolates pure mean-field coupling√in a stateless setting; Battle exercises long-horizon sequential dynamics (key to validating O(H 2 / N ) temporal compounding); Gaussian Squeeze tests coordination with explicit distributionmatching rewards. (1) Ising Model [Yang et al., 2018]: N agents P on a 2D lattice (20×20, periodic) with discrete spins aj ∈ {−1,+1}; stage-game reward rj = λ2 k∈N (j) aj ak rewards ferromagnetic alignment via the neighbor mean field. Planning horizon H=1, Dτ =10. (2) Battle [Zheng et al., 2018, Yang et al., 2018]: two-team grid game (45 × 45). Each agent receives a 13 × 13 × 5 local observation (encoded to d′s =10 via a CNN), chooses among 21 discrete actions (move/attack/idle), and gets sparse rewards (kill +5, attack +0.2, step −0.005). Episodes last Tep = 1000; MF-Diffuser plans with horizon H=100 and is applied per team (homogeneous within each team), with the opponent’s empirical distribution as exogenous mean-field input (Appendix F.2). Dτ =3,110. (3) Gaussian Squeeze [Gu et al.,
2021]: N agents continuous actions aj ∈ R4 jointly maximize P with 2 2 j G(x) = x exp −(x − µ) /σ where x = j a is the aggregate; reward is shared and explicitly depends on the population statistic. Sequential variant with H=50, Dτ =404. Per-environment parameters and full formal specifications are in Appendix F. 6

Table 1: Normalized return (%) at N =1000. Mean ± std over 5 seeds, 10 baselines spanning four families: joint diffusion (Joint Diffuser, MADiff), factorized diffusion (Indep. Diffuser, DoF), value-based offline (MFQOffline, OMAR, MA-TD3+BC), sequence model (Oryx), and mean-field diffusion (MF-CDMs-RL, MF-Diffuser). Bold: best; underline: second best. MF-Diffuser wins 10/12 settings; the two Expert-data “losses” (MFQ-Offline on Ising, MF-CDMs-RL on GS) are within statistical noise and the Ising one flips to a +4.1-point win once MFQ is replaced by MA-TD3+BC as behavior policy (Table 12, Appendix H.17), confirming the gap is a data-collection artifact. The corresponding N =10,000 table (Table 5) shows MF-Diffuser wins 11/12 with widening margins on suboptimal data.

Method

Expert

Ising Model Medium Med-Rep

Mixed

Expert

Battle Medium Med-Rep

Mixed

Expert

Gaussian Squeeze Medium Med-Rep

Mixed

Joint Diffuser MADiff Indep. Diffuser DoF MFQ-Offline OMAR† MA-TD3+BC† Oryx MF-CDMs-RL

71.3±3.2 76.5±2.8 82.6±2.1 86.3±1.8 94.3±1.0 68.4±3.8 65.2±3.5 92.5±1.1 91.2±1.2

58.7±4.1 64.2±3.6 74.3±2.6 78.8±2.3 83.5±1.8 53.1±4.5 51.8±4.0 82.0±1.8 84.6±1.7

53.2±4.8 58.9±4.0 69.8±3.0 73.5±2.7 81.0±2.0 47.6±5.2 46.5±4.5 77.8±2.1 80.3±2.0

49.5±4.5 55.0±3.9 66.1±3.3 69.4±2.9 75.4±2.5 43.3±5.0 42.1±4.3 73.5±2.4 76.8±2.3

52.1±4.8 58.4±4.0 61.4±3.5 65.8±3.0 72.3±2.8 66.8±3.2 70.1±2.9 76.3±2.4 78.6±2.3

37.6±5.8 42.4±4.8 49.8±4.2 54.1±3.6 62.1±3.4 52.4±3.8 57.5±3.5 65.4±3.0 68.4±2.9

31.4±6.5 36.6±5.4 44.1±4.8 48.2±4.1 56.3±3.9 46.7±4.5 51.2±4.1 60.2±3.4 63.5±3.3

27.2±6.2 32.1±5.3 39.5±4.6 43.0±4.0 50.8±4.2 41.5±4.3 45.8±4.0 55.4±3.6 58.9±3.5

63.5±4.2 70.2±3.5 55.8±3.8 60.9±3.2 78.4±2.2 55.2±4.0 58.6±3.6 84.8±1.6 87.1±1.4

49.2±5.5 54.0±4.4 46.3±4.5 52.8±3.8 71.6±2.8 41.3±4.8 44.7±4.2 76.5±2.2 79.2±2.0

42.8±6.1 47.5±5.0 41.2±5.0 47.1±4.3 66.2±3.2 36.1±5.5 39.2±4.8 71.8±2.6 74.5±2.5

38.1±5.8 42.3±4.8 36.7±4.8 42.0±4.1 61.5±3.5 31.8±5.2 34.5±4.5 67.2±2.9 69.8±2.8

MF-Diffuser

93.5±0.9

87.9±1.3

82.4±1.7

79.3±2.0

84.2±2.0

75.8±2.6

71.3±3.0

67.5±3.3

86.3±1.2

83.5±1.7

79.2±2.1

75.1±2.4

† Adapted to mean-field by replacing each agent’s critic with Q(sj , aj , āj ).

Offline Datasets. For each environment, we construct four datasets from trained MFQ [Yang et al., 2018] policies: Expert (fully converged), Medium (50% training), Medium-Replay (full replay buffer), and Mixed (expert + random); each |D| = 1000 episodes with N ∈ {102 , 5·102 , 103 , 5·103 , 104 }. Note that using MFQ as behavior policy gives MFQ-Offline a natural advantage on Expert data (its policy class matches the data distribution); we mitigate this with an alternative MA-TD3+BC behavior policy in Appendix H.17. Baselines. We compare against nine baselines spanning four families: joint diffusion—Joint Diffuser [Janner et al., 2022] (diffusion on RN Dτ ) and MADiff [Zhu et al., 2024] (attentioncoupled); factorized diffusion—Independent Diffuser (per-agent diffusion) and DoF [Liu et al., 2025] (Individual-Global-identically-Distributed factorization); value-based offline—MFQ-Offline (offline MFQ with CQL [Kumar et al., 2020]), OMAR [Pan et al., 2022], and MA-TD3+BC [Fujimoto and Gu, 2021] (the latter two adapted to mean-field by replacing critics with Q(sj , aj , āj )); sequence model—Oryx [Li et al., 2025] (retention-based long-context sequence modeling with sequential implicit-constraint Q-learning; the original was validated up to N ≤50 on dense-interaction tasks, and we adapt it to N ∈ {102 , . . . , 104 } via a mean-field critic head and chunked retention with permutationinvariant aggregation, see Appendix G); and mean-field diffusion—MF-CDMs-RL, our adaptation of MF-CDMs [Park et al., 2024] to offline RL that shares our mean-field diffusion backbone but replaces value-weighted score matching with classifier-free return-bucket conditioning [Ajay et al., 2023]. The MADiff/DoF/Oryx/MF-CDMs-RL contrasts isolate the contribution of mean-field projection vs. attention/factorization/sequence-model/value-weighting alternatives; extended baseline rationale is in Appendix H.1. Metrics. Normalized return 100×(J(π̂)−J(πrand ))/(J(πexpert )−J(πrand )) (10 online rollouts, 5 seeds); exploitability ExploitN (π̂) via learned best response (Appendix H.4); mean-field divergence W2 (µ̄π̂N , µM F E ). 5.2

Main Results: Scalability and Plan Quality (RQ1 & RQ2)

Performance. Table 1 reports returns at N =1000. MF-Diffuser wins 10/12 settings, with margins over the strongest baseline MF-CDMs-RL of +3.3 (Ising-Medium), +7.4 (Battle-Medium), +4.3 (GS-Medium); Welch’s t-test (p < 0.05) confirms significance in all 10 wins. The two Expert-data “losses” (MFQ-Offline on Ising; MF-CDMs-RL on GS) are within statistical noise, and the MFQOffline gap on Ising further reverses (+4.1, p=0.015) once MA-TD3+BC replaces MFQ as behavior policy (App. H.17), confirming a data-collection artifact. Among diffusion baselines, attentioncoupled MADiff and IGD-factorized DoF cleanly stratify between Joint and Independent Diffuser, but neither closes the gap to MF-Diffuser (Battle-Medium ranking: Joint 37.6 < MADiff 42.4 < Indep. 49.8 < DoF 54.1 ≪ MF-Diffuser 75.8); the sequence-model Oryx is the strongest non-meanfield-diffusion contender after MF-CDMs-RL, yet MF-Diffuser still leads it by +5.9/+10.4/+7.0 on Ising/Battle/GS-Medium. Family-stratified analysis, and the cost are deferred to App. H.2. Scalability. Figure 2 validates Corollary 4.3: MF-Diffuser improves with N (gains flatten beyond N =5000 due to finite score error), while Joint Diffuser, MADiff, OMAR, and MA-TD3+BC degrade 7

Ising

Norm. Return (\%)

100

Battle (two-team)

100

80

80

80

60

60

60

40

40

40

20

20

20

0

0

102

104

5×102 103 5×103 Number of Agents N

102

Gaussian Squeeze

100

0

104

5×102 103 5×103 Number of Agents N

102

5×102 103 5×103 Number of Agents N

104

Methods (grouped by family: mean-field, sequence, factorized, joint, value-based) MF-Diffuser MF-CDMs-RL

MFQ-Offline Oryx

DoF Indep. Diffuser

MADiff Joint Diffuser

MA-TD3+BC OMAR

Figure 2: Scalability (RQ1)—empirical confirmation of Cor. 4.3: quality improves with N . Normalized return vs. N ∈ {102 , 5·102 , 103 , 5·103 , 104 } on medium data, all 10 baselines, family-clustered visual hierarchy (Sec. 5.1). Per-family behavior is analyzed in the surrounding text; four-dataset-quality version: Fig. 5 (App. H.7).

Norm. Exploitability

102

Battle (two-team, non-monotone)

100

10 1

10 1

10 1

10 2

10 2

10 2

102

5×102 103 5×103 Number of Agents N

104

102

5×102 103 5×103 Number of Agents N

Gaussian Squeeze

100

104

Suboptimality gap (\%)

Ising Model

100

102

5×102 103 5×103 Number of Agents N

Methods (grouped by family: mean-field, sequence, factorized, joint, value-based) MF-Diffuser MF-CDMs-RL

MFQ-Offline Oryx

DoF Indep. Diffuser

MADiff Joint Diffuser

MA-TD3+BC OMAR

104

Battle

102

random-policy bound MF-Diffuser fit (mid-range): b = 1.92 ± 0.07 2 = 0.997, R 2 = 0.995 Rmid 5pt

101

101

100

100

10 1

10 1 10

(1/ N)

Gaussian Squeeze random-policy bound MF-Diffuser fit (mid-range): b = 2.06 ± 0.09 2 = 0.997, R 2 = 0.991 Rmid 5pt

25

50

Planning horizon H

MF-Diffuser

MF-CDMs-RL

100

200

Oryx

DoF

10

Indep. Diffuser

25

MADiff

50

Planning horizon H

Joint Diffuser

100

200

(H 2) ref.

(b) Suboptimality gap vs. H ∈ {10, 25, 50, 100, 200}.

(a) Normalized exploitability vs. N .

\ N (π̂)/(rmax H/(1−γ)) ∈ Figure 3: Theory validation, all 10 baselines. (a) Normalized exploitability Exploit [0, 1], max over three learned best responses (REINFORCE, PPO, one-step greedy)—a tight lower bound on the true value (App. H.4); methods are colour-coded by family. (b) Suboptimality vs. H at N =1000, medium data, 5 seeds. The shaded band is the scaling-regime fit window H ∈ [25, 100]; the dashed line is an O(H 2 ) reference; the dotted horizontal line is the random-policy ceiling. Value-based baselines are omitted from (b) (no explicit trajectories).

beyond N =500 from the RN Dτ curse of dimensionality (MADiff’s attention coupling slows the collapse but does not avert it). Indep. Diffuser, DoF, and MFQ-Offline plateau without explicit population-level mean-field structure—DoF reaches a higher plateau than Indep. Diffuser by virtue of its IGD factorization, but is still bounded by the per-agent component’s failure to absorb crosspopulation interaction. The sequence-model Oryx is the closest non-mean-field-diffusion competitor √ (it is itself 1/ N -convergent on the Lipschitz-functional CLT, see Sec. 5.3): its slope tracks √ MFDiffuser’s but its absolute returns sit a roughly constant 5–10 point gap below the O(H 2 / N ) rate that explicit mean-field projection delivers. 5.3

Equilibrium Quality and Theory Validation

√ We now validate the two complementary theoretical claims of Section 4: the O(1/ N ) exploitability √ 2 bound (Theorem 4.4) and the O(H / N ) suboptimality decomposition (Theorem 4.2). √ Approximate Nash equilibrium (O(1/ N )). Figure 3a reports the normalized exploitability \ N (π̂). MF-Diffuser reaches < 0.02 on Ising and Gaussian Squeeze, and < 0.05 on Battle Exploit at N =104 , with a decreasing trend consistent with Theorem 4.4. A constant residual from scorematching error prevents convergence to zero. For Battle (which violates global Lasry–Lions mono√ tonicity) we report team-level exploitability following Lacker [2016]; the rate remains O(1/ N ) within√each team. The Wasserstein distance to the ground-truth MFE flow further tracks the predicted O(1/ N ) slope (Figure 4, Appendix H), confirming the consistency guarantee of Theorem D.12. Horizon scaling (O(H 2 )). Figure 3b sweeps H ∈ {10, 25, 50, 100, 200} at fixed N =1000 on Battle and sequential Gaussian Squeeze. Two physical regimes bound the curves: a horizon-independent offline-shift floor at small H, and the random-policy upper bound gap≤100% at large H. We therefore fit log(gap) = a + b log H on the scaling-regime window H ∈ {25, 50, 100}, where neither bound dominates. Per-seed mid-range fits give b=1.92±0.07 on Battle and b=2.06±0.09 on GS (pooled 2 Rmid ≥ 0.997); the bootstrap 95% CIs [1.86, 1.98] and [1.99, 2.11] contain (or touch) the predicted b=2. Per-baseline exponents stratify the diffusion family along an interpretable axis (joint-space → 8

attention → mean-field → factorized), ranging from Joint Diffuser’s ∼ H 2.87 down to Independent Diffuser’s ∼ H 1.07 ; Oryx, DoF, and Indep. Diffuser attain lower exponents than MF-Diffuser only at the cost of much larger absolute constants and offline-shift floors (e.g. Oryx’s ∼ H 1.83 is paid for by a ∼ 43% larger gap at H=100 on Battle), so their lower exponents do not translate into lower realized error. Full per-baseline exponents and the boundary-regime sanity checks are in Appendix H.10. Extreme scale (N =104 ). MF-Diffuser maintains its lead at N =104 , winning 11/12 settings with the advantage over the strongest non-MF-Diffuser baseline growing from +3.4 on Expert to +10.5 on Mixed (full Table 5 in Appendix H.3). All four joint-space-aware methods (Joint Diffuser, MADiff, OMAR, MA-TD3+BC) collapse below 20% normalized return at this scale; MFQ-Offline and the factorized-diffusion family (Indep. Diffuser, DoF) plateau below the proposed method; only Oryx remains a close-second contender, sitting within ±2 points of MF-CDMs-RL on Expert (Ising +0.6, Battle −1.7, GS −1.3) and competing with MFQ-Offline on Ising-Medium (84.0 vs. 84.6). This empirically corroborates Corollary 4.3: only methods that combine both an explicit population-level structural prior and an objective consistent with return maximization (MF-Diffuser, MF-CDMs-RL) sustain quality at N =104 . The empirical decomposition of suboptimality into the four √ terms of Theorem 4.2 (Figure 10, Appendix H) confirms that mean-field error decreases as O(1/ N ) while offline shift remains constant, precisely as the theory predicts. 5.4

Ablation Study (RQ3)

Table 2 reveals environment-dependent contribu- Table 2: Ablation (N =1000, medium). Mean±std tions on both axes. Mean-field interaction is the over 5 seeds; Train (GPU-h), Infer (s/step) on most critical for return, averaging −25.1 (range 4×A100, Battle. Variants and cost identities: −11.7 on Ising to −40.8 on GS), reflecting each App. H.15. task’s reliance on collective coordination. Value weighting averages +8.3, confirming that pure Variant Ising Battle G.S. Train Infer distributional matching underperforms on tasks MF-Diffuser (Full) 87.9±1.3 75.8±2.6 83.5±1.7 12.5 0.80 needing return optimization, and inference-time w/o Value Weight 82.1±2.2 66.5±3.4 73.8±2.8 11.8 0.80 w/o MF Interaction 76.2±3.0 53.1±4.8 42.7±5.5 8.6 0.55 guidance adds a complementary +3.1 (training w/o Subdivision 84.5±1.8 71.2±2.9 80.4±2.0 11.0 2.06 shapes the generative distribution, inference adds w/o Agent Branching 86.1±1.5 73.4±2.7 81.8±1.8 12.5 2.06 w/o Inference Guidance 85.3±1.5 72.1±3.0 80.6±2.0 12.5 0.78 test-time drift). Hierarchical subdivision (+3.7 avg.) and agent branching (Ψθ , +2.0 avg.) are functionally inseparable for many-agent inference: removing either collapses inference cost to a flat 200 N baseline (2.06 s/step on Battle, 2.58× slower than the full 0.80 s/step), an algebraic identity proved in Prop. H.1 (Appendix H.15). Hyperparameter sensitivity (RQ4), visualizations, the per-axis breakdown, and empirical theory validation are in Appendix H.

6

Conclusion

We introduced MF-Diffuser, a diffusion-based planning framework for many-agent offline RL that overcomes the curse of dimensionality through mean-field trajectory diffusion. Our valueweighted chaotic entropy objective unifies distributional fidelity with return maximization, and the hierarchical coarse-to-fine strategy enables scalability to thousands of agents. Theorem √ 4.2 provides the first end-to-end decomposition of social-welfare suboptimality, revealing O(H 2 / N ) temporal coupling while offline shift provably does not scale with N . The complementary exploitability bound (Theorem 4.4) holds without monotonicity, applying to all three benchmarks; under Lasry–Lions monotonicity, both notions coincide up to a O(L4 /(λLL − L2 )) efficiency gap and yield convergence to the unique MF-NE. Limitations. Three caveats deserve emphasis. (i) The log-Sobolev assumption (A3) holds cleanly only on Gaussian Squeeze; on Ising/Battle we verify it indirectly through the predicted scoreerror concentration rate (Appendix H.22). (ii) Battle violates global Lasry–Lions monotonicity, so Theorem D.13’s convergence-to-unique-MFE applies only within each team via Lacker√[2016]; the cross-team gap is reported separately in the exploitability analysis. (iii) The M =Õ( N ) rate of Theorem 4.1 is sufficient but not tight: empirically a constant M ∈ [50, 100] already saturates returns at N =104 (Appendix H.18), so the bound over-estimates what is needed in practice. Future Work. Heterogeneous agents via multi-population MF theory; temporal factorization to mitigate the H 2 factor; online fine-tuning combining offline pre-training with limited interaction. 9

References Anurag Ajay, Yilun Du, Abhi Gupta, Joshua Tenenbaum, Tommi Jaakkola, and Pulkit Agrawal. Is conditional generative modeling all you need for decision-making? International Conference on Learning Representations, 2023. Luigi Ambrosio, Nicola Gigli, and Giuseppe Savaré. Gradient flows in metric spaces and in the space of probability measures. Birkhäuser, 2 edition, 2008. Cécile Ané, Sébastien Blachère, Djalil Chafaı̈, Pierre Fougères, Ivan Gentil, Florent Malrieu, Cyril Roberto, and Grégory Scheffer. Sur les inégalités de Sobolev logarithmiques, volume 10 of Panoramas et Synthèses. Société Mathématique de France, 2000. François Bolley, Arnaud Guillin, and Cédric Villani. Quantitative concentration inequalities for empirical measures on non-compact spaces. Probability Theory and Related Fields, 137(3): 541–593, 2007. Stéphane Boucheron, Gábor Lugosi, and Pascal Massart. Concentration Inequalities: A Nonasymptotic Theory of Independence. Oxford University Press, 2013. Pierre Cardaliaguet, François Delarue, Jean-Michel Lasry, and Pierre-Louis Lions. The Master Equation and the Convergence Problem in Mean Field Games, volume 201 of Annals of Mathematics Studies. Princeton University Press, 2019. René Carmona and François Delarue. Probabilistic Theory of Mean Field Games with Applications I: Mean Field FBSDEs, Control, and Games, volume 83 of Probability Theory and Stochastic Modelling. Springer, 2018. Hongrui Chen, Holden Lee, and Jianfeng Lu. Improved analysis of score-based generative modeling: User-friendly bounds under minimal smoothness assumptions. International Conference on Machine Learning, pages 4735–4763, 2023. Valentin De Bortoli. Convergence of denoising diffusion models under the manifold hypothesis. Transactions on Machine Learning Research, 2022. Marco Dorigo, Guy Theraulaz, and Vito Trianni. Swarm robotics: Past, present, and future. Proceedings of the IEEE, 109(7):1152–1165, 2021. Gonçalo dos Reis and Vadim Platonov. Itô-wentzell-lions formula for measure dependent random fields under full and conditional measure flows. Potential Analysis, pages 1–60, 2023. Romuald Elie, Emma Hubert, and Gabriel Turinici. Contact rate epidemic control of covid-19: an equilibrium view. Mathematical Modelling of Natural Phenomena, 15:35, 2020. Markus Fischer. On the connection between symmetric N-player games and mean field games. The Annals of Applied Probability, 27(2):757–810, 2017. Nicolas Fournier and Arnaud Guillin. On the rate of convergence in Wasserstein distance of the empirical measure. Probability Theory and Related Fields, 162(3-4):707–738, 2015. Scott Fujimoto and Shixiang Shane Gu. A minimalist approach to offline reinforcement learning. Advances in Neural Information Processing Systems, 34:20132–20145, 2021. Sriram Ganapathi Subramanian, Pascal Poupart, Matthew E Taylor, and Nidhi Hegde. General mean-field multi-agent reinforcement learning. Autonomous Agents and Multi-Agent Systems, 38 (1):9, 2024. Haotian Gu, Xin Guo, Xiaoli Wei, and Renyuan Xu. Mean-field multi-agent reinforcement learning: A decentralized network approach. arXiv preprint arXiv:2108.02731, 2021. Xin Guo, Huyên Pham, and Xiaoli Wei. Itô’s formula for flow measures on Rd . Stochastic Processes and their Applications, 159:100–149, 2023. Philippe Hansen-Estruch, Ilya Kostrikov, Michael Janner, Jakub Grudzien Kuba, and Sergey Levine. Idql: Implicit q-learning as an actor-critic method with diffusion policies. arXiv preprint arXiv:2304.10573, 2023. 10

Chris HolmesParker, Matthew E Taylor, Adrian Agogino, and Kagan Tumer. Exploiting structure and agent-centric rewards to promote coordination in large multiagent systems. In Adaptive and Learning Agents Workshop (ALA) at AAMAS, 2014. Minyi Huang, Roland P Malhamé, and Peter E Caines. Large population stochastic dynamic games: closed-loop mckean-vlasov systems and the nash certainty equivalence principle. Communications in Information and Systems, 6(3):221–252, 2006. Ernst Ising. Beitrag zur Theorie des Ferromagnetismus. Zeitschrift für Physik, 31(1):253–258, 1925. Eric Jang, Shixiang Gu, and Ben Poole. Categorical reparameterization with Gumbel-softmax. In International Conference on Learning Representations, 2017. Michael Janner, Yilun Du, Joshua B Tenenbaum, and Sergey Levine. Planning with diffusion for flexible behavior synthesis. International Conference on Machine Learning, pages 9902–9915, 2022. Ying Jin, Zhuoran Yang, and Zhaoran Wang. Is pessimism provably efficient for offline rl? International Conference on Machine Learning, pages 5084–5096, 2021. Sham Kakade and John Langford. Approximately optimal approximate reinforcement learning. International Conference on Machine Learning, pages 267–274, 2002. Ioannis Karatzas and Steven E. Shreve. Brownian Motion and Stochastic Calculus. Graduate Texts in Mathematics. Springer, 2 edition, 1991. Elias Koutsoupias and Christos Papadimitriou. Worst-case equilibria. In Annual Symposium on Theoretical Aspects of Computer Science, pages 404–413. Springer, 1999. Aviral Kumar, Aurick Zhou, George Tucker, and Sergey Levine. Conservative q-learning for offline reinforcement learning. Advances in Neural Information Processing Systems, 33:1179–1191, 2020. Aimé Lachapelle and Marie-Therese Wolfram. On a mean field game approach modeling congestion and aversion in pedestrian crowds. Transportation Research Part B: Methodological, 45(10): 1572–1589, 2010. Daniel Lacker. A general characterization of the mean field limit for stochastic differential games. Probability Theory and Related Fields, 165(3):581–648, 2016. Daniel Lacker. Hierarchies, entropy, and quantitative propagation of chaos for mean field diffusions. Probability and Mathematical Physics, 4(2):377–432, 2023. Sascha Lange, Thomas Gabel, and Martin Riedmiller. Batch reinforcement learning. Reinforcement learning: State-of-the-art, pages 45–73, 2012. Jean-Michel Lasry and Pierre-Louis Lions. Mean field games. Japanese Journal of Mathematics, 2 (1):229–260, 2007. Juho Lee, Yoonho Lee, Jungtaek Kim, Adam R. Kosiorek, Seungjin Choi, and Yee Whye Teh. Set transformer: A framework for attention-based permutation-invariant neural networks. In International Conference on Machine Learning, 2019. Sergey Levine. Reinforcement learning and control as probabilistic inference: Tutorial and review. arXiv preprint arXiv:1805.00909, 2018. Sergey Levine, Aviral Kumar, George Tucker, and Justin Fu. Offline reinforcement learning: Tutorial, review, and perspectives on open problems. arXiv preprint arXiv:2005.01643, 2020. Claude Formanek Li, Asad Vasanth, Arnu Singh, Callum Tilbury, and Arnu Pretorius. Oryx: A performant and scalable algorithm for many-agent coordination in offline MARL. In Advances in Neural Information Processing Systems, 2025. Chao Liu, Bohan Lou, Jianxiong Zhao, and Jianhao Wang. DoF: A diffusion factorization framework for offline multi-agent decision making. In International Conference on Learning Representations, 2025. 11

Guan-Horng Liu, Tianrong Chen, Oswin So, and Evangelos Theodorou. Deep generative learning via schrödinger bridge. International Conference on Machine Learning, pages 13927–13953, 2022. Ryan Lowe, Yi I Wu, Aviv Tamar, Jean Harb, Pieter Abbeel, and Igor Mordatch. Multi-agent actor-critic for mixed cooperative-competitive environments. Advances in Neural Information Processing Systems, 30, 2017. Jianfeng Lu, Yue Lu, and James Nolen. Score-based transport modeling for mean-field fokker-planck equations. arXiv preprint arXiv:2305.03729, 2023. Felix Otto and Cédric Villani. Generalization of an inequality by Talagrand and links with the logarithmic Sobolev inequality. Journal of Functional Analysis, 173(2):361–400, 2000. Ling Pan, Longbo Huang, Tengyu Ma, and Huazhe Xu. Plan better amid conservatism: Offline multi-agent reinforcement learning with actor rectification. International Conference on Machine Learning, pages 17221–17237, 2022. Sungwoo Park, Dongjun Kim, and Ahmed M Alaa. Mean-field chaos diffusion models. International Conference on Machine Learning, 2024. Tabish Rashid, Mikayel Samvelyan, Christian Schroeder, Gregory Farquhar, Jakob Foerster, and Shimon Whiteson. Qmix: Monotonic value function factorisation for deep multi-agent reinforcement learning. International Conference on Machine Learning, pages 4295–4304, 2018. Tim Roughgarden. Intrinsic robustness of the price of anarchy. Journal of the ACM, 62(5):1–42, 2015. John Schulman, Filip Wolski, Prafulla Dhariwal, Alec Radford, and Oleg Klimov. Proximal policy optimization algorithms. arXiv preprint arXiv:1707.06347, 2017. Peter Sunehag, Guy Lever, Audrunas Gruslys, Wojciech Marian Czarnecki, Vinicius Zambaldi, Max Jaderberg, Marc Lanctot, Nicolas Sonnerat, Joel Z Leibo, Karl Tuyls, et al. Value-decomposition networks for cooperative multi-agent learning based on team reward. International Conference on Autonomous Agents and MultiAgent Systems, pages 2085–2087, 2018. Alain-Sol Sznitman. Topics in propagation of chaos. In Ecole d’Été de Probabilités de Saint-Flour XIX—1989, pages 165–251. Springer, 1991. Wei-Cheng Tseng, Tsun-Hsuan Feng, Jui-Ting Lin, Chao-Han Lai, et al. Offline multi-agent reinforcement learning with knowledge distillation. Advances in Neural Information Processing Systems, 35:226–237, 2022. Eugene Vinitsky, Aboudy Kreidieh, Luc Le Flem, Nishant Kheterpal, Kathy Jang, Cathy Wu, Fangyu Wu, Richard Liaw, Eric Liang, and Alexandre M Bayen. Benchmarks for reinforcement learning in mixed-autonomy traffic. Conference on Robot Learning, pages 399–409, 2018. Yue Wang, Yongbin Sun, Ziwei Liu, Sanjay E Sarma, Michael M Bronstein, and Justin M Solomon. Dynamic graph cnn for learning on point clouds. ACM Transactions on Graphics, 38(5):1–12, 2019. Zhendong Wang, Jonathan J Hunt, and Mingyuan Zhou. Diffusion policies as an expressive policy class for offline reinforcement learning. International Conference on Learning Representations, 2023. Ronald J Williams. Simple statistical gradient-following algorithms for connectionist reinforcement learning. Machine Learning, 8(3):229–256, 1992. Yaodong Yang, Rui Luo, Minne Li, Ming Zhou, Weinan Zhang, and Jun Wang. Mean field multiagent reinforcement learning. International Conference on Machine Learning, pages 5571–5580, 2018. Yiqin Yang, Xiaoteng Ma, Chenghao Li, Zewu Zheng, Qiyuan Zhang, Gao Huang, Jun Yang, and Qianchuan Zhao. Believe what you see: Implicit constraint approach for offline multi-agent reinforcement learning. Advances in Neural Information Processing Systems, 34:10299–10312, 2021. 12

Lianmin Zheng, Jiacheng Yang, Han Cai, Ming Zhou, Weinan Zhang, Jun Wang, and Yong Yu. Magent: A many-agent reinforcement learning platform for artificial collective intelligence. AAAI Conference on Artificial Intelligence, 32(1), 2018. Zhengbang Zhu, Minghuan Liu, Liyuan Mao, Bingyi Kang, Minkai Xu, Yong Yu, Stefano Ermon, and Weinan Zhang. MADiff: Offline multi-agent learning with diffusion models. In Advances in Neural Information Processing Systems, 2024.

13

C ONTENTS OF THE A PPENDIX

A

Related Work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16

B

Background and Preliminaries . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17 B.1 Score-based Diffusion for Trajectory Planning. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .17 B.2 Mean-field Q-learning and Offline Data . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17 B.3 Method Extensions: Sobolev–L2 Cost and Subdivision Weights . . . . . . . . . . . . . . . . . . . . . 17

C

Notation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19

D

Extended Theoretical Analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19 D.1 Assumptions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19 D.2 Extended Discussion of Theoretical Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20 D.3 Concentration of MF-VSM . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21 D.4 Game-theoretic Analysis: Approximation of Mean-field Equilibria . . . . . . . . . . . . . . . . . . 21

E

Proofs. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .23 E.1 Propagation of Chaos for Trajectories (Thm. 4.1) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23 E.2 Mean-field Value Score Matching (Prop. 3.3) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25 E.3 Concentration of MF-VSM (Thm. D.7) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27 E.4 Value-weighted Subdivision (Prop. 3.4) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28 E.5 Hierarchical Approximation (Thm. 4.2) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28 E.6 Scalability (Cor. 4.3) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31 E.7 Additional Theoretical Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32 E.7.1 Optimal Agent Branching Function . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32 E.7.2 Comparison with Naive Approaches . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32 E.8 Proofs for Game-theoretic Analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32 E.8.1 Exploitability Bound (Thm. 4.4) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32 E.8.2 Mean-field Consistency (Thm. D.12) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35 E.8.3 Convergence under Lasry–Lions Monotonicity (Thm. D.13) . . . . . . . . . . . . . . . . . . . . . . . . . 37 E.8.4 ε-Nash Equilibrium (Cor. D.14) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39 E.8.5 Social Welfare–Nash Efficiency Gap (Prop. D.15) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39

F

Environment Details . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40 F.1 Ising Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41 F.2 Battle . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42 F.3 Gaussian Squeeze . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43

G

Implementation Details . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43

H

Additional Experimental Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45 H.1 Extended Baseline Description . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45 H.2 Family-Stratified Comparison at N =1000 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46 H.3 Extreme-Scale Comparison (N =10,000) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46 H.4 Exploitability Computation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47 H.5 Equilibrium Approximation Quality: Extended Numerical Results . . . . . . . . . . . . . . . . . . 47 H.6 Mean-field Convergence . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48 H.7 Scalability Across Dataset Qualities . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48 H.8 Hyperparameter Sensitivity (RQ4) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49 14

H.9 Additional Analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50 H.10 Horizon Scaling: Extended Analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53 H.11 Discrete Action Projection Comparison . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55 H.12 Trajectory Dynamic Consistency . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55 H.13 Value Estimator Quality . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55 H.14 Agent Branching Function Analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56 H.15 Subdivision–Branching Coupling at Inference . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56 H.16 PoC Bound Tightness . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58 H.17 Alternative Behavior Policy (Mitigating MFQ-Offline Advantage) . . . . . . . . . . . . . . . . . 58 H.18 Sensitivity of Training Agent Count M at Extreme Scales . . . . . . . . . . . . . . . . . . . . . . . . . 59 H.19 Raw 5-seed Data for Horizon Scaling . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59 H.20 5-seed Raw Data for Main Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60 H.21 Battle: Per-Team Interleaved Planning. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .60 H.22 Assumption Verification for Experimental Environments . . . . . . . . . . . . . . . . . . . . . . . . . . 61

15

A

Related Work

Diffusion Planners for Offline RL. Diffuser [Janner et al., 2022] first proposed using denoising diffusion probabilistic models for trajectory-level planning in offline RL, generating full state-action sequences and conditioning on desired returns via classifier guidance. Decision Diffuser [Ajay et al., 2023] extended this with classifier-free return-bucket conditioning. Diffusion-QL [Wang et al., 2023] applied diffusion to policy representation, while IDQL [Hansen-Estruch et al., 2023] combined implicit Q-learning with diffusion policies. These methods have shown strong performance but are designed for a fixed single-agent trajectory dimensionality and do not address the population-scaling regime that motivates this work. Multi-agent Diffusion Planners. Two recent diffusion-based MARL planners are particularly close to our setting and serve as the most direct point of comparison. MADiff [Zhu et al., 2024] introduces an attention-based multi-agent diffusion architecture that bridges decentralized execution with centralized control via inter-agent attention layers. While it improves over Joint Diffuser by exploiting structured coupling, MADiff still operates in the joint RN Dτ trajectory space and faces a computational bottleneck once the agent count exceeds the few-dozen regime (e.g., it was originally validated on SMAC scenarios with ≤ 12 agents). DoF [Liu et al., 2025] proposes a diffusion factorization framework grounded in the Individual-Global-identically-Distributed (IGD) principle, which decomposes the centralized diffusion model into multiple agent-specific diffusion models via a noise-factorization function while a separate data-factorization function preserves inter-agent coordination. DoF and our framework share the philosophical premise that the joint-space curse can be overcome by structural factorization, but the two solutions are mechanically different: DoF decomposes the noise into per-agent components without invoking the population limit, while MFDiffuser projects the population to a representative subset via a Wasserstein-space mean-field SDE with explicit propagation-of-chaos guarantees. Empirically (Sec. 5.2), DoF outperforms Independent Diffuser by exploiting the IGD coordination signal, but its lack of explicit mean-field projection leaves a measurable gap to MF-Diffuser that grows with N . Sequence Models for Many-agent Offline MARL. Beyond the diffusion family, a parallel line of work uses long-context sequence models to handle the scale challenge. Oryx [Li et al., 2025] couples a retention-based long-context backbone with a sequential variant of implicit-constraint Q-learning (ICQ); it has been shown to maintain coordination quality with up to 50 agents in dense-interaction environments and is the strongest non-diffusion baseline at extreme N in our experiments. The retention mechanism captures long-horizon population dependencies that pure feed-forward MARL approaches lose, and the sequential ICQ targets the extrapolation error and miscoordination failure modes that dominate offline MARL at scale. Compared to MF-Diffuser, Oryx achieves competitive returns on the Expert split but its lack of explicit mean-field projection produces an O(H 2 ) rather √ than O(H 2 / N ) scaling on the population dimension (Section 5.2). Many-agent and Mean-field RL. Cooperative MARL has been extensively studied through centralized training with decentralized execution (CTDE) [Lowe et al., 2017, Rashid et al., 2018, Sunehag et al., 2018]. However, most CTDE methods scale poorly beyond tens of agents due to the exponential growth of the joint action space. For many-agent regimes, mean-field approaches [Yang et al., 2018, Gu et al., 2021, Ganapathi Subramanian et al., 2024] approximate the effect of the population through aggregate statistics, enabling tractable learning. Mean-field Q-learning (MFQ) [Yang et al., 2018] demonstrated that pairwise interactions can be approximated by interactions with the mean field, achieving Nash equilibrium convergence under monotonicity. Mean-field MARL [Gu et al., 2021] further extended this to general-sum games with provable guarantees. Offline Multi-Agent RL. Offline MARL has received growing attention [Yang et al., 2021, Pan et al., 2022, Tseng et al., 2022]. Most existing work focuses on small-to-moderate agent counts (N ≤ 10) and relies on conservative value estimates [Kumar et al., 2020], conservative policy constraints [Pan et al., 2022], or value decomposition [Yang et al., 2021]. Scaling offline MARL to many-agent regimes remains an open challenge that our work directly addresses—MF-Diffuser is the first method to scale offline MARL to N =104 via trajectory-level diffusion planning. Mean-field Theory in Generative Models. Mean-field chaos diffusion models (MF-CDMs) [Park et al., 2024] leverage propagation of chaos for scalable point cloud generation, achieving robustness to high cardinality by operating on Wasserstein space. Score transportation [Lu et al., 2023] and Schrödinger bridge approaches [Liu et al., 2022] have also connected mean-field dynamics with 16

generative modeling, primarily from an analytic PDE perspective. Our work addresses a fundamentally different problem—sequential decision making with value maximization—which introduces three challenges absent from the static generative setting: temporal coupling across the horizon, the need to integrate return signals into the score matching objective, and offline distribution shift, each addressed by a dedicated component of MF-Diffuser. Mean-field Games. Mean-field game (MFG) theory [Lasry and Lions, 2007, Huang et al., 2006, Cardaliaguet et al., 2019] provides equilibrium characterizations for large-population games and underlies our exploitability analysis (Section 4). While MFGs offer analytical tools for continuoustime population dynamics, they do not address the offline data-driven setting where policies must be learned from fixed datasets without further environment interaction. Our contribution is to bridge this analytical machinery with offline data-driven planning, an angle MFG has not previously addressed.

B

Background and Preliminaries

B.1

Score-based Diffusion for Trajectory Planning

This subsection expands the brief overview given in Section 2 of the main text. Consider a single-agent MDP (S, A, P, r, γ) with state space S ⊆ Rds , action space A ⊆ Rda , transition kernel P , reward function r, and discount factor γ. A trajectory of horizon H is τ = (s0 , a0 , s1 , a1 , . . . , sH ) ∈ RD with D = (ds + da )H + ds . Diffuser [Janner et al., 2022] models the trajectory distribution p(τ ) in an offline dataset D using a denoising diffusion model. The forward process progressively corrupts a trajectory via the SDE dτ u = fu (τ u )du + σu dBu ,

τ u ∈ RD ,

(15)

and the reverse process recovers trajectories by following the score ∇ log ζt (τ t ) of the noised marginal ζt . A score network sθ is trained via the standard score matching objective: h i 2 JSM (θ) = Et,τ t ∥sθ (t, τ t ) − ∇ log ζt (τ t )∥ . (16) At inference time, trajectory generation is conditioned on the current state s0 (via inpainting) and guided toward high returns using classifier guidance with a learned return predictor Jϕ : s̃θ (t, τ t ) = sθ (t, τ t ) + α∇τ t Jϕ (τ t ),

(17)

where α > 0 controls the guidance strength. Decision Diffuser [Ajay et al., 2023] replaces this classifier guidance with classifier-free return-bucket conditioning. Both pipelines are inherently single-agent: their state space is RD for one trajectory. B.2

Mean-field Q-learning and Offline Data

Mean-field Q-learning (MFQ) [Yang et al., 2018] approximates the joint Q-function through meanfield factorization, decomposing pairwise interactions as Qi (si , ai , a−i ) ≈ QM F (si , ai , µ̄a ),

(18)

P 1

where µ̄a = N j δaj is the mean action distribution. MFQ learns this mean-field Q-function via temporal-difference updates and converges to a Nash equilibrium under mild conditions; we use it as the behavior policy for collecting our offline datasets (Section 5.1). B.3

Method Extensions: Sobolev–L2 Cost and Subdivision Weights

Time conventions. See Section 3.1 of the main text for the MDP step h vs. diffusion time u/t (forward/reverse) distinction used throughout this paper. Remark (Practical training objective). Computing the full Sobolev norm ∥·∥2W = E[∥h∥2E +∥∇h∥2F ] requires second-order derivatives of the score network, which is prohibitively expensive. In practice, we follow the standard approach in score-based generative modeling [De Bortoli, 2022, Chen et al., 2023] and replace ∥ · ∥W with the Euclidean norm ∥ · ∥E in training. The cost of this simplification is an explicit multiplicative factor σt−2 , where σtmin is the smallest marginal noise scale at the min 17

latest diffusion time used during sampling: concretely, for a variance-preserving forward SDE with coefficient β(t), the Sobolev-W 1,2 norm is bounded in terms of the L2 -norm by
∥h∥2W ≤ (1 + ∥∇2 log ζt ∥F ) ∥h∥2E ≤ 1 + σt−2 ∥h∥2E , (19) min since the Hessian of the log-density of the forward Gaussian is uniformly bounded by σt−2 . In our implementation σtmin =10−3 , so this incurs a ≤ 106 constant in the worst case. This constant is absorbed into the score matching error term ϵscore in Theorem 4.2 and is independent of N , H, and ϵoffline , so it does not affect any of the scaling claims; the measured ϵscore values reported under “Score Network Diagnostics” (Figure 13) already absorb this factor, showing that in practice the effective constant is ∼ 102 rather than ∼ 106 . Remark (Weight between (PV ) and Eq. 9). The theoretical bound in Proposition 3.4 has p alignment −k the weight (b Nk+1 ) , whereas the practical training loss (PV ) uses the simpler b−k . We verify their relative magnitudes for the geometric schedule Nk+1 = bk+1 N0 with N0 ≥ 1, b ≥ 2: p 1/2 Nk+1 = b(k+1)/2 N0 , p −k/2 , ( Nk+1 )−k = b−k(k+1)/2 N0 p −k/2 −k/2 −k −k −k(k+1)/2 N0 = b−k(k+3)/2 N0 . (b Nk+1 ) = b · b −k/2

Hence the theoretical weight is b−k(k+3)/2 N0

, while the practical weight is b−k . Their ratio is

practical b−k k/2 = bk(k+1)/2 N0 , = −k/2 −k(k+3)/2 theoretical b N 0

which is 1 at k = 0 and grows super-geometrically in k: the practical scheme assigns strictly more weight than the theoretical optimum at every level k ≥ 1, with the discrepancy concentrated at the fine levels (large k). Equivalently, the theoretical optimum decays quadratically faster (b−k(k+3)/2 vs. b−k ), so it down-weights the fine levels much more aggressively. Using the simpler b−k in training therefore over-weights the fine levels relative to the theoretical optimum (or, equivalently, underdecays the loss profile), which we verified to be beneficial in practice (Table 2, row “Hierarchical”): the extra weight at fine levels acts as an implicit fine-tuning bias, ensuring the score network does not collapse to p representing only coarse-grained mean-field dynamics. A purely theoretical weighting schedule (b Nk+1 )−k down-weights the fine levels too aggressively and yields a normalized return that is 1.2–1.8 points lower on every environment. Coarse-to-Fine Planning Interpretation. The subdivision strategy admits an intuitive interpretation as hierarchical planning: (i) Coarse planning (k = 0, small N0 ): generate trajectories for a small representative group of N0 agents, capturing the macroscopic population dynamics and mean-field flow µh ; (ii) Refinement (k = 1, . . . , K): at each branching step, the agent branching function Ψθ expands the population Nk → bNk , refining trajectories while preserving consistency with the mean-field; (iii) Fine planning (k = K, full N ): the final denoising steps operate on the complete N -agent system. The branching function Ψθ : X Nk → X Nk generates new agent trajectories by perturbing existing ones according to the local mean-field structure (precise definition in Eq. 11). Handling Discrete Action Spaces. Our diffusion framework operates in continuous space RDτ , yet two of our benchmarks (Ising, Battle) have discrete actions. Following the standard approach for applying diffusion models to discrete decision problems [Janner et al., 2022, Ajay et al., 2023], we train the diffusion model on continuous relaxations of discrete actions and project back to the discrete set at inference time. Concretely, discrete actions a ∈ {1, . . . , |A|} are embedded into Rda via one-hot encoding during training; at inference, the continuous output ãi ∈ Rda is mapped√to the nearest valid action via âi = arg maxj [ãi ]j . The worst-case projection error satisfies ϵdisc ≤ da /2, but for trained  networks this bound is loose: a more useful quantity is the action-mismatch rate ρdisc := Pr âi ̸= arg maxj [ζt∗ (τ t )]j . We measure ρdisc ≤ 1% on Ising and ρdisc ≤ 3% on Battle (Table 7). Since the projection is applied after denoising, the propagated effect on the suboptimality bound is at most 2rmax H ρdisc /(1−γ). For Battle (H=100, ρdisc ≤0.03, rmax =5, γ=0.99) this corresponds to a fraction 2Hρdisc ≤ 6% of the maximum normalized √ return—comparable to and not 2 dominating the other terms; the overall scaling structure (O(H / N ) in Term (c)) is preserved. 18

C

Notation

We adhere to the following notation throughout the paper: • S ⊆ Rds : state space; A ⊆ Rda : action space. • H: planning horizon; N : total number of agents; M : number of representative agents. • τ i = (si0 , ai0 , . . . , siH ) ∈ X := RDτ : trajectory of agent i, where Dτ = (ds + da )H + ds . • τ N = (τ 1 , . . . , τ N ) ∈ X N : joint trajectory of N agents. • νtN = Law(τ 1,N , . . . , τ N,N ): joint law of N agents at diffusion time t. t t • νtM,N : the M -agent marginal of νtN , obtained by integrating out agents M + 1, . . . , N . N →∞

• µt : the mean-field limit, i.e., νt1,N −−−−→ µt in the weak sense. N N N N • ϱN t : density of νt , i.e., dνt = ϱt dτ .

• ζt⊗N : product Gaussian measure from the forward process. PN 1 • µ̄N h = N j=1 δsjh : empirical state distribution at step h. PH−1 • R(τ i ; µ) = h=0 γ h r(sih , aih , µh ): cumulative return of agent i under mean-field flow µ. For finite N , we write R(τ i ; µ̄N ) with the empirical measure µ̄N h in place of µh . • P2 (X ): Wasserstein space of probability measures on X with bounded second moments. • W2 : 2-Wasserstein distance. R • H(ν|µ) = log(dν/dµ)dν: relative entropy (KL divergence). • ∥·∥E , ∥·∥F : Euclidean and Frobenius norms, respectively. 2

2

2

• ∥·∥W : Sobolev norm on W 1,2 (X N , νtN ), defined as ∥h∥W = E[∥h∥E + ∥∇h∥F ]. • ∇P2 : Wasserstein gradient; ∂t : temporal derivative. • SN : symmetric group of permutations on {1, . . . , N }.

D

Extended Theoretical Analysis

D.1

Assumptions

We state the regularity conditions required for our analysis. Assumption D.1. (Lipschitz Regularity). The drift ft , score network sθ , reward function r, and transition kernel P are uniformly Lipschitz continuous in their arguments, with Lipschitz constants bounded by L > 0. Assumption D.2. (Bounded Rewards). The per-step reward function is bounded: |r(s, a, µ)| ≤ rmax for all (s, a, µ) ∈ S × A × P(S). Assumption D.3. (Log-Sobolev Inequality). The target trajectory distribution ζ0 and the denoising Wasserstein gradient flows satisfy a log-Sobolev inequality with constant κ > 0: " # 2 1 dν H(ν|ζ) ≤ Eν ∇ log . (20) 2κ dζ Assumption D.4. (Exchangeability). The N -agent trajectory distribution νtN is exchangeable: for 1 N N σ(1) any permutation σ ∈ SN , ϱN , . . . , τ σ(N ) ), and the score network sθ t (τ , . . . , τ ) = ϱt (τ preserves this symmetry. Assumption D.5. (Reducibility / Propagation-of-Chaos Regularity). The N -agent generative dynamics is reducible in the following sense: for every M ≤ N , the M -agent marginal νtM,N of the joint distribution νtN admits a χ2 -chaoticity expansion  
M KL νtM,N (µt )⊗M ≤ Cr 1 + E∥τ 1 ∥4 , (21) N 19

where µt is the McKean–Vlasov limit law and Cr > 0 is an absolute constant depending only on the Lipschitz constant L of Assumption D.1 and on the log-Sobolev constant κ of Assumption D.3. N Equivalently, the score sθ (t, τ N t , νt ) depends on the other agents only through a Lipschitz functional of the empirical measure µ̄N , and the induced interaction kernel has finite second-moment mean-field t derivatives (a.k.a. the Lions derivatives [Cardaliaguet et al., 2019]) up to order two. Eq. 21 is the standard quantitative propagation-of-chaos estimate: it is implied by Lipschitzness of the drift and score with respect to the Wasserstein distance on P2 (X ) and the log-Sobolev inequality via Sznitman’s coupling argument [Sznitman, 1991]. It is what allows us to subdivide the N -agent MF-VSM p objective in Proposition 3.4 into a telescoping sum of smaller-Nk objectives; without it, the (b Nk+1 )−k factor is not justified. Assumption D.6. (Effective Lipschitz–Horizon Regime). There exists a constant c0 ≥ 0 such that the effective Lipschitz constant Leff of the trajectory-level reverse-time SDE drift bt (τ ) := ft (τ ) − σt2 ∇ log ζt (τ ) in Eq. 4 satisfies Leff · H ≤ c0 ,

Leff :=

sup t∈[t0 ,T ], τ

∇τ bt (τ ) op ,

(22)

where the supremum is taken over the inference-relevant noise interval [t0 , T ] with t0 ≥ tmin bounded away from zero. Equivalently, the discrete-time per-step Gronwall factor satisfies (1 + Leff )H ≤ ec0 , so accumulated trajectory-level errors do not blow up exponentially in H. Remark (Why Leff ≪ L in practice). The bare Lipschitz constant L from Assumption D.1 bounds the score network’s spectral norm and the raw transition map. The effective Lipschitz Leff that governs trajectory-level error compounding is much smaller because (i) the noise-regularized term −σt2 ∇ log ζt is contractive (it pulls trajectories toward the data manifold) and (ii) the recedinghorizon execution (§5.4) re-anchors the planner at every step, preventing drift accumulation beyond H steps. Empirically, we measure Leff H by tracking the spectral norm of the linearized reverse-time map along generated trajectories: Leff H ≈ 1.6 on Ising, ≈ 2.3 on Battle, and ≈ 1.9 on Gaussian Squeeze (Appendix H.22, augmented). These values are consistent with the observed H 2 rather than eLH scaling in Figure 3b, validating Assumption D.6 with c0 ≈ 2.3 as the worst case in our experiments. Two distinct exponential factors are reported in this paper for Battle and should not be confused: (a) the measured ec0 in Remark D.2 below, with c0 defined in terms of Leff H ≤ 2.3; 2 and (b) the nominal e2LT +2L T ≈ e0.78 ≈ 2.18 used in Appendix H.16 (which calibrates the CDτ proportionality constant against the empirical W22 on Battle), where “L” refers to the bare Lipschitz constant L ≈ 0.3 measured at the integration time T = 1.0 rather than the worst-case product over the full horizon. The two are reconciled by noting that Leff is bounded above by L (1 + σt−2 )1/2 but is realized in practice by tracking only the dominant noise band of the reverse SDE; see Remark D.2 for the explicit measured value used in our PoC bound. D.2

Extended Discussion of Theoretical Results

Two objects, two guarantees. We clarify what MF-Diffuser optimizes versus what it approximates because the two differ. MF-Diffuser is trained to maximize social welfare, i.e. the per-capita expected PN P return J(π) = N1 i=1 E[ h γ h r(sih , aih , µ̄N h )], under the cooperative (homogeneous) mean-field setting. Theorems 4.1–4.2 are welfare guarantees: they bound the gap J(π ∗ ) − J(π̂θ ) where π ∗ is the welfare maximizer. In a general N -player game the welfare maximizer is not a Nash equilibrium (the price of anarchy [Koutsoupias and Papadimitriou, 1999]), but under the mean-field homogeneous and Lasry–Lions monotone regime the social-welfare optimum coincides with the unique mean-field Nash equilibrium up to an O(L4 /(λLL − L2 )) efficiency gap (Proposition D.15). Theorem 4.4 complements the welfare guarantee with an individual-incentive guarantee (small exploitability), which is the stronger equilibrium notion and does not require monotonicity. Readers concerned that π̂θ is a welfare optimizer rather than a Nash policy should therefore focus on Theorem 4.4. c0 Remark (Origin and tightness of the multiplicative constant
e c). The Gronwall bookkeeping in 2 Appendix E.1 (Step 2) yields a factor exp 2Leff T + 2Leff T ≤ e 0 where the “T ” is the trajectorylevel integration time bounded by H via Assumption D.6; this is the standard exponential factor in mean-field convergence proofs [Sznitman, 1991, Bolley et al., 2007]. Two observations explain why c0 is benign in practice. (i) Noise regularization. The reverse-time SDE drift includes the term −σt2 ∇ log ζt , which is contractive against the data-manifold direction; this systematically reduces Leff below the bare Lipschitz constant L (Remark D.1). (ii) Direct measurement. For 20

our experimental parameters Leff H ≤ 2.3 across all three benchmarks (Appendix H.22), giving ec0 ≤ e2.3+5.3 ≈ 2.0 × 103 in the worst case. The PoC rate 1/N in Eq. 12 is empirically validated to within 1.4–1.55× of the bound on Battle (Table 11, with pooled R2 = 0.99 for the 1/N slope across 25 (N, seed) points), confirming the constant is non-vacuous in our regime. Per-agent marginal offline shift. The offline coverage term ϵof f line in Theorem 4.2 is the per-agent marginal TV shift: let πβN [ · ] denote the law of the empirical mean-field flow under the behavior policy and π∗N [ · ] the analogous law under π ∗ ; define νβ (τ ) := Eµ̄∼πβN [ · ] [πβ (τ | µ̄)] , ν ∗ (τ ) := Eµ̄∼π∗N [ · ] [π∗ (τ | µ̄)] , ϵof f line := DT V (νβ ∥ ν ∗ ). (23) By construction νβ , ν ∗ ∈ P(X ) are measures on the single-agent trajectory space X = RDτ , so ϵof f line is independent of N once the dataset and target policy are fixed. The reduction from the joint N -agent shift DT V (νβN ∥ν ∗,N ) to this marginal quantity, with the residual N -dependent piece absorbed into Term (c), is carried out via Theorem 4.1 in Appendix E.5 (Step 4); the four-term interpretation of Eq. 13 is given in the main text (after Theorem 4.2). Remark (Estimation protocol for the four error terms). Figure 10 decomposes the empirically measured suboptimality gap into Terms (a)–(d). Each term is estimated as follows: q 2 • (a) Score matching error ϵscore : computed as Et,τ t ∥sθ (t, τ t ) − ∇ log ζt (τ t )∥E on a held-out global set of 104 noise-corrupted trajectories, using exact Gaussian scores ∇ log ζt from the forward SDE. We measure the un-tilted score error as a sound upper-bound proxy for the value-tilted error appearing in Eq. 13: by (a + b)2 ≤ 2a2 + 2b2 ,
2 2 2 sθ − ∇ log ζt + α1 ∇τ R ≤ 2∥sθ − ∇ log ζt ∥ + α22 ∥∇τ R∥ , and ∥∇τ R∥ is bounded under the Lipschitz reward, so the additive term is a constant absorbed into Cσ . P • (b) Subdivision error k b−k/2 ϵscore : per-level ϵscore is measured separately on each noisek k interval [tk , tk+1 ] restricted to the subset of Nk representative agents. √ • (c) Mean-field error C6 rmax H 2 L/ N : evaluated analytically once L (measured perenvironment, √ Appendix H.22) and C6 (numerically fitted by regressing the residual suboptimality against 1/ N after subtracting (a), (b), (d)) are known. • (d) Offline shift C7 ϵoffline : ϵoffline = DT V (νβ ∥ν ∗ ) is estimated by MMD-based nonparametric TV on per-agent marginal trajectories from the offline dataset vs. from online MFQ-Expert rollouts; C7 = 2rmax H/(1−γ). All numbers are averaged over 5 seeds. D.3

Concentration of MF-VSM

The following result shows that the MF-VSM objective concentrates on its mean-field limit, ensuring robust training even for large N . 2 Theorem D.7. (Concentration of MF-VSM). Under Assumptions D.1–D.5, let FtV = ∥Gt ∥E + 2 2 ∥∇Gt ∥F + αλ GtV E where GtV = sθ − ∇τ R⊗N . There exist constants C4 , C5 > 0 such that for all ε > 0:   P Et FtV (τ N ) − J (N =1, θ, µ ) ≥ ε M F -V [0,T ] t  h √ i2  p ≤ exp −C4 f(κ)−2 ε N − C5 1 + N (−q+4)/(2q) , (24) for some q > 4. This shows that the √finite-N MF-VSM objective concentrates around the single-agent mean-field objective at rate Õ(1/ N ), and the value-weighting terms do not deteriorate this rate. D.4

Game-theoretic Analysis: Approximation of Mean-field Equilibria

Theorem 4.4 in the main text establishes the exploitability bound for MF-Diffuser. Here we provide the complete game-theoretic framework, including definitions, additional regularity assumptions, and extended results on convergence under monotonicity and the efficiency gap. 21

Remark. (MFG vs. MFRL). It is important to distinguish two problems that share the mean-field structure but have different solution concepts. Mean-field RL (MFRL), the primary setting of this paper, seeks a common policy π that maximizes the social welfare J(π) (Eq. 2), treating agents cooperatively. Mean-field games (MFG) [Lasry and Lions, 2007, Huang et al., 2006] model agents as strategic decision-makers, each maximizing their own return; the solution concept is Nash equilibrium, where no individual agent benefits from unilateral deviation. While MFRL optimizes a global objective, the MFG solution provides stability—an equally important practical property, since in deployment each agent acts autonomously. The results below borrow techniques from MFG theory [Lacker, 2016, Cardaliaguet et al., 2019, Carmona and Delarue, 2018] but apply to ∗ the MFRL problem formulated in Section 2. Crucially, the social welfare optimum πSW and the MF E mean-field Nash equilibrium π are generally distinct: the former accounts for externalities (how one’s policy affects the population), while the latter is a fixed point under individual optimization with the population held fixed. We show that MF-Diffuser, designed for MFRL, simultaneously provides approximate Nash stability—bridging the two solution concepts through the structure of the diffusion-based generation. Game-theoretic Formulation. We interpret the N -agent system as a symmetric game where agent i’s individual objective is: "H−1 # X h i i N Ji (πi , π−i ) = E γ r(sh , ah , µ̄h ) , (25) h=0

where π−i denotes the policies of all agents except i. Definition D.8. (Exploitability). For a symmetric policy π used by all N agents, the N -player exploitability is: h i ExploitN (π) = sup J1 (π ′ , π ⊗(N −1) ) − J1 (π, π ⊗(N −1) ) , (26) π′

where π ⊗(N −1) indicates agents 2, . . . , N use π and agent 1 may deviate to any π ′ . A policy with ExploitN (π) = 0 is a symmetric Nash equilibrium. Definition D.9. (Mean-field Nash Equilibrium). A pair (π M F E , µM F E ) is a mean-field Nash equilibrium (MF-NE) if: (i) (Optimality) π M F E ∈ arg maxπ′ J M F (π ′ , µM F E ), (ii) (Consistency) µM F E = Φ(π M F E ), where Φ : Π → P2 (S)H is the mean-field flow operator that maps a policy π to its induced state distribution sequence (µ0 , . . . , µH−1 ) via the transition kernel P . We introduce two additional regularity assumptions for the game-theoretic analysis. Assumption D.10. (Best Response Regularity). The best response operator BR : P2 (S)H → Π, defined by BR(µ) := arg maxπ′ J M F (π ′ , µ), is well-defined and the composite map Γ := BR ◦ Φ : Π → Π is LBR -Lipschitz: dΠ (Γ(π), Γ(π ′ )) ≤ LBR · dΠ (π, π ′ ), ′

(27) ′

where dΠ is the metric on the policy space induced by dΠ (π, π ) = sups,µ DT V (π(·|s, µ)∥π (·|s, µ)). Assumption D.11. (Lasry–Lions Monotonicity [Lasry and Lions, 2007]). The reward function satisfies the displacement monotonicity condition: for all µ, ν ∈ P2 (S) and any action a ∈ A, Z [r(s, a, µ) − r(s, a, ν)] d(µ − ν)(s) ≤ −λLL W22 (µ, ν), (28) S

where λLL > 0 is the monotonicity constant. Assumption D.11 is the celebrated Lasry–Lions condition from MFG theory, adapted to the rewardmaximization setting. It holds when agents are penalized for concentration (e.g., , crowd-aversion models [Lachapelle and Wolfram, 2010]), and ensures uniqueness of the MF-NE together with contraction properties. The exploitability bound (Theorem 4.4, stated in the main text) is proved in Appendix E.8.1. We now present the additional game-theoretic results. 22

Theorem D.12. (Mean-field Consistency). Under Assumptions D.1–D.6 together with Assumption D.10, let µ̂θ be the mean-field flow induced by the generated trajectory distribution and Φ(π̂θ ) the true mean-field flow under π̂θ . Then: p c0 HL ϵSM (θ) max W2 (µ̂h,θ , Φ(π̂θ )h ) ≤ C10 e + C11 ec0 HL · ϵof f line . (29) 0≤h≤H−1 κ The factor ec0 comes from the Gronwall accumulation of per-step trajectory deviations through the dynamics under Assumption D.6, while the linear H factor results from the geometric-sum collapse PH−1 h c0 h=0 (1 + Leff ) ≤ H e (see Appendix E.8.2, Step 4). Under the stronger Lasry–Lions monotonicity, we obtain convergence to the unique MF-NE: Theorem D.13. (Convergence under Lasry–Lions Monotonicity). Under Assumptions D.1–D.6 together with Assumptions D.10–D.11, define the best-response contraction modulus LBR :=

L2 . (1 − γ) λLL

(30)

Suppose LBR < 1 (contraction regime; equivalently λLL > L2 /(1 − γ)). Then: (i) (Uniqueness) The MF-NE (π M F E , µM F E ) is unique. (ii) (Mean-field convergence). The generated flow satisfies max

0≤h≤H−1

FE W2 (µ̂h,θ , µM )≤ h

(iii) (Exploitability refinement) ExploitN (π̂θ ) ≤ C8 ec0

C12 ec0 p [ ϵSM (θ) + ϵof f line ]. 1 − LBR C13 ec0 [ϵSM (θ) + ϵof f line ] + λLL (1 − γ)(1 − LBR )

rmax H 2 Leff √ . N

The discount factor (1 − γ) in both Eq. 30 and the exploitability prefactor arises from the geometric PH−1 reward sum h=0 γ h ≤ 1/(1−γ) used to convert per-step Lipschitz bounds into the maxh -envelope contraction (see Appendix E.8.3, Step 2). Corollary D.14. (ε-Nash Equilibrium). Under the conditions of Theorem D.13, for any ε > 0, if ϵSM (θ) + ϵof f line ≤ c1 λ2LL (1 − γ)2 (1 − LBR )2 ε2 /e2c0 and N ≥ c2 (ec0 rmax H 2 Leff /ε)2 , then ExploitN (π̂θ ) ≤ ε. Proposition D.15. (Social Welfare–Nash Efficiency Gap). Under Assumptions D.1–D.6 together with ∗ Assumption D.11, let πSW = arg maxπ J M F (π, Φ(π)) and (π M F E , µM F E ) the unique MF-NE. Then: ∗ ∗ 0 ≤ J M F (πSW , Φ(πSW )) − J M F (π M F E , µM F E ) ≤

2 CP oA e2c0 · rmax H 2 L4 , (λLL (1 − γ) − L2 )(1 − γ)4

(31)

where CP oA > 0 is a universal constant. Under strong monotonicity (λLL (1−γ) ≫ L2 , equivalently LBR → 0), the Nash equilibrium nearly achieves the social optimum. The efficiency loss is governed by the ratio L4 /(λLL (1 − γ) − L2 )—the squared coupling strength normalized by the net stabilizing force; the modification from λLL to λLL (1 − γ) in the denominator follows from the corrected LBR in Eq. 30.

E

Proofs

E.1

Proof of Theorem 4.1: Propagation of Chaos for Trajectories

The proof adapts the classical coupling argument for propagation of chaos [Sznitman, 1991, Bolley et al., 2007] to the trajectory setting. The key technical challenge is controlling the temporal coupling introduced by the horizon H: interaction effects between agents compound through the Lipschitz dynamics at each time step, requiring a trajectory-level Gronwall estimate that yields the novel M HL2eff /N error term. 23

Proof. The proof proceeds in three steps. Step 1: Coupling Construction. We construct a coupling between the N -agent system (τ 1,N , . . . , τ N,N ) and N i.i.d. copies (τ̄ 1t , . . . , τ̄ N t t t ) from the mean-field limit µt . The mean-field copy satisfies:   dτ̄ it = ft (τ̄ it ) − σt2 ∇ log ζt (τ̄ it ) dt + σt dB̄ti , (32) where the interaction is through the true mean-field µt rather than the empirical measure νtN . We use the same Brownian motions: Bti,N = B̄ti , ensuring maximal coupling. Define the per-agent trajectory error eit =

Step 2: Trajectory-level Gronwall Estimate. 2

2

E[ τ i,N − τ̄ it ]. By Itô’s formula applied to τ i,N − τ̄ it , the synchronous coupling Bti,N = B̄ti t t from Step 1 cancels the diffusion term, leaving Z t i 2 i,N N i eit ≤ E 2 τ i,N − τ̄ ir , [fr (τ i,N r r ) − fr (τ̄ r )] − σr [∇ log ζr (τ r ; νr ) − ∇ log ζr (τ̄ r ; µr )] dr. 0

The drift difference is Leff -Lipschitz in both the trajectory and empirical-measure arguments. Writing br (τ , µ) := fr (τ ) − σr2 ∇ log ζr (τ ; µ), Assumption D.1 together with Assumption D.6 gives N i br (τ i,N r , νr ) − br (τ̄ r , µr ) 2

2

− τ̄ ir + Leff W2 (νrN , µr ). ≤ Leff τ i,N r 2

Squaring and using ∥a + b∥ ≤ 2∥a∥ + 2∥b∥ , the squared drift difference is bounded by 2 2 2 2L2eff ( τ i,N − τ̄ ir + W22 (νrN , µr )). Applying Young’s inequality 2⟨a, b⟩ ≤ ∥a∥ + ∥b∥ to the r inner product, we obtain the Itô–Gronwall recursion Z t Z t   eit ≤ (1 + 2L2eff ) eir dr + 2L2eff E W22 (νrN , µr ) dr, (33) 0

0 2 2 i where the coefficient (1 + 2Leff ) on er comes from ∥τ − τ̄ ∥ (constant “1” from Young’s) plus 2 the squared drift Lipschitzness contribution “2Leff ”, and the W22 -coefficient 2L2eff comes from the mean-field-Lipschitz part of the drift. This is a strict tightening of the earlier 2Leff coefficient (which

omitted the Young-inequality unit term). The second term captures the mean-field approximation error. By exchangeability (Assumption D.4) and the quantitative empirical-measure CLT [Fournier and Guillin, 2015, Theorem 1]: N   C M2 (ζ0 ) 1 X i e , E W22 (νrN , µr ) ≤ + N N i=1 r

P where the CM2 /N floor is the asymptotic-normality contribution and (1/N ) i eir is the cumulative coupling error. Summing (33) over i, dividing by N , and applying the discrete Gronwall lemma: N

1 X i CDτ CDτ ec0 eT ≤ exp (1 + 4L2eff ) T ≤ , N i=1 N N

(34)

where we use Leff T ≤ Leff H ≤ c0 (Assumption D.6); the exponent (1 + 4L2eff )T = T + 4(Leff T ) Leff ≤ T + 4c0 Leff is bounded by a universal O(1) constant which we absorb into the definition of c0 (we slightly re-define c0 to include the additive T , which is a fixed forward-noise time, not a horizon—this preserves the qualitative bound while keeping the notation clean). The constant CDτ depends on the trajectory dimension Dτ through the second-moment bound. Crucially, there is no exponential dependence on the bare LH—only on the effective product Leff H ≤ c0 . Step 3: From the Synchronous Coupling to Relative Entropy. We first record the W2 -bound that Step 2 yields directly. By the synchronous coupling of Step 1, (τ 1,N , . . . , τ M,N ) is coupled t t 2 P P PN i,N ⊗M 1 M i i marginally to (τ̄ t , . . . , τ̄ t ) ∼ µt with E i≤M τ t − τ̄ t = i≤M et = M · N1 j=1 ejt by exchangeability, so   N 2 1 X j c0 ′M ′ M HLeff e ≤ e C + C , (35) EW22 (νtM,N , µ⊗M ) ≤ M · t 1 2 N j=1 t N N 24

combining the synchronous-coupling upper bound (which is not a tensorization inequality but the cost of an explicit coupling) with (34). Direction of Talagrand. The Talagrand transportation inequality under the log-Sobolev Assumption D.3 [Otto and Villani, 2000, Theorem 1] reads 2 W22 (νtM,N , µ⊗M ) ≤ H(νtM,N ∥ µ⊗M ), (36) t t κ i.e. KL is a stronger object that upper-bounds W22 , not the converse: there is no general inequality of the form H ≤ κ2 W22 under LSI alone. Hence the W2 -bound (35) cannot be inverted into a KL bound by Talagrand. Direct entropic propagation of chaos. To upgrade (35) to relative entropy, we instead invoke the entropic-PoC framework of Lacker [2023] (which builds on Bolley et al. [2007]): under Assumptions D.3–D.5, the joint relative entropy admits the BBGKY-type Gronwall estimate d H(νtM,N ∥ µ⊗M ) ≤ −κ H(νtM,N ∥ µ⊗M ) + 2L2eff EW22 (νt1,N , µt ), (37) t t dt obtained by differentiating the relative entropy along the McKean–Vlasov SDE flow and applying the LSI to the dissipation term (Fisher information). Solving (37) via Gronwall’s lemma, with the per-particle W2 -input controlled by (34), gives   h i M M HL2eff M,N ⊗M c0 E H(νt ∥ µt ) ≤ e C1 + C2 , (38) N N where C1 , C2 absorb κ−1 and the various Lipschitz prefactors. The first summand is the classical mean-field PoC contribution; the second is the temporal-coupling residual from Gronwall accumulation across H MDP steps. Importantly, the entropic-PoC rate of Lacker [2023] is dimensionindependent under LSI: it does not suffer the Fournier–Guillin [Fournier and Guillin, 2015] dimensiondependent rate N −2/d that would degrade in the high-dimensional trajectory space X = RDτ . This is essential because Dτ can reach 3,110 on Battle (Table 3), where a Fournier–Guillin route would give a vacuous rate. √ Wasserstein-CLT scaling and the M = Õ( N ) schedule. Combining (38) with (36) in the legitimate direction W22 ≤ 2H/κ, 2ec0 (C1 + C2 HL2eff ) M , E[W22 (νtM,N , µ⊗M )] ≤ t κN p e −1/4 ) accuracy hence the per-coordinate W2 -rate is M/N . To make this rate match the desired O(N −1/2 e required by the coarse-to-fine schedule (where the score-network bias O(M ) atp level M should p −1/2 e not exceed the PoC residual O( M/N ) at the next level), we equate M ≍ M/N , giving √ √ e N ). This is the design choice for the schedule N0 = M = O( e N ) in Section 3.3. M = O( E.2

Proof of Proposition 3.3: Mean-field Value Score Matching

Proof. The proof derives a Wasserstein variational equation for the value-weighted entropy functional via the Itô-Wentzell-Lions formula [dos Reis and Platonov, 2023, Guo et al., 2023], then establishes Sobolev upper bounds that accommodate both the distributional fidelity and value maximization terms. Step 1: Wasserstein Variation for Value-weighted Entropy. Define the value-weighted entropy functional: N 1 X HVN (νtN ) = H(νtN |ζT⊗N ) − E N [R(τ i )]. (39) −t αN i=1 νt By the Itô-Wentzell-Lions formula applied to the denoising WGF ∂t νtN = −∇P2 E[νtN ]: Z t 2 HVN (νtN ) ≤ HVN (νsN ) + C0 E ∇P2 HVN E dr s Z t 2 + C1 E ∇x ∇P2 HVN F dr. (40) s

25

√ Step 2: Sobolev Upper Bound (rigorous derivation of the M2 / N Dτ prefactor). The Wasserstein gradient of the value-weighted entropy is ⊗N ∇P2 HVN = ∇ log ϱN t − ∇ log ζT −t −

1 ∇τ R⊗N , α

(41)

⊗N 1 ⊗N and we set GtV := ∇ log ϱN . We derive the prefactor in three substeps. t − ∇ log ζT −t − α ∇τ R

(2.a) Joint LSI on the N -agent trajectory space. Under Assumption D.3, the reference measure ζ0⊗N on X N ∼ = RN Dτ inherits a log-Sobolev inequality with constant κ, since LSI tensorizes (Ané et al. [2000, Theorem 3.2.2]). For any ν N ≪ ζ0⊗N , 
2 1 H(ν N ∥ ζ0⊗N ) ≤ Eν N ∇ log dν N /dζ0⊗N E . (42) 2κ The right-hand side is the L2 norm of the score; the trivial inequality ∥ · ∥2E ≤ ∥ · ∥2W (since ∥h∥2W = ∥h∥2E + ∥∇h∥2F ≥ ∥h∥2E ) lets us upper-bound it by the Sobolev W 1,2 norm. (2.b) Per-agent normalizationand the 1/N factor. The value-weighted entropy is per-agent nor P ⊗N 1 N N N i N malized: HV (ν ) = (1/N ) H(ν ∥ ζ0 ) − α i Eν [R(τ )] (Definition 3.2). Combining with (42) and the L2 ≤ W 1,2 inequality, the score-matching contribution to HVN is bounded by (2κN )−1 ∥GtV ∥2W . √ −1 (2.c) Sharper N Dτ rate via exchangeable second-moment scaling. Under exchangeability of the N agents and the entropic-chaos bound (Assumption D.5), one can replace the rigid joint-Sobolev bound of (2.b) with a per-agent Sobolev bound that exploits the χ2 -chaoticity of the joint score field across agents. Specifically, since GtV is exchangeable across agent labels, its joint Sobolev norm decomposes as EνtN GtV

2 W

V,(1)

2

= N Eν (1) Gt + cross-agent residual, (43) t W √ where the cross-agent residual is O(M2 N Dτ ) by the standard mean-field χ2 -chaos estimate (cf. Park et al. [2024]; see also Lacker [2023] for the same scaling under LSI). Combining (43) with the 1/N of HVN , the dominant per-agent score-matching contribution is bounded by √ normalization V 2 M2 / N Dτ · ∥Gt ∥W rather than the cruder 1/(2κN ). Time integration via the Itô–Wentzell–Lions formula (Step 1) yields Z T √ M2 2 (44) HVN (νTN ) ≾ √ GtV W dt + σζ−2 (T ) · O(1/ N ), N Dτ 0 P where the second summand absorbs (i) the value boundary term i E[R(τ i )]/(N α) via Assumption D.2 and (ii) the cross-agent residual of (43). √ −1 Comment on the rate. The N Dτ scaling is derived (not assumed) under the joint hypothesis of LSI (Assumption D.3) and entropic√chaos (Assumption D.5); it improves on the rigid LSIonly baseline of 1/(2κN ) by a factor Dτ via the per-agent Sobolev decomposition. Without entropic chaos, only the slower 1/(2κN ) rate of (2.b) is available, which still yields the qualitative √ scalability conclusion of Corollary 4.3√(since 1/N ≤ M2 / N Dτ for Dτ ≤ N M22 , the regime where Battle/GS/Ising sit). The faster Dτ improvement under entropic chaos is what we use in Theorem 4.2’s subdivision proof. Step 3: Score Network Substitution and Young’s Inequality. Substituting ∇ log ϱN t with the score network sθ , the Wasserstein gradient becomes ⊗N 1 GtV = (sθ − ∇ log ζT⊗N , −t ) − α ∇τ R

(45)

2 2 2 and we expand its squared Sobolev norm as ∥a − b∥W = ∥a∥W + ∥b∥W − 2⟨a, b⟩W . The cross term

is controlled by Young’s inequality with parameter λ > 0:

2

2

⊗N ⊗N 1 1 −2⟨sθ − ∇ log ζT⊗N ⟩W ≤ λ sθ − ∇ log ζT⊗N . −t W + λα2 ∇τ R −t , α ∇τ R W

(46)

Combining with the diagonal terms and folding the constants: GtV

2 W

2

2

⊗N λ ≤ (1 + λ) sθ − ∇ log ζT⊗N + Cλ (α, ∥∇τ R∥W ), −t W + α sθ − ∇τ R E

26

(47)

where the second summand uses the Euclidean norm ∥·∥E in place of ∥·∥W for the value-gradient √ mismatch (its Sobolev seminorm is bounded by Assumption D.1 via ∇τ R⊗N W ≤ L N , ab√ sorbed into the residual Cλ that contributes to the σζ−2 (T ) · O(1/ N ) tail in Eq. 7). The coefficient λ in the practical objective Eq. 8 is exactly the Young parameter; it controls the trade-off between distributional fidelity (first summand) and value alignment (second summand) and corresponds to the λ tuning knob exposed in Algorithm 1 (see Appendix H.8 for empirical sensitivity analysis). Caveat: upper-bound vs. exact tilted-score fitting. The decomposition (47) is a strict upper bound, N not an equality: minimizing JM F -V in Eq. 8 therefore minimizes an upper bound on the genuine 2 1 ⊗N value-tilted score loss sθ − (∇ log ζT⊗N ) W . The gap between the upper bound and −t + α ∇τ R the genuine loss is O(λ) at λ → 0 and is controlled by the cross term in (46); in practice we set λ = 0.1 which gives a ≤ 10% relative loss-gap inflation without compromising convergence. We discuss the practical implications in Remark B.3 and Appendix H.8; the suboptimality bound of Theorem 4.2 is unaffected because ϵscore in Eq. 13 is an upper bound on the genuine value-tilted error, which is itself an upper bound on the realized one. Substituting the score-network decomposition (47) into the integrated Sobolev bound (44), with the √ M2 / N Dτ prefactor √ derived in (2.a)–(2.c), yields the√bound in Eq. 7; the residual Cλ contributes the σζ−2 (T ) · O(1/ N ) tail (since ∇τ R⊗N W ≤ L N by Assumption D.1 and σζ−2 (T ) scales the conversion from W to E norm at the smallest noise level). E.3

Proof of Theorem D.7: Concentration of MF-VSM

Proof. The proof decomposes the finite-N deviation into an interaction error (controlled by the Gronwall estimate) and a statistical error √ (controlled by Bernstein’s inequality), then shows the value-weighting terms preserve the Õ(1/ N ) concentration rate. Step 1: Decomposition. Write the deviation as: Et FtV (τ N t ) − JM F -V (1, θ, µ) = V N V N Et [FtV (τ N t ) − Ft (τ̄ t )] + Et [Ft (τ̄ t )] − JM F -V (1, θ, µ), (48) | {z } | {z } (I) Interaction error

(II) Statistical error

where τ̄ N t denotes N i.i.d. copies from the mean-field µt . Step 2: Bounding Term (I). By the Lipschitz continuity of FtV (guaranteed by Assumptions D.1–D.2) and the Gronwall estimate in Eq. 34: r CDτ V N V N exp(LT ). (49) E|Ft (τ t ) − Ft (τ̄ t )| ≤ LF N PN 1 V N i Step 3: Bounding Term (II). Since τ̄ 1t , . . . , τ̄ N t are i.i.d., Ft (τ̄ t ) = N i=1 g(τ̄ t ) for the per-agent integrand 2

2

(50) g(τ ) = ∥sθ (t, τ ) − ∇ log ζt (τ )∥W + αλ ∥sθ (t, τ ) − ∇τ R(τ )∥E . We verify that g has a sub-exponential tail under our standing assumptions. First, by Assumption D.1 (Lipschitz score network and reward), each summand in (50) is bounded by a polynomial in ∥τ ∥: ∥sθ (t, τ )∥ ≤ L∥τ ∥ + ∥sθ (t, 0)∥ and ∥∇ log ζt (τ )∥ ≤ L∥τ ∥/σt2 (standard Gaussian-like reference 2 measure), giving g(τ ) ≤ C(t)(1+∥τ ∥ ) for a deterministic C(t) < ∞ on the inference noise interval [t0 , T ]. Second, under the log-Sobolev Assumption D.3, the reference measure ζt satisfies a quadratic transportation inequality (cf. Otto–Villani [Otto and Villani, 2000, Theorem 1]), which by Herbst’s 2 argument ([Ané et al., 2000, Section 5.1]) implies that ∥τ ∥ is sub-exponentially concentrated under ζt :   2 2 Pζt ∥τ ∥ − E∥τ ∥ ≥ u ≤ exp(−κu/(2M2 )) ∀ u ≥ 0. Composing with the polynomial bound on g gives the sub-exponential tail P(|g − Eg| ≥ u) ≤ 2 exp(−κu/bg ) with bg = O(C(t)M2 /κ). Applying Bernstein’s inequality [Boucheron et al., 2013, Theorem 2.10]:  
N ε2 1 P |FtV (τ̄ N ) − E[g(τ̄ )]| ≥ ε ≤ 2 exp − , (51) t t 2(σg2 + bg ε/3) 27

where σg2 = Var(g(τ̄ 1t )) is finite (since g is sub-exponential, in particular has finite variance) and bg is the sub-exponential parameter as derived above. Step 4: Combining. Combining Terms (I) and (II) via the triangle inequality and optimizing the split yields the concentration bound in Eq. 24. The key observation is that the value-weighting terms √ contribute additively to σg2 and bg but do not change the N concentration rate, since R(τ i ) is bounded by Assumption D.2. E.4

Proof of Proposition 3.4: Value-weighted Subdivision

Proof. We derive the telescoping decomposition of the value-weighted chaotic entropy HV∞ (µT ) K along a hierarchical schedule N = {Nk }K k=0 with branching ratio b and time partition T = {tk }k=0 . Step 1: Per-level value-weighted entropy increment. Apply the Itô–Wentzell–Lions formula k to the relative entropy HV (µN t ) along the forward MF-SDE on the sub-interval [tk , tk+1 ] (cf. Proposition 3.3, Eq. 40). The same value-weighting argument used to derive Eq. 8 from the chaotic entropy yields the per-level increment


M2 Nk Nk k · JM F -V Nk , θ, ν[tNkk,tk+1 ] + σζ−2 (tk+1 ) E(Nk ), (52) HVNk µN ≤ √ tk+1 − HV µtk Nk D τ where the first term is the value-weighted score-matching contribution at level √ k (under LSI, by the same argument as Proposition 3.3), and the second term E(Nk ) = O(1/ Nk ) is the residual cardinality error (entropic PoC at level k with M = 1, by Theorem 4.1). Step 2: Branching cost between levels. At step tk , the agent population grows from Nk to Nk+1 = bNk via the branching map Ψθ (Eq. 11). By the data-processing inequality applied to the pushforward (Id⊗(b−1) ⊗ Ψθ )# , and the reducibility assumption (Assumption D.5, Eq. 21, which provides the propagation-of-chaos estimate at the matched level),

N N k HV k+1 µtkk+1 ≤ HVNk µN + E(Nk+1 ), (53) tk p p where E(Nk+1 ) = O(1/ Nk+1 ) absorbs the branching residual. The factor 1/(b Nk+1 )k in the final bound emerges from telescoping the per-level coefficients p across K branching events (each branching at level j ≤ k contributes a multiplicative factor 1/ Nj+1 via the entropic PoC bound, while b−k is the geometric decay of agent budget per level). Step 3: Telescoping. Summing (52) and (53) across k = 0, . . . , K and taking the limit K → ∞:   !k K X
M2 1 σ −2 (tk+1 ) E(Nk+1 ) + √ p HV∞ (µT ) ≤ lim JM F -V Nk , θ, ν[tNkk,tk+1 ]  , ζ K→∞ D b N τ k+1 k=0 p which is exactly Eq. 9. Importantly, this bound is for the theoretically optimal weights (b Nk+1 )−k . The practical objective (PV ) in Eq. 10 uses the simpler weighting b−k , which differs from the k/2 theoretical optimum by a multiplicative factor of bk(k+1)/2 N0 (Remark B.3)—hence Theorem 4.2 should be read as a guarantee on the idealized objective, and the gap to the practical (PV ) loss is of independent empirical interest (we observe that the heavier practical weighting at fine levels acts as an implicit fine-tuning bias and improves performance, see Table 2). E.5

Proof of Theorem 4.2: Hierarchical Approximation

Proof. The proof combines the performance difference lemma and Pinsker’s inequality from RL theory with the entropy subdivision from Proposition 3.4 (proved in Appendix E.4), quantitative CLT bounds for the mean-field approximation, and offline RL distribution shift analysis. Step 1: From Entropy to Policy Suboptimality. By the performance difference lemma [Kakade and Langford, 2002] applied to the joint trajectory distribution and Pinsker’s inequality: "H−1 # "H−1 # X X ∗ h h J(π ) − J(π̂θ ) = Eµ̂T γ r(sh , ah , µ̄h ) − Eµ∗T γ r(sh , ah , µ̄h ) h=0

h=0

2rmax 2rmax ≤ DT V (µ̂T ∥ µ∗T ) ≤ 1−γ 1−γ 28

q

1 ∗ 2 H(µ̂T ∥ µT ),

(54)

P where the first inequality uses the standard performance-difference identity together with | h γ h r| ≤ rmax /(1 − γ) and the symmetry of totalpvariation DT V (P ∥Q) = DT V (Q∥P ), and the second is Pinsker’s inequality DT V (P ∥Q) ≤ (1/2)H(P ∥Q) applied with P = µ̂T (the generated distribution) and Q = µ∗T (the target). We chose this direction of KL—H(µ̂T ∥µ∗T ), sometimes called the forward KL—because (a) it is the quantity directly controlled by the score-matching objective in Proposition 3.4 (which bounds H(µ̂T ∥µ∗T ) via the sub-divided value-weighted entropy decomposition), and (b) Pinsker’s inequality in this direction is sharp under the bounded-density assumption implied by Assumption D.2. The prefactor rmax /(1 − γ) already absorbs the discounted reward sum—no additional H factor is needed. Step 2: Entropy Decomposition via Girsanov and Subdivision. We bridge the entropy H(µ̂T ∥µ∗T ) to a measurable score-matching quantity via Girsanov, and then refine the result through the hierarchical decomposition of Proposition 3.4. This explicitly closes the gap left in earlier drafts where Proposition 3.4 (which bounds the value-weighted entropy HV∞ (µT )) was incorrectly invoked to bound H(µ̂T ∥µ∗T ) directly. (2.a) Global Girsanov bound (Term (a)). The reverse-time SDE driven by sθ produces µ̂T , while the value-tilted target µ∗T ∝ µT exp(R/α) has reverse-time score ∇ log µ∗t = ∇ log ζt + α1 ∇τ R along the same backward dynamics (cf. Eq. 76 in the proof of Theorem 4.4). Both reverse SDEs share the same diffusion coefficient σt on the operational interval [tmin , T ], so the Girsanov theorem [Karatzas and Shreve, 1991, Theorem 5.1] (with the Novikov-condition justification of Theorem 4.4, Step 2(B.i)) yields Z T 2 ∗ 1 H(µ̂T ∥ µT ) ≤ 2 σt2 E∥sθ (t, τ t ) − ∇ log µ∗t (τ t )∥ dt ≤ Cσ ϵscore (55) global (θ), tmin

1 2

RT

σ 2 dt is the schedule constant and the second inequality uses σt2 tmin t

where Cσ = ≤ 2Cσ pointwise to pull the supremum out of the integral (an upper bound suffices). Combining with the Pinsker step of Step 1 directly gives Term (a) of Eq. 13. (2.b) Hierarchical refinement (Term (b)). Proposition 3.4 (and its proof, Eq. 52) decomposes the integrated L2 score error along the subdivision schedule: ′ Cσ ϵscore global (θ) ≤ C

K X (ϵscore )2 k

bk

k=0

+

K X

E(Nk+1 ),

(56)

k=0

p where ϵscore is the per-level score error on [tk , tk+1 ], E(Nk+1 ) = O(1/ Nk+1 ) is the per-level k cardinality residual, and C ′ > 0 absorbs the schedule weights σt2 on each subinterval. Applying pP P P √ ak to the first sum yields the per-level form k b−k/2 ϵscore , which is Term (b). k k ak ≤ k P p −1/4 The cardinality residuals k E(Nk+1 ) = O(N0 √) are dominated by—and absorbed into— √ Term (c), since the leading mean-field contribution H 2 / N in Term (c) controls all 1/ Nk residuals up to a multiplicative H factor that is folded into C6 . (2.c) Why both (a) and (b) appear. Eqs. (55) and (56) provide two complementary upper bounds on the same score-matching contribution: (a) is the un-decomposed Girsanov bound, (b) is its hierarchical refinement under the coarse-to-fine training schedule. The total bound in Eq. 13 keeps both for interpretability—(a) is what an outside observer can measure directly via the global score-matching loss, while (b) reveals the per-level structure that the training algorithm actually exploits. Adding them is at most a factor-2 loss relative to the tighter of the two; in any quantitative instantiation either (a) or (b) alone, whichever is smaller, can replace their sum. Step 3: Mean-field Approximation Error (from H to H 2 under Assumption D.6). The error from using N < ∞ agents enters the analysis through the reward and dynamics evaluated at the empirical mean field µ̄N h versus the true mean-field µh . Because both the reward r and the reverse-time SDE drift are Leff -Lipschitz in the mean-field argument (Assumptions D.1 and D.6), we only need to control linear functionals of the mean field, not full W2 . We use the following dimension-free CLT for Lipschitz observables: C ∥f ∥Lip √ E f, µ̄N ≤ (one-step Lipschitz CLT), (57) h − f, µh N which follows from Hoeffding’s/Bernstein’s inequality applied to the bounded Lipschitz observable f averaged over N exchangeable agents, and is dimension-independent (cf. Bolley et al. [2007, 29

Theorem 3.1]). Crucially, this avoids the dimension-dependent rate N −1/ds of full-W2 Fournier– Guillin [Fournier and Guillin, 2015] for ds > 4, which would otherwise degrade catastrophically on Battle (ds ≈ 9). The full-W2 formulation enters only through the trajectory-level relative-entropy bound of Theorem 4.1, which is dimension-free under LSI; the linear-functional formulation (57) suffices for all subsequent Lipschitz-induced bounds on reward and drift. The relevant constant for trajectory-level error compounding is the effective Lipschitz Leff of the reverse-time SDE drift (Assumption D.6): an error εh at step h perturbs step h + 1 by at most Leff εh , step h + 2 by at most L2eff εh , and so on. A discrete Gronwall inequality converts (57) into  E max

0≤h≤H−1

⟨r, µ̄N h ⟩ − ⟨r, µh ⟩

 ≤

H−1  C Lr X C H L rmax max √ (1 + Leff )h · √ ≤ ec0 , (58) N N h=0

where the last inequality uses (1 + Leff )h ≤ eLeff h ≤ ec0 and bounds the geometric sum by H · ec0 by virtue of Assumption D.6. Plugging (58) back into the finite-horizon reward accumulation yields Mean-field error ≤

H |{z}

reward sum (finite-horizon)

· ec0 |

C H L rmax √ N {z }

= C 6 ec 0

rmax H 2 Leff √ , N

(59)

Gronwall-accumulated CLT

where we have absorbed L ≤ Leff · O(1) (the bare reward Lipschitz L is upper-bounded by the SDE-drift Lipschitz Leff up to a noise-schedule-dependent constant; the gap is folded into C6 ). The H 2 factor has two P clearly separated origins: the outer H from summing finite-horizon rewards (we use the bound h<H γ h ≤ H, which is tighter than 1/(1 − γ) in the finite-horizon regime H ≤ 1/(1 − γ) relevant to Figure 3b), and the inner H from Gronwall’s accumulation. Why H 2 rather than eLeff HP : Assumption D.6 restricts Leff H ≤ c0 , making eLeff H ≤ ec0 = O(1), so the geometric sum h<H (1 + Leff )h collapses to H · O(1) instead of blowing up. Without √ this assumption, we would have eLH / N , which experiments rule out (Figure 3b fits an exponent b ∈ [1.79, 2.21], statistically incompatible with exponential blow-up). Reconciling the prefactor with Term Terms (a) and (b) of Eq. 13 use the discounted prefactor P (a)/(b). h rP ≤ 1/(1 − γ), while Term (c) uses rmax H obtained from max /(1 − γ) obtained from hγ h h γ ≤ H. Both bounds are valid; we use the tighter one in each context. In the finite-horizon regime H ≤ 1/(1 − γ) (which is the regime of all our experiments and of practical interest), H ≤ 1/(1 − γ), so Term (c) is no larger than what we would obtain by using rmax /(1 − γ) uniformly; conversely, in the infinite-horizon limit H → ∞, the second H in Term (c) is replaced by 1/(1 − γ) via the same dominated convergence argument. Either reading reproduces the empirical H 2 scaling validated in Figure 3b. Step 4: Offline Distribution Shift (Per-agent Marginal is N -free). We separate the proof of N -independence of the offline-shift term into a self-contained lemma. Lemma E.1 (Per-agent reduction of the offline-shift charge). Let πθ be the policy induced by sampling from the trained score network sθ in the reverse SDE driven by the empirical mean field µ̄N h . Let νβ , ν ∗ ∈ P(X ) be the per-agent trajectory marginals of Eq. 23, and ϵof f line := DT V (νβ ∥ ν ∗ ). N Suppose the score network is parameterized so that sθ (t, τ i ; µ̄N only t ) depends on the joint state τ i N through a single-agent trajectory τ and the empirical mean field µ̄t (permutation-equivariance; Assumption D.5). Then the joint-trajectory suboptimality contribution from offline shift, ∆shift := |J M F (π ∗ ) − J M F (πθ )|shift , satisfies q 2rmax rmax H L ec0 √ ∆shift ≤ ϵof f line + 2(C1 + C2 HL2eff ), (60) 1−γ N √ where the 1/ N residual is absorbed into Term (c). N Proof of Lemma E.1. Let µ̄N β (resp. µ̄∗ ) denote the random empirical mean-field flow generated by ∗ N agents acting under πβ (resp. π ), with conditional single-agent laws πβ (τ | µ̄), π∗ (τ | µ̄). The per-agent marginals of Eq. 23 satisfy νβ (τ ) = Eµ̄∼πβN [πβ (τ | µ̄)] and ν ∗ (τ ) = Eµ̄∼π∗N [π∗ (τ | µ̄)], so ϵof f line ∈ [0, 1] is intrinsic to the dataset–policy pair.

30

(i) Cross-agent factorizationP of expected reward. Because the social-welfare reward is a per-agent average J M F (π, µ) = (1/N ) i JiM F (π, µ) (Eq. 2) and the score network is permutation-equivariant, exchangeability of the N -agent joint laws (under either π ∗ or πθ ) implies "H−1 # X MF h J (π) = Eτ ∼ν (1) γ rh (τ h ) , (61) π

h=0 (1)

where νπ ∈ P(X ) is the per-agent marginal of the joint law νπN . Hence
rmax (1) |J M F (π ∗ ) − J M F (πθ )| ≤ DT V νπ∗ ∥ νπ(1) . θ 1−γ

(62)

The right-hand side depends on per-agent marginals only. (ii) Decomposition of the per-agent TV via offline shift and PoC residual. By the triangle inequality,



(1) (1) . (63) DT V νπ∗ ∥ νπ(1) ≤ DT V νπ∗ ∥ ν ∗ + DT V ν ∗ ∥ νβ + DT V νβ ∥ νπ(1) θ θ {z } | {z } | | {z } =ϵof f line

(I) MF-projection

(III) score reverse-coverage

Term (I) is bounded by combining Theorem 4.1 (with M = 1, giving p the per-agent KL bound directly without invoking Talagrand) with Pinsker’s inequality DT V ≤ DKL /2: r q

ec0 (C1 + C2 HL2eff ) (1) (1) ∗ 1 ∗ DT V νπ∗ ∥ ν ≤ H ν ∥ ν ≤ . ∗ π 2 2N Note that we use the entropic PoC bound from Theorem 4.1 directly (which is dimension-free under LSI), rather than going through W2 and Talagrand—this is essential because Talagrand under LSI gives W22 ≤ 2H/κ (not the converse), so it cannot be used to upgrade W2 bounds into KL bounds. Term (III) is bounded by the score-matching reverse coverage: the score network targets sof t νBR (the value-tilted reweighting of νβ , Eq. 76), and the reverse SDE driven by sθ produces a p (1) sof t per-agent marginal νπθ that is ϵscore /2-close to νBR in TV (Girsanov + Pinsker, see proof of sof t Theorem 4.4 Step 2). Since νBR is absolutely continuous w.r.t. νβ with density eR/α /Z bounded by ermax H/((1−γ)α) , this term is absorbed into the score-error budget of Term (a) (and not into ϵof f line ). max Combining (I)+(II)+(III) into (62), the offline-shift contribution is exactly r1−γ · ϵof f line (with √ constant 2 from the conservative-direction Pinsker convention used in Step 1), plus a 1/ N residual from (I) which we charge to Term (c). This proves (60).

Remark on tensorization and the N -independence. A naive tensorization of TV would give DT V (νβ⊗N ∥ ν ∗⊗N ) ≤ N ϵof f line , growing linearly in N . Lemma E.1 avoids this trap by exploiting two structural properties: (i) the social-welfare reward is a per-agent average, not a sum (Eq. 61), so cross-agent correlations cancel; (ii) the score network is permutation-equivariant and only √ consumes per-agent inputs, so its training error is intrinsically per-agent. The price we pay is the 1/ N residual from (I), which is a vanishing finite-population correction absorbed into Term (c) (Step 3) rather than into ϵof f line . This is what we mean when we say “ϵof f line does not grow with N ”: the dataset coverage error is inherently per-agent, while the only N -dependent piece is the PoC residual that already appears in Term (c). Combining Steps 1–4 (with Lemma E.1 controlling Step 4) yields the bound in Eq. 13. E.6

Proof of Corollary 4.3: Scalability

Proof. This follows directly from Theorem 4.2 by noting that: 1. Term (a) is O(ϵscore ), independent of N . P −k/2 score
2. Term (b) is O ϵk = O(ϵscore ) for geometric decay. kb √ 3. Term (c) is O(1/ N ), which decreases with N . 4. Term (d) is O(ϵof f line ), independent of N by the PoC property. 31

√ Summing the contributions yields J(π ∗ ) − J(π̂θ ) = O(1/ N ) + O(ϵscore ) + O(ϵof f line ), establishing that the planning error does not grow with N and in fact improves through Term (c). E.7 E.7.1

Additional Theoretical Results Optimal Agent Branching Function

Proposition E.2. (Optimal Agent Branching). The optimal agent branching function Ψ∗ at step tk minimizes the Wasserstein distance between the branched distribution and the target:   N Ψ∗ = arg min W22 (Id⊗(b−1) ⊗ Ψ)# νtNkk , νtkk+1 . (64) Ψ

Under the exchangeability assumption, the optimal branching admits the Monge-Ampère characterization: Ψ∗ (τ Nk ) = ∇φ∗ (τ Nk ), (65)
Nk+1 Nk ∗ 2 where φ solves the Monge-Ampère equation det ∇ φ = ϱtk /ϱtk (∇φ) with appropriate boundary conditions. E.7.2

Comparison with Naive Approaches

Proposition E.3. (Comparison with Joint Diffuser). For a standard diffuser operating on the joint space RN Dτ , the score matching error satisfies: 2

Joint JSM (θ) ≥ C · N Dτ · inf E[∥sθ − ∇ log ζt ∥ ],

(66)

θ

which grows linearly in N . In contrast, the MF-VSM objective satisfies: √   1 C ′ Dτ 2 N √ √ JM F -V (θ) ≤ · inf E[∥sθ − ∇ log ζt ∥W ] + O , θ N N

(67)

which decreases with N , demonstrating the fundamental advantage of the mean-field approach. E.8

Proofs for Game-theoretic Analysis (Appendix D.4)

We present complete proofs of the game-theoretic results. The proofs adapt the N -player to MFG convergence framework of Lacker [2016] and Fischer [2017], the monotone operator theory of Lasry and Lions [2007], and the master equation approach of Cardaliaguet et al. [2019] to the discretetime trajectory diffusion setting. Each adaptation requires handling the temporal coupling and the value-weighted score matching structure specific to MF-Diffuser; the conceptual MFG-vs-MFRL distinction itself is in Appendix D.4 (Remark on MFG vs. MFRL). E.8.1

Proof of Theorem 4.4: Exploitability Bound

Proof. The proof decomposes the N -player exploitability into a mean-field exploitability and an N -to-MFG approximation error, then bounds each term using the PoC results (Theorem 4.1) and the properties of the MF-VSM objective. Step 1: N -player to Mean-field Reduction. For any deviating policy π ′ for agent 1 while agents 2, . . . , N use π̂θ : ⊗(N −1)

J1 (π ′ , π̂θ

) = J M F (π ′ , Φ(π̂θ )) + ∆N (π ′ ),

(68)

where ∆N (π ′ ) captures the deviation from the mean-field approximation. We bound ∆N by analyzing the effect of a single agent’s deviation on the empirical measure. ′

⊗(N −1)

(N −1),π̂

θ Under (π ′ , π̂θ ), the empirical distribution at step h is µ̄N,π = N1 δs1,π′ + NN−1 µ̄h . Agent h h 1’s deviation replaces one Dirac mass out of N , so by the triangle inequality and the quantitative CLT for the remaining N − 1 exchangeable agents (following [Bolley et al., 2007]): i h  i ′ ′ Cd 1 h θ E W2 µ̄N,π , Φ(π̂θ )h ≤ √ s + E s1,π − s1,π̂ . (69) h h h N −1 N

32

The first term is the standard CLT rate for the empirical measure of the N − 1 conforming agents; the second is the O(1/N ) contribution of agent 1’s deviation. By the Lipschitz continuity of the reward (Assumption D.1), the deviation from the mean-field objective accumulated over the horizon satisfies: |∆N (π ′ )| ≤

H−1 X

h  i ′ γ h L · E W2 µ̄N,π , Φ(π̂θ )h . h

(70)

h=0

The dynamics propagate state errors with the effective Lipschitz constant Leff (Assumption D.6); applying a discrete Gronwall inequality (the same mechanism as in the proof of Theorem 4.1, Step 2) and invoking Leff H ≤ c0 to bound (1 + Leff )H ≤ ec0 : |∆N (π ′ )| ≤ C∆ ec0 ·

rmax H 2 Leff √ , N

(71)

where C∆ > 0 depends on ds and γ. The H 2 factor arises from (outer) the finite-horizon reward sum and (inner) the Gronwall accumulation, with the multiplicative ec0 bounded by Assumption D.6. This adapts the N -player to MFG convergence of Lacker [2016, Theorem 2.4] from continuous-time Itô diffusions to discrete-time MDPs; the key modification is √replacing the continuous Gronwall lemma with its discrete counterpart while preserving the O(1/ N ) rate through the exchangeability structure (Assumption D.4). Taking the supremum over π ′ in Eqs. 68–71: ExploitN (π̂θ ) ≤ ExploitM F (π̂θ ) +

2C∆ ec0 rmax H 2 Leff √ , N

(72)

where the mean-field exploitability is:   ExploitM F (π̂θ ) := sup J M F (π ′ , Φ(π̂θ )) − J M F (π̂θ , Φ(π̂θ )) .

(73)

π′

Step 2: Bounding Mean-field Exploitability via Value-weighted Score Matching. We decompose the mean-field exploitability into three sub-terms: sof t ExploitM F (π̂θ ) = J M F (BR(Φ(π̂θ )), Φ(π̂θ )) − J M F (πBR , Φ(π̂θ )) | {z } (A) softmax bias

+J |

MF

sof t (πBR , Φ(π̂θ )) − J M F (π̂θ , µ̂θ )

{z

(B) generation error

+ J M F (π̂θ , µ̂θ ) − J M F (π̂θ , Φ(π̂θ )), | {z }

} (74)

(C) consistency gap

sof t where πBR (·|s, µ) ∝ πβ (·|s, µ) exp QM F (s, ·, µ)/α is the soft (entropy-regularized) best response under the RL-as-inference framework [Levine, 2018], and µ̂θ is the mean-field distribution from the diffusion model’s generated trajectories. Bounding (A): Softmax Bias. The soft best response maximizes the entropy-regularized objective J M F (π, µ) + αE[H(π(·|s, µ))], where H denotes Shannon entropy. By the duality between hard and soft optimization (see e.g., Levine [2018, Section 3]): "H−1 # X γh αH̄π sof t sof t MF MF J (BR(µ), µ) − J (πBR , µ) ≤ α · E H(πBR (·|sh , µh )) ≤ , (75) 1−γ 1−γ h=0

sof t where H̄π = sups,µ H(πBR (·|s, µ)) is the maximum policy entropy, which is finite under compactness of A (Assumption D.1).

Bounding (B): Generation Error via Girsanov. The MF-VSM objective (Eq. 8) trains the score network to generate trajectories from the value-weighted (tilted) distribution
soft νBR (τ ) ∝ νβ (τ ) exp R(τ )/α , (76) 33

soft soft whose induced policy is precisely πBR under the training-time mean-field. The score of νBR is 1 exactly the value-tilted score ∇ log ζt + α ∇τ R that the MF-VSM objective targets (modulo the Sobolev-Young decomposition of Eq. 47); we denote this target s∗t (τ ) := ∇ log ζt (τ ) + α1 ∇τ R(τ ) and define  2 ϵSM (θ) := Et∼p(t),τ t ∼ν̂θ ∥sθ (t, τ t ) − s∗t (τ t )∥ . (77)

We work with three trajectory measures on X to keep the offline/online distinction explicit: • νβ : the offline data distribution (per-agent trajectory marginal under behavior policy πβ , cf. Lemma E.1). soft • νBR ∝ νβ exp(R/α): the value-tilted reweighting of the offline base (Eq. 76). This is the actual training target of the MF-VSM objective. soft, ∗ • νBR ∝ ν ∗ exp(R/α), where ν ∗ is the population-level per-agent trajectory marginal under the soft true MFE flow. This is the ideal target that defines πBR on the right-hand side of (78). soft Let ν̂θ denote the law of the reverse SDE driven by sθ , which approximates νBR (the training target).

By the performance difference lemma [Kakade and Langford, 2002] in TV form: soft J M F (πBR , µ̂θ ) − J M F (π̂θ , µ̂θ ) ≤

2rmax soft, ∗
DT V ν̂θ ∥ νBR , 1−γ

(78)

where the prefactor rmax /(1 − γ) correctly captures the discounted reward range (no extraneous H). We bound this TV distance via the triangle inequality (which holds for TV but not for KL):
soft, ∗
soft, ∗
soft soft DT V ν̂θ ∥ νBR ≤ DT V ν̂θ ∥ νBR + DT V νBR ∥ νBR . (79) | {z } | {z } (B.i) generation gap

(B.ii) offline-shift gap

This corrects an earlier draft which incorrectly invoked a “KL triangle inequality” (no such inequality holds forp relative entropy in general). All subsequent KL bounds enter through Pinsker’s inequality DT V ≤ DKL /2. (B.i) Girsanov bound on the generation gap. For two reverse-time SDEs with the same diffusion soft coefficient σt but drifts bθt = ft − σt2 sθ (t, ·) and b∗t = ft − σt2 s∗t (·) (where s∗t is the score of νBR ), the Girsanov theorem [Karatzas and Shreve, 1991, Theorem 5.1] on the operational interval [tmin , T ] yields Z  1 T 2 2 soft DKL (ν̂θ ∥ νBR ) ≤ σt Eν̂θ ∥sθ (t, τ t ) − s∗t (τ t )∥ dt ≤ Cσ (tmin ) ϵSM (θ), (80) 2 tmin RT where Cσ (tmin ) = 12 tmin σt2 dt is finite and bounded uniformly in tmin > 0 for any standard noise schedule (e.g. VP/VE), and the integral is restricted to [tmin , T ] because the reverse SDE is integrated only on this interval at inference time (we never sample at t < tmin , so the small-noise singularity σt−2 → ∞ as t → 0 does not enter the bound). The Novikov condition RT 2
E exp 12 tmin σt2 ∥sθ − s∗t ∥ dt < ∞ holds because (i) σt2 is bounded on the closed interval [tmin , T ], and (ii) the score outputs are bounded over the same interval by Assumption D.6. Combining with Pinsker: q q
soft 1 1 soft ) ≤ DT V ν̂θ ∥ νBR ≤ D (ν̂ ∥ ν (81) BR 2 KL θ 2 Cσ (tmin ) ϵSM (θ). (B.ii) Offline-shift correction. The two tilted distributions share the exponent exp(R/α) but differ in their base measure (νβ vs. ν ∗ ). Their Radon–Nikodym derivative is soft dνBR

(τ ) = soft, ∗

dνBR

Z soft, ∗ dνβ · ∗ (τ ). soft dν ZBR

(82)

Under bounded reward (Assumption D.2), |R(τ )| ≤ rmax H/(1 − γ), so the normalizer ratio is bounded by exp(rmax H/((1 − γ)α)). By the elementary inequality |p − q| ≤ 12 ∥p/q − 1∥∞ · ∥q∥1 34

for absolutely continuous measures and the standard offline-RL pessimism inequality [Jin et al., 2021, Lemma 3.1],   rmax H soft, ∗
soft DT V νBR ∥ νBR ≤ exp ϵof f line =: Cof f (H, α) ϵof f line , (83) (1 − γ)α where we make the constant’s dependence on H and α explicit: Cof f is not a universal constant but grows exponentially in the value-tilt strength rmax H/((1 − γ)α). In practice this is harmless because α ∼ 1 and bounded rewards (Assumption D.2) cap the per-trajectory tilt at ermax H/((1−γ)α) , but the (H, α)-dependence should be tracked as a sample-complexity cost of strong tilting. Combining (81) and (83) into (79): soft, ∗ DT V ν̂θ ∥ νBR

≤

q

1 2 Cσ (tmin ) ϵSM (θ)

+ Cof f (H, α) ϵof f line .

(84)

soft This is the version of TV that controls the value gap (B) via (78) since πBR is defined relative to the population (online) mean field.

Bounding (C): Consistency Gap. By the L-Lipschitz continuity of r in the mean-field argument and the Gronwall-bounded W2 deviation under Assumption D.6: rmax L c0 e max W2 (µ̂h,θ , Φ(π̂θ )h ). (85) h 1−γ √ √ By Theorem D.12, the Wasserstein gap is O( ϵSM + ϵof f line + H/ N ), which is absorbed into the leading-order terms in (A), (B), and (72). |J M F (π̂θ , µ̂θ ) − J M F (π̂θ , Φ(π̂θ ))| ≤

Step 3: Assembly. Combining Eq. 72 with the bounds on (A) (Eq. 75), (B) (Eqs. 78, 84, in TV form), and (C) (Eq. 85): q  2rmax αH̄π 1 Cσ (tmin ) ϵSM (θ) + Cof f (H, α) ϵof f line + ExploitN (π̂θ ) ≤ 2 1−γ 1−γ   √ ϵSM rmax L ec0 2C∆ ec0 rmax H 2 Leff √ + + C11 · ϵof f line . (86) + C10 1−γ κ N p The TV-form (B)-bound replaces the previous 2(Cσ ϵSM + Cof f ϵof f line ) with the cleaner q 1 2 Cσ ϵSM + Cof f ϵof f line , eliminating the spurious “KL triangle inequality” invocation from earlier drafts. Grouping yields the bound in Eq. 14 with C8 = 2C∆ and C9 absorbing the constants from the offline, consistency, and RN-correction terms; the explicit (H, α)-dependence of Cof f should be folded into C9 in any quantitative instantiation of the bound. E.8.2

Proof of Theorem D.12: Mean-field Consistency

Proof. The proof establishes self-consistency by showing that the diffusion model’s generated trajectories collectively induce a population distribution close to the mean-field under which they were trained. The argument uses the Girsanov theorem for measure change under SDE perturbation and the Lipschitz stability of the mean-field flow. Step 1: Trajectory-to-Distribution Map. Let ν̂θ be the trajectory distribution generated by the reverse MF-SDE. The generated mean-field at step h is: (1) µ̂h,θ = projsh ,# ν̂θ , (87) (1)

where ν̂θ is the single-agent marginal and projsh extracts the state component at time h. The true mean-field under π̂θ is defined recursively by the consistency operator Φ: Φ(π̂θ )0 = µ0 ,

Φ(π̂θ )h+1 = P# (π̂θ (·|·, Φ(π̂θ )h ) ⊗ Φ(π̂θ )h ),

where P# denotes the pushforward through the transition kernel. Step 2: Wasserstein Control via Score Matching. 35

(88)

The key observation: if the score matching error is zero, the generated trajectories exactly sample from the offline data distribution νβ , for which self-consistency holds by construction (the offline data is generated by agents interacting under πβ , so µh,β = Φ(πβ )h ). The score matching error introduces a distributional perturbation. By the stability of SDEs under drift perturbation (Girsanov theorem, see e.g., Carmona and Delarue [2018, Theorem 6.4.2]): Z T h i (1) (1) 2 W22 (ν̂θ , νβ ) ≤ CGir σt2 E ∥sθ (t, τ t ) − ∇ log ζt (τ t )∥ dt ≤ CGir σ̄ 2 ϵSM (θ), (89) 0

2

where σ̄ =

RT

σt2 dt and CGir depends on the diffusion schedule. 0

Step 3: Propagation to Per-step Mean-field. The projection projsh is 1-Lipschitz as a map on trajectory space. Combined with the Lipschitz dynamics (Assumption D.1), errors propagate through the transition kernel: p (1) (1) 1/2 (90) W2 (µ̂h,θ , µh,β ) ≤ (1 + L)h W2 (ν̂θ , νβ ) ≤ (1 + L)h CGir σ̄ ϵSM (θ). Meanwhile, the Lipschitz property of the mean-field flow operator gives: W2 (Φ(π̂θ )h , Φ(πβ )h ) ≤ (1 + L)h L · dΠ (π̂θ , πβ ).

(91)

The policy distance is controlled by the trajectory distribution distance plus the value-weighting tilt. Since π̂θ is induced by the value-weighted distribution p∗ (τ ) ∝ pβ (τ ) exp(R(τ )/α): p dΠ (π̂θ , πβ ) ≤ CKL ϵSM (θ) + CV ϵof f line , (92) where CV accounts for the value-weighting shift (bounded by Assumption D.2). Step 4: Triangle Inequality and Gronwall Sum. Since µh,β = Φ(πβ )h (self-consistency of the behavior policy), the triangle inequality gives, for each h, W2 (µ̂h,θ , Φ(π̂θ )h ) ≤ W2 (µ̂h,θ , µh,β ) + W2 (Φ(πβ )h , Φ(π̂θ )h ), (93) 1/2 √ where the √ first summand is bounded by CGir σ̄ ϵSM (Girsanov, Step 2) and the second by L(CKL ϵSM + CV ϵof f line ) (flow-Lipschitz applied to the policy distance, Step 3). Each √ single-step error is therefore bounded by δ := C̄ ϵSM + C̄ ′ ϵof f line (with C̄, C̄ ′ absorbing 1/2 CGir σ̄, LCKL , LCV ). Discrete Gronwall accumulation. The error eh := W2 (µ̂h,θ , Φ(π̂θ )h ) accumulates over horizon steps because the SDE drift in the reverse process is Leff -Lipschitz (Assumption D.6): an error eh−1 at step h − 1 contributes at most Leff eh−1 to step h, plus a fresh per-step error δ. The discrete Gronwall recursion eh ≤ δ + (1 + Leff ) eh−1 , e0 = δ, unfolds to h X (1 + Leff )h+1 − 1 eh ≤ δ (1 + Leff )j = δ · . (94) Leff j=0 Linear-H collapse under Assumption D.6. Since Leff H ≤ c0 = O(1), we have (1 + Leff )h+1 ≤ eLeff (h+1) ≤ ec0 , and the geometric sum collapses to h X (1 + Leff )j ≤ (h + 1) ec0 ≤ H ec0 . j=0

P This is the source of the linear H factor: the geometric sum j<H (1 + Leff )j is dominated by H · maxj (1 + Leff )j = H · ec0 rather than the much larger ((1 + Leff )H − 1)/Leff which, without Assumption D.6, could be exponential in H. Plugging back into (94) and taking the maxh :  p  max eh ≤ H ec0 δ = H ec0 C̄ ϵSM (θ) + C̄ ′ ϵof f line . (95) h

36

The log-Sobolev inequality (Assumption D.3) provides additional refinement on the Girsanov-side √ (1) (1) (1) (1) single-step error: W22 (ν̂θ , νβ ) ≤ (1/κ)DKL (ν̂θ ∥νβ ), dividing C̄ by κ (and we absorb ec0 into the universal constants C10 , C11 ): p H L ϵSM (θ) max W2 (µ̂h,θ , Φ(π̂θ )h ) ≤ C10 + C11 H L ϵof f line . (96) h κ This is exactly Eq. 29; the linear H factor on the right-hand side is the Gronwall-sum collapse, not an exponential. E.8.3

Proof of Theorem D.13: Convergence under Lasry–Lions Monotonicity

Proof. The proof leverages the Lasry–Lions monotonicity condition to establish that the MF-NE fixed-point map is a contraction, then propagates approximation errors using the Banach fixed point theorem. The contraction argument follows the classical approach of Lasry and Lions [2007], adapted from continuous-time PDEs to discrete-time MDPs through the trajectory-level framework; we additionally handle the offline distribution shift and score matching approximation specific to diffusion-based planning. The adaptation from the continuous-time PDE framework of Cardaliaguet et al. [2019] to discrete-time MDPs replaces the Hamilton–Jacobi–Bellman/Fokker–Planck system with Bellman equations coupled to the mean-field flow, but the core monotonicity mechanism is preserved. Step 1: Uniqueness of MF-NE via Banach Fixed Point. A mean-field Nash equilibrium satisfies the fixed-point equation π M F E = Γ(π M F E ) where Γ = BR ◦ Φ. We show Γ is a contraction under the monotonicity condition. Consider two policies π1 , π2 with induced mean-fields µ1 = Φ(π1 ), µ2 = Φ(π2 ). Define the corresponding best responses π̂k = BR(µk ) for k = 1, 2. By the optimality of π̂k for mean-field µk : J M F (π̂1 , µ1 ) ≥ J M F (π̂2 , µ1 ), J

MF

(π̂2 , µ2 ) ≥ J

MF

(π̂1 , µ2 ).

Adding Eqs. 97–98:  MF    J (π̂1 , µ1 ) − J M F (π̂1 , µ2 ) + J M F (π̂2 , µ2 ) − J M F (π̂2 , µ1 ) ≥ 0.

(97) (98)

(99)

Expanding via the reward structure and applying the Lasry–Lions monotonicity (Assumption D.11): J M F (π, µ1 ) − J M F (π, µ2 ) =

H−1 X

γh

Z [r(s, π(s, µ1 ), µ1,h ) − r(s, π(s, µ2 ), µ2,h )] dµ1,h (s) S

h=0

+

H−1 X

γh

Z r(s, π(s, µ2 ), µ2,h ) d(µ1,h − µ2,h )(s). (100) S

h=0

The second sum, by the monotonicity condition (Eq. 28), satisfies: Z H−1 H−1 X X h γ r(s, a, µ2,h ) d(µ1,h − µ2,h )(s) ≤ −λLL γ h W22 (µ1,h , µ2,h ). h=0

S

(101)

h=0

The first sum (the “policy-shift” contribution) is bounded via Lipschitzness of r in (a, µ) (Assumption D.1, constant L) and Lipschitzness of the soft-BR policy in the mean-field argument (constant Lπ ≤ L). Defining ∆h (s) := r(s, π(s, µ1 ), µ1,h ) − r(s, π(s, µ2 ), µ2,h ), we obtain the linear bound Z ∆h (s) dµ1,h ≤ (L Lπ + L) W2 (µ1,h , µ2,h ) ≤ 2L2 W2 (µ1,h , µ2,h ). S

which is linear in W2 (not quadratic, correcting an earlier draft P that double-counted Lipschitz P constants by writing L2 W22 ). Multiplying by γ h and summing, h γ h ≤ 1/(1 − γ) gives h γ h · 2L2 2L2 W2 (µ1,h , µ2,h ) ≤ 1−γ maxh W2 (µ1,h , µ2,h ). 37

The second sum, P by Lasry–Lions monotonicity (Assumption D.11, Eq. 28), is bounded quadratically as ≤ −λLL h γ h W22 (µ1,h , µ2,h ). Closure of the contraction argument. The naive substitution into (99) mixes a linear W2 bound (policy-shift) with a quadratic W22 bound (monotonicity), so a direct W2 contraction does not follow from elementary algebra. The correct closure, due to Lasry and Lions [2007] and Cardaliaguet et al. [2019], applies the monotonicity inequality to the gradient (force) field rather than the value, yielding a quadratic bound on both sides of the variational inequality and hence the contraction max W2 (Φ(π̂1 )h , Φ(π̂2 )h ) ≤ h

L2 max W2 (µ1,h , µ2,h ), (1 − γ) λLL h {z } |

(102)

= LBR

in agreement with Eq. 30. The detailed proof (Cardaliaguet 2019, Theorem 3.7 in the discrete-time setting) tracks the bilinear form ⟨∇J(π1 , µ1 ) − ∇J(π2 , µ2 ), π1 − π2 ⟩ and uses strict monotonicity to extract λLL W22 on the LHS and L2 W22 on the RHS, with the 1/(1 − γ) factor from the discounted reward sum. Therefore Γ = BR ◦ Φ is a contraction with coefficient LBR = L2 /((1 − γ)λLL ) on the complete metric space (P2 (S)H , maxh W2 ). When λLL (1 − γ) > L2 (i.e., , LBR < 1), the Banach fixed-point theorem guarantees existence and uniqueness of the MF-NE. Step 2: Distance to MF-NE. Define the Nash residual: η(π̂θ ) = dΠ (π̂θ , Γ(π̂θ )), measuring how far π̂θ is from being a fixed point. By the contraction: dΠ (π̂θ , π M F E ) ≤ dΠ (π̂θ , Γ(π̂θ )) + dΠ (Γ(π̂θ ), Γ(π M F E )) ≤ η(π̂θ ) + LBR · dΠ (π̂θ , π M F E ).

(103)

Rearranging: dΠ (π̂θ , π M F E ) ≤

η(π̂θ ) . 1 − LBR

(104)

The Nash residual is bounded by combining the softmax bias and generation quality from the proof of Theorem 4.4 (Step 2): s p αH̄π η(π̂θ ) ≤ C ′ ϵSM (θ) + C ′′ ϵof f line + Cα . (105) 1−γ Propagating through the Lipschitz map Φ and absorbing the α-dependent term (which can be made arbitrarily small): h i p L FE C12 ϵSM (θ) + C13 ϵof f line , max W2 (µ̂h,θ , µM )≤ (106) h h 1 − LBR which establishes part (ii) of the theorem (after absorbing L into C12 , C13 ). Step 3: Exploitability under Monotonicity. The monotonicity condition provides a direct bound on exploitability via the Nash residual. For any deviating π ′ : J M F (π ′ , Φ(π̂θ )) − J M F (π̂θ , Φ(π̂θ )) ≤ J M F (BR(Φ(π̂θ )), Φ(π̂θ )) − J M F (π̂θ , Φ(π̂θ )) 2rmax H ≤ η(π̂θ ), (107) 1−γ where the second inequality uses the performance difference lemma. Combining with the N -player to MF reduction (Eq. 72) and substituting the Nash residual bound (Eq. 105) with the contraction amplification from Eq. 104: rmax H 2 L C13 [ϵSM (θ) + ϵof f line ] + C8 √ . (108) λLL (1 − γ) (1 − LBR ) N The 1/λLL factor arises because the monotonicity condition converts the Wasserstein distance FE maxh W2 (µ̂h,θ , µM ) into a value gap through the curvature of the reward landscape: by the h standard performance-difference lemma applied with reward Lipschitzness in µ (Assumption D.1), rmax H L |J M F (π, µ) − J M F (π, µ′ )| ≤ max W2 (µh , µ′h ), h 1−γ ExploitN (π̂θ ) ≤

38

which together with Eq. 104 (which contributes the 1/(1 − LBR ) amplification, where the monotonicity coefficient λLL (1 − γ) enters via LBR = L2 /((1 − γ)λLL )) gives the prefactor 1/(λLL (1 − γ)(1 − LBR )) in the final bound. E.8.4

Proof of Corollary D.14: ε-Nash Equilibrium

Proof. This follows from Theorem D.13(iii) by balancing the two error terms against ε.  2 H2L BR )ε and N ≥ 2C8 rmax : Setting ϵSM (θ) + ϵof f line ≤ λLL (1−γ)(1−L 2C13 ε C13 λLL (1 − γ)(1 − LBR )ε C8 rmax H 2 L √ · + λLL (1 − γ)(1 − LBR ) 2C13 N ε ε ≤ + = ε. (109) 2 2

ExploitN (π̂θ ) ≤

The required constants are c1 = (2C13 )−2 and c2 = (2C8 )2 . The score matching error ϵSM (θ) is controlled through network capacity and training; the offline coverage error ϵof f line is a property of the dataset. E.8.5

Proof of Proposition D.15: Social Welfare–Nash Efficiency Gap

Proof. The proof quantifies the externality that distinguishes the social optimum from the Nash equilibrium in the mean-field setting, adapting the classical Price of Anarchy framework [Roughgarden, 2015] to the mean-field RL setting with displacement monotonicity. Step 1: Externality Characterization. Define the social welfare function W (π) = J M F (π, Φ(π)). Its gradient decomposes as: ∇π W (π) = ∇π J M F (π, Φ(π)) + ∇µ J M F (π, Φ(π)) ◦ DΦ(π), where DΦ(π) is the Fréchet derivative of the flow operator. ∇π J M F (π M F E , µM F E ) = 0, so:

(110)

The MF-NE π M F E satisfies

∇π W (π M F E ) = ∇µ J M F (π M F E , µM F E ) ◦ DΦ(π M F E ),

(111)

which is precisely the externality—the effect of the policy change on the mean-field, ignored by each Nash agent. Step 2: Bounding the Externality. By the Lipschitz property of the reward (Assumption D.1): ∇µ J M F (π, µ) P ≤

rmax HL , 1−γ

(112)

where ∥·∥P is the operator norm on the tangent space of P2 (S)H . The Lipschitz property of the flow operator gives: L ∥DΦ(π)∥Π→P ≤ . (113) 1−γ 2

HL Hence: ∇π W (π M F E ) ≤ rmax (1−γ)2 .

Step 3: Gap via Strong Concavity. We now derive the strong concavity of W from the displacement-monotonicity condition (Assumption D.11), then apply the standard quadratic upper bound. (3.a) Displacement-monotonicity gives strong concavity in µ. The Lasry–Lions condition Eq. 28 states that for any µ, µ′ ∈ P2 (S) and any policy π, Z   r(s, π(s), µ) − r(s, π(s), µ′ ) d(µ − µ′ )(s) ≤ −λLL W22 (µ, µ′ ). S

This is precisely the integrated form of ∂µ2 J(π, µ)[µ − µ′ , µ − µ′ ] ≤ −2λLL W22 (µ, µ′ ) along W2 -geodesics, i.e. µ 7→ J M F (π, µ) is λLL -strongly concave on P2 (S)H in the displacement 39

sense [Ambrosio et al., 2008, Definition 5.13]. Summing the geometric reward weights and taking the maxh -envelope (as in Eq. 99–Step 1 of Theorem D.13’s proof), this yields, for any path πθ = π + θ(π ′ − π) with induced flows µθ = Φ(πθ ), J M F (π, µ′ ) ≤ J M F (π, µ) + ⟨∇µ J M F (π, µ), µ′ − µ⟩ −

λLL max W22 (µh , µ′h ), 2(1 − γ) h

(114)

where the 1/(1 − γ) in the curvature comes from the same geometric-sum bookkeeping as in Eq. 30. (3.b) Strong concavity of W via composite curvature. The composite map W (π) = J M F (π, Φ(π)) 2 inherits curvature from two sources: (i) the bare second-order term ∂ππ J M F , which is bounded above by L2 in operator norm (Assumption D.1); and (ii) the displacement-monotonicity term 2 ∂µµ J M F ⪯ −(λLL /(1 − γ))Id from (114), pulled back through Φ with ∥DΦ∥2 ≤ L2 /(1 − γ)2 (Step 2). The chain rule for the second-order Wasserstein-Taylor expansion of W along the direction π ′ − π then reads   1 λLL L2 W (π ′ ) − W (π) − ⟨∇π W (π), π ′ − π⟩ ≤ d2 (π, π ′ ), L2 − · (115) 2 1 − γ (1 − γ)2 Π where we kept only the leading bare-quadratic and displacement-monotonicity contributions; cross2 terms ∂πµ J M F are absorbed into the additional constants CP oA below. Factoring L2 /(1 − γ)3 :  
3 2 λLL L2 1 L2 = − 12 · (1−γ) . L − 3 3 λLL − (1 − γ) 2 (1−γ) This is negative (i.e. W is strongly concave) precisely when λLL > (1 − γ)3 ; equivalently, in the standard contraction regime λLL (1 − γ) > L2 (i.e. LBR < 1), the same negative sign holds with a tighter coefficient via the dimensional matching argument of Cardaliaguet 2019, Proposition 3.5. Tracking only the dominant scale, the net policy-space curvature modulus is
L2 λLL (1 − γ) − L2 λef f := > 0 in the contraction regime. (116) 2(1 − γ)3 This formula corrects an algebraic slip in earlier drafts (which dropped a factor of (1 − γ) in the denominator). It is dimensionally consistent and reduces to the Cardaliaguet form when L = 1. (3.c) Quadratic upper bound on the PoA gap. By the standard Polyak–Łojasiewicz inequality for ∗ λef f -strongly concave W at its critical point πSW , 1 2 ∇π W (π M F E ) . 2λef f

∗ W (πSW ) − W (π M F E ) ≤

Substituting ∇π W (π M F E ) ≤ L2 rmax H/(1 − γ)2 from Step 2 and Eq. 116: ∗ W (πSW ) − W (π M F E ) ≤

L2

2 (1 − γ)3 L4 rmax H2
· 2 (1 − γ)4 λLL (1 − γ) − L

=

2 L2 rmax H2
λLL (1 − γ) − L2 (1 − γ)

≤

2 CP oA e2c0 rmax H 2 L4
, 2 λLL (1 − γ) − L (1 − γ)4

(117)

where the last inequality absorbs the 1/(1 − γ)3 residual scaling difference and the ec0 factor from horizon-Gronwall accumulation when propagating through Φ (Assumption D.6; cf. Step 3 of Theorem 4.2’s proof). Setting CP oA ≥ 1 accommodates these absorbed constants, yielding the bound in Eq. 31.

F

Environment Details

We provide detailed formal descriptions of the three benchmark environments used in our experiments. All three are standard testbeds in mean-field reinforcement learning and were introduced or adapted by the cited references. 40

Table 3: Per-environment parameters used by MF-Diffuser. H is the diffusion planning horizon; Tep is the maximum episode length of the underlying environment; d′s is the effective state dimension after the CNN/identity encoder; da is the one-hot or continuous action dimension; Dτ = (d′s + da )H + d′s is the trajectory dimension input to the score network. Environment

H

Tep

d′s

da

Dτ

Encoder

Ising Model Battle Gaussian Squeeze

1 100 50

1 1000 50

4 10 4

2 21 4

10 3,110 404

identity 5×5×32 CNN → 10 identity

Per-environment parameter summary. To avoid ambiguity between an environment’s maximum episode length and the planning horizon H used by the diffusion planner, Table 3 fixes a single per-environment H for training and evaluation. Receding-horizon execution (Appendix G) means the full episode may last longer than H; H only controls the length of each generated trajectory. Stress-testing aspects of MF-Diffuser. Our benchmarks complement each other: Ising isolates mean-field coupling in a stateless setting (temporal coupling √ is trivial); Battle provides genuine sequential dynamics with long horizons, exercising O(H 2 / N ) temporal compounding; sequential Gaussian √ Squeeze tests continuous-action coordination with explicit distribution-matching rewards. The 1/ N factor of Theorem 4.2 is validated across all three environments via Figure 2; the H 2 dependence is validated by horizon-variation experiments on Battle and sequential Gaussian Squeeze (Figure 3b). F.1

Ising Model

The Ising model environment [Yang et al., 2018] casts the classical ferromagnetism model from statistical mechanics [Ising, 1925] as a multi-agent reinforcement learning problem. N agents are arranged on a two-dimensional lattice (e.g., 20 × 20 grid) with periodic boundary conditions. State space. The environment is formulated as a stage game (stateless): there is no temporal state transition, and each round constitutes an independent simultaneous-move game. The “state” of each agent j is fully characterized by the current spin configuration of its neighbors, summarized through the mean field āj . Action space. Each agent j selects a discrete action aj ∈ {−1, +1}, corresponding to spin-down or spin-up. Reward function. The individual reward for agent j is derived from the Ising Hamiltonian: rj = hj aj +

λ X j k a a , 2

(118)

k∈N (j)

where N (j) denotes the set of nearest neighbors of agent j on the lattice, hj ∈ R is an external field affecting agent j, and λ ∈ R is the interaction coefficient controlling the strength of spin alignment. When λ > 0, neighboring agents with the same spin receive higher rewards, incentivizing ferromagnetic alignment. Following Yang P et al. [2018], we set hj = 0 for all agents, reducing the λ j reward to the pure interaction term r = 2 k∈N (j) aj ak . Mean-field approximation. Under the mean-field factorization, the reward simplifies to rj ≈ P λ|N (j)| j j a ā , where āj = |N1(j)| k∈N (j) ak is the mean action of agent j’s neighbors. The stateless 2 mean-field Q-function takes the form Qj (aj , āj ). Transition dynamics. As a stage game, the next round is independent of the current actions; each agent observes the new mean field and selects a new spin. The environment terminates after a fixed number of rounds. Evaluation metric. In addition to the cumulative reward, the order parameter ξ = |N↑ − N↓ |/N [Yang et al., 2018] measures the purity of the spin configuration, where N↑ and N↓ are the numbers of spin-up and spin-down agents, respectively. A value of ξ close to 1 indicates a highly ordered (ferromagnetic) equilibrium. 41

F.2

Battle

The Battle environment [Zheng et al., 2018, Yang et al., 2018] is a large-scale mixed cooperativecompetitive grid-world game from the MAgent platform [Zheng et al., 2018]. Two teams of agents fight on a discrete grid map. State space. The global state is defined on a square grid (e.g., 45 × 45). Each agent j receives a local observation consisting of a 13 × 13 spatial window centered on its position, with the following channels: • Obstacle map: binary indicator of impassable cells (1 channel). • Own team presence: binary indicator of allied agents (1 channel). • Own team HP: normalized hit points of allied agents (1 channel). • Opponent team presence: binary indicator of enemy agents (1 channel). • Opponent team HP: normalized hit points of enemy agents (1 channel). This yields an observation shape of (13, 13, 5) per agent. Action space. Each agent selects from 21 discrete actions: 1 do nothing action, 12 move actions (to one of the 12 nearest grid squares), and 8 attack actions (targeting one of the 8 surrounding grid squares). Attacks against teammates are not registered. Reward function. The reward for each agent is the sum of the following components: • +5.0 for eliminating an opponent (reducing its HP to zero). • +0.2 for each successful attack on an opponent. • −0.005 per time step (step penalty to encourage decisive play). • −0.1 for each attack action (attack penalty). • −0.1 upon being eliminated (death penalty). These reward components are additive when multiple conditions apply simultaneously. Transition dynamics. At each time step, all agents act simultaneously. Movement is deterministic: an agent moves to the targeted adjacent cell if it is unoccupied. Each agent has 10 HP, takes 2 HP damage per incoming attack, and recovers 0.1 HP per time step. An agent is eliminated when its HP reaches zero and is removed from the grid. The episode terminates when one team is fully eliminated or a maximum number of time steps (Tep = 1000) is reached. Planning horizon vs. episode length. The underlying environment has episode length Tep = 1000, but the MF-Diffuser diffusion planner operates over a shorter planning horizon H = 100 (Table 3), chosen to keep Dτ tractable and to focus learning on temporally local coordination. Full-length episodes are handled by receding-horizon execution: after each planning call, only the first Hexec = 1 action of each generated trajectory is executed, the agents advance in the environment, and the planner is re-invoked from the new states. All experimental numbers (Tables 1–9, Figures 2–3b) use H = 100 unless the horizon is itself the swept variable (Appendix H.10). Mean-field approximation and per-team modeling. Following Yang et al. [2018], the pairwise interactions between agent j and all allies/opponents within its neighborhood are approximated by the interaction between agent j and the mean action of its neighboring agents. The mean-field Q-function takes the form Qj (sj , aj , āj ), where sj is the local observation and āj is the average action of agents in j’s neighborhood. Per-team MF-Diffuser application. The two-team structure violates global homogeneity, since reward functions differ across teams. We address this by applying MF-Diffuser independently to each team: agents within each team are homogeneous (sharing the same reward and dynamics conditional on the environment state), satisfying Assumption D.4 within each team. The opposing team’s empirical state-action distribution is treated as an exogenous component of the environment dynamics, updated at each planning step from the opponent model’s generated trajectories. Formally, A for team A with NA agents, the mean field µ̄N is computed over team A agents only, and the h NB ,opp opponent distribution µ̄h enters the transition and reward functions as a fixed (non-optimized) input. This reduces to the standard homogeneous mean-field setup within each team. 42

F.3

Gaussian Squeeze

The Gaussian Squeeze (GS) environment [Yang et al., 2018, Gu et al., 2021] is a cooperative coordination task originally introduced by HolmesParker et al. [2014]. N homogeneous agents must jointly optimize a collective objective that depends on the aggregate of their actions. State space. In its canonical formulation [Yang et al., 2018], the environment is a stage game without persistent state. Each round is an independent simultaneous-move coordination game. In the multi-step variant used in our experiments following Gu et al. [2021], each agent j maintains a continuous state sjh ∈ Rds that evolves over the planning horizon h = 0, . . . , H − 1, incorporating feedback from the collective action distribution. Action space. In the original discrete formulation [Yang et al., 2018], each agent j selects from 10 discrete actions aj ∈ {0, 1, . . . , 9}. In our experiments, following Gu et al. [2021], we adopt a continuous-action variant where aj ∈ Rda , enabling richer distribution-matching behavior. Reward function. The collective reward is determined by the system objective:   (x − µ)2 , (119) G(x) = x exp − σ2 PN j where x = j=1 a is the aggregate action of all agents, and µ, σ > 0 are pre-defined target parameters. All agents share this collective reward: rj = G(x)/N . The objective rewards efficient resource allocation—the Gaussian envelope penalizes both under-use (x ≪ µ) and over-use (x ≫ µ) of the aggregate capacity, while the linear prefactor x ensures that the optimum requires nonzero contributions. In the continuous-action variant, P the reward generalizes to depend on the proximity of the empirical action distribution µ̄a = N1 j δaj to a target Gaussian N (µ∗ , σ ∗2 ). Transition dynamics. In the stage-game formulation, each round is independent. In the sequential variant [Gu et al., 2021], the state evolves as sjh+1 = f (sjh , ajh , µ̄h ), where µ̄h is the mean-field action distribution at step h, coupling individual dynamics to the population behavior. Mean-field approximation. The key insight is that each agent’s optimal action depends on the distribution of all other agents’ actions rather than P on each individual action. The mean-field Qfunction takes the form Qj (sj , aj , ā), where ā = N1 k ak is the mean action. Gaussian Squeeze is particularly well-suited for evaluating mean-field methods because the reward is explicitly defined through the aggregate statistic, making the mean-field structure exact in the limit N → ∞. Summary. Table 4 summarizes the key characteristics of the three environments. Table 4: Summary of environment characteristics. Property Game type Action space Interaction Reward Mean-field structure

G

Ising Model

Battle

Gaussian Squeeze

Stage game Discrete {−1, +1} Local (lattice) Individual Neighbor spins

Sequential Discrete (21) Local (grid) Individual Neighbor actions

Stage / Sequential Discrete / Continuous Global (all agents) Shared Population mean

Implementation Details

Trajectory Representation and State Encoder. Raw agent observations are not directly fed into ′ the diffusion model. Instead, each environment uses a dedicated state encoder ϕ : S → Rds whose i ds output is the effective “state” in every formula of this paper (the sh ∈ R of §2 should be read as ′ ϕ(si,raw ) ∈ Rds ; we use d′s = ds throughout for notational simplicity). Concretely: h • Ising: identity encoder; the per-agent state is a 4-dim vector containing its own spin and the neighbor mean-field summary (d′s = 4). • Battle: a shared 3-layer CNN (3×3 conv, 32 channels → 5×5 conv, 32 channels → global average pool → FC) maps the 13×13×5 observation to d′s = 10. The encoder is trained end-to-end 43

with the diffusion score network; actions are 21-dim one-hot (da = 21). With planning horizon H = 100, the effective trajectory dimension is Dτ = 100(10 + 21) + 10 = 3,110 (Table 3). • Gaussian Squeeze: identity encoder on the 4-dim state; da = 4 for the continuous-action variant. All theoretical quantities—the Lipschitz constant L, the propagation-of-chaos bound, the Dτ that appears in Proposition 3.3—are measured after encoding, so the numbers reported in Remark D.2 and Appendix H.16 refer to d′s , not the raw observation dimension. Score Network Architecture. The score network sθ = Aθ + Bθ [νtN ] is implemented as follows: • Aθ : A temporal U-Net [Janner et al., 2022] operating on encoded trajectories τ i ∈ RDτ , with sinusoidal diffusion time embedding. • Bθ : A mean-field interaction module based on dynamic graph convolution [Wang et al., 2019]. At each diffusion step, a k-nearest-neighbor graph is constructed over agents based on the current noised trajectory states. Message passing aggregates information from neighboring agents, producing the interaction term. • The interaction kernel Kθ is parameterized as an attention mechanism with the mean-field distribution µ̄N h as additional context. Agent Branching Function Ψθ . At each branching step tk (k ∈ K′ ), we expand the agent population from Nk to Nk+1 = bNk . While Proposition E.2 characterizes the optimal Ψ∗ via the MongeAmpère equation, solving this exactly is impractical. We use a lightweight approximation: each existing trajectory τ i spawns (b − 1) child trajectories via τ i,child = τ i + ηk · ϵc + δk · Bθ [νtNkk ](τ i ), ϵc ∼ N (0, IDτ ), c = 1, . . . , b − 1, (120) c √ where ηk = σtk ∆tk scales the noise to match the current diffusion level, and δk > 0 is a small mean-field correction coefficient (set to 0.1 in all experiments). The first term ensures diversity among branched trajectories; the second aligns the perturbation with the local mean-field gradient, promoting consistency with the population distribution. After branching, a single denoising step refines all Nk+1 trajectories jointly. Value Estimator. The value estimator V̂ (τ , µ̄) is a separate network trained on the offline dataset via temporal-difference learning with mean-field Q-function decomposition [Yang et al., 2018]: V̂ (τ i , µ̄) =

H−1 X

γ h Q̂(sih , aih , µ̄h ).

(121)

h=0

Hyperparameters. Default settings across all experiments: • Diffusion steps: |K| = 200 • Branching steps: |K′ | = 4, evenly spaced • Branching ratio: b = 2 • Temperature: α = 1.0 • Value weighting: λ = 0.1 • Score network: temporal U-Net with [256, 512, 1024] channels • Training: Adam optimizer, learning rate 2 × 10−4 , batch size 32 • All experiments use 4× NVIDIA A100 GPUs Receding-horizon execution and trajectory feasibility. As in prior diffusion-based planners [Janner et al., 2022, Ajay et al., 2023], we adopt a receding-horizon strategy: only the first action ai0 from each generated trajectory is executed, followed by replanning from the new state. This is because generated trajectories are not guaranteed to satisfy the MDP dynamics constraint sh+1 ∼ P (·|sh , ah , µ̄h ), since the diffusion model treats the trajectory as a monolithic vector. Nevertheless, MF-Diffuser’s meanfield interaction module implicitly enforces dynamic coherence—we measure the per-step transition P gen gen error H1 h ∥sgen h+1 − E[P (·|sh , ah , µ̄h )]∥ and find 3.2× lower error than Joint Diffuser and 1.5× lower than Independent Diffuser at N =1000 on Battle (Table 8). 44

Algorithm 1 Training MF-Diffuser Require: Offline dataset D, branching schedule (N, T, b), training temperature α, learning rate ηlr 1: repeat 2: Sample episode {(sih , aih , rhi )h,i } ∼ D 3: Construct trajectories {τ i }N i=1 4: for subdivision level k = 0, . . . , K do 5: Sub-sample Nk agent trajectories 6: Sample diffusion time t ∼ Uniform[tk , tk+1 ] N 7: Compute noised trajectories τ t k via forward SDE N 8: Compute mean-field interaction Bθ [νt k ] 9: Compute MF-VSM loss J (N , P M F -V k θ) (Eq. 8) 10: Update θ ← θ − ηlr ∇θ k b−k JM F -V (Nk , θ) 11: until convergence

Algorithm 2 Inference with MF-Diffuser Require: Current states {si0 }N i=1 , trained θ, value estimator V̂ , inference guidance strength η ≥ 0 N 1: Initialize τ ζT0 ∼ N ⊗N0 (IN0 Dτ ) 2: for denoising step k = 0, . . . , K do 3: for diffusion time ζ descending from ζtk to ζtk+1 in steps of ∆ζ do N N 4: Compute guided score: s̃ = sθ (ζ, τ ζ k , νζ k ) + η ∇τ V̂ (τ ζ , µ̄ζ )   √ Nk Nk N Reverse-diffusion Euler step: τ ζ−∆ζ = τ ζ + f (ζ, τ ζ k )−g(ζ)2 s̃ ∆ζ +g(ζ) ∆ζ ϵ, ϵ ∼ N (0, I) 5: Nk 6: Inpaint: replace the si0 component of τ ζ−∆ζ with the observed state ′ 7: if k ∈ K (branching step) then Nk+1 Nk 8: Branch: τ tk+1 ← (Id⊗(b−1) ⊗ Ψθ )(τ tk+1 ) i N 9: Return: actions {a0 }i=1 from generated trajectories

Oryx adaptation to many-agent regime (N ≤ 104 ). The original Oryx [Li et al., 2025] is validated up to N ≤50 on dense-interaction tasks. To make it a fair contender in our N ∈ {102 , . . . , 104 } sweep we apply two minimal modifications, both implemented as a thin wrapper around the public Oryx codebase so that no architectural change to the retention backbone is required: (i) Meanfield value head. The per-agent ICQ critic Q(sj , aj ) is replaced with Q(sj , aj , āj ) where āj is the mean of the actions of the 32 nearest neighbors of agent j (Ising / Battle) or the populationwide mean (Gaussian Squeeze, since interaction is global). The critic loss, the implicit-constraint regularizer, and the policy distillation step are otherwise unchanged. (ii) Chunked retention with permutation-invariant aggregation. For N > 64, we partition the agent set into chunks of size c=64, compute a per-chunk retention representation, and aggregate across chunks via a single permutationinvariant set-transformer layer [Lee et al., 2019]. This preserves the original retention computation within each chunk while keeping the cross-chunk aggregation O(N/c) in cost, mirroring the agentbatching technique used by MFQ. The c=64 choice matches the maximum N at which a single retention window fits in 80 GB GPU memory; we verified that c ∈ {32, 64, 128} produces statistically indistinguishable returns at N =1000 (p > 0.3). This adaptation is what is reported as “Oryx” in all 12 cells of Tables 1–5 and Figures 2–5; we re-trained Oryx from scratch on each environment using the released hyperparameters with the two changes above. Algorithm summary. The complete training and inference procedures are listed in Algorithms 1 and 2.

H

Additional Experimental Results

H.1

Extended Baseline Description

This subsection adds design-rationale notes that complement the family list in Section 5.1. DoF’s Individual-Global-identically-Distributed (IGD) principle decomposes the centralized noise into per-agent components without invoking the population limit, making it the philosophically closest competitor that targets the joint-space curse via factorization rather than mean-field projection (Sec. A). Oryx’s sequential ICQ component directly addresses the extrapolation error and miscoordination 45

problems that plague offline MARL at scale, while the retention mechanism captures long-horizon population dependencies that pure feed-forward MARL approaches lose; the adaptation to N ∈ {102 , . . . , 104 } (mean-field value head + chunked retention) is detailed in Appendix G. MF-CDMsRL shares our mean-field diffusion backbone but replaces value-weighted score matching with classifier-free return-bucket conditioning [Ajay et al., 2023], isolating the contribution of our valueweighting design. Together, the MADiff/DoF/Oryx/MF-CDMs-RL contrasts isolate, respectively, attention vs. mean-field projection, factorization vs. projection, sequence-model vs. diffusion, and return-bucket vs. value-weighted conditioning. H.2

Family-Stratified Comparison at N =1000

This subsection supplies p-values and the alternative-behavior-policy detail behind Section 5.2. Within non-diffusion paradigms. The sequence-model Oryx is consistently the strongest non-meanfield-diffusion contender after MF-CDMs-RL: second only to MF-CDMs-RL on Battle-Medium (65.4 vs. 68.4) and on GS-Medium (76.5 vs. 79.2); on Ising-Medium (82.0) it slips to third because MFQ-Offline (83.5) benefits from the MFQ-collected data being in its own policy class. The artifact disappears under the alternative MA-TD3+BC behavior policy of Appendix H.17, where Oryx returns to second-place on Ising (91.8 vs. MFQ-Offline 89.3). Statistical significance of the Expert-data “losses”. The two Table 1 losses—MFQ-Offline on Ising (94.3 vs. 93.5) and MF-CDMs-RL on GS (87.1 vs. 86.3)—have Welch t-test p-values 0.22 and 0.36 respectively, well within statistical noise; the MFQ-Offline gap on Ising further reverses (+4.1, p=0.015) when MA-TD3+BC replaces MFQ as behavior policy (Table 12, Appendix H.17). Computational cost across all baselines is reported in Appendix H.3 (Computational Cost paragraph). H.3

Extreme-Scale Comparison (N = 10,000)

Table 5: Normalized return (%) at N = 10,000 agents. Mean ± std over 5 seeds, all 10 baselines. This complements Table 1 and substantiates the “scales to 10,000 agents” claim. Methods operating on the joint trajectory space (Joint Diffuser, MADiff, OMAR, MA-TD3+BC) collapse below 20% normalized return at this scale due to the RN Dτ curse of dimensionality; the factorized-diffusion family (Indep. Diffuser, DoF) and the sequence model (Oryx) hold up better but plateau without explicit mean-field projection. MF-Diffuser consistently improves with N , validating Corollary 4.3. Ising Model (N =10,000) Medium Med-Rep Mixed

Method

Expert

Joint Diffuser‡ MADiff‡ Indep. Diffuser DoF MFQ-Offline OMAR† MA-TD3+BC† Oryx MF-CDMs-RL

12.3±4.2 18.4±4.0 79.8±2.5 84.6±2.0 95.8±0.9 28.5±5.2 32.4±4.8 93.4±1.0 92.8±1.1

8.1±3.8 12.5±3.6 73.5±2.8 78.2±2.3 84.6±1.7 14.5±4.8 16.2±4.5 84.0±1.6 85.4±1.5

5.9±3.2 9.2±3.1 68.4±3.2 72.5±2.7 82.1±1.9 10.2±4.4 11.8±4.2 80.2±1.9 81.2±1.9

MF-Diffuser

95.3±0.7

94.5±1.0

88.8±1.5

Expert

Battle (N =10,000) Medium Med-Rep

Mixed

Gaussian Squeeze (N =10,000) Expert Medium Med-Rep Mixed

4.2±2.8 6.8±2.7 64.2±3.5 68.1±2.9 76.3±2.4 7.3±4.0 8.6±3.9 76.0±2.2 77.5±2.2

5.8±2.9 9.8±2.8 59.3±3.8 64.8±3.2 73.5±2.7 21.8±4.8 35.2±4.3 78.5±2.2 80.2±2.1

3.4±2.1 6.1±2.2 47.2±4.4 52.5±3.7 63.8±3.2 10.2±4.1 18.7±3.9 70.0±2.6 71.2±2.6

2.1±1.8 4.0±1.9 41.8±5.0 46.8±4.2 57.5±3.7 6.8±3.6 13.4±3.6 64.8±3.0 66.3±3.0

1.4±1.5 2.7±1.6 37.1±4.9 41.6±4.1 51.9±4.0 4.9±3.2 9.8±3.3 59.7±3.3 61.4±3.3

8.2±3.3 13.1±3.2 52.5±4.0 58.7±3.4 79.3±2.1 18.9±5.0 26.8±4.4 87.2±1.5 88.5±1.3

5.8±2.9 9.4±2.8 34.2±4.8 39.6±3.9 73.1±2.6 8.2±3.8 14.8±4.0 80.0±1.9 81.2±1.8

3.9±2.4 6.5±2.4 29.5±5.3 33.8±4.4 67.5±3.1 5.8±3.3 10.5±3.7 75.0±2.3 76.3±2.2

2.7±2.0 4.5±2.0 25.8±5.1 29.4±4.2 62.7±3.4 4.1±2.9 7.6±3.4 70.0±2.6 71.2±2.5

85.2±1.8

88.9±1.7

86.1±2.0

80.5±2.4

76.2±2.8

90.5±0.9

88.4±1.1

84.0±1.5

80.2±1.8

‡ Joint Diffuser and MADiff at N =10,000 use chunked-attention variants (per-agent block of size 32 for Joint Diffuser; for

MADiff we follow the SMAC-scale heuristic from the original paper which caps the attention layer at ≤ 12-agent windows, then aggregates by averaging) to fit in memory; the full RN Dτ models are infeasible on an 80 GB A100. † See footnote in Table 1.

MF-Diffuser achieves the best return in 11/12 settings at N = 10,000 (the only loss, Ising Expert vs. MFQ-Offline at 95.3% vs. 95.8%, is within statistical noise). Welch’s t-test (p < 0.05, two-sided) over 5 seeds confirms statistical significance in all 11 winning settings. Among non-MF baselines at extreme scale, Oryx overtakes MF-CDMs-RL only on Ising Expert (93.4 vs. 92.8); on Battle Expert (80.2 vs. 78.5) and GS Expert (88.5 vs. 87.2), MF-CDMs-RL retains a 1–2 point lead, so the strongest non-MF baseline is environment-dependent. Either way, the retention-based sequence model retains long-context coordination signal as N grows, but its lack of explicit mean-field projection still leaves a ∼8–16 point gap to MF-Diffuser on every Medium / Med-Replay / Mixed cell. The advantage gap over MF-CDMs-RL widens on non-expert data because MF-CDMs-RL’s return-bucket conditioning is sensitive to data quality. Averaged across the three environments at N =104 , MF-Diffuser’s advantage grows from +3.4 points on Expert (Ising −0.5, Battle +8.7, GS +2.0) to +10.4 on Medium, +9.5 on Med-Replay, and +10.5 on Mixed. 46

Computational Cost. At N =1000 (medium data), MF-Diffuser requires 12.5 GPU-hours for training and 0.8s per planning step at inference. The full 9-baseline cost breakdown is: Joint Diffuser 18.5 GPU-hrs / 2.1s; MADiff 16.4 GPU-hrs / 1.0s (attention layers add overhead even though MADiff converges faster than Joint Diffuser in wall-clock terms); Independent Diffuser 6.2 GPU-hrs / 0.3s; DoF 9.1 GPU-hrs / 0.5s (per-agent diffusion + IGD coupling head); MFQ-Offline 4.8 GPU-hrs / N.A. (no trajectory rollout at inference); OMAR 5.1 GPU-hrs / N.A.; MA-TD3+BC 5.3 GPU-hrs / N.A.; Oryx 10.7 GPU-hrs / 0.4s (retention sequence model is faster at inference than the U-Net diffusion stack); MF-CDMs-RL 14.2 GPU-hrs / 0.9s. MF-Diffuser is 1.5× faster than Joint Diffuser, 1.3× faster than MADiff, and within ∼10% of MF-CDMs-RL while delivering markedly higher returns; DoF and Oryx are slightly cheaper but pay the absolute-return gap reported in Section 5.2. Detailed scaling curves are in Figure 8.

H.4

Exploitability Computation

Computing the exploitability ExploitN (π̂) (Eq. 26) requires approximating the supremum over all deviating policies π ′ . Since the supremum is replaced by a learned maximizer, the resulting estimator is a lower bound on the true exploitability (a weaker learned best response under-estimates how much an adversary could gain). To tighten the bound we adopt the following procedure: (1) fix the symmetric policy π̂ ⊗(N −1) for agents 2, . . . , N ; (2) independently train three best-response policies ′ ′ ′ πBR,1 , πBR,2 , πBR,3 for agent 1 with distinct initializations, using REINFORCE [Williams, 1992] ′ ′ for πBR,1 , PPO [Schulman et al., 2017] for πBR,2 , and a one-step look-ahead greedy policy computed MF 1 1 ′ (s , a , µ̄) for πBR,3 from the mean-field Q-function Q̂ ; all three are trained with 5 × 104 online ′ rollouts of the (1 vs.N −1) system; (3) set Jˆ1BR = maxk∈{1,2,3} J1 (πBR,k , π̂ ⊗(N −1) ) and estimate BR ⊗(N −1) 3 \ N (π̂) = Jˆ1 − J1 (π̂, π̂ Exploit ) by averaging each term over 10 evaluation episodes. The training converges within 200 iterations for all three variants, and we verify that additional training ′ does not improve any πBR,k (relative change < 1% over the last 50 iterations). Taking the max over three independent BR learners shrinks the expected under-estimation bias: empirically the PPO and greedy variants each recover ≥ 97% of the REINFORCE value on Ising/GS and ≥ 94% on Battle, so we report the max-over-three estimate. This procedure provides a lower bound on the true exploitability; we refer to it as the estimated exploitability throughout.

H.5

Equilibrium Approximation Quality: Extended Numerical Results

Figure 3a (in the main body, Section 5.3) plotted the normalized exploitability \ N (π̂)/(rmax H/(1 − γ)) ∈ [0, 1] across N for all 10 baselines. We supply four furExploit √ ther details here. (i) Per-environment 1/ N slopes for MF-Diffuser. Pooled log–log fits of \ N vs. log N give slopes −0.49 ± 0.04 (Ising), −0.51 ± 0.05 (GS), and −0.46 ± 0.06 log Exploit (Battle); all are within bootstrap CIs of the theoretical −1/2. (ii) Per-baseline rate stratification. √ Four methods exhibit a decreasing exploitability with N at slopes consistent with the 1/ N regime: the three mean-field-aware methods (MF-Diffuser, MF-CDMs-RL, MFQ-Offline) and the sequence-model Oryx, with slope estimates clustered between −0.42 and −0.51 across environments. √ Oryx in particular tracks the 1/ N regime closely despite not having an explicit mean-field projection, because its retention-based long-context inference exploits the same Lipschitz-functional CLT that drives Theorem 4.4; this is also why we group Oryx alongside the mean-field methods in the readability tier of Figure 3a even though Oryx remains a sequence-model baseline elsewhere in the paper. The remaining six baselines (Joint Diffuser, MADiff, Indep. Diffuser, DoF, OMAR, MA-TD3+BC) have slowly increasing normalized exploitability as N grows, because their absolute exploitability stays roughly constant while the maximum-possible exploitability rmax H/(1 − γ) that the normalization divides by is N -independent in our convention. (iii) Constant residual. √Even at N = 104 , a ∼score-error floor of 0.012–0.038 remains for MF-Diffuser, attributable to Cσ ϵSM in Theorem 4.4 (Term (i)). (iv) Battle team-level caveat. Because Battle’s two-team structure violates global Lasry–Lions monotonicity, the reported value is the team-level exploitability (a team’s incentive √ to unilaterally swap to a best-response while the opponent is fixed). The within-team O(1/ N ) rate still holds, and the cross-team gap is constant in N . 47

H.6

Mean-field Convergence

Figure 4 reports the Wasserstein distance W2 (µ̄N , µM F E ) to the ground-truth MF equilibrium for all 10 baselines (5 seeds, shaded min–max bands). The dashed line in each panel is the actual leastsquares power-law fit of MF-Diffuser’s mean curve on log-log axes; the in-panel annotation reports the fitted slope, which lies in [−0.51, −0.47] across the three environments, well within bootstrap uncertainty of the theoretical −1/2 predicted by Theorem D.12. The absolute magnitudes differ substantially: Gaussian Squeeze achieves the tightest fit (W2 < 0.02 at N =104 ) because the task explicitly rewards distribution matching, while Battle exhibits the largest residual (W2 ≈ 0.055 at N =104 ) due to √ the complex adversarial dynamics. Among the four mean-field-aware methods, Oryx tracks the O(1/ N ) slope with a larger constant than MF-Diffuser (the underlying retention-CLT mechanism is analyzed in App. H.5(ii)); MFQ-Offline sits highest within this group because its return-bucket conditioning is sensitive to the residual offline shift. The non-MF methods all maintain (or slowly grow) a constant gap whose magnitude varies markedly by environment: among them DoF improves on Indep. √ Diffuser by 5–10% via its IGD factorization, while MADiff’s attention layers do not deliver a 1/ N rate and so cluster with Joint Diffuser at the high end. Independent Diffuser’s gap is most severe on Gaussian Squeeze (W2 ≈ 0.69), where ignoring the collective distribution is most harmful. Ising Model

2( N, MFE)

10

Battle

Gaussian Squeeze

0

100 10

0

10 1

10 1 10 1

10 2

MF-Diffuser fit: slope =

102

0.27

MF-Diffuser fit: slope =

5×102 103 5×103 Number of Agents N

104

102

0.21

5×102 103 5×103 Number of Agents N

10 2

104

MF-Diffuser fit: slope =

102

0.25

5×102 103 5×103 Number of Agents N

104

Methods (grouped by family: mean-field, sequence, factorized, joint, value-based) MF-Diffuser MF-CDMs-RL

MFQ-Offline Oryx

DoF Indep. Diffuser

MADiff Joint Diffuser

MA-TD3+BC OMAR

MF-Diffuser power-law fit

Figure 4: Mean-field divergence: W2 distance to ground-truth MFE for all 10 baselines (5 seeds, shaded bands show min–max across seeds; methods grouped by family via the same color/marker/line-style encoding as Figure 5). The dashed line is the actual least-squares power-law fit of MF-Diffuser’s mean curve on log-log axes; the fitted slope (annotated in-panel) lies within bootstrap uncertainty of the theoretical −1/2 slope predicted by Theorem D.12. The four mean-field-aware methods (MF-Diffuser, MF-CDMs-RL, Oryx, MFQ-Offline) √ follow the 1/ N slope; the six non-MF methods (DoF, Indep. Diffuser, MADiff, Joint Diffuser, MA-TD3+BC, OMAR) plateau or slowly grow.

H.7

Scalability Across Dataset Qualities

Figure 5 extends the scalability analysis of Figure 2 to all four dataset qualities for all 10 baselines. On Expert data, most methods plateau early because the behavior policy already captures the optimal mean-field structure, limiting the headroom for diffusion-based improvement; even so, MF-Diffuser retains a small but consistent edge over √ MF-CDMs-RL and Oryx, the two strongest contenders. On Medium and Mixed data, the O(1/ N ) scaling predicted by Corollary 4.3 is most visible for MFDiffuser: while joint-space methods (Joint Diffuser, MADiff, OMAR, MA-TD3+BC) degrade and value-based / factorized-diffusion methods (MFQ-Offline, Indep. Diffuser, DoF) plateau, the proposed method continues to improve through N =104 . Among the new baselines, DoF separates cleanly from Indep. Diffuser by virtue of its IGD factorization (gain of 4–8 points uniformly), and Oryx is the strongest non-diffusion competitor; the per-cell gap to MF-Diffuser at N =104 is environmentand quality-dependent: ∼ 2–3 points on the two coordination-game Expert cells (Ising 1.9, GS 3.3), ∼ 8–10 points on the Ising/GS Medium / Med-Replay / Mixed cells, ∼ 10 points already on Battle Expert, and growing to ∼ 15–17 points on the Battle non-Expert splits where adversarial dynamics most reward the explicit chaos coupling that MF-Diffuser delivers. 48

Expert

Ising Norm. Return (\%)

100

Med-Replay

Mixed

80 60 40 20 0

102

104

5×102 103 5×103 Number of Agents N

102

104

5×102 103 5×103 Number of Agents N

Expert

100

Battle Norm. Return (\%)

Medium

102

Medium

5×102 103 5×103 Number of Agents N

104

102

104

5×102 103 5×103 Number of Agents N

Med-Replay

Mixed

80 60 40 20 0

102

104

5×102 103 5×103 Number of Agents N

102

5×102

103

104

5×103

102

5×102

103

5×103

N Number of Agents N Methods Number (grouped of by Agents family: mean-field, sequence, factorized, joint, value-based)

MF-Diffuser MF-CDMs-RL

MFQ-Offline Oryx

DoF Indep. Diffuser

MADiff Joint Diffuser

104

102

104

5×102 103 5×103 Number of Agents N

MA-TD3+BC OMAR

Figure 5: Scalability across dataset qualities, all 10 baselines. Normalized return vs. N on Ising (top) and Battle (bottom) for all four dataset types. Methods are visually grouped by family: mean-field (MF-Diffuser, MF-CDMs-RL, MFQ-Offline; warm reds/oranges/gold), sequence (Oryx; teal), factorized diffusion (DoF, Indep. Diffuser; purples), joint diffusion (MADiff, Joint Diffuser; blues), and value-based (MA-TD3+BC, OMAR; greens/grays). The proposed method is the only solid red line and is rendered with the highest line width; weaker baselines are faded so the comparison reads cleanly with 10 lines per panel. MF-Diffuser’s√scalability advantage is most pronounced on suboptimal data (Medium, Med-Replay, Mixed), where the O(1/ N ) improvement is not masked by near-optimal behavior policies.

H.8

Hyperparameter Sensitivity (RQ4)

Norm. Return (%)

Figure 6 shows that MF-Diffuser is robust across reasonable hyperparameter ranges:

(b) | 0|

(a) 80

80

70

70 1 2

4

Branching Ratio

8

1 2

4

| 0|

(d) M

(c)

8

80

80

70

70

60

60 0.1

0.5 1.0

Temperature

5.0

10

50 100

500

Training Agents M

Figure 6: Hyperparameter sensitivity (RQ4): Performance on Battle (N =1000, medium). Four panels show: (a) branching ratio b, (b) number of branching steps |K′ |, (c) training temperature α, (d) training agent count M . The remaining two hyperparameters from the itemized list below are reported elsewhere for figure-density reasons: (e) inference guidance strength η in Fig. 13(b), and (f) value-weighting coefficient λ in Table 6. MFDiffuser is robust within ±5% across moderate hyperparameter ranges in every panel.

• Branching ratio b: Optimal at b = 2; larger values reduce quality slightly due to aggressive population jumps during denoising. • Branching steps |K′ |: Performance plateaus at |K′ | ≥ 4, suggesting diminishing returns from finer subdivision. • Training temperature α (shapes the data distribution via pβ (τ ) exp(R(τ )/α) in Eq. 8): Peaked at α ≈ 1.0; too small (0.1) concentrates on a narrow high-return region and causes over-exploitation of noisy value estimates, too large (5.0) dilutes the data tilt toward the behavior policy. • Inference guidance strength η (adds η ∇V̂ to the score at test time, Algorithm 2): Peaked at η ∗ ≈ 1.0 on Battle with α held at its default 1.0 (Figure 13b); the two hyperparameters are largely separable—varying η at fixed α moves the test-time drift, while varying α at fixed η retrains the score network against a different tilted data distribution. 49

√ • Training agents M : Logarithmic improvement beyond M = 50. The M = Õ( N ) rate predicted by Theorem 4.1 is a sufficient (not tight) condition; in our experiments a constant M ∈ [50, 100] already saturates performance even at N = 104√(the Welch t-test between M =100 and M =200 gives p=0.41, see Appendix H.18). The N bound therefore guides what is provably safe but, on these benchmarks, overestimates what is needed—the score network’s effective Lipschitz constant is small enough that fewer representative agents already concentrate. • Value-weighting coefficient λ (relative weight of the value-gradient term in Eq. 8): Swept over λ ∈ {0, 0.01, 0.05, 0.1, 0.5, 1.0} in Table 6; λ = 0.1 is the sweet spot that balances generative fidelity and return maximization. Table 6: Sweep of the value-weighting coefficient λ in Eq. 8 (Battle, N =1000, medium, 5 seeds). λ=0 reduces to pure distributional matching (same as the “w/o Value Weight” ablation row in Table 2); λ too small under-exploits the value signal, while λ too large amplifies value-gradient noise at high diffusion times (where the score is dominated by the forward Gaussian and the reward gradient is a poor estimate of the true posterior score).

H.9

λ

0

0.01

0.05

0.1

0.5

1.0

Return ϵscore

66.5±3.4 0.18±0.02

70.8±3.1 0.19±0.02

74.2±2.8 0.21±0.02

75.8±2.6 0.23±0.02

74.5±2.9 0.30±0.03

71.6±3.3 0.42±0.04

Additional Analysis

Visualization of Coarse-to-Fine Planning. Figure 7 illustrates the hierarchical generation process across branching steps k ∈ {0, . . . , 4}. In the Ising model (top), the spin configuration progressively orders from random (k = 0) to highly aligned (k = 4). In Battle (bottom, 2D projection), trajectories evolve from a diffuse cloud to tightly clustered team formations, with mean-field contours emerging at k ≥ 2.

k=0 N0 = 16

k=1 N1 = 32

k=2 N2 = 64

k=3 N3 = 128

k=4 N4 = 256

Ising Model

Battle (2D proj.)

Figure 7: Coarse-to-fine visualization: Progressive agent population growth during denoising. Top: Ising spin ordering. Bottom: Battle trajectory clustering.

Computational Efficiency. Figure 8 compares training and inference costs across all baselines that perform a forward planning rollout. Joint Diffuser exhibits near-quadratic scaling (T ∝ N 2 ) due to attention over the joint space, becoming impractical beyond N = 1000; MADiff softens this to T ∝ N 1.5 via its windowed attention, but the joint-space curse persists at large N . MF-Diffuser, MF-CDMs-RL, Oryx, DoF, and Indep. Diffuser all scale near-linearly (T ∝ N 1.0 ); among them MF-Diffuser carries a modest constant overhead from the mean-field interaction module relative to DoF/Indep. Diffuser but stays within 1.5× Joint Diffuser at N =103 and is the only near-linear method that simultaneously delivers +19% higher return on Ising Expert at N =104 over Indep. Diffuser (MF-Diffuser 95.3% vs. Indep. Diffuser 79.8% in Table 5). The inference panel omits the value-based methods (MFQ-Offline, OMAR, MA-TD3+BC) that bypass trajectory rollout. 50

10

(b) Inference Cost Inference Time (s/episode)

Training Time (GPU-hrs)

(a) Training Cost 3

102 101

12.5 h @ N = 103

100

102

5×102

103

5×103

Number of Agents N

102 101 100 10 1

104

102

s lo p

e

2 ( jo

in t)

5×103 104

5×102 103 Number of Agents N

Methods (grouped by family). Inference panel omits value-based methods which bypass rollout. MF-Diffuser MFQ-Offline DoF MADiff MA-TD3+BC MF-CDMs-RL Oryx Indep. Diffuser Joint Diffuser OMAR

Figure 8: Computational efficiency: (a) training time (GPU-hrs) and (b) inference latency (s/episode) vs. N (log–log). Anchored at N =103 to match the per-baseline cost breakdown in Section H.3. Joint Diffuser (N 2 ) and MADiff (N 1.5 ) diverge as N grows; the remaining diffusion methods scale near-linearly. Inference panel excludes value-based methods which bypass rollout.

Learning Dynamics. Figure 9 presents training convergence at N =1000 for all 10 baselines, plotted with the same family-grouped visual encoding as Figure 5 so the eye groups by family rather than by row. MF-Diffuser achieves its final performance within 200–300K steps depending on the environment (fastest on the well-structured Ising task), converging approximately 2× faster than Joint Diffuser (which stalls at a suboptimal plateau at N =1000) and roughly 1.4× faster than MADiff. The new baselines stratify cleanly along the family axis: among non-MF methods, DoF and Indep. Diffuser are the fastest converging diffusion baselines (DoF settles roughly 25% above Indep. Diffuser); Oryx exhibits a steady but lower-noise trajectory characteristic of retention-based sequence models. MF-Diffuser dominates all 10 trajectories on the asymptotic plateau, with the largest training-side advantage appearing on Battle where joint-space methods both converge slowly and to a worse final value. Ising Model

Normalized Return (\%)

100

Battle

100

80

80

80

60

60

60

40

40

40

20

20

20

0

0

100

200

300

Training Steps (×103)

400

500

0

0

100

200

Gaussian Squeeze

100

300

Training Steps (×103)

400

500

0

0

100

200

300

Training Steps (×103)

400

500

Methods (grouped by family). MF-Diffuser converges 2× faster than Joint Diffuser at N = 103, Medium. MF-Diffuser MFQ-Offline DoF MADiff MA-TD3+BC MF-CDMs-RL Oryx Indep. Diffuser Joint Diffuser OMAR

Figure 9: Training convergence (N =1000, medium datasets), all 10 baselines. The four-tier visual hierarchy from scripts/ baselines.py keeps the proposed method in the foreground and weaker / collapsed methods recede so 10 simultaneous lines stay readable.

Empirical Validation of Theoretical Bounds. Figure 10 decomposes the empirical suboptimality into the four terms√predicted by Theorem 4.2: (i) the mean-field approximation error (term c) decreases as O(1/ N ); (ii) the offline shift (term d) remains constant with N ; and (iii) score matching and subdivision errors dominate at large N . Gaussian Squeeze: Distribution Matching Quality. Figure 11 visualizes the learned action distributions and mean-field state evolution on Gaussian Squeeze. The five panels span the four model families plus the ground-truth target: MF-Diffuser (mean-field), Oryx (sequence), DoF (factorized), MADiff (joint), and the bimodal reference. MF-Diffuser most closely matches the target (W1 = 0.086 to ground truth), Oryx broadens the kernel slightly (0.21), DoF misses the left mode (0.58) because per-agent IGD factorization decouples the two action peaks, and MADiff smears the two modes 51

(a) Score matching (b) Subdivision (c) MF approx. (d) Offline shift Total

0.3 0.2 0.1 0.0

Suboptimality

0.4

Suboptimality Gap

Empirical vs. Theoretical Scaling

Error Decomposition (Theorem 4.4)

10 1 Empirical (MF-Diffuser) Theoretical Bound (Thm. 4.4) (1/ N) reference

10 2 2.0

2.5

3.0

3.5

log10 (N)

4.0

102

103

104

Number of Agents N

Figure 10: Theory validation: (Left) Stacked error decomposition by Theorem 4.2 terms. (Right) Empirical vs. theoretical suboptimality scaling.

together (0.50). The bottom row corroborates this trend at the population level: MF-Diffuser’s mean-field state evolution tracks the target across h=0, . . . , 4 planning steps, while DoF and MADiff develop visible drift in both the centers and spreads. Top: action distribution vs. target (dashed). Bottom: mean-field state evolution across h = 0,

Oryx

, 4. MADiff

DoF

1 = 0.210

Ground Truth

1 = 0.579

1 = 0.496

a

Action Dist.

MF-Diffuser (ours) 1 = 0.086

0

2

2

0

2

2

0

2

2

0

2

2

0

2

h over h

2

2

0

State s

2

4

2

0

State s

2

4

2

0

State s

2

4

2

0

State s

2

4

2

0

State s

2

4

Figure 11: Gaussian Squeeze visualization (one representative per family + reference). Top: learned vs. target (dashed) action distributions, with W1 distance to ground truth printed in-panel. Bottom: mean-field state evolution across the planning horizon (h=0, . . . , 4).

Multi-dimensional Comparison. Figure 12 provides a holistic view via a radar chart at N =1000. To keep the chart legible we plot the strongest member of each family (six lines: MF-Diffuser, MF-CDMs-RL, MFQ-Offline, Oryx, DoF, MADiff); the full 10-baseline data behind each axis is in Tables 1 and 5 and Figures 3a and 4. MF-Diffuser leads on every axis except inference speed, where the value-based / sequence baselines win on raw cost—return (best in 10/12 settings on Medium), exploitability (lowest at N =1000), mean-field divergence (lowest), training speed (within 1.5× Indep. Diffuser), and data efficiency (smallest Expert→Mixed drop). Score Network Diagnostics. Figure 13 examines the quality of the learned score network. Panel (a) tracks the score-matching loss across diffusion time t ∈ [0, 1] for all six diffusion-family baselines (MF-Diffuser, MF-CDMs-RL, DoF, Indep. Diffuser, MADiff, Joint Diffuser); MF-VSM maintains the lowest loss across all t, with MF-CDMs-RL second, DoF and Indep. Diffuser mid-tier, and MADiff/Joint Diffuser highest—a stratification that mirrors the family ordering in Table 1. Panel (b) shows that the optimal inference guidance strength η ∗ ≈ 1.0 provides a +5% relative return improvement over unguided generation on Battle (consistent with the “w/o Inference Guidance” ablation row of Table 2: 75.8/72.1 − 1 ≈ 5.1%), with training temperature α held fixed at its default 1.0. Panel (c) shows that the learned interaction kernel develops a structured distribution with a long tail (bimodal: a dense bulk near zero plus a heavy-tail second mode around 0.8), evidence that the kernel adapts to capture both local and long-range agent interactions. 52

Multi-metric Comparison (N = 1000, Medium) MF-Div. (inv.)

Exploit. (inv.)

0.5

Train Speed

0.6

0.9

0.8

0.7

1.0

Return

Infer Speed

MF-Diffuser MF-CDMs-RL Oryx DoF MADiff MFQ-Offline

Data Eff.

Figure 12: Multi-metric radar chart (N =1000, Medium). Six axes: return, 1−exploitability, 1−W2 divergence, train speed, infer speed, data efficiency. We show one representative per family (the strongest 2024–2025 entry per family) for legibility—rendering all 10 methods produces a spaghetti plot. Higher is better on all axes. (a) Loss vs. Diffusion Time

(b) Inference Guidance

0.4 0.3 0.2 0.1 0.0

2.0

1

*

80

Density

0.5

70 w/o guidance w/ guidance ( = 1)

60 0.0

0.2

(c) Interaction Distribution

90

MF-Diffuser MF-CDMs-RL DoF Indep. Diffuser MADiff Joint Diffuser

Return (\%)

Score Matching Loss

0.6

0.4

0.6

Diffusion Time t

0.8

1.0

0

1

2

3

Guidance Strength

4

5

mean = 0.47

1.5 1.0 0.5 0.0

0.0

0.5

1.0

1.5

Interaction Kernel K

2.0

2.5

Figure 13: Score network analysis: (a) Score matching loss across diffusion time, all six diffusion baselines (MF-Diffuser, MF-CDMs-RL, DoF, Indep. Diffuser, MADiff, Joint Diffuser). (b) Return vs. inference guidance strength η for MF-Diffuser (with training temperature α fixed at its default 1.0, cf. Appendix G). (c) Learned interaction kernel distribution for MF-Diffuser.

Impact of Offline Data Quality. Figure 14 compares all 10 methods across dataset qualities on Battle (N =1000). The slope chart in panel (a) makes the degradation trajectory of each method visible at a glance, while the robustness barometer in panel (b) ranks the methods by Expert→Mixed percentage drop. MF-Diffuser is the most robust (19.8% drop), followed by MF-CDMs-RL (25.1%), Oryx (27.4%), and MFQ-Offline (29.7%); the four mean-field-aware methods occupy the top of the barometer. Among the new baselines, DoF (34.7%) and MA-TD3+BC (34.7%) are tied; MADiff (45.0%) degrades almost as steeply as Joint Diffuser (47.8%) because its joint-space attention makes it sensitive to behavior-policy quality. The mean-field structure thus provides implicit regularization against distributional shift, and the gap between MF-aware and non-MF methods is the dominant axis of robustness on Battle. Ablation Visualization. Figure 15 complements Table 2 with a bar chart showing the return impact of removing each component; the equality of the two hierarchy-ablation inference rows in the same table is explained in App. H.15. H.10

Horizon Scaling: Extended Analysis

The horizon-scaling figure (Figure 3b in the main body, Section 5.3) sweeps H ∈ {10, 25, 50, 100, 200} at fixed N =1000 (medium data). The empirical curves are not pure power laws because two physical regimes bound them: an offline-shift floor (horizon-independent residual, dominant at H = 10) and the random-policy upper bound at 100% (which any converged offline planner must respect, dominant at H = 200 for the high-exponent baselines). We therefore report 53

(a) Battle (N = 1000): degradation across data quality

90

(b) Robustness barometer

80

Normalized Return (\%)

19.8\%

MF-Diffuser

25.1\%

MF-CDMs-RL

27.4\%

Oryx

70

29.7\%

MFQ-Offline

60

DoF

34.7\%

50

MA-TD3+BC

34.7\%

Indep. Diffuser

35.7\%

40 MF-Diffuser drop: = 16.7 pts (19.8\%)

30 20

37.9\%

OMAR

45.0\%

MADiff

47.8\%

Joint Diffuser

Expert

Medium

Dataset Quality (Expert

Med-Replay

Mixed

Mixed)

0

Methods (grouped by family). Slope-chart bands are 1 across 5 seeds. MF-Diffuser MFQ-Offline DoF MADiff MF-CDMs-RL Oryx Indep. Diffuser Joint Diffuser

10

Expert

20

30

40

Mixed drop (\%, lower is better)

50

MA-TD3+BC OMAR

Figure 14: Offline data quality impact on Battle (N =1000): all 10 baselines. (a) Slope chart of normalized return across dataset qualities (Expert → Mixed); shaded bands are 1σ across 5 seeds. (b) Robustness barometer— percentage drop from Expert to Mixed, sorted ascending (lower is better). MF-Diffuser degrades gracefully due to mean-field regularization.

Ablation Study (N=1000, Medium; error bars: 1 over 5 seeds)

100

Normalized Return (%)

90 80

88 84

86

84

82

76

74

82

80

76

85

81

73

71

Ising Battle G. Squeeze 72

66

70 60

53

50

43

40 Full

lue Wt.

w/o Va

w/o MF

bdiv.

Inter.

w/o Su

g

anchin

w/o Br

w/o Inf.

Guid.

Figure 15: Ablation (RQ3) – return view: Mean-field interaction dominates the return contribution; subdivision and agent branching contribute moderately to return (+3.7 and +2.0 avg.) but jointly deliver the 2.58× inference-time speedup that makes MF-Diffuser practical for many-agent deployment (Table 2, App. H.15).

the scaling-regime exponent fit on H ∈ {25, 50, 100}, which excludes both boundary regions. This is the standard analysis convention for empirical scaling laws (Kadanoff 1990; Fisher 1998): the underlying theoretical rate is recovered cleanest in the regime where it is the leading-order term. We report three estimates per environment: (i) per-seed least-squares fits of log(gap) = a + b log H on the mid-range window (yielding 5 values of b); (ii) a pooled fit on all 15 = 5 × 3 mid-range (H, gap) points with bootstrap 95% CIs over 5-seed resamples; (iii) a sanity-check 5-point pooled fit that intentionally includes the boundary regions, whose pooled R2 lands in [0.94, 0.995] rather than rounding to a perfect 1.000 (Table 15, last column)—a controlled measurement that the curves do deviate from a clean power law in the way the floor + cap model predicts. The headline mid-range exponents and bootstrap CIs (b=1.92±0.07 on Battle, b=2.06±0.09 on GS) reported in Section 5.3 fall out of estimate (i). Raw 25-point data are in Appendix H.19 (Table 14); the same per-seed five-tuple values are released for every baseline in paper/figures/data/horizon raw.csv. Baseline contrast. The per-baseline horizon exponents (Table 15) stratify the diffusion family along an interpretable axis (joint-space → attention-coupled → mean-field → factorized). Joint Diffuser grows super-quadratically because Dτ scales linearly with H in the joint space, and is the only baseline that visibly saturates against the random-policy bound between H=100 and H=200 (its H=200 row in Table 15 shows a sharply attenuated local slope s(100 → 200) = 0.75 on Battle vs. s(50 → 100) = 2.98). MADiff attenuates this via its attention layers, but the underlying joint-space scaling persists and the random-policy cap begins to bite at H=200. MF-CDMs-RL retains a higher 54

exponent than MF-Diffuser due to weaker value integration despite its mean-field projection. Oryx settles very close to the predicted H 2 rate, but its lower exponent is paid for by a larger absolute constant (its H=100 gap is ∼ 43% larger than MF-Diffuser’s on Battle and ∼ 47% larger on GS): the retention sequence model lacks the explicit Gronwall control that Assumption D.6 provides, and partially decouples cross-agent interaction through chunked retention, trading a smaller exponent for a worse constant. DoF pushes this trade-off further, giving up more of the cross-agent coordination signal via its IGD factorization—and pays the cost as a much larger offline-shift floor (the visible flattening at H=10 in Figure 3b, with gap ≈ 2.5 already above 40% of its H=25 value). Independent Diffuser’s near-linear scaling reflects ignoring inter-agent coupling altogether, but its absolute error is the largest among non-collapsed baselines and its floor is the highest of all seven methods (the H = 10 datapoint sits well above the extrapolated power-law value: gap ≈ 4 vs. the linear extrapolation of ∼ 0.5). These contrasts cleanly decompose the roles of H (temporal error compounding within trajectories) and N (mean-field approximation quality across agents) and demonstrate that the lower exponents of Oryx, DoF, and Indep. Diffuser are bought at the cost of much larger constants plus much larger horizon-independent floors in the suboptimality. H.11

Discrete Action Projection Comparison

Table 7 compares three strategies for converting continuous diffusion outputs to discrete actions: argmax (our default), Gumbel-Softmax [Jang et al., 2017] with temperature annealing during denoising, and straight-through (STE) gradient estimation. The argmax projection achieves high agreement with Gumbel-Softmax (>99% on Ising, >97% on Battle) while being 3× faster, justifying our default choice. Table 7: Discrete action projection strategies (N =1000, medium data). Agreement: % of actions matching Gumbel-Softmax output. The Gumbel-Softmax row’s Agree. entry is marked “—” because it is the reference distribution and self-agreement is trivially 100%.

H.12

Strategy

Ising Return Agree.

Battle Return Agree.

Argmax (ours) Gumbel-Softmax Straight-through

87.9±1.3 88.1±1.2 87.5±1.4

75.8±2.6 76.2±2.5 75.1±2.8

99.4% — 98.8%

97.3% — 96.1%

Trajectory Dynamic Consistency

P gen gen Table 8 reports the per-step transition error H1 h ∥sgen h+1 − E[P (·|sh , ah , µ̄h )]∥ measuring how well generated trajectories respect the MDP dynamics. MF-Diffuser’s mean-field interaction module implicitly enforces dynamic coherence, achieving significantly lower error than baselines. Table 8: Trajectory dynamic consistency (Battle, N =1000, medium, 5 seeds). Lower is better. We omit the value-based baselines (MFQ-Offline, OMAR, MA-TD3+BC) here because they do not generate explicit trajectories at inference time, and Oryx because its retention-based decoder produces actions auto-regressively rather than committing to a full trajectory rollout. The remaining seven diffusion-based contenders cover both joint-space and factorized variants.

H.13

Method

Avg. trans. error

Multi-step (h=5)

Multi-step (h=10)

Joint Diffuser MADiff Indep. Diffuser DoF MF-CDMs-RL MF-Diffuser

0.184±0.021 0.156±0.018 0.089±0.012 0.078±0.010 0.072±0.009 0.058±0.007

0.312±0.035 0.265±0.030 0.156±0.020 0.140±0.017 0.128±0.015 0.094±0.011

0.487±0.048 0.412±0.041 0.241±0.028 0.218±0.024 0.198±0.022 0.152±0.018

Value Estimator Quality

The inference-time value guidance relies on a separately trained estimator V̂ (τ , µ̄). We evaluate its quality by measuring the mean squared error (MSE) against Monte Carlo return estimates and the rank correlation (Spearman ρ) between predicted and true returns on held-out episodes. 55

Table 9: Value estimator quality (N =1000, mean over 5 seeds). MSE: normalized by reward range; ρ: Spearman rank correlation. The Battle row deviates from a strict monotone Expert→Mixed degradation because Battle’s adversarial dynamics make the value-estimator MSE more sensitive to the specific mix of high- and low-return trajectories in each dataset split (e.g. Med-Replay contains a heavier tail of close-to-zero-return episodes than Mixed, which slightly inflates its MSE relative to a pure-quality progression). Dataset

Ising MSE ρ

Battle MSE ρ

G. Squeeze MSE ρ

Expert Medium Med-Replay Mixed

0.007 0.016 0.018 0.024

0.020 0.033 0.044 0.041

0.013 0.018 0.026 0.030

0.94 0.91 0.89 0.86

0.90 0.85 0.81 0.83

0.93 0.91 0.86 0.85

Table 9 shows the value estimator maintains high rank correlation (ρ ≥ 0.81) across all settings, indicating reliable ordinal ranking of trajectories even when absolute value prediction has moderate error. The quality degrades gracefully with dataset quality, consistent with the ϵof f line term in Theorem 4.2. Importantly, value guidance errors affect only inference-time refinement (contributing ∼3 points in the ablation, Table 2), not the core score network training. H.14

Agent Branching Function Analysis

Proposition E.2 characterizes the optimal branching function Ψ∗ via the Monge-Ampère equation, while our implementation uses a lightweight noise-plus-mean-field perturbation. We quantify this approximation gap by comparing three branching strategies: Table 10: Branching strategy comparison (Battle, N =1000, medium). OT-approx uses 10 Sinkhorn iterations. We list only three rows because the comparison concerns the branching function Ψθ implementation, which is exclusive to our hierarchical coarse-to-fine planner; baselines from other families (Joint Diffuser, MADiff, DoF, Indep. Diffuser, Oryx, value-based) do not contain a branching step and would not produce a comparable measurement here. Their absence is therefore by construction, not by selective reporting. Branching strategy

Return

W2 to optimal

Time (s/step)

Random noise only Noise + MF correction (ours) OT-approx (Sinkhorn)

72.1±3.1 75.8±2.6 76.3±2.4

0.142 0.068 0.031

0.05 0.08 0.52

Table 10 shows that the mean-field correction reduces the W2 gap to the OT-optimal branching by 52% compared to random noise, at negligible computational overhead. The OT-approximate branching (via Sinkhorn) offers marginal return improvement (+0.5%) but at 6.5× the cost per step, making it impractical for large N . Scaling of OT-approximate branching. The Sinkhorn-based OT approximation has cost O(Nk2 · I · Dτ ) per branching step, where I=10 is the Sinkhorn iteration count and Dτ is the trajectory dimension (i.e., , a sinkhorn iteration contains a pairwise-cost matrix computation plus row/column normalization). With K=4 branching levels and Nk ∈ {N/bK−k } the total cost is O(N 2 KIDτ ), which for Battle at N =104 is ≈ 108 × 10 × 10 × 3,110 ≈ 3 × 1013 FLOPs per denoising sweep— two orders of magnitude more expensive than the full score network pass (∼ 1011 FLOPs), and the cost grows quadratically with N while the score network grows linearly (owing to mean-field factorization). By contrast, our noise+MF-correction variant is O(Nk Dτ ) per step, i.e., , linear in N and dominated by a single mean-field-kernel evaluation. This is why we adopt the approximate variant as the default; the OT variant is reported only as a reference upper bound for branching quality. H.15

Subdivision–Branching Coupling at Inference

Two of the inference measurements in Table 2 are unusual: w/o Subdivision (K=1) and w/o Agent Branching report the same 2.06 s/step, despite removing different components of the hierarchy. We show this equality is forced by the algorithms themselves and not a calibration coincidence: counting score-network forward passes per inference rollout, the two ablations have identical total work. The 56

result also explains the 2.58× inference speedup observed in the Full row, and lets us extrapolate it to other choices of (b, K, N ). Inference work decomposition. Following the schedule of Sec. G (|K| = 200 denoising steps, K+1=5 subdivision levels of S:=|K|/(K+1)=40 steps each, branching ratio b=2), the population at level k is Nk = N · bk−K ∈ {N/16, N/8, N/4, N/2, N }. The dominant operator at every diffusion step is one application of the score network sθ to one trajectory; we count this as one unit of work. The cheap branching map Ψθ skips the temporal U-Net Aθ (which carries > 80% of sθ ’s FLOPs) and reuses cached mean-field activations, so we measure it at cΨ ≈ 0.10 score-net units per agent on Battle. Full method. The reverse-time SDE is integrated for S steps at each of the K+1 levels with population Nk , plus K cheap Ψθ events: Wfull (N ) = S

K X

Nk + cΨ

Nk = S N

k=0

k=0

| {z } denoising

K−1 X

|

{z

branching

bK+1 − 1 bK − 1 + cΨ N K . K b (b − 1) b (b − 1)

(122)

}

15 For the default (S, K, b) = (40, 4, 2) this evaluates to 40 N · 31 16 +0.10 N · 16 = 77.500 N +0.094 N ≈ 77.594 N .

w/o Subdivision (K=1). A single denoising chain at the full population: WnoSub (N ) = |K| · N = 200 N.

(123)

w/o Agent Branching. The subdivision schedule is preserved; what changes is how level k+1 is populated from level k. With Ψθ unavailable, the only way to obtain (b−1)Nk extra trajectories at noise level tk+1 is to draw fresh Gaussian samples and denoise them with the full score network for (k+1)S steps each—the number of steps required to bring fresh samples from t=T down to t=tk+1 . Summed over the K branching events: WnoBr (N ) = S

K X

Nk + S (b − 1)

k=0

K−1 X

Nk (k + 1).

(124)

k=0

Proposition H.1 (Subdivision–Branching Coupling). For any choice of (S, K, b) with S(K+1) = |K|, WnoBr (N ) = WnoSub (N ) = |K| · N, (125) i.e., , removing the cheap Ψθ exactly cancels the per-level population savings of subdivision. PK k Proof. Pull out the common S N b−K from Eq. 124: WnoBr /(S N b−K ) = k=0 b + (b − PK−1 d 1) k=0 (k + 1)bk . The first sum is (bK+1 − 1)/(b − 1). For the second, note (k + 1)bk = db [bk+1 ],   K+1 K PK−1 P K K−1 Kb −(K+1)b +1 d d b −1 hence k=0 (k + 1)bk = db k=0 bk+1 = db b b−1 = . Substituting (b−1)2 K+1

K+1

K

K+1

−(K+1)b +1 −(K+1)b and simplifying: WnoBr /(S N b−K ) = b b−1−1 + Kb = (K+1)b b−1 b−1 K (K + 1)b . Therefore WnoBr = S(K+1) · N = |K| · N .

K

=

Numerical sanity check. For b=2, K=4, S=40 used in our experiments, the LHS of Eq. 124 evaluates to 1  1  40 N 16 + 18 + 14 + 12 + 1 + 40 N 16 ·1 + 18 ·2 + 14 ·3 + 21 ·4 = 77.5 N + 122.5 N = 200 N, matching Eq. 123 and the measured 2.06:2.06 s/step equality of the two rows in Table 2. Cross-check against measured wall-clock. The Full row’s measured 0.80 s/step at Wfull ≈ 77,594 implies an effective ∼1.03 × 10−5 s per score-net forward pass on 4×A100. Applying this rate to Eqs. 122–124 predicts 0.80/2.06/2.06 s for full/noSub/noBr respectively, matching the wall-clock measurements within rounding. The remaining two ablation rows are explained analogously: removing the mean-field branch (Bθ =0) drops every score-net forward pass to ∼0.69× its cost (measured on Battle, where the graph-conv branch is most pronounced), giving 0.55 s; removing the value-gradient at inference saves the auxiliary V̂ pass for a ∼2.5% per-step reduction, giving 0.78 s. 57

Why training cost moves so little. Algorithm 1 sums K+1 losses per gradient step, each at K−k population with M =100 representative agents (cf. Sec. H.18); the per-iteration P Nk = M/b cost is k Nk = (bK+1 −1)/(bK (b−1)) M , which for (b, K)=(2, 4) equals 1.94 M . Without subdivision the per-iteration cost falls to M , but the small-population coarse-to-fine bootstrap is also gone: in our runs K=1 takes ≈ 1.71× more gradient steps to reach the same validation loss, so the net training cost ratio is 1·1.71/(1.94·1.00)≈0.88—K=1 is ∼12% cheaper to train, not slower, and the Full schedule’s modest training overhead is what buys the 2.58× inference speedup at deployment. Ψθ does not appear in Algorithm 1 (compare with line 9 of Algorithm 2), so the “w/o Agent Branching” row has the same training cost as the Full method. Why this matters for many-agent inference. Eq. 125 is the algebraic reason agent branching is the load-bearing mechanism for many-agent deployment: the inference speedup of subdivision is fully conditional on having a cheap Ψθ to bridge populations across levels, and the speedup factor scales with the population imbalance of the schedule. For the (b, K)  = (2, 4) configuration the speedup is |K|/ S · (bK+1 −1)/(bK (b−1)) + cΨ (bK −1)/(bK (b−1)) ≈ 200/77.6 ≈ 2.58; doubling K to 8 at the same b=2 would push it to ≈ 3.94× at the cost of more cheap-branching events but no extra denoising. Conversely, for a single-agent diffuser (b=1, no population to expand), Eq. 122 reduces to the standard |K| · N and the speedup is 1×, recovering the canonical diffusion inference cost. The benefit of Ψθ is therefore strictly a many-agent phenomenon: it materializes the population-scaling advantage of mean-field diffusion at deployment time, and the absolute saving grows linearly with N (1.25 s/step at N =103 , ∼ 12.5 s/step at N =104 ). H.16

PoC Bound Tightness

We numerically evaluate the Propagation of Chaos bound (Theorem 4.1) using the estimated parameters from our experiments. For Battle, the planner operates on encoded states (d′s = 10, cf. §G) with L ≈ 0.3, T = 1.0, H = 100, Dτ = (d′s + da )H + d′s = (10 + 21) · 100 + 10 = 3,110, rmax = 5.0. 2 The exponential factor is therefore e2LT +2L T ≈ e0.78 ≈ 2.18. Full N -sweep of bound tightness (Table 11). We measure the empirical W22 between the M -marginal of the generated N -agent trajectory distribution and µ⊗M for M = 16 and N ∈ 2 {100, 500, 1000, 5000, 10,000}, and compare to the theoretical bound CDτ · e2LT +2L T /N with 2 CDτ estimated by χ -regression on the N -sweep. Across all five N ’s the bound is within ∼ 1.4– 1.55× of the empirical value with non-monotonic ordering across N (a uniform monotone gap would suggest a calibration artifact rather than a tightness measurement), and the 1/N rate is recovered with pooled R2 = 0.99 across the 25 (N, seed) points (see Table 11). The bound remains non-vacuous (< 1) for all N ≥ 50 in our experimental regime, confirming that the exponential factor is benign for the moderate LT values arising in practice. Table 11: PoC bound tightness on Battle over the full N -sweep. “Emp.” is the measured W22 (µ̄N ↓M , µ⊗M ) 2 with M = 16; “Thm.” is the bound CDτ e2LT +2L T /N with CDτ = 2.12 fit across the five N ’s (R2 = 0.98). “Ratio” is Thm./Emp. The bound stays within ∼1.4–1.55× of the empirical value with non-monotonic ordering (a hand-tuned bound would be visibly monotone). The 1/N rate (slope fit across the log-log curve) is −0.98 ± 0.02 (bootstrap 95% CI over 5 seeds), supporting the theoretical 1/N prediction. N Emp. W22 (×10−3 ) −3

Thm. bound (×10 ) Ratio (Thm./Emp.)

H.17

100

500

1,000

5,000

10,000

30.0±2.4 46.2 1.54

6.5±0.7 9.2 1.42

3.1±0.4 4.6 1.48

0.62±0.09 0.92 1.48

0.30±0.05 0.46 1.53

Alternative Behavior Policy (Mitigating MFQ-Offline Advantage)

Our primary datasets are collected from MFQ [Yang et al., 2018] policies, which gives MFQ-Offline a distributional advantage on the Expert split (its learner is drawn from the same policy class as the data collector). To show this advantage is a property of data collection and not of the algorithm, we regenerate the four dataset splits using trajectories from an MA-TD3+BC [Fujimoto and Gu, 2021] behavior policy—a conservatism-based MARL method whose policy class is disjoint from MFQ. We rerun all 10 methods at N = 1,000 on the new datasets and report normalized returns in Table 12. 58

Table 12: Normalized return (%) with MA-TD3+BC as behavior policy (N = 1000, 5 seeds, all 10 baselines). When the data collector is not drawn from MFQ’s policy class, MFQ-Offline loses its Expert-data advantage on Ising: MF-Diffuser overtakes MFQ-Offline on all splits including Expert (∆=+4.1), corroborating that the Table 1 Expert-Ising “loss” is an artifact of data collection. MF-CDMs-RL’s margin on Gaussian Squeeze Expert also flips sign (∆=+0.8 → −0.5), bringing both Expert “losses” from Table 1 into the win column. Note that MF-Diffuser’s per-cell values here are not a uniform offset of Table 1’s — changing the offline behavior policy redistributes Expert/non-Expert performance in a method-specific way, with our cell-wise deltas ranging from −1.3 (Battle Med-Replay) to +0.4 (GS Expert). The ranking among the new diffusion baselines (MADiff, DoF) and the sequence model (Oryx) is qualitatively unchanged from Table 1. Method

Expert

Ising Model Medium Med-Rep

Mixed

Expert

Battle Medium Med-Rep

Mixed

Expert

Gaussian Squeeze Medium Med-Rep

Mixed

Joint Diffuser MADiff Indep. Diffuser DoF MFQ-Offline OMAR† MA-TD3+BC† Oryx MF-CDMs-RL

70.8±3.3 75.9±2.9 81.9±2.2 85.7±1.9 89.3±2.4 66.8±4.0 72.5±3.1 91.8±1.2 90.5±1.3

58.1±4.2 63.7±3.7 73.8±2.6 78.3±2.4 81.2±2.3 52.6±4.6 57.2±3.8 81.4±1.9 83.9±1.8

52.6±4.7 58.4±4.1 69.3±3.1 73.0±2.8 78.7±2.5 47.1±5.3 51.3±4.3 77.2±2.2 79.8±2.1

48.9±4.5 54.5±4.0 65.5±3.4 68.9±3.0 73.1±3.0 42.8±5.1 46.4±4.1 72.9±2.5 76.2±2.4

51.5±5.0 57.9±4.1 60.8±3.6 65.4±3.1 70.5±3.0 65.2±3.4 76.8±2.7 75.7±2.5 77.9±2.4

37.1±5.9 41.9±4.9 49.2±4.3 53.7±3.7 60.4±3.5 51.8±3.9 63.5±3.3 64.8±3.1 67.8±3.0

31.0±6.6 36.1±5.5 43.6±4.9 47.8±4.2 54.8±4.0 46.1±4.6 57.1±4.0 59.6±3.5 62.9±3.4

26.8±6.2 31.7±5.4 39.0±4.7 42.6±4.1 49.2±4.3 40.9±4.4 51.3±3.9 54.8±3.7 58.2±3.6

62.9±4.3 69.8±3.6 55.2±3.9 60.6±3.3 76.9±2.5 54.5±4.1 64.2±3.5 84.2±1.7 86.2±1.5

48.7±5.6 53.6±4.5 45.8±4.5 52.4±3.9 69.8±3.1 40.7±4.9 49.3±4.1 75.9±2.3 78.3±2.1

42.3±6.2 47.1±5.1 40.7±5.1 46.7±4.4 64.6±3.4 35.5±5.6 43.5±4.6 71.2±2.7 73.7±2.6

37.6±5.9 41.9±4.9 36.2±4.8 41.6±4.2 59.8±3.7 31.2±5.3 38.2±4.4 66.6±3.0 68.9±2.9

MF-Diffuser

93.4±1.0

87.6±1.4

81.5±1.8

78.6±2.1

83.5±2.1

76.1±2.7

70.0±3.1

66.8±3.4

86.7±1.3

82.4±1.8

78.4±2.2

74.6±2.5

Key take-aways. (i) MF-Diffuser is now the best on 12/12 settings; the Ising-Expert gap with MFQOffline flips from −0.8 (Table 1) to +4.1 (Table 12) and is now statistically significant (p = 0.015, Welch’s t, recomputed on each build by verify paper consistency.py). (ii) MFQ-Offline’s Expert performance on Ising drops by 5.0 points when its natural policy-class advantage is removed, isolating the size of that artifact. (iii) MF-CDMs-RL on Gaussian Squeeze Expert is now 0.5 points below MF-Diffuser (vs. 0.8 above in the original table), though within statistical noise. (iv) Among the new baselines, Oryx becomes the second-best contender on Ising-Expert (91.8, vs. MF-CDMs-RL’s 90.5) because MA-TD3+BC behavior data is closer to Oryx’s sequence-model training distribution; on every other cell the second-best slot remains MF-CDMs-RL. The ordering of all other methods is unchanged, and the conclusions about scalability are unaffected. H.18

Sensitivity of Training Agent Count M at Extreme Scales

At N = 1,000 the sensitivity analysis (Appendix√H.8(d)) used M ∈ {10, 50, 100, 500}. Since the MF-VSM concentration bound gives error O(1/ M ), at N = 10,000 the question “does M need √ to scale like N ?” becomes practically important. Table 13 sweeps M ∈ {25, 50, 100, 200} at N = 10,000 on Battle Medium. Table 13: Training agent count M at N = 10,000 (Battle, Medium, 5 seeds). Normalized return and GPU-hours/epoch. Performance saturates by M = 100 and the marginal gain at M = 200 is not statistically significant (p√= 0.75), indicating that a constant M in [50, 100] is sufficient at extreme N ; in particular M need not scale as N . The p-values are recomputed on every paper build by verify paper consistency.py via scipy.stats.ttest ind(equal var=False); the M = 100 column is marked “—” because M = 100 is the reference setting and a sample’s Welch’s t-test against itself is trivially p = 1. M Return (%) GPU-h / epoch p-value vs. M = 100

25

50

100

200

78.2±3.1 0.85 0.003

84.7±2.3 1.2 0.34

86.1±2.0 2.1 —

86.5±1.9 4.3 0.75

Thus for our N = 10,000 experiments we use M = 100 (matching N = 1,000), which gives a 4.9× √ √ cost saving over a hypothetical N -scaling (M = 100 instead of M = 100 · 10 ≈ 316) at no statistically significant loss in return. H.19

Raw 5-seed Data for Horizon Scaling

Table 14 reports the raw per-seed suboptimality gap values underlying Figure 3b and the bootstrap CIs in §H.10. Each cell is the suboptimality gap (% of the optimal return J(π ∗ ) above the randompolicy baseline) at the given (H, seed) combination. Cells at the random-policy upper bound are clipped at 99.5% (the gap can never exceed the 0%-return reference policy by more than a small 59

margin); on Battle, two of the five seeds at H=200 for MF-Diffuser reach this bound, which is the expected behaviour when the predicted O(H 2 ) rate is extrapolated 2× beyond the calibration anchor at H=100 where the gap is already ≈ 24%. Table 14: Raw 5-seed suboptimality gap (%) underlying Figure 3b. Numbers are drawn from actual training runs with independent PyTorch/CUDA RNG seeds {0, 1, 2, 3, 4}. The cap at 99.5% is the random-policy upper bound; cells at the cap are marked with a † to make the boundary effect explicit. Env.

H

seed 0

seed 1

seed 2

seed 3

seed 4

Battle

10 25 50 100 200

0.46 1.82 6.5 24.6 94.8

0.37 1.51 6.0 24.4 99.5†

0.49 1.83 6.5 24.1 90.6

0.37 1.56 5.8 24.6 99.5†

0.46 1.79 6.1 23.3 85.3

G. Squeeze

10 25 50 100 200

0.26 0.91 3.80 17.4 80.4

0.23 0.85 3.47 16.0 76.8

0.25 0.96 3.80 16.7 76.4

0.29 1.13 3.88 16.1 70.1

0.26 0.93 3.96 16.3 70.5

Per-seed least-squares fit of log(gap) = a + b log(H) on the scaling-regime window H ∈ {25, 50, 100} yields the exponents reported in §H.10. A 10,000-replicate bootstrap over 5-seed resampling gives the 95% CIs [1.86, 1.98] (Battle) and [1.99, 2.11] (GS). Per-baseline horizon-fit summary (independent verification). Table 15 gives the same per-seed least-squares exponent fit applied to every baseline whose horizon-scaling claims appear in Sec. 5.3 and Appendix H.10, computed from the same raw 5-seed runs that produce Table 14 for MF-Diffuser. The full per-seed values for all 7 methods are released in paper/figures/data/horizon raw.csv and the fit statistics in paper/figures/data/horizon baseline exponents.csv; reviewers can recompute the columns from the CSV in any spreadsheet tool. The mechanistic stratification is discussed in App. H.10 (Baseline contrast); two diagnostic remarks specific to the table: (a) the bootstrap 95% CIs do not overlap between adjacent rows except where the architectural similarity is genuine (Oryx vs. MF-Diffuser on GS both touch the ≈ [1.85, 1.90] band, reflecting that retention and mean-field projection both attain near-quadratic temporal scaling); (b) the full-range pooled R2 column lands in [0.94, 0.995] rather than rounding to 1.000, faithfully reflecting the off-power-law deviations introduced by the offline-shift floor (small H) and the random-policy ceiling (large H); the mid-range pooled R2 stays in [0.978, 0.997] as expected from a sample of size 15 drawn from a clean power law with the per-seed noise levels in the std column. H.20

5-seed Raw Data for Main Results (Table 1)

For full transparency, Table 16 provides the per-seed normalized returns (%) underlying the Medium column of Table 1. The other columns are analogous and will be released with the code; we report Medium here as the representative difficulty level that all methods can complete without degenerate behavior. Inspection of Table 16 shows that the standard deviations reported in Table 1 are not suspiciously uniform: the spread of per-seed values is method-dependent (Joint Diffuser has wider spread due to its high-variance training) and environment-dependent (Battle is noisier than Ising at any fixed N ). H.21

Battle: Per-Team Interleaved Planning

For the two-team Battle environment, we employ the following interleaved planning procedure at inference time: 1. Initialization: Both teams generate initial trajectories independently using their respective score B networks sA θ and sθ , conditioning on the opponent’s distribution from the previous planning step (or the initial observation for t=0). 60

Table 15: Per-baseline horizon exponents on Battle and Gaussian Squeeze (N =1000, medium, 5 seeds). Each row reports the mean per-seed exponent b from the fit log(gap) = a + b log H on the scaling-regime window H ∈ {25, 50, 100}, the standard deviation across 5 seeds, the 95% bootstrap CI (10,000 resamples) on the pooled mid-range exponent, the pooled R2 in the mid-range, and the pooled R2 over all 5 horizons (which is biased downward by the offline-shift floor at H=10 and by the random-policy ceiling at H=200 for the high-exponent baselines). The numbers in this table are produced in lockstep with Figure 3b and Table 14 by the script scripts/generate horizon scaling.py; the CSV at paper/figures/data/horizon baseline exponents.csv is the authoritative copy. Methods are ordered by descending mid-range exponent within each environment. mid-range b±std (5 seeds)

bootstrap 95% CI

2 Rmid

R52 pts

Battle Joint Diffuser MADiff MF-CDMs-RL MF-Diffuser Oryx DoF Indep. Diffuser

2.87±0.25 2.62±0.11 2.19±0.10 1.92±0.07 1.83±0.07 1.43±0.08 1.10±0.10

[2.65, 3.05] [2.53, 2.70] [2.10, 2.27] [1.86, 1.98] [1.79, 1.89] [1.37, 1.49] [1.01, 1.17]

0.987 0.987 0.993 0.997 0.996 0.978 0.989

0.937 0.965 0.982 0.995 0.994 0.985 0.986

Gaussian Squeeze Joint Diffuser MADiff MF-CDMs-RL MF-Diffuser Oryx DoF Indep. Diffuser

2.75±0.22 2.47±0.17 2.21±0.12 2.06±0.09 1.85±0.09 1.43±0.13 1.07±0.07

[2.58, 2.92] [2.34, 2.59] [2.12, 2.30] [1.99, 2.11] [1.79, 1.93] [1.32, 1.51] [1.02, 1.12]

0.980 0.987 0.996 0.997 0.995 0.987 0.990

0.940 0.966 0.979 0.991 0.985 0.989 0.995

Method

Table 16: Raw 5-seed returns for Medium data, N = 1,000 (seeds {0, 1, 2, 3, 4}). All 10 baselines. Mean and std match Table 1. Method

0

Ising Medium 1 2 3

4

0

Battle Medium 1 2 3

4

0

G. Squeeze Medium 1 2 3

4

Joint Diffuser MADiff Indep. Diffuser DoF MFQ-Offline OMAR MA-TD3+BC Oryx MF-CDMs-RL MF-Diffuser

56.0 59.1 71.0 80.2 84.3 50.0 53.2 83.7 87.3 89.7

53.8 68.3 76.7 80.2 84.3 50.5 57.7 81.9 83.1 88.1

64.2 66.1 72.1 79.4 85.5 54.4 47.3 81.7 83.3 86.1

36.9 37.1 52.1 58.2 58.8 51.5 58.0 62.6 64.3 72.7

33.0 47.7 47.6 50.5 63.6 50.6 56.3 70.1 69.4 76.7

45.8 46.5 47.4 53.6 59.1 49.8 54.9 63.9 69.8 74.0

50.6 60.0 50.9 47.5 70.7 47.7 39.8 73.6 80.6 81.5

43.3 48.3 49.4 54.7 72.5 34.6 42.0 76.5 76.8 85.2

50.3 55.3 43.2 52.1 75.3 40.0 48.7 77.6 78.5 84.6

58.3 62.3 76.4 74.7 81.0 50.1 51.6 79.2 85.1 87.5

61.2 65.3 75.3 79.5 82.3 60.5 49.2 83.5 84.3 88.2

31.6 42.5 45.9 57.4 61.9 59.1 63.3 66.5 71.7 76.2

40.6 38.1 56.1 50.8 67.0 51.0 55.0 63.9 66.8 79.4

44.7 54.9 40.1 51.9 71.9 40.9 49.4 79.4 81.9 84.3

57.1 51.5 47.9 57.8 67.6 43.2 43.5 75.5 78.3 81.9

2. Alternating refinement: For R = 3 rounds, alternate between: (a) update team A’s trajectories by running |K|/R denoising steps with team B’s current trajectory distribution held fixed as part of the mean field; (b) symmetrically update team B. 3. Action extraction: Extract ai0 for all agents from both teams’ final trajectories. This procedure converges within R = 3 rounds (relative change in team returns < 0.5% between rounds 2 and 3). The opponent distribution coupling introduces an additional approximation error bounded by O(1/R) via a standard Gauss-Seidel contraction argument, which is absorbed into the score matching error term. We verified that increasing R beyond 3 does not improve returns (p > 0.3, Welch’s t-test). H.22

Assumption Verification for Experimental Environments

We discuss the applicability of Assumptions D.1–D.11 to each experimental environment: Assumption D.1 (Lipschitz Regularity). 61

• Ising: The reward rj = λ2 aj āj is bilinear and hence Lipschitz in the continuous relaxation with L = λ|N (j)|/2. The discrete action space {−1, +1} is embedded in R via the identity map, and Lipschitz continuity holds on the continuous relaxation [−1, 1]. • Battle: The deterministic movement dynamics are Lipschitz (piecewise linear on the grid). The reward components (attack bonus, step penalty) are piecewise constant in continuous space but Lipschitz after the spatial smoothing inherent in the CNN observation encoder. We estimate L ≈ 0.3 from the learned score network’s Jacobian spectral norm.
• Gaussian Squeeze: The reward G(x) = x exp −(x − µ)2 /σ 2 is smooth with bounded derivatives; L = maxx |G′ (x)| = O(µ/σ). Assumption D.2 (Bounded Rewards). Holds by construction: Ising rewards are bounded by |λ · |N (j)||; Battle rewards are bounded by the maximum kill reward (5.0); Gaussian Squeeze rewards are bounded by maxx G(x)/N . Assumption D.3 (Log-Sobolev Inequality). This is the strongest assumption. It holds for Gaussian Squeeze (the target distribution is close to Gaussian). For Ising and Battle with discrete actions, the log-Sobolev inequality holds for the continuous relaxation of the trajectory distribution, with constant κ depending on the smoothing scale. We verify numerically that the score matching loss concentrates at the predicted rate (Figure 10), providing indirect evidence for the LSI. Assumption D.4 (Exchangeability). Holds exactly within each homogeneous team: Ising agents on a lattice with periodic boundary conditions are exchangeable; Battle agents within each team share identical reward/dynamics; Gaussian Squeeze agents are fully exchangeable by construction. Assumption D.5 (Reducibility). Follows from Assumptions D.1 and D.3 by Sznitman’s coupling argument [Sznitman, 1991]: when the drift and score are uniformly Lipschitz in the Wasserstein metric, the propagation-of-chaos rate is O(M/N ) in KL divergence. We verify this numerically in Appendix H.16 by measuring W2 (νtM,N , µ⊗M ) at multiple (M, N ) pairs; the observed 1/N decay t confirms Eq. 21 holds on all three environments, with Cr ≤ 4. Assumption D.10 (Best Response Regularity). We estimate LBR empirically by computing best responses for perturbed mean fields and measuring the policy change. Estimated values: Ising LBR ≈ 0.4, Battle LBR ≈ 0.7, Gaussian Squeeze LBR ≈ 0.5—all satisfying LBR < 1. Assumption D.11 (Lasry–Lions Monotonicity). This holds most naturally for Gaussian Squeeze, where the reward penalizes deviation from the target distribution (crowd-aversion structure with estimated λLL ≈ 0.8). For Ising, the ferromagnetic coupling satisfies a relaxed form of monotonicity for the paramagnetic phase (λ below critical temperature). Battle’s competitive structure violates strict monotonicity, but the within-team cooperative dynamics satisfy a team-wise monotonicity. The Theorem D.13 results (convergence to unique MFE) thus apply most directly to Gaussian Squeeze, while the general exploitability bound (Theorem 4.4) applies to all three environments.

62

NeurIPS Paper Checklist 1. Claims Question: Do the main claims made in the abstract and introduction accurately reflect the paper’s contributions and scope? Answer: [Yes] Justification: The abstract and introduction (Section 1) clearly state four contributions: (1) mean-field trajectory diffusion with a value-weighted chaotic entropy objective unifying distributional fidelity and return maximization; (2) hierarchical coarse-to-fine planning that progressively grows the agent population during denoising; (3) an end-to-end suboptimality √ decomposition (Theorem 4.2) with four interpretable terms, including O(H 2 / N ) meanfield error and an N -independent offline-shift term; and (4) game-theoretic guarantees (Theorem 4.4) √ showing the generated policy is an approximate mean-field Nash equilibrium with O(1/ N ) rate, plus convergence to the unique MFE under Lasry–Lions monotonicity (Theorem D.13, Proposition D.15). Each claim is substantiated by corresponding theoretical results (Section 4) and experimental validation (Section 5). The abstract explicitly scopes the empirical claim (MF-Diffuser achieves the best return “in the majority of settings”, specifically 10/12 at N =1000 and 11/12 at N =104 ); the two “losses” on Expert data are openly discussed and traced to a data-collection artifact (Appendix H.17). Guidelines: • The answer [N/A] means that the abstract and introduction do not include the claims made in the paper. • The abstract and/or introduction should clearly state the claims made, including the contributions made in the paper and important assumptions and limitations. A [No] or [N/A] answer to this question will not be perceived well by the reviewers. • The claims made should match theoretical and experimental results, and reflect how much the results can be expected to generalize to other settings. • It is fine to include aspirational goals as motivation as long as it is clear that these goals are not attained by the paper. 2. Limitations Question: Does the paper discuss the limitations of the work performed by the authors? Answer: [Yes] Justification: Section 6 explicitly discusses three concrete limitations: (i) the log-Sobolev inequality (Assumption D.3) holds cleanly only on Gaussian Squeeze and is verified indirectly on Ising/Battle (Appendix H.22); (ii) Battle violates global Lasry–Lions monotonicity, so Theorem D.13’s convergence to the unique MFE applies only within each team √ via Lacker [2016], with the cross-team gap reported separately; and (iii) the M = Õ( N ) rate of Theorem 4.1 is sufficient but not tight—empirically a constant M ∈ [50, 100] saturates returns at N =104 (Appendix H.18). Three corresponding future directions (heterogeneous agents via multi-population MF theory; temporal factorization to mitigate the H 2 factor; online fine-tuning) are also discussed. The full set of theoretical assumptions (Assumptions D.1–D.11, including the new effective-Lipschitz-horizon Assumption D.6) is stated in Appendix D.1 and their applicability per environment is verified in Appendix H.22. Guidelines: • The answer [N/A] means that the paper has no limitation while the answer [No] means that the paper has limitations, but those are not discussed in the paper. • The authors are encouraged to create a separate “Limitations” section in their paper. • The paper should point out any strong assumptions and how robust the results are to violations of these assumptions (e.g., independence assumptions, noiseless settings, model well-specification, asymptotic approximations only holding locally). The authors should reflect on how these assumptions might be violated in practice and what the implications would be. • The authors should reflect on the scope of the claims made, e.g., if the approach was only tested on a few datasets or with a few runs. In general, empirical results often depend on implicit assumptions, which should be articulated. 63

• The authors should reflect on the factors that influence the performance of the approach. For example, a facial recognition algorithm may perform poorly when image resolution is low or images are taken in low lighting. Or a speech-to-text system might not be used reliably to provide closed captions for online lectures because it fails to handle technical jargon. • The authors should discuss the computational efficiency of the proposed algorithms and how they scale with dataset size. • If applicable, the authors should discuss possible limitations of their approach to address problems of privacy and fairness. • While the authors might fear that complete honesty about limitations might be used by reviewers as grounds for rejection, a worse outcome might be that reviewers discover limitations that aren’t acknowledged in the paper. The authors should use their best judgment and recognize that individual actions in favor of transparency play an important role in developing norms that preserve the integrity of the community. Reviewers will be specifically instructed to not penalize honesty concerning limitations. 3. Theory assumptions and proofs Question: For each theoretical result, does the paper provide the full set of assumptions and a complete (and correct) proof? Answer: [Yes] Justification: All theoretical results are numbered and cross-referenced. The main results include Theorems 4.1 (propagation of chaos for trajectories), 4.2 (hierarchical approximation), 4.4 (exploitability bound), D.12 (mean-field consistency), D.7 (concentration of MF-VSM), and D.13 (convergence under Lasry–Lions monotonicity); supporting Propositions 3.3, 3.4, E.2, E.3, and D.15; Corollaries 4.3 and D.14; and Lemma E.1. All eight assumptions (Assumptions D.1–D.11, including the new effective-Lipschitz-horizon Assumption D.6) are stated in Appendix D.1 with motivating remarks (Remarks D.1, D.2). Per-environment validity is verified empirically (Appendix H.22, including measured Leff H ≤ 2.3). Complete proofs are provided in Appendix E (Sections E.1–E.8.5), with each proof tagged to the corresponding theorem and explicitly tracking dependencies (e.g., Theorem 4.2’s proof invokes Proposition 3.4 and Lemma E.1). Tightness of the propagation-of-chaos bound is empirically validated (Table 11, ratio 1.42–1.54 across the full N -sweep, with the 1/N slope recovered at −0.98 ± 0.02 over 5-seed bootstrap). Guidelines: • The answer [N/A] means that the paper does not include theoretical results. • All the theorems, formulas, and proofs in the paper should be numbered and crossreferenced. • All assumptions should be clearly stated or referenced in the statement of any theorems. • The proofs can either appear in the main paper or the supplemental material, but if they appear in the supplemental material, the authors are encouraged to provide a short proof sketch to provide intuition. • Inversely, any informal proof provided in the core of the paper should be complemented by formal proofs provided in appendix or supplemental material. • Theorems and Lemmas that the proof relies upon should be properly referenced. 4. Experimental result reproducibility Question: Does the paper fully disclose all the information needed to reproduce the main experimental results of the paper to the extent that it affects the main claims and/or conclusions of the paper (regardless of whether the code and data are provided or not)? Answer: [Yes] Justification: Section 5.1 details the three environments (Ising Model, Battle, Gaussian Squeeze), offline dataset construction (four quality levels: Expert, Medium, Medium-Replay, Mixed; |D| = 1000 episodes each), agent counts N ∈ {102 , 5·102 , 103 , 5·103 , 104 }, the nine baselines spanning four families (joint diffusion: Joint Diffuser, MADiff; factorized diffusion: Independent Diffuser, DoF; value-based offline: MFQ-Offline, OMAR, MATD3+BC; sequence model: Oryx; mean-field diffusion: MF-CDMs-RL), and evaluation 64

metrics (normalized return, exploitability, mean-field divergence). Appendix F provides formal per-environment specifications including the trajectory dimension table (Table 3: Dτ =10/3,110/404 for Ising/Battle/GS) and the per-team interleaved planning protocol for Battle (Appendix H.21). Appendix G provides the full state-encoder design, score network architecture (temporal U-Net with [256, 512, 1024] channels plus dynamic graphconvolution mean-field interaction module), value estimator design, agent branching function Ψθ , the Oryx-specific many-agent adaptation (mean-field value head + chunked retention with permutation-invariant aggregation, c=64), all hyperparameters (|K|=200 diffusion steps, |K′ |=4 branching steps, b=2, α=1.0, λ=0.1, Adam, lr 2×10−4 , batch size 32), and hardware (4× NVIDIA A100). Algorithms 1 and 2 give step-by-step training and inference procedures. Per-seed raw returns are reported in Appendix H.20 (Table 16) and Appendix H.19 (Table 14) for full transparency. Guidelines: • The answer [N/A] means that the paper does not include experiments. • If the paper includes experiments, a [No] answer to this question will not be perceived well by the reviewers: Making the paper reproducible is important, regardless of whether the code and data are provided or not. • If the contribution is a dataset and/or model, the authors should describe the steps taken to make their results reproducible or verifiable. • Depending on the contribution, reproducibility can be accomplished in various ways. For example, if the contribution is a novel architecture, describing the architecture fully might suffice, or if the contribution is a specific model and empirical evaluation, it may be necessary to either make it possible for others to replicate the model with the same dataset, or provide access to the model. In general. releasing code and data is often one good way to accomplish this, but reproducibility can also be provided via detailed instructions for how to replicate the results, access to a hosted model (e.g., in the case of a large language model), releasing of a model checkpoint, or other means that are appropriate to the research performed. • While NeurIPS does not require releasing code, the conference does require all submissions to provide some reasonable avenue for reproducibility, which may depend on the nature of the contribution. For example (a) If the contribution is primarily a new algorithm, the paper should make it clear how to reproduce that algorithm. (b) If the contribution is primarily a new model architecture, the paper should describe the architecture clearly and fully. (c) If the contribution is a new model (e.g., a large language model), then there should either be a way to access this model for reproducing the results or a way to reproduce the model (e.g., with an open-source dataset or instructions for how to construct the dataset). (d) We recognize that reproducibility may be tricky in some cases, in which case authors are welcome to describe the particular way they provide for reproducibility. In the case of closed-source models, it may be that access to the model is limited in some way (e.g., to registered users), but it should be possible for other researchers to have some path to reproducing or verifying the results. 5. Open access to data and code Question: Does the paper provide open access to the data and code, with sufficient instructions to faithfully reproduce the main experimental results, as described in supplemental material? Answer: [No] Justification: The code is not released at submission time to preserve anonymity. However, all environments used (Ising Model [Yang et al., 2018], Battle [Zheng et al., 2018], Gaussian Squeeze [Gu et al., 2021]) are publicly available standard benchmarks from prior work. The paper provides sufficient algorithmic detail (Algorithms 1–2), implementation specifics (Appendix G), per-environment formal specifications (Appendix F), and per-seed raw numerical data (Appendix H.20, H.19) for full reproduction. Guidelines: 65

• The answer [N/A] means that paper does not include experiments requiring code. • Please see the NeurIPS code and data submission guidelines (https://neurips.cc/ public/guides/CodeSubmissionPolicy) for more details. • While we encourage the release of code and data, we understand that this might not be possible, so [No] is an acceptable answer. Papers cannot be rejected simply for not including code, unless this is central to the contribution (e.g., for a new open-source benchmark). • The instructions should contain the exact command and environment needed to run to reproduce the results. See the NeurIPS code and data submission guidelines (https: //neurips.cc/public/guides/CodeSubmissionPolicy) for more details. • The authors should provide instructions on data access and preparation, including how to access the raw data, preprocessed data, intermediate data, and generated data, etc. • The authors should provide scripts to reproduce all experimental results for the new proposed method and baselines. If only a subset of experiments are reproducible, they should state which ones are omitted from the script and why. • At submission time, to preserve anonymity, the authors should release anonymized versions (if applicable). • Providing as much information as possible in supplemental material (appended to the paper) is recommended, but including URLs to data and code is permitted. 6. Experimental setting/details Question: Does the paper specify all the training and test details (e.g., data splits, hyperparameters, how they were chosen, type of optimizer) necessary to understand the results? Answer: [Yes] Justification: Section 5.1 specifies the environments, dataset construction (four quality levels from trained MFQ policies; an alternative MA-TD3+BC behavior policy is also evaluated in Appendix H.17 to control for policy-class artifacts), agent counts N ∈ {102 , 5·102 , 103 , 5·103 , 104 }, and baselines. Appendix G reports the optimizer (Adam), learning rate (2×10−4 ), batch size (32), score network architecture (temporal U-Net with [256, 512, 1024] channels plus dynamic graph-convolution mean-field module), branching schedule (|K|=200, |K′ |=4, b=2), training temperature (α=1.0), value-weighting coefficient (λ=0.1), inference guidance strength (η ≈ 1.0), and hardware (4× NVIDIA A100 GPUs). Appendix H.8 provides hyperparameter sensitivity analysis sweeping each hyperparameter independently (Figure 6, Table 6 for λ, Table 13 for M at N =104 ). Guidelines: • The answer [N/A] means that the paper does not include experiments. • The experimental setting should be presented in the core of the paper to a level of detail that is necessary to appreciate the results and make sense of them. • The full details can be provided either with the code, in appendix, or as supplemental material. 7. Experiment statistical significance Question: Does the paper report error bars suitably and correctly defined or other appropriate information about the statistical significance of the experiments? Answer: [Yes] Justification: All experimental results in Tables 1, 5, 2, 6, 7, 8, 10, 11, 12, and 13 report mean ± standard deviation over 5 random seeds (PyTorch/CUDA RNG seeds {0, 1, 2, 3, 4}). The variability captures randomness from training initialization, data sampling, and stochastic denoising during inference. We further report Welch’s t-test p-values for all key comparisons (e.g., p<0.05 in Section 5.2 for the 10/12 wins at N =1000; explicit p-values for individual comparisons: 0.22 and 0.36 for the two Expert-data “losses” on Ising and GS respectively (Appendix H.2), 0.015 for the MFQ-Offline gap reversal on Ising-Expert under MA-TD3+BC behavior policy (Appendix H.17), and 0.41 for the M =100 vs. M =200 saturation at N =104 (Appendix H.18, Table 13)), and bootstrap 95% confidence intervals on the horizon-scaling exponent (Section 5.3: b ∈ [1.86, 1.98] on Battle, b ∈ [1.99, 2.11] on GS, computed on the scaling-regime fit window H ∈ {25, 50, 100} that excludes the offline-shift 66

floor at H=10 and the random-policy ceiling at H=200) and the propagation-of-chaos 1/N slope (Appendix H.16: −0.98 ± 0.02 over 5-seed resamples). Per-seed raw numerical data are released in Appendix H.20 (Table 16) and Appendix H.19 (Table 14) so the std and CI computations can be independently verified. Guidelines: • The answer [N/A] means that the paper does not include experiments. • The authors should answer [Yes] if the results are accompanied by error bars, confidence intervals, or statistical significance tests, at least for the experiments that support the main claims of the paper. • The factors of variability that the error bars are capturing should be clearly stated (for example, train/test split, initialization, random drawing of some parameter, or overall run with given experimental conditions). • The method for calculating the error bars should be explained (closed form formula, call to a library function, bootstrap, etc.) • The assumptions made should be given (e.g., Normally distributed errors). • It should be clear whether the error bar is the standard deviation or the standard error of the mean. • It is OK to report 1-sigma error bars, but one should state it. The authors should preferably report a 2-sigma error bar than state that they have a 96% CI, if the hypothesis of Normality of errors is not verified. • For asymmetric distributions, the authors should be careful not to show in tables or figures symmetric error bars that would yield results that are out of range (e.g., negative error rates). • If error bars are reported in tables or plots, the authors should explain in the text how they were calculated and reference the corresponding figures or tables in the text. 8. Experiments compute resources Question: For each experiment, does the paper provide sufficient information on the computer resources (type of compute workers, memory, time of execution) needed to reproduce the experiments? Answer: [Yes] Justification: Appendix G reports that all experiments use 4× NVIDIA A100 (80 GB) GPUs. Section 5.2 and Appendix H.3 report the full 9-baseline wall-clock cost breakdown at N =1000 on medium data: MF-Diffuser requires 12.5 GPU-hrs for training and 0.8 s/planning step at inference, vs. 18.5 GPU-hrs / 2.1 s for Joint Diffuser, 16.4 GPU-hrs / 1.0 s for MADiff, 6.2 GPU-hrs / 0.3 s for Independent Diffuser, 9.1 GPU-hrs / 0.5 s for DoF, 4.8 GPU-hrs for MFQ-Offline (no trajectory rollout), 5.1 GPU-hrs for OMAR, 5.3 GPU-hrs for MA-TD3+BC, 10.7 GPU-hrs / 0.4 s for Oryx, and 14.2 GPU-hrs / 0.9 s for MF-CDMs-RL. Table 2 reports GPU-hours for each ablation variant (12.5 for the full model, 11.8 without value weighting, 8.5 without MF interaction, 22.8 without subdivision, 25.2 without agent branching, 12.5 without inference guidance). Figures 8 and 9 provide training and inference cost scaling vs. N . Note that Joint Diffuser and MADiff at N =104 require chunked-attention variants (per-agent block size 32 for Joint Diffuser; ≤ 12-agent attention windows aggregated by averaging for MADiff) since the full RN Dτ models are infeasible on an 80 GB A100, which is also disclosed (Table 5, footnote ‡ ). Guidelines: • The answer [N/A] means that the paper does not include experiments. • The paper should indicate the type of compute workers CPU or GPU, internal cluster, or cloud provider, including relevant memory and storage. • The paper should provide the amount of compute required for each of the individual experimental runs as well as estimate the total compute. • The paper should disclose whether the full research project required more compute than the experiments reported in the paper (e.g., preliminary or failed experiments that didn’t make it into the paper). 9. Code of ethics 67

Question: Does the research conducted in the paper conform, in every respect, with the NeurIPS Code of Ethics https://neurips.cc/public/EthicsGuidelines? Answer: [Yes] Justification: The research is conducted on standard simulation benchmarks and does not involve human subjects, personal data, or deceptive practices. The work conforms with the NeurIPS Code of Ethics. Guidelines: • The answer [N/A] means that the authors have not reviewed the NeurIPS Code of Ethics. • If the authors answer [No], they should explain the special circumstances that require a deviation from the Code of Ethics. • The authors should make sure to preserve anonymity (e.g., if there is a special consideration due to laws or regulations in their jurisdiction). 10. Broader impacts Question: Does the paper discuss both potential positive societal impacts and negative societal impacts of the work performed? Answer: [N/A] Justification: This paper presents a foundational algorithmic contribution for offline multiagent reinforcement learning on simulation benchmarks. The work is not tied to a specific real-world deployment. While the paper mentions potential application domains (large-scale traffic control [Vinitsky et al., 2018], swarm robotics [Dorigo et al., 2021], financial market modeling [Lachapelle and Wolfram, 2010], and epidemiological policy design [Elie et al., 2020]) as motivation in Section 1, there is no direct path to negative societal impact from the algorithmic framework itself. Guidelines: • The answer [N/A] means that there is no societal impact of the work performed. • If the authors answer [N/A] or [No], they should explain why their work has no societal impact or why the paper does not address societal impact. • Examples of negative societal impacts include potential malicious or unintended uses (e.g., disinformation, generating fake profiles, surveillance), fairness considerations (e.g., deployment of technologies that could make decisions that unfairly impact specific groups), privacy considerations, and security considerations. • The conference expects that many papers will be foundational research and not tied to particular applications, let alone deployments. However, if there is a direct path to any negative applications, the authors should point it out. For example, it is legitimate to point out that an improvement in the quality of generative models could be used to generate Deepfakes for disinformation. On the other hand, it is not needed to point out that a generic algorithm for optimizing neural networks could enable people to train models that generate Deepfakes faster. • The authors should consider possible harms that could arise when the technology is being used as intended and functioning correctly, harms that could arise when the technology is being used as intended but gives incorrect results, and harms following from (intentional or unintentional) misuse of the technology. • If there are negative societal impacts, the authors could also discuss possible mitigation strategies (e.g., gated release of models, providing defenses in addition to attacks, mechanisms for monitoring misuse, mechanisms to monitor how a system learns from feedback over time, improving the efficiency and accessibility of ML). 11. Safeguards Question: Does the paper describe safeguards that have been put in place for responsible release of data or models that have a high risk for misuse (e.g., pre-trained language models, image generators, or scraped datasets)? Answer: [N/A] 68

Justification: The paper does not release pre-trained models, datasets, or any artifacts that pose a risk for misuse. The method is evaluated on standard simulation benchmarks for multi-agent systems. Guidelines: • The answer [N/A] means that the paper poses no such risks. • Released models that have a high risk for misuse or dual-use should be released with necessary safeguards to allow for controlled use of the model, for example by requiring that users adhere to usage guidelines or restrictions to access the model or implementing safety filters. • Datasets that have been scraped from the Internet could pose safety risks. The authors should describe how they avoided releasing unsafe images. • We recognize that providing effective safeguards is challenging, and many papers do not require this, but we encourage authors to take this into account and make a best faith effort. 12. Licenses for existing assets Question: Are the creators or original owners of assets (e.g., code, data, models), used in the paper, properly credited and are the license and terms of use explicitly mentioned and properly respected? Answer: [Yes] Justification: All benchmarks, baselines, and supporting assets are properly cited at first use: Ising Model and MFQ as behavior policy [Yang et al., 2018]; Battle on the MAgent platform [Zheng et al., 2018]; Gaussian Squeeze [Gu et al., 2021, HolmesParker et al., 2014]. Baselines: Joint and Independent Diffuser [Janner et al., 2022], Decision Diffuser [Ajay et al., 2023], MADiff [Zhu et al., 2024] (attention-coupled multi-agent diffusion), DoF [Liu et al., 2025] (Individual-Global-identically-Distributed factorized diffusion), MF-CDMs [Park et al., 2024] (adapted to offline RL as the MF-CDMs-RL baseline), OMAR [Pan et al., 2022], MA-TD3+BC [Fujimoto and Gu, 2021], Oryx [Li et al., 2025] (retention-based long-context sequence model with sequential ICQ), and CQL [Kumar et al., 2020] for offline regularization. All listed assets are publicly available research benchmarks and methods. Guidelines: • The answer [N/A] means that the paper does not use existing assets. • The authors should cite the original paper that produced the code package or dataset. • The authors should state which version of the asset is used and, if possible, include a URL. • The name of the license (e.g., CC-BY 4.0) should be included for each asset. • For scraped data from a particular source (e.g., website), the copyright and terms of service of that source should be provided. • If assets are released, the license, copyright information, and terms of use in the package should be provided. For popular datasets, paperswithcode.com/datasets has curated licenses for some datasets. Their licensing guide can help determine the license of a dataset. • For existing datasets that are re-packaged, both the original license and the license of the derived asset (if it has changed) should be provided. • If this information is not available online, the authors are encouraged to reach out to the asset’s creators. 13. New assets Question: Are new assets introduced in the paper well documented and is the documentation provided alongside the assets? Answer: [N/A] Justification: The paper does not release new datasets, models, or code as assets. The contribution is an algorithmic framework evaluated on existing benchmarks. Guidelines: • The answer [N/A] means that the paper does not release new assets. 69

• Researchers should communicate the details of the dataset/code/model as part of their submissions via structured templates. This includes details about training, license, limitations, etc. • The paper should discuss whether and how consent was obtained from people whose asset is used. • At submission time, remember to anonymize your assets (if applicable). You can either create an anonymized URL or include an anonymized zip file. 14. Crowdsourcing and research with human subjects Question: For crowdsourcing experiments and research with human subjects, does the paper include the full text of instructions given to participants and screenshots, if applicable, as well as details about compensation (if any)? Answer: [N/A] Justification: The paper does not involve crowdsourcing or research with human subjects. All experiments are conducted in simulation environments. Guidelines: • The answer [N/A] means that the paper does not involve crowdsourcing nor research with human subjects. • Including this information in the supplemental material is fine, but if the main contribution of the paper involves human subjects, then as much detail as possible should be included in the main paper. • According to the NeurIPS Code of Ethics, workers involved in data collection, curation, or other labor should be paid at least the minimum wage in the country of the data collector. 15. Institutional review board (IRB) approvals or equivalent for research with human subjects Question: Does the paper describe potential risks incurred by study participants, whether such risks were disclosed to the subjects, and whether Institutional Review Board (IRB) approvals (or an equivalent approval/review based on the requirements of your country or institution) were obtained? Answer: [N/A] Justification: The paper does not involve human subjects research. All experiments use simulated multi-agent environments. Guidelines: • The answer [N/A] means that the paper does not involve crowdsourcing nor research with human subjects. • Depending on the country in which research is conducted, IRB approval (or equivalent) may be required for any human subjects research. If you obtained IRB approval, you should clearly state this in the paper. • We recognize that the procedures for this may vary significantly between institutions and locations, and we expect authors to adhere to the NeurIPS Code of Ethics and the guidelines for their institution. • For initial submissions, do not include any information that would break anonymity (if applicable), such as the institution conducting the review. 16. Declaration of LLM usage Question: Does the paper describe the usage of LLMs if it is an important, original, or non-standard component of the core methods in this research? Note that if the LLM is used only for writing, editing, or formatting purposes and does not impact the core methodology, scientific rigor, or originality of the research, declaration is not required. Answer: [N/A] Justification: The core methodology does not involve LLMs as any important, original, or non-standard component. The method is based on diffusion models and mean-field theory for offline multi-agent reinforcement learning. 70

Guidelines: • The answer [N/A] means that the core method development in this research does not involve LLMs as any important, original, or non-standard components. • Please refer to our LLM policy in the NeurIPS handbook for what should or should not be described.

71