Proof of Weighted Second-Moment Evolution of the State Fluctuation Process

Write $\mathbf{1}=\mathbf{1}_{\Omega_0}$ (equal to $1$ on the regular event $\Omega_0$ and $0$ off it). Since $\Omega_0$ has probability $1$, expectations are unchanged when integrands are modified off $\Omega_0$, and we use this silently. All uses of the Tonelli and Fubini theorems are on the product of $[0,T]$ (trace Borel $\sigma$-algebra, restricted Lebesgue measure, total mass $T$ by the toolkit) with the probability space $(\Omega,\mathcal{F},P)$, both finite, hence $\sigma$-finite.

Step 0 (measurability and bounds). Every $\Sigma_t$ lies in the probability simplex (each agent occupies exactly one state, by the derived notation of the solution definition), so $|\mathfrak{s}_t|\le2\sqrt{N}$ everywhere, any two points of the simplex having Euclidean norm at most $1$. By the joint measurability lemma and measurability of sequentially continuous functions of measurable maps, the maps $\mathbf{1}\mathfrak{s}^\gamma_t$, $\mathbf{1}b^\gamma(\Sigma_t,\alpha_t)$, and $\mathbf{1}\Theta^{\gamma\delta}(\Sigma_t,\alpha_t)$ are product-measurable: for the latter two, replace $(\Sigma_t,\alpha_t)$ off $\Omega_0$ by a fixed point of $\Delta^l\times\mathcal{A}$ (which is nonempty) and compose the componentwise-measurable modified map with the sequentially continuous functions $b^\gamma$ and $\Theta^{\gamma\delta}$, whose sequential continuity follows from the joint continuity clause of the transition-rate family and continuity of the coordinate factors in their defining formulas; for $\mathbf{1}\mathfrak{s}^\gamma_t$, subtract the product-measurable $(t,\omega)\mapsto\mathbf{1}S^\gamma_t$ (continuous in $t$, via the composition lemma applied to $(t,\omega)\mapsto t$). On $\Omega_0$, $|b^\gamma|\le2(l-1)B$ and $|\Theta^{\gamma\delta}|\le2(l-1)B$ by part (a) of the martingale decomposition theorem, so $|g^\gamma_s|\le4\sqrt{N}(l-1)B$ there. Since each $\dot{z}^{\gamma\delta}$ is continuous on $[0,T]$, so is each $Z^{\gamma\delta}$ (claim 4 of the componentwise toolkit, the indefinite Riemann integral of a continuous integrand being continuous), and both are bounded there by the extreme value theorem, every expectation named in part (a) of the statement is finite and bounded in $s$, and its measurability in $s$ follows from the Fubini theorem, the ($\mathbf{1}$-modified) integrands being bounded and product-measurable on a finite product measure. This proves (a).

Step 1 (integral representation and cross moments). By part (b) of the martingale decomposition theorem and condition 2 of the mean-field trajectory pair (with claim 3 of the integral toolkit for the Riemann-Lebesgue agreement; the indicator $\mathbf{1}_{\Omega_0}$ in the decomposition's integral equals $1$ on $\Omega_0$), almost surely, for all $t$ and $\gamma$,

$$\mathfrak{s}^\gamma_t=\mathfrak{s}^\gamma_0+F^\gamma_t+\mathfrak{m}^\gamma_t,\qquad F^\gamma_t=\int_{[0,t]}g^\gamma_s\,ds,\qquad \mathfrak{m}^\gamma_t=\sqrt{N}\,M^\gamma_t,$$

where each $M^\gamma$ is a square-integrable martingale with $M^\gamma_0=0$. We record four facts, for all $\gamma,\delta$ and $0\le s\le t\le T$.

(1a) If $X$ is a bounded random variable that is measurable with respect to the system filtration entry $\mathcal{F}^{\mathrm{sys}}_s$, then $\mathbb{E}[X\,(\mathfrak{m}^\delta_t-\mathfrak{m}^\delta_s)]=0$. Indeed, with $|X|\le c$, the dyadic truncations $X_n=2^{-n}\lfloor2^nX\rfloor$ (a finite sum $\sum_k k2^{-n}\mathbf{1}_{D_{n,k}}$ over the finitely many levels with $|k|2^{-n}\le c+1$, each $D_{n,k}\in\mathcal{F}^{\mathrm{sys}}_s$) satisfy $|X_n-X|\le2^{-n}$; the defining property of the square-integrable martingale $M^\delta$ gives $\mathbb{E}[\mathbf{1}_{D_{n,k}}(M^\delta_t-M^\delta_s)]=0$ for each level set, hence $\mathbb{E}[X_n(\mathfrak{m}^\delta_t-\mathfrak{m}^\delta_s)]=0$ by linearity, and $|\mathbb{E}[(X-X_n)(\mathfrak{m}^\delta_t-\mathfrak{m}^\delta_s)]|\le2^{-n}(\mathbb{E}|\mathfrak{m}^\delta_t|+\mathbb{E}|\mathfrak{m}^\delta_s|)\to0$, using that square-integrable variables are integrable (Cauchy-Schwarz for the mean-square norm against the constant $1$).

(1b) $\mathbb{E}[\mathfrak{s}^\gamma_0\,\mathfrak{m}^\delta_t]=0$: apply (1a) with $s=0$, $X=\mathfrak{s}^\gamma_0$ (bounded by $2\sqrt{N}$; $\Sigma_0$ is $\mathcal{F}^{\mathrm{sys}}_0$-measurable by part (iv) of the existence theorem) and $\mathfrak{m}^\delta_0=0$.

(1c) $\mathbb{E}[g^\gamma_s\,\mathfrak{m}^\delta_t]=\mathbb{E}[g^\gamma_s\,\mathfrak{m}^\delta_s]$: apply (1a) with $X=g^\gamma_s$, which is bounded and $\mathcal{F}^{\mathrm{sys}}_s$-measurable ($\Sigma_s$ and $\alpha_s$ are adapted by part (iv) of the existence theorem, and $b^\gamma$ composed with them is measurable by the composition lemma; $b^\gamma(S_s,A_s)$ is a constant).

(1d) $\mathbb{E}[\mathfrak{m}^\gamma_t\,\mathfrak{m}^\delta_t]=\int_{[0,t]}\mathbb{E}[\Theta^{\gamma\delta}(\Sigma_s,\alpha_s)]\,ds$: this is the covariation identity, part (c) of the decomposition theorem, with $r=0$, $D=\Omega$, multiplied by $N$ (its time integrand carrying the factor $\mathbf{1}_{\Omega_0}$, which changes nothing under the expectation, $\Omega_0$ having probability $1$), together with the Fubini theorem to exchange $\mathbb{E}$ and the time integral of the (after $\mathbf{1}$-modification) product-measurable bounded integrand.

Step 2 (evolution of the second-moment matrix). Fix $\gamma,\delta$ and set $\Psi^{\gamma\delta}(t)=\mathbb{E}[\mathfrak{s}^\gamma_t\mathfrak{s}^\delta_t]$. Expanding the product of the two three-term representations of Step 1 and taking expectations termwise (all nine terms are integrable, every factor being bounded except the martingales, which are square-integrable):

$$\Psi^{\gamma\delta}(t)=\mathbb{E}[\mathfrak{s}^\gamma_0\mathfrak{s}^\delta_0]+\mathbb{E}[\mathfrak{s}^\gamma_0F^\delta_t]+\mathbb{E}[F^\gamma_t\mathfrak{s}^\delta_0]+\mathbb{E}[F^\gamma_tF^\delta_t]+\mathbb{E}[F^\gamma_t\mathfrak{m}^\delta_t]+\mathbb{E}[\mathfrak{m}^\gamma_tF^\delta_t]+\mathbb{E}[\mathfrak{m}^\gamma_t\mathfrak{m}^\delta_t],$$

the terms $\mathbb{E}[\mathfrak{s}^\gamma_0\mathfrak{m}^\delta_t]$ and $\mathbb{E}[\mathfrak{m}^\gamma_t\mathfrak{s}^\delta_0]$ vanishing by (1b). Now: $\mathbb{E}[\mathfrak{s}^\gamma_0F^\delta_t]=\int_{[0,t]}\mathbb{E}[\mathfrak{s}^\gamma_0g^\delta_s]ds$ by the Fubini theorem (bounded integrand). Pathwise, the integration by parts lemma on $[0,t]$ with $u_0=v_0=0$ (the case $t=0$ being trivial, both sides vanishing) gives $F^\gamma_tF^\delta_t=\int_{[0,t]}(g^\gamma_sF^\delta_s+F^\gamma_sg^\delta_s)ds$ almost surely, so $\mathbb{E}[F^\gamma_tF^\delta_t]=\int_{[0,t]}\mathbb{E}[g^\gamma_sF^\delta_s+F^\gamma_sg^\delta_s]ds$ (Fubini; $|F^\gamma_s|\le4\sqrt{N}(l-1)BT$ on $\Omega_0$). Also $F^\gamma_t\mathfrak{m}^\delta_t=\int_{[0,t]}g^\gamma_s\,ds\cdot \mathfrak{m}^\delta_t=\int_{[0,t]}g^\gamma_s \mathfrak{m}^\delta_t\,ds$ pathwise, so by the Fubini theorem (the integrand is dominated by $4\sqrt{N}(l-1)B\,|\mathfrak{m}^\delta_t|$, which is integrable on the product) and (1c),

$$\mathbb{E}[F^\gamma_t\mathfrak{m}^\delta_t]=\int_{[0,t]}\mathbb{E}[g^\gamma_s\mathfrak{m}^\delta_t]\,ds=\int_{[0,t]}\mathbb{E}[g^\gamma_s\mathfrak{m}^\delta_s]\,ds,$$

and symmetrically for $\mathbb{E}[\mathfrak{m}^\gamma_tF^\delta_t]$. Combining with (1d) and collecting the integrands via $\mathfrak{s}^\delta_s=\mathfrak{s}^\delta_0+F^\delta_s+\mathfrak{m}^\delta_s$ (almost surely):

$$\Psi^{\gamma\delta}(t)=\Psi^{\gamma\delta}(0)+\int_{[0,t]}\psi^{\gamma\delta}(s)\,ds,\qquad \psi^{\gamma\delta}(s)=\mathbb{E}[g^\gamma_s\mathfrak{s}^\delta_s]+\mathbb{E}[g^\delta_s\mathfrak{s}^\gamma_s]+\mathbb{E}[\Theta^{\gamma\delta}(\Sigma_s,\alpha_s)],$$

with each $\psi^{\gamma\delta}$ bounded, and measurable by the same Fubini argument as in Step 0 applied to the bounded product-measurable integrands $\mathbf{1}g^\gamma\mathfrak{s}^\delta$ and $\mathbf{1}\Theta^{\gamma\delta}$.

Step 3 (weighting by $Z$). Fix $t\in[0,T]$; the case $t=0$ is trivial, so let $t>0$. Apply the integration by parts lemma on $[0,t]$ to $u=Z^{\gamma\delta}$ (with density $\dot{z}^{\gamma\delta}$, continuous, Riemann and Lebesgue integrals agreeing) and $v=\Psi^{\gamma\delta}$ (with density $\psi^{\gamma\delta}$):

$$Z^{\gamma\delta}_t\Psi^{\gamma\delta}(t)=Z^{\gamma\delta}_0\Psi^{\gamma\delta}(0)+\int_{[0,t]}\big(\dot{z}^{\gamma\delta}(s)\Psi^{\gamma\delta}(s)+Z^{\gamma\delta}_s\psi^{\gamma\delta}(s)\big)ds .$$

Sum over $\gamma,\delta\in\{1,\dots,l\}$. On the left, $\sum_{\gamma,\delta}Z^{\gamma\delta}_t\Psi^{\gamma\delta}(t)=\mathbb{E}[\mathfrak{s}_t\cdot Z_t\mathfrak{s}_t]$ by linearity of the expectation, and likewise at $0$. In the integrand, $\sum_{\gamma,\delta}\dot{z}^{\gamma\delta}(s)\Psi^{\gamma\delta}(s)=\mathbb{E}[\mathfrak{s}_s\cdot\dot{z}(s)\mathfrak{s}_s]$, while by the symmetry of $Z_s$ and relabeling of the summation indices,

$$\sum_{\gamma,\delta}Z^{\gamma\delta}_s\big(\mathbb{E}[g^\gamma_s\mathfrak{s}^\delta_s]+\mathbb{E}[g^\delta_s\mathfrak{s}^\gamma_s]\big)=2\sum_{\gamma,\delta}Z^{\gamma\delta}_s\,\mathbb{E}[\mathfrak{s}^\gamma_sg^\delta_s]=2\,\mathbb{E}[\mathfrak{s}_s\cdot Z_sg_s].$$

This yields exactly the displayed identity of part (b).