Proof of 4. Let \(V^\pi(s,\mu)\) be the value function of the coupled representative-population problem in [37], associated with the same fixed policy \(\pi\). By [37] Theorem 3.3, under 2, the coupled value function satisfies \(V^\pi\in \mathcal{C}^{2,2}_{\mathrm{poly}} (\mathbb{R}^d\times\mathcal{P}_2(\mathbb{R}^d))\) (the space of functions with polynomial growth that are twice differentiable in \(s\) and \(\mu\); see [37]) and solves \[\label{eq:coupled-HJB-used-mkv} (\mathscr L_{b,\Sigma}^{\pi}-\beta)V^\pi(s,\mu) +\ell_\lambda^\pi(s,\mu)=0,\tag{35}\] where \(\mathscr L_{b,\Sigma}^{\pi}\) being the generator on \(\mathbb{R}^d\times\mathcal{P}_2(\mathbb{R}^d)\) defined as \[\label{eq:L-pi} \begin{align} (\mathscr{L}_{b,\Sigma}^{\pi}F)(s,\mu) &:= b^{\pi}(s,\mu)\cdot\nabla_s F(s,\mu) + \tfrac{1}{2}\,\Sigma^{\pi}(s,\mu) : \nabla^2_{ss} F(s,\mu)\\ &\quad + \int_{\mathbb{R}^d} \Bigl[ b^{\pi}(\xi,\mu)\cdot \partial_{\mu}F(s,\mu)(\xi) + \tfrac{1}{2}\,\Sigma^{\pi}(\xi,\mu) : D_\xi\partial_{\mu}F(s,\mu)(\xi) \Bigr]\,\mu(d\xi), \end{align}\tag{36}\] On the other hand, since the initial condition in 5 satisfies \(s_0\sim\mu\), conditioning on \(s_0=s\), using Fubini’s theorem, we have \[\begin{align} \label{eq:mkv-value-as-integrated-coupled-value} \begin{aligned} \mathcal{V}^\pi(\mu) &= \mathbb{E}\left[ \int_0^\infty e^{-\beta t} \ell_\lambda^\pi(s_t,\mu_t)\,dt \,\bigg|\,\mu_0=\mu \right] = \int_{\mathbb{R}^d} \mathbb{E}\left[ \int_0^\infty e^{-\beta t} \ell_\lambda^\pi(s_t,\mu_t)\,dt \,\bigg|\,s_0=s,\;\mu_0=\mu \right]\mu(ds) \\ &= \int_{\mathbb{R}^d}V^\pi(s,\mu)\,\mu(ds). \end{aligned} \end{align}\tag{37}\] Since \(V^\pi\in\mathcal{C}^{2,2}_{\mathrm{poly}}\), applying the Lions derivative formula for maps of the form \(\mu\mapsto\int_{\mathbb{R}^d} V^\pi(s,\mu)\,\mu(ds)\) [7] yields \[\label{eq:mkv-value-lions-derivative} \partial_\mu\mathcal{V}^\pi(\mu)(x) = \nabla_sV^\pi(x,\mu) + \int_{\mathbb{R}^d} \partial_\mu V^\pi(s,\mu)(x)\,\mu(ds),\tag{38}\] and \[\label{eq:mkv-value-lions-second-derivative} D_x\partial_\mu\mathcal{V}^\pi(\mu)(x) = \nabla^2_{ss}V^\pi(x,\mu) + \int_{\mathbb{R}^d} D_x\partial_\mu V^\pi(s,\mu)(x)\,\mu(ds).\tag{39}\] The polynomial bounds in the definition of \(\mathcal{C}^2_{\mathrm{poly}}(\mathcal{P}_2(\mathbb{R}^d))\) follow directly from the corresponding bounds for \(V^\pi\). Hence \(\mathcal{V}^\pi\in \mathcal{C}^2_{\mathrm{poly}}(\mathcal{P}_2(\mathbb{R}^d))\). Using 38 and 39 , we compute \[\begin{align} (\mathcal{L}_{b,\Sigma}^{\pi}\mathcal{V}^\pi)(\mu) &= \int_{\mathbb{R}^d} \left[ b^\pi(x,\mu)\cdot\nabla_sV^\pi(x,\mu) + \frac{1}{2}\Sigma^\pi(x,\mu):\nabla^2_{ss}V^\pi(x,\mu) \right]\mu(dx) \\ &\quad+ \int_{\mathbb{R}^d}\int_{\mathbb{R}^d} \left[ b^\pi(x,\mu)\cdot\partial_\mu V^\pi(s,\mu)(x) + \frac{1}{2}\Sigma^\pi(x,\mu): D_x\partial_\mu V^\pi(s,\mu)(x) \right]\mu(dx)\mu(ds) \\ &= \int_{\mathbb{R}^d} (\mathscr L_{b,\Sigma}^{\pi}V^\pi)(s,\mu)\,\mu(ds). \end{align}\] Integrating 35 with respect to \(\mu(ds)\) and using 37 gives \[(\mathcal{L}_{b,\Sigma}^{\pi}\mathcal{V}^\pi)(\mu) - \beta\mathcal{V}^\pi(\mu) + r_\lambda^\pi(\mu) = 0.\] Thus \(\mathcal{V}^\pi\) solves ?? .

We now prove uniqueness. Let \(U\in\mathcal{C}^2_{\mathrm{poly}}(\mathcal{P}_2(\mathbb{R}^d))\) be another classical solution of ?? . Applying the Itô-Lions formula to \(e^{-\beta t}U(\mu_t)\) gives \[\frac{d}{dt} \left(e^{-\beta t}U(\mu_t)\right) = e^{-\beta t} \left[ (\mathcal{L}_{b,\Sigma}^{\pi}U)(\mu_t) - \beta U(\mu_t) \right] = -e^{-\beta t}r_\lambda^\pi(\mu_t).\] Therefore, for every \(T>0\), \[U(\mu) = e^{-\beta T}U(\mu_T) + \int_0^T e^{-\beta t}r_\lambda^\pi(\mu_t)\,dt.\] By the polynomial growth of \(U\), there exists \(C_U>0\) such that \[|U(\nu)|\le C_U\bigl(1+m_2(\nu)\bigr), \qquad \nu\in\mathcal{P}_2(\mathbb{R}^d).\] Moreover, by [37], there exists \(C>0\) such that \[m_2(\mu_T)\le C e^{\beta_0(2)T}\bigl(1+m_2(\mu)\bigr), \qquad T\ge0.\] Since \(\beta>\bar\beta_\pi\) and, by [37], \(\bar\beta_\pi\ge \beta_0(2)\), we obtain \[e^{-\beta T}|U(\mu_T)| \le C e^{-\beta T} + C e^{-(\beta-\beta_0(2))T}\bigl(1+m_2(\mu)\bigr) \longrightarrow 0, \qquad T\to\infty.\] Hence \[U(\mu) = \int_0^\infty e^{-\beta t}r_\lambda^\pi(\mu_t)\,dt = \mathcal{V}^\pi(\mu).\] It remains to prove ?? . In [37], the finite-horizon derivatives of the coupled value functions are represented through the variational flows. The estimates [37] therein imply, after integration against \(e^{-\beta t}\), that \[\begin{align} |\nabla_sV^\pi(s,\nu)| + |\partial_\mu V^\pi(s,\nu)(x)| &\le \frac{C_\pi}{\beta-\bar\beta_\pi} \bigl(1+|s|+|x|+m_2(\nu)\bigr), \\ \|\nabla^2_{ss}V^\pi(s,\nu)\| + \|D_x\partial_\mu V^\pi(s,\nu)(x)\| &\le \frac{C_\pi}{\beta-\bar\beta_\pi} \bigl(1+m_2(\nu)\bigr). \end{align}\] The passage from the finite-horizon bounds to the infinite-horizon bounds follows from the locally uniform convergence of the derivatives proved in [37] Step 4 of the proof of Theorem 3.3.

Combining these estimates with 38 gives \[\begin{align} |\partial_\mu\mathcal{V}^\pi(\nu)(x)| &\le |\nabla_sV^\pi(x,\nu)| + \int_{\mathbb{R}^d} |\partial_\mu V^\pi(s,\nu)(x)|\,\nu(ds)\le \frac{C_\pi}{\beta-\bar\beta_\pi} \bigl(1+|x|^2+m_2(\nu)\bigr). \end{align}\] Similarly, using 39 , \[\begin{align} \|D_x\partial_\mu\mathcal{V}^\pi(\nu)(x)\| &\le \|\nabla^2_{ss}V^\pi(x,\nu)\| + \int_{\mathbb{R}^d} \|D_x\partial_\mu V^\pi(s,\nu)(x)\|\,\nu(ds) \le \frac{C_\pi}{\beta-\bar\beta_\pi} \bigl(1+|x|^2+m_2(\nu)\bigr). \end{align}\] This proves ?? . ◻

Proof of 8. Let \(V^{\pi^\varepsilon}\) denote the value function of the coupled state-law problem studied in [37]. Set \[J(\varepsilon) := \int_{\mathbb{R}^d} V^{\pi^\varepsilon}(s,\mu)\,\mu(ds).\] By 37 , \(J(\varepsilon)=\mathcal{V}^{\pi^\varepsilon}(\mu)=\mathcal{J}(\varepsilon)\). It is therefore enough to compute the derivative of \(J\) at \(\varepsilon=0\).

We apply [37] with initial population law \(\mu\) and initial representative law also equal to \(\mu_0=\mu\). Denote by \(\widehat\rho^\pi\) the corresponding discounted occupancy measure of the coupled state-law process. Then \(J\) is differentiable at \(\varepsilon=0\), and \[\begin{align} \frac{d}{d\varepsilon}J(\varepsilon)\bigg|_{\varepsilon=0} &=\frac{1}{\beta} \mathbb{E}_{(s,\eta)\sim\widehat\rho^\pi}\!\left[ \int_{\mathcal{A}} q_{\mathrm{rep}}^\pi(s,\eta,a)\, \varphi(da\mid s,\eta) + \int_{\mathbb{R}^d}\!\int_{\mathcal{A}} q_{\mathrm{pop}}^\pi(s,\eta,\xi,a)\, \varphi(da\mid \xi,\eta)\,\eta(d\xi) \right], \end{align}\] where, according to [37], \[\begin{align} q_{\mathrm{rep}}^\pi(s,\eta,a) &= r(s,\eta,a)-\lambda\log p_\pi(a\mid s,\eta) +b(s,\eta,a)\cdot\nabla_s V^\pi(s,\eta) \\ &\qquad\qquad\qquad+\frac{1}{2}\,\Sigma(s,\eta,a):D^2_{ss}V^\pi(s,\eta)-\beta V^\pi(s,\eta), \\ q_{\mathrm{pop}}^\pi(s,\eta,\xi,a) &= b(\xi,\eta,a)\cdot\partial_\mu V^\pi(s,\eta)(\xi) +\frac{1}{2}\,\Sigma(\xi,\eta,a):D_\xi\partial_\mu V^\pi(s,\eta)(\xi). \end{align}\]

Since the initial representative law is \(\mu\), the representative component in the coupled state-law system has law \(\mu_t^\pi\) at time \(t\). Hence, for every measurable \(F\) for which the following quantities are well defined, \[\begin{align} \label{eq:rho-transfer-mkv} \mathbb{E}_{(s,\eta)\sim\widehat\rho^\pi}[F(\eta,s)] = \int_0^\infty e^{-\beta t} \int_{\mathbb{R}^d} F(\mu_t^\pi,\xi)\,\mu_t^\pi(d\xi)\,dt = \mathbb{E}_{(\eta,\xi)\sim\rho^\pi}[F(\eta,\xi)]. \end{align}\tag{40}\] Applying this identity with \[F(\eta,\xi) := \int_{\mathcal{A}} \left[ q_{\mathrm{rep}}^\pi(\xi,\eta,a) + \int_{\mathbb{R}^d} q_{\mathrm{pop}}^\pi(s,\eta,\xi,a)\,\eta(ds) \right] \varphi(da\mid \xi,\eta),\] we obtain \[\begin{align} \label{eq:derivative95j95almost95finished} \begin{aligned} \frac{d}{d\varepsilon}\mathcal{J}(\varepsilon)\bigg|_{\varepsilon=0} &=\frac{1}{\beta} \mathbb{E}_{(\eta,\xi)\sim\rho^\pi}\!\left[ \int_{\mathcal{A}} \left[ q_{\mathrm{rep}}^\pi(\xi,\eta,a) + \int_{\mathbb{R}^d} q_{\mathrm{pop}}^\pi(s,\eta,\xi,a)\,\eta(ds) \right] \varphi(da\mid \xi,\eta) \right]. \end{aligned} \end{align}\tag{41}\] It remains to identify the term in brackets. By 38 , 39 , and [37], the representative and population kernels satisfy \[\begin{align} \label{eq:realtion95q95paper95195neq95q} q_{\mathrm{rep}}^\pi(\xi,\eta,a) + \int_{\mathbb{R}^d} q_{\mathrm{pop}}^\pi(s,\eta,\xi,a)\,\eta(ds) = q^\pi(\eta,\xi,a), \end{align}\tag{42}\] where \(q^\pi\) is the McKean-Vlasov advantage function defined in ?? . Substituting 42 into 41 gives ?? . ◻

Proof of 13. Fix \(\omega\in\mathbb{R}^p\) and set \(\pi:=\pi_\omega\). Define \(J(\omega):=\int_{\mathbb{R}^d}V^{\pi_\omega}(s,\mu)\,\mu(ds).\) As in the proof of 8, one has \(J(\omega)=\mathcal{V}^{\pi_\omega}(\mu)=\mathcal{J}(\omega)\), so it is enough to compute the gradient of \(J\).

By the parametric policy gradient theorem for \(J\) in the first paper, \(J\) is differentiable at \(\omega\) and \[\begin{align} \nabla_\omega J(\omega) &= \frac{1}{\beta}\mathbb{E}_{(s,\eta)\sim\widehat\rho^\pi}\!\left[ \int_{\mathcal{A}} q_{\mathrm{rep}}^\pi(s,\eta,a)\, \nabla_\omega\log p_\pi(a\mid s,\eta)\, \pi(da\mid s,\eta) \right]\\ &\qquad+ \frac{1}{\beta}\mathbb{E}_{(s,\eta)\sim\widehat\rho^\pi}\!\left[ \int_{\mathbb{R}^d}\!\int_{\mathcal{A}} q_{\mathrm{pop}}^\pi(s,\eta,\xi,a)\, \nabla_\omega\log p_\pi(a\mid \xi,\eta)\, \pi(da\mid \xi,\eta)\, \eta(d\xi) \right]. \end{align}\] Now repeat verbatim the reduction carried out in the proof of 8: use the transfer identity 40 to pass from \(\widehat\rho^\pi\) to \(\rho^\pi\). Then use the reduced expression 41 and the relation between the representative and population kernels with the single advantage function introduced in 42 to conclude ?? . The equivalent expectation form with \(a\sim\pi(\cdot\mid\xi,\eta)\) follows immediately from the definition of conditional expectation. ◻

Proof. Since \((\mu_t)_{t\in[0,\Delta t]}\) is deterministic, the Itô-Lions formula applied to \(F(s_t,\mu_t)\) gives \[F(s_t,\mu_t) = F(s,\mu) + \int_0^t (\mathscr L_{b,\Sigma}^aF)(s_\tau,\mu_\tau)\,d\tau + M_t^F,\] with \[M_t^F = \int_0^t \nabla_xF(s_\tau,\mu_\tau) \sigma(s_\tau,\mu_\tau,a)\,dB_\tau .\] The polynomial growth of \(\nabla_xF\), the boundedness of the coefficients in 16, and the moment estimates for the frozen-action dynamics imply that \(M^F\) is a square-integrable martingale. Hence \[\mathbb{E}[F(s_t,\mu_t)]-F(s,\mu) = \int_0^t \mathbb{E}[ (\mathscr L_{b,\Sigma}^aF)(s_\tau,\mu_\tau) ]\,d\tau .\] Applying the same argument to \(\mathscr L_{b,\Sigma}^aF\) yields, for \(\tau\in[0,\Delta t]\), \[\mathbb{E}[ (\mathscr L_{b,\Sigma}^aF)(s_\tau,\mu_\tau) ] = (\mathscr L_{b,\Sigma}^aF)(s,\mu) + \int_0^\tau \mathbb{E}[ (\mathscr L_{b,\Sigma}^a)^2F(s_u,\mu_u) ]\,du .\] Combining the two identities gives \[\mathbb{E}[F(s_t,\mu_t)]-F(s,\mu) = t(\mathscr L_{b,\Sigma}^aF)(s,\mu) + \int_0^t\int_0^\tau \mathbb{E}[ (\mathscr L_{b,\Sigma}^a)^2F(s_u,\mu_u) ]\,du\,d\tau .\] Taking \(t=\Delta t\) and dividing by \(\Delta t\) proves ?? , with \[R_F^{a,\Delta t}(s,\mu) := \frac{1}{\Delta t} \int_0^{\Delta t}\int_0^\tau \mathbb{E}[ (\mathscr L_{b,\Sigma}^a)^2F(s_u,\mu_u) ]\,du\,d\tau .\] The estimate ?? follows from \[\frac{1}{\Delta t}\int_0^{\Delta t}\int_0^\tau du\,d\tau = \frac{\Delta t}{2}.\]

For \(F_k(x,\nu):=x_k\), \[\mathbb{E}[F_k(s_{\Delta t},\mu_{\Delta t})]-F_k(s,\mu) = \mathbb{E}[(s_{\Delta t})_k-s_k].\] Thus the first identity in ?? follows from the definition of \(\hat{b}_k\). Since \(\nabla_xF_k=e_k\), \(D_{xx}^2F_k=0\), and \(\partial_\nu F_k=0\), we also have \[(\mathscr L_{b,\Sigma}^aF_k)(s,\mu) = b(s,\mu,a)\cdot e_k = b_k(s,\mu,a).\]

For \(G_{ij}^s(x,\nu):=(x_i-s_i)(x_j-s_j)\), \(G_{ij}^s(s,\mu)=0\) and \[\mathbb{E}[G_{ij}^s(s_{\Delta t},\mu_{\Delta t})] = \mathbb{E}[ ((s_{\Delta t})_i-s_i)((s_{\Delta t})_j-s_j) ].\] Hence the first identity in ?? follows from the definition of \(\hat{\Sigma}_{ij}\). Moreover, \[\nabla_xG_{ij}^s(x,\nu) = e_i(x_j-s_j)+e_j(x_i-s_i), \qquad D_{xx}^2G_{ij}^s(x,\nu) = e_ie_j^\top+e_je_i^\top, \qquad \partial_\nu G_{ij}^s=0.\] Therefore \[(\mathscr L_{b,\Sigma}^aG_{ij}^s)(x,\nu) = b_i(x,\nu,a)(x_j-s_j) + b_j(x,\nu,a)(x_i-s_i) + \Sigma_{ij}(x,\nu,a).\] Evaluating at \((x,\nu)=(s,\mu)\) gives the second identity in ?? .

The final identities follow by applying ?? to \(F_k\) and \(G_{ij}^s\), and then using ?? and ?? . ◻

Proof. Step 1: bounds for the short-time operator \(P_u^a\). For \(u\in[0,1]\), define \[P_u^a\Phi(s,\mu):=\mathbb{E}[\Phi(s_u,\mu_u)],\] where \((s_u,\mu_u)\) is the frozen-action dynamics introduced in 1. We first prove that \[\label{eq:semigroup-uniform-bound-mkv} \left\vert\kern-0.25ex\left\vert\kern-0.25ex\left\vert P_u^a\Phi(s,\mu) \right\vert\kern-0.25ex\right\vert\kern-0.25ex\right\vert_{\xi}\le C C_\Phi, \qquad u\in[0,1],\tag{43}\] whenever \[\left\vert\kern-0.25ex\left\vert\kern-0.25ex\left\vert \Phi(x,\nu) \right\vert\kern-0.25ex\right\vert\kern-0.25ex\right\vert_{\eta}\le C_\Phi \qquad \text{for all }(x,\nu,\eta)\in\mathbb{R}^d\times\mathcal{P}_2(\mathbb{R}^d)\times\mathbb{R}^d.\] By [37], applied to the frozen coefficients \(b(\cdot,\cdot,a)\) and \(\sigma(\cdot,\cdot,a)\), the derivatives of the frozen-action flow with respect to the initial state and the initial law satisfy \[\sup_{u\in[0,1]}\mathbb{E}\!\left[ \|\nabla_s s_u\|^2+\|D_{ss}^2s_u\|^2+\|\partial_\mu s_u(\xi)\|^2+\|D_\xi\partial_\mu s_u(\xi)\|^2 \right]\le C,\] where \(C\) is independent of \(a\), \(u\), \(s\), \(\mu\), \(\xi\), and \(\Delta t\). Therefore, differentiating \(P_u^a\Phi(s,\mu)=\mathbb{E}[\Phi(s_u,\mu_u)]\) and using the chain rule gives \[\|\nabla_s P_u^a\Phi(s,\mu)\|+\|D_{ss}^2P_u^a\Phi(s,\mu)\|\le C C_\Phi.\] The same argument, using the Lions chain rule and the estimates for \(\partial_\mu s_u(\xi)\) and \(D_\xi\partial_\mu s_u(\xi)\) in [37], gives \[|\partial_\mu P_u^a\Phi(s,\mu)(\xi)|+\|D_\xi\partial_\mu P_u^a\Phi(s,\mu)(\xi)\|\le C C_\Phi.\] Together with \(|P_u^a\Phi(s,\mu)|\le C_\Phi\), this proves 43 . We also use the following linear-growth variant of 43 . If \[\left\vert\kern-0.25ex\left\vert\kern-0.25ex\left\vert \Phi_z(x,\nu) \right\vert\kern-0.25ex\right\vert\kern-0.25ex\right\vert_{\eta}\le C_\Phi(1+|x-z|) \qquad \text{for all }(x,\nu,\eta,z),\] then \[\label{eq:semigroup-linear-growth-bound-mkv} \left\vert\kern-0.25ex\left\vert\kern-0.25ex\left\vert P_u^a\Phi_z(s,\mu) \right\vert\kern-0.25ex\right\vert\kern-0.25ex\right\vert_{\xi} \le C C_\Phi(1+|s-z|+u^{1/2}), \qquad u\in[0,1].\tag{44}\] Indeed, the proof of 43 gives the same estimate with \(1+\mathbb{E}[|s_u-z|]\) on the right-hand side. Moreover, \[\mathbb{E}[|s_u-z|]\le |s-z|+\mathbb{E}[|s_u-s|]\le |s-z|+Cu^{1/2},\] where the last estimate follows from the boundedness of \(b\) and \(\sigma\) and the Burkholder-Davis-Gundy inequality. This proves 44 . Step 2: the remainder for \(F_k\). Let \(F_k(x,\nu)=x_k\) and set \[\Phi_k^a(x,\nu):=(\mathscr L_{b,\Sigma}^a)^2F_k(x,\nu).\] Since \(F_k\) is linear and the coefficients have bounded derivatives up to order four, uniformly in \(a\), applying \(\mathscr L_{b,\Sigma}^a\) twice gives \[\left\vert\kern-0.25ex\left\vert\kern-0.25ex\left\vert \Phi_k^a(x,\nu) \right\vert\kern-0.25ex\right\vert\kern-0.25ex\right\vert_{\xi}\le C.\] By Step 1, \[\left\vert\kern-0.25ex\left\vert\kern-0.25ex\left\vert P_u^a\Phi_k^a(s,\mu) \right\vert\kern-0.25ex\right\vert\kern-0.25ex\right\vert_{\xi}\le C, \qquad u\in[0,1].\] Using ?? , \[R_{F_k}^{a,\Delta t}(s,\mu) = \frac{1}{\Delta t} \int_0^{\Delta t}\int_0^\tau P_u^a\Phi_k^a(s,\mu)\,du\,d\tau .\] Therefore, \[\left\vert\kern-0.25ex\left\vert\kern-0.25ex\left\vert R_{F_k}^{a,\Delta t}(s,\mu) \right\vert\kern-0.25ex\right\vert\kern-0.25ex\right\vert_{\xi} \le \frac{1}{\Delta t} \int_0^{\Delta t}\int_0^\tau C\,du\,d\tau \le C\Delta t .\]

Step 3: the structure of the remainder for \(G_{ij}^s\). Fix first \(z\in\mathbb{R}^d\) and define \[G_{ij}^z(x,\nu):=(x_i-z_i)(x_j-z_j).\] A direct computation gives \[(\mathscr L_{b,\Sigma}^aG_{ij}^z)(x,\nu) = b_i(x,\nu,a)(x_j-z_j) + b_j(x,\nu,a)(x_i-z_i) + \Sigma_{ij}(x,\nu,a).\] This formula is the key point. The only terms that can grow in \(x\) or \(z\) are the factors \(x_j-z_j\) and \(x_i-z_i\). All coefficients and all their derivatives needed below are bounded by 16.

Set \[H_{ij}^{a,z}:=\mathscr L_{b,\Sigma}^aG_{ij}^z, \qquad \Phi_{ij}^{a,z}:=\mathscr L_{b,\Sigma}^aH_{ij}^{a,z}.\] From the displayed formula for \(H_{ij}^{a,z}\), and since the operator \(\mathscr L_{b,\Sigma}^a\) does not depend on \(z\), differentiating with respect to \(z\) gives, for every multi-index \(\ell\) with \(|\ell|\le2\), \[\partial_z^\ell\Phi_{ij}^{a,z} = \mathscr L_{b,\Sigma}^a(\partial_z^\ell H_{ij}^{a,z}).\] Now \(\partial_z^\ell H_{ij}^{a,z}\) is a sum of bounded coefficients multiplied either by \(1\) or by one of the factors \(x_r-z_r\). Therefore, using again the boundedness of the coefficients and of their derivatives, \[\left\vert\kern-0.25ex\left\vert\kern-0.25ex\left\vert \partial_z^\ell\Phi_{ij}^{a,z}(x,\nu) \right\vert\kern-0.25ex\right\vert\kern-0.25ex\right\vert_{\xi} \le C(1+|x-z|), \qquad |\ell|\le2.\] This is the desired bound: it follows from the explicit expression of \(\mathscr L_{b,\Sigma}^aG_{ij}^z\) and from the fact that applying \(\mathscr L_{b,\Sigma}^a\) adds only bounded coefficients and bounded derivatives. Step 4: integration of the bound for \(G_{ij}^z\). By 44 applied to \(\partial_z^\ell\Phi_{ij}^{a,z}\), and using the bound proved in Step 3, we obtain \[\label{eq:semigroup-bound-Gij-uniform-mkv} \sum_{|\ell|\le2} \left\vert\kern-0.25ex\left\vert\kern-0.25ex\left\vert P_u^a(\partial_z^\ell\Phi_{ij}^{a,z})(s,\mu) \right\vert\kern-0.25ex\right\vert\kern-0.25ex\right\vert_{\xi} \le C(1+|s-z|+u^{1/2}), \qquad u\in[0,1].\tag{45}\] Define \[R_{ij}^{a,\Delta t}(s,\mu,z) := \frac{1}{\Delta t} \int_0^{\Delta t}\int_0^\tau P_u^a\Phi_{ij}^{a,z}(s,\mu)\,du\,d\tau .\] Differentiating under the integral with respect to \(z\) and using 45 , we get \[\sum_{|\ell|\le2} \left\vert\kern-0.25ex\left\vert\kern-0.25ex\left\vert \partial_z^\ell R_{ij}^{a,\Delta t}(s,\mu,z) \right\vert\kern-0.25ex\right\vert\kern-0.25ex\right\vert_{\xi} \le \frac{1}{\Delta t} \int_0^{\Delta t}\int_0^\tau C(1+|s-z|+u^{1/2})\,du\,d\tau .\] Since \(\Delta t\le1\), \[\label{eq:RGij-z-bound-uniform-mkv} \sum_{|\ell|\le2} \left\vert\kern-0.25ex\left\vert\kern-0.25ex\left\vert \partial_z^\ell R_{ij}^{a,\Delta t}(s,\mu,z) \right\vert\kern-0.25ex\right\vert\kern-0.25ex\right\vert_{\xi} \le C(1+|s-z|)\Delta t .\tag{46}\]

Step 5: set \(z=s\). By definition, \[R_{G_{ij}^s}^{a,\Delta t}(s,\mu) = R_{ij}^{a,\Delta t}(s,\mu,s).\] For the spatial derivatives, the chain rule gives \[\nabla_s R_{G_{ij}^s}^{a,\Delta t}(s,\mu) = \left. (\nabla_s+\nabla_z)R_{ij}^{a,\Delta t}(s,\mu,z) \right|_{z=s},\] and \[D_{ss}^2R_{G_{ij}^s}^{a,\Delta t}(s,\mu) = \left. (D_{ss}^2+D_{sz}^2+D_{zs}^2+D_{zz}^2) R_{ij}^{a,\Delta t}(s,\mu,z) \right|_{z=s}.\] For the Lions derivatives, \(z=s\) is independent of \(\mu\), so \[\partial_\mu R_{G_{ij}^s}^{a,\Delta t}(s,\mu)(\xi) = \left. \partial_\mu R_{ij}^{a,\Delta t}(s,\mu,z)(\xi) \right|_{z=s},\] and the same identity holds for \(D_\xi\partial_\mu\). Evaluating 46 at \(z=s\) gives \[\left\vert\kern-0.25ex\left\vert\kern-0.25ex\left\vert R_{G_{ij}^s}^{a,\Delta t}(s,\mu) \right\vert\kern-0.25ex\right\vert\kern-0.25ex\right\vert_{\xi} \le C\Delta t .\] Together with the estimate for \(F_k\), this proves ?? . ◻

Proof. Fix \(a\in\mathcal{A}\), \(\Delta t\in(0,1]\), and \((s,\mu,\xi)\in\mathbb{R}^d\times\mathcal{P}_2(\mathbb{R}^d)\times\mathbb{R}^d\).

For the drift, let \(F_k(x,\nu):=x_k\). By ?? , ?? , and the definition of the remainder, \[\hat{b}_k(s,\mu,a)-b_k(s,\mu,a) = R_{F_k}^{a,\Delta t}(s,\mu), \qquad k=1,\dots,d.\] Hence, by 2, \[\left\vert\kern-0.25ex\left\vert\kern-0.25ex\left\vert \hat{b}_k(s,\mu,a)-b_k(s,\mu,a) \right\vert\kern-0.25ex\right\vert\kern-0.25ex\right\vert_{\xi} = \left\vert\kern-0.25ex\left\vert\kern-0.25ex\left\vert R_{F_k}^{a,\Delta t}(s,\mu) \right\vert\kern-0.25ex\right\vert\kern-0.25ex\right\vert_{\xi} \le C\Delta t.\] Summing over \(k=1,\dots,d\) and using equivalence of norms in finite dimension gives \[\left\vert\kern-0.25ex\left\vert\kern-0.25ex\left\vert \hat{b}(s,\mu,a)-b(s,\mu,a) \right\vert\kern-0.25ex\right\vert\kern-0.25ex\right\vert_{\xi}\le C\Delta t.\]

For the diffusion matrix, let \(G_{ij}^s(x,\nu):=(x_i-s_i)(x_j-s_j)\). By ?? , ?? , and the definition of the remainder, \[\hat{\Sigma}_{ij}(s,\mu,a)-\Sigma_{ij}(s,\mu,a) = R_{G_{ij}^s}^{a,\Delta t}(s,\mu), \qquad i,j=1,\dots,d.\] Therefore, by 2, \[\left\vert\kern-0.25ex\left\vert\kern-0.25ex\left\vert \hat{\Sigma}_{ij}(s,\mu,a)-\Sigma_{ij}(s,\mu,a) \right\vert\kern-0.25ex\right\vert\kern-0.25ex\right\vert_{\xi} = \left\vert\kern-0.25ex\left\vert\kern-0.25ex\left\vert R_{G_{ij}^s}^{a,\Delta t}(s,\mu) \right\vert\kern-0.25ex\right\vert\kern-0.25ex\right\vert_{\xi} \le C\Delta t.\] Summing over \(i,j=1,\dots,d\) and using equivalence of norms in finite dimension gives \[\left\vert\kern-0.25ex\left\vert\kern-0.25ex\left\vert \hat{\Sigma}(s,\mu,a)-\Sigma(s,\mu,a) \right\vert\kern-0.25ex\right\vert\kern-0.25ex\right\vert_{\xi}\le C\Delta t.\] This proves ?? . ◻

Proof. Fix \(\pi\in\Pi_*\). By 18, we may write \(\pi(da\mid s,\mu)=p_\pi(a\mid s,\mu)\nu(da)\). Moreover, there exists \(\ell_\pi\in L^1(\nu)\) such that \(p_\pi\) and all the derivatives of \(p_\pi\) appearing below are bounded by \(\ell_\pi(a)\), uniformly in \((s,\mu,\xi)\). More precisely, for \(\nu\)-a.e. \(a\in\mathcal{A}\), \[|p_\pi(a\mid s,\mu)| +\|\nabla_s p_\pi(a\mid s,\mu)\| +\|D_{ss}^2p_\pi(a\mid s,\mu)\| +|\partial_\mu p_\pi(a\mid s,\mu)(\xi)| +\|D_\xi\partial_\mu p_\pi(a\mid s,\mu)(\xi)\| \le \ell_\pi(a).\] In addition, for every \(R>0\), the same derivatives are locally Lipschitz in \((s,\mu,\xi)\) on the set \(|s|+\sqrt{m_2(\mu)}+|\xi|\le R\), with a Lipschitz constant bounded by some function \(\ell_{\pi,R}\in L^1(\nu)\).

We first prove an auxiliary statement. Let \(h=h(s,\mu,a)\) be scalar valued and assume that \(h\), \(\nabla_s h\), \(D_{ss}^2h\), \(\partial_\mu h\), and \(D_\xi\partial_\mu h\) are bounded and locally Lipschitz in their variables, uniformly with respect to \(a\). Define \[h^\pi(s,\mu):=\int_{\mathcal{A}}h(s,\mu,a)p_\pi(a\mid s,\mu)\nu(da).\] Since the integrands obtained by differentiating with respect to \(s\) and \(\mu\) are bounded by an \(L^1(\nu)\) function, the differentiation under the integral sign is justified. Thus \[\nabla_s h^\pi = \int_{\mathcal{A}} \bigl[ (\nabla_s h)p_\pi+h\nabla_s p_\pi \bigr]\,\nu(da),\] \[D_{ss}^2h^\pi = \int_{\mathcal{A}} \bigl[ (D_{ss}^2h)p_\pi +(\nabla_s h)\otimes(\nabla_s p_\pi) +(\nabla_s p_\pi)\otimes(\nabla_s h) +hD_{ss}^2p_\pi \bigr]\,\nu(da),\] and \[\partial_\mu h^\pi(\xi) = \int_{\mathcal{A}} \bigl[ \partial_\mu h(\xi)p_\pi +h\partial_\mu p_\pi(\xi) \bigr]\,\nu(da), \qquad D_\xi\partial_\mu h^\pi(\xi) = \int_{\mathcal{A}} \bigl[ D_\xi\partial_\mu h(\xi)p_\pi +hD_\xi\partial_\mu p_\pi(\xi) \bigr]\,\nu(da).\] Here all terms are evaluated at \((s,\mu,a)\), except the Lions derivatives, which are evaluated at \((s,\mu,a,\xi)\). The uniform bounds on \(h\) and its derivatives, together with the bound by \(\ell_\pi\), give \[\|\nabla_s h^\pi(s,\mu)\| +\|D_{ss}^2h^\pi(s,\mu)\| +|\partial_\mu h^\pi(s,\mu)(\xi)| +\|D_\xi\partial_\mu h^\pi(s,\mu)(\xi)\| \le K_\pi .\]

We now check local Lipschitz continuity. Let \(Z=(s,\mu,\xi)\) and \(\bar Z=(\bar s,\bar\mu,\bar\xi)\) belong to \(|s|+\sqrt{m_2(\mu)}+|\xi|\le R\) and \(|\bar s|+\sqrt{m_2(\bar\mu)}+|\bar\xi|\le R\). Set \[d(Z,\bar Z):= |s-\bar s|+\mathcal{W}_2(\mu,\bar\mu)+|\xi-\bar\xi|.\] For each term in the four displayed formulas for \(\nabla_s h^\pi\), \(D_{ss}^2h^\pi\), \(\partial_\mu h^\pi\), and \(D_\xi\partial_\mu h^\pi\), we use \[|AB-\bar A\bar B| \le |A-\bar A|\,|B|+|\bar A|\,|B-\bar B|.\] For instance, in the term \(\int_{\mathcal{A}}(\nabla_s h)p_\pi\,\nu(da)\), we have \[\begin{align} &\left| (\nabla_s h)(s,\mu,a)p_\pi(a\mid s,\mu) - (\nabla_s h)(\bar s,\bar\mu,a)p_\pi(a\mid \bar s,\bar\mu) \right|\le C_R d(Z,\bar Z)\ell_\pi(a) + C\,\ell_{\pi,R}(a)d(Z,\bar Z). \end{align}\] The right-hand side is integrable with respect to \(\nu(da)\). The same argument applies to each product appearing in the four formulas above, using the local Lipschitz bounds for \(h\) and its derivatives and the \(L^1(\nu)\) local Lipschitz bounds for \(p_\pi\) and its derivatives. Therefore, \[\begin{align} &\|\nabla_s h^\pi(s,\mu)-\nabla_s h^\pi(\bar s,\bar\mu)\| +\|D_{ss}^2h^\pi(s,\mu)-D_{ss}^2h^\pi(\bar s,\bar\mu)\| \\ &\quad +|\partial_\mu h^\pi(s,\mu)(\xi) -\partial_\mu h^\pi(\bar s,\bar\mu)(\bar\xi)| +\|D_\xi\partial_\mu h^\pi(s,\mu)(\xi) -D_\xi\partial_\mu h^\pi(\bar s,\bar\mu)(\bar\xi)\|\le C_{\pi,R}d(Z,\bar Z). \end{align}\] This proves the auxiliary statement.

Applying the auxiliary statement to \(h=b_i\), for \(i=1,\dots,d\), gives the required boundedness and local Lipschitz regularity of \(b^\pi\). Applying the same statement to \(h=\Sigma_{ij}\), for \(i,j=1,\dots,d\), gives the corresponding properties for \(\Sigma^\pi\). Indeed, \(\Sigma=\sigma\sigma^\top\), and 16 gives the required bounds and local Lipschitz regularity for \(\sigma\) and its derivatives.

It remains to pass from \(\Sigma^\pi\) to \(\sigma^\pi=(\Sigma^\pi)^{1/2}\). By 16, \[v^\top\Sigma^\pi(s,\mu)v = \int_{\mathcal{A}} v^\top\Sigma(s,\mu,a)v\,p_\pi(a\mid s,\mu)\nu(da) \ge \alpha |v|^2.\] Since \(\sigma\) is bounded, there exists \(M>0\) such that \(\|\Sigma^\pi(s,\mu)\|\le M\) for all \((s,\mu)\). Hence the image of \(\Sigma^\pi\) is contained in \[\mathcal{S}_{\alpha,M} := \{A\in\mathbb{S}^d:\alpha I_d\preceq A,\;\|A\|\le M\}.\] The principal square-root map \(\Psi(A)=A^{1/2}\) is smooth on the positive definite cone. Hence \(D\Psi\) and \(D^2\Psi\) are bounded and locally Lipschitz on \(\mathcal{S}_{\alpha,M}\).

Since \(\sigma^\pi=\Psi(\Sigma^\pi)\), the chain rule gives \[\nabla_s\sigma^\pi = D\Psi(\Sigma^\pi)[\nabla_s\Sigma^\pi],\] \[D_{ss}^2\sigma^\pi = D\Psi(\Sigma^\pi)[D_{ss}^2\Sigma^\pi] + D^2\Psi(\Sigma^\pi)[\nabla_s\Sigma^\pi,\nabla_s\Sigma^\pi],\] and \[\partial_\mu\sigma^\pi(\xi) = D\Psi(\Sigma^\pi)[\partial_\mu\Sigma^\pi(\xi)], \qquad D_\xi\partial_\mu\sigma^\pi(\xi) = D\Psi(\Sigma^\pi)[D_\xi\partial_\mu\Sigma^\pi(\xi)].\] The bounds and local Lipschitz regularity already proved for \(\Sigma^\pi\), together with the boundedness and local Lipschitz regularity of \(D\Psi\) and \(D^2\Psi\) on \(\mathcal{S}_{\alpha,M}\), imply the same coefficient regularity for \(\sigma^\pi\). Therefore \((b^\pi,\sigma^\pi)\) satisfies the coefficient part of 2. ◻

Proof. Fix \(\pi\in\Pi_*\). By 18, we may write \(\pi(da\mid s,\mu)=p_\pi(a\mid s,\mu)\nu(da)\). Moreover, there exists \(\ell_\pi\in L^1(\nu)\) such that \(p_\pi\) and all derivatives of \(p_\pi\) appearing below are bounded by \(\ell_\pi(a)\), uniformly in \((s,\mu,\xi)\).

Set \[\delta b(s,\mu,a):=\hat{b}(s,\mu,a)-b(s,\mu,a), \qquad \delta\Sigma(s,\mu,a):=\hat{\Sigma}(s,\mu,a)-\Sigma(s,\mu,a).\] By 3, for every component \(i,j\), \[\left\vert\kern-0.25ex\left\vert\kern-0.25ex\left\vert \delta b_i(s,\mu,a) \right\vert\kern-0.25ex\right\vert\kern-0.25ex\right\vert_\xi + \left\vert\kern-0.25ex\left\vert\kern-0.25ex\left\vert \delta\Sigma_{ij}(s,\mu,a) \right\vert\kern-0.25ex\right\vert\kern-0.25ex\right\vert_\xi \le C\Delta t,\] uniformly in \((s,\mu,\xi,a)\).

For each component of the averaged drift error, \[\hat{b}_i^\pi(s,\mu)-b_i^\pi(s,\mu) = \int_{\mathcal{A}}\delta b_i(s,\mu,a)p_\pi(a\mid s,\mu)\nu(da).\] Differentiating this identity with respect to \(s\) and \(\mu\) gives the same product formulas used in the proof of 4, with \(h\) replaced by \(\delta b_i\). Since \(\delta b_i\) and the derivatives entering \(\left\vert\kern-0.25ex\left\vert\kern-0.25ex\left\vert \cdot \right\vert\kern-0.25ex\right\vert\kern-0.25ex\right\vert_\xi\) are bounded by \(C\Delta t\), and since \(p_\pi\) and its corresponding derivatives are bounded by \(\ell_\pi\in L^1(\nu)\), we obtain \[\left\vert\kern-0.25ex\left\vert\kern-0.25ex\left\vert \hat{b}_i^\pi(s,\mu)-b_i^\pi(s,\mu) \right\vert\kern-0.25ex\right\vert\kern-0.25ex\right\vert_\xi \le C_\pi\Delta t.\] Summing over \(i=1,\dots,d\) gives ?? .

The proof of ?? is the same. Indeed, \[\hat{\Sigma}_{ij}^\pi(s,\mu)-\Sigma_{ij}^\pi(s,\mu) = \int_{\mathcal{A}} \delta\Sigma_{ij}(s,\mu,a)p_\pi(a\mid s,\mu)\nu(da),\] and the same differentiation formulas, now with \(\delta\Sigma_{ij}\) in place of \(\delta b_i\), give \[\left\vert\kern-0.25ex\left\vert\kern-0.25ex\left\vert \hat{\Sigma}_{ij}^\pi(s,\mu)-\Sigma_{ij}^\pi(s,\mu) \right\vert\kern-0.25ex\right\vert\kern-0.25ex\right\vert_\xi \le C_\pi\Delta t.\] Summing over \(i,j=1,\dots,d\) proves ?? .

We now verify the coefficient regularity of \(\hat{b}^\pi\). By 4, \(b^\pi\) satisfies the coefficient part of 2. Hence the derivatives of \(b^\pi\) entering \(\left\vert\kern-0.25ex\left\vert\kern-0.25ex\left\vert \cdot \right\vert\kern-0.25ex\right\vert\kern-0.25ex\right\vert_\xi\) are bounded and locally Lipschitz. The estimate ?? gives boundedness of the corresponding derivatives of \(\hat{b}^\pi-b^\pi\). Their local Lipschitz continuity follows from the local Lipschitz part of 2, together with the \(L^1(\nu)\) local Lipschitz bounds for the derivatives of \(p_\pi\) in 18, using the same product estimate as in the proof of 4. Therefore \(\hat{b}^\pi\) has the boundedness and local Lipschitz regularity required in the coefficient part of 2.

It remains to treat \(\hat{\sigma}^\pi\). From ?? , \[\|\hat{\Sigma}^\pi(s,\mu)-\Sigma^\pi(s,\mu)\|\le C_\pi\Delta t.\] By 4, \(\Sigma^\pi(s,\mu)\succeq \alpha I_d\). Choose \(\Delta_\pi\in(0,1]\) such that \(C_\pi\Delta_\pi\le\alpha/2\). Then, for \(\Delta t\in(0,\Delta_\pi]\), \[\hat{\Sigma}^\pi(s,\mu)\succeq \frac{\alpha}{2}I_d.\] Moreover, since \(\Sigma^\pi\) is uniformly bounded and ?? holds, there exists \(M_\pi>0\) such that \(\|\hat{\Sigma}^\pi(s,\mu)\|\le M_\pi\) for all \((s,\mu)\). Thus the image of \(\hat{\Sigma}^\pi\) is contained in \[\mathcal{S}_{\alpha/2,M_\pi} := \{A\in\mathbb{S}^d:\tfrac{\alpha}{2}I_d\preceq A,\;\|A\|\le M_\pi\}.\] The principal square-root map \(\Psi(A)=A^{1/2}\) is smooth on the positive definite cone. Hence \(D\Psi\) and \(D^2\Psi\) are bounded and locally Lipschitz on \(\mathcal{S}_{\alpha/2,M_\pi}\).

Since \(\hat{\sigma}^\pi=\Psi(\hat{\Sigma}^\pi)\), the chain rule gives \[\nabla_s\hat{\sigma}^\pi = D\Psi(\hat{\Sigma}^\pi)[\nabla_s\hat{\Sigma}^\pi],\] \[D_{ss}^2\hat{\sigma}^\pi = D\Psi(\hat{\Sigma}^\pi)[D_{ss}^2\hat{\Sigma}^\pi] + D^2\Psi(\hat{\Sigma}^\pi)[\nabla_s\hat{\Sigma}^\pi,\nabla_s\hat{\Sigma}^\pi],\] and \[\partial_\mu\hat{\sigma}^\pi(\xi) = D\Psi(\hat{\Sigma}^\pi)[\partial_\mu\hat{\Sigma}^\pi(\xi)], \qquad D_\xi\partial_\mu\hat{\sigma}^\pi(\xi) = D\Psi(\hat{\Sigma}^\pi)[D_\xi\partial_\mu\hat{\Sigma}^\pi(\xi)].\] The bounds and local Lipschitz regularity of the derivatives of \(\hat{\Sigma}^\pi\), together with the boundedness and local Lipschitz regularity of \(D\Psi\) and \(D^2\Psi\) on \(\mathcal{S}_{\alpha/2,M_\pi}\), imply that \(\hat{\sigma}^\pi\) satisfies the coefficient regularity required in 2. Therefore \((\hat{b}^\pi,\hat{\sigma}^\pi)\) satisfies the coefficient part of 2. ◻

Proof of 20. Let \(\Pi_*:=\{\pi^*,\hat{\pi}^*\}\). Choose \(\Delta_0\in(0,1]\) and \(\bar\beta_*>0\) so that, for every \(\pi\in\Pi_*\) and every \(\Delta t\in(0,\Delta_0]\), both coefficient pairs \((b^\pi,\sigma^\pi)\) and \((\hat{b}^\pi,\hat{\sigma}^\pi)\) satisfy the coefficient part of 2, and 4 applies to the corresponding evaluation equations. Choose \(\tilde{\beta}>\max\{\bar\beta_*,\beta_0(2)\}\), and fix \(\beta>\tilde{\beta}\) and \(\Delta t\in(0,\Delta_0]\).

Step 1: fixed-policy estimate. Fix \(\pi\in\Pi_*\). By the choice of \(\bar\beta_*\), the functions \(\mathcal{V}^\pi\) and \(\hat{\mathcal{V}}^\pi\) are the unique classical solutions of \[(\mathcal{L}_{b,\Sigma}^{\pi}-\beta)\mathcal{V}^\pi+r_\lambda^\pi=0\] and \[(\mathcal{L}_{\hat{b},\hat{\Sigma}}^{\pi}-\beta)\hat{\mathcal{V}}^\pi+r_\lambda^\pi=0,\] respectively. Set \(w^\pi:=\mathcal{V}^\pi-\hat{\mathcal{V}}^\pi\). Subtracting the two equations gives \[\label{eq:w-linear-HJB-fixed-policy-mkv} (\mathcal{L}_{b,\Sigma}^{\pi}-\beta)w^\pi+ (\mathcal{L}_{b,\Sigma}^{\pi}-\mathcal{L}_{\hat{b},\hat{\Sigma}}^{\pi}) \hat{\mathcal{V}}^\pi=0.\tag{47}\]

Let \((\mu_t)_{t\ge0}\) be the law flow generated by 5 under the policy \(\pi\), with \(\mu_0=\mu\). Applying the Itô-Lions formula to \(e^{-\beta t}w^\pi(\mu_t)\) and using 47 , we obtain, for every \(T>0\), \[w^\pi(\mu)=e^{-\beta T}w^\pi(\mu_T)+\int_0^T e^{-\beta t}\Bigl[(\mathcal{L}_{b,\Sigma}^{\pi}-\mathcal{L}_{\hat{b},\hat{\Sigma}}^{\pi})\hat{\mathcal{V}}^\pi\Bigr](\mu_t)\,dt.\] By the polynomial growth of \(\mathcal{V}^\pi\) and \(\hat{\mathcal{V}}^\pi\), together with the second-moment estimate associated with \(\beta_0(2)\), \[e^{-\beta T}|w^\pi(\mu_T)|\le C e^{-(\beta-\beta_0(2))T}(1+m_2(\mu))\longrightarrow0, \qquad T\to\infty.\] Hence \[\label{eq:w-representation-fixed-policy-mkv} w^\pi(\mu)=\int_0^\infty e^{-\beta t}\Bigl[(\mathcal{L}_{b,\Sigma}^{\pi}-\mathcal{L}_{\hat{b},\hat{\Sigma}}^{\pi})\hat{\mathcal{V}}^\pi\Bigr](\mu_t)\,dt.\tag{48}\]

For \(\nu\in\mathcal{P}_2(\mathbb{R}^d)\), by 6 and 12 , \[\begin{align} \Bigl[(\mathcal{L}_{b,\Sigma}^{\pi}-\mathcal{L}_{\hat{b},\hat{\Sigma}}^{\pi})\hat{\mathcal{V}}^\pi\Bigr](\nu) &=\int_{\mathbb{R}^d}\bigl(b^\pi(x,\nu)-\hat{b}^\pi(x,\nu)\bigr)\cdot\partial_\mu\hat{\mathcal{V}}^\pi(\nu)(x)\,\nu(dx) \\ &\quad+\frac{1}{2}\int_{\mathbb{R}^d}\bigl(\Sigma^\pi(x,\nu)-\hat{\Sigma}^\pi(x,\nu)\bigr):D_x\partial_\mu\hat{\mathcal{V}}^\pi(\nu)(x)\,\nu(dx). \end{align}\] By 5, \[\|b^\pi(x,\nu)-\hat{b}^\pi(x,\nu)\|+\|\Sigma^\pi(x,\nu)-\hat{\Sigma}^\pi(x,\nu)\|\le C_\pi\Delta t.\] Moreover, the derivative estimate in 4, applied to the MF-PhiBE, gives \[|\partial_\mu\hat{\mathcal{V}}^\pi(\nu)(x)|+\|D_x\partial_\mu\hat{\mathcal{V}}^\pi(\nu)(x)\|\le\frac{C_\pi}{\beta-\bar\beta_*}(1+|x|^2+m_2(\nu)).\] Combining the previous two estimates yields \[\begin{align} \left|\Bigl[(\mathcal{L}_{b,\Sigma}^{\pi}-\mathcal{L}_{\hat{b},\hat{\Sigma}}^{\pi})\hat{\mathcal{V}}^\pi\Bigr](\nu)\right| &\le\frac{C_\pi\Delta t}{\beta-\bar\beta_*}\int_{\mathbb{R}^d}(1+|x|^2+m_2(\nu))\,\nu(dx)\le\frac{C_\pi}{\beta-\bar\beta_*}(1+m_2(\nu))\Delta t. \end{align}\] Using this bound in 48 , we get \[|w^\pi(\mu)|\le\frac{C_\pi\Delta t}{\beta-\bar\beta_*}\int_0^\infty e^{-\beta t}(1+m_2(\mu_t))\,dt.\] The moment estimate associated with \(\beta_0(2)\) gives \[1+m_2(\mu_t)\le C e^{\beta_0(2)t}(1+m_2(\mu)).\] Since \(\beta>\beta_0(2)\), we obtain \[\label{eq:fixed-policy-error-final-mkv} |\mathcal{V}^\pi(\mu)-\hat{\mathcal{V}}^\pi(\mu)| \le \frac{C_\pi}{(\beta-\bar\beta_*)(\beta-\beta_0(2))}(1+m_2(\mu))\Delta t.\tag{49}\]

Step 2: optimality argument. By the definition of \(\pi^*\), \(\mathcal{V}^*(\mu)=\mathcal{V}^{\pi^*}(\mu).\) By the optimality of \(\hat{\pi}^*\) for the MF-PhiBE, \[\label{eq:estimated-optimality-proof-mkv} \hat{\mathcal{V}}^{\hat{\pi}^*}(\mu)\ge\hat{\mathcal{V}}^{\pi^*}(\mu).\tag{50}\] Adding and subtracting \(\hat{\mathcal{V}}^{\pi^*}(\mu)\) and \(\hat{\mathcal{V}}^{\hat{\pi}^*}(\mu)\) gives \[\begin{align} \mathcal{V}^*(\mu)-\mathcal{V}^{\hat{\pi}^*}(\mu) &=\mathcal{V}^{\pi^*}(\mu)-\mathcal{V}^{\hat{\pi}^*}(\mu) \\ &=\bigl(\mathcal{V}^{\pi^*}-\hat{\mathcal{V}}^{\pi^*}\bigr)(\mu) +\bigl(\hat{\mathcal{V}}^{\pi^*}-\hat{\mathcal{V}}^{\hat{\pi}^*}\bigr)(\mu) +\bigl(\hat{\mathcal{V}}^{\hat{\pi}^*}-\mathcal{V}^{\hat{\pi}^*}\bigr)(\mu). \end{align}\] By 50 , the middle term is nonpositive. Therefore, \[\mathcal{V}^*(\mu)-\mathcal{V}^{\hat{\pi}^*}(\mu) \le |\mathcal{V}^{\pi^*}(\mu)-\hat{\mathcal{V}}^{\pi^*}(\mu)| + |\hat{\mathcal{V}}^{\hat{\pi}^*}(\mu)-\mathcal{V}^{\hat{\pi}^*}(\mu)|.\] Applying 49 with \(\pi=\pi^*\) and with \(\pi=\hat{\pi}^*\) yields \[\mathcal{V}^*(\mu)-\mathcal{V}^{\hat{\pi}^*}(\mu) \le \frac{C}{(\beta-\bar\beta_*)(\beta-\beta_0(2))}(1+m_2(\mu))\Delta t.\] Finally, since \(\hat{\pi}^*\in\Pi_{\mathrm{add}}\), the definition 1 gives \(\mathcal{V}^*(\mu)-\mathcal{V}^{\hat{\pi}^*}(\mu)\ge0\), which proves ?? . ◻

Proof. For \(\hat{A}\), we have \[\begin{align} \label{eq:estimate95A95hat95A} \hat{A}-A = \frac{e^{A\Delta t}-I_d-A\Delta t}{\Delta t} = \frac{1}{\Delta t}\int_0^{\Delta t}(e^{A\tau}-I_d)A\,d\tau. \end{align}\tag{51}\] Taking norm in 51 , and using that \(\|e^{A\tau}-I_d\| \le \tau c_A e^{c_A\tau},\) we obtain \[\|\hat{A}-A\| \le \frac{c_A}{\Delta t}\int_0^{\Delta t}\tau c_A e^{c_A\tau}\,d\tau \le \frac{c_A^2}{\Delta t}e^{c_A\Delta t}\int_0^{\Delta t}\tau\,d\tau = \frac{1}{2} c_A^2 e^{c_A\Delta t}\Delta t.\]

For \(\hat{\bar A}\), we write \[\hat{\bar A}-\bar A = \frac{1}{\Delta t}\int_0^{\Delta t} \Bigl( e^{A(\Delta t-\tau)}\bar A e^{\tilde{A}\tau}-\bar A \Bigr)\,d\tau.\] Adding and subtracting \(\bar A e^{\tilde{A}\tau}\) and taking norm, we have \[\begin{align} \label{eq:estimate95A95hat95A95tilde} \|\hat{\bar A}-\bar A\| \le \frac{c_{\bar A}}{\Delta t}\int_0^{\Delta t} \Bigl[ \|e^{A(\Delta t-\tau)}-I_d\|\,\|e^{\tilde{A}\tau}\| + \|e^{\tilde{A}\tau}-I_d\| \Bigr]\,d\tau. \end{align}\tag{52}\] Combining \[\|e^{A(\Delta t-\tau)}-I_d\| \le (\Delta t-\tau)c_A e^{c_A(\Delta t-\tau)}, \qquad \|e^{\tilde{A}\tau}\|\le e^{c_{\tilde{A}}\tau},\quad \text{and}\quad\|e^{\tilde{A}\tau}-I_d\| \le \tau c_{\tilde{A}}e^{c_{\tilde{A}}\tau},\] with 52 , we conclude that \[\|\hat{\bar A}-\bar A\| \le \frac{c_{\bar A}}{\Delta t} \int_0^{\Delta t} \Bigl[ c_A(\Delta t-\tau) + c_{\tilde{A}}\tau \Bigr] e^{(c_A+c_{\tilde{A}})\Delta t}\,d\tau\le \frac{1}{2} c_{\bar A}(c_A+c_{\tilde{A}}) e^{(c_A+c_{\tilde{A}})\Delta t}\Delta t.\]

Finally, for \(B\) we have \[\begin{align} \label{eq:error95B95bhat} \hat{B}-B = \frac{1}{\Delta t}\int_0^{\Delta t} (e^{A(\Delta t-\tau)}-I_d)B\,d\tau + \frac{1}{\Delta t}\int_0^{\Delta t} e^{A(\Delta t-\tau)}\bar A M_{\tilde{A},\tau}B\,d\tau. \end{align}\tag{53}\] For the first term on the right-hand side of 53 we have, \[\begin{align} \label{eq:estimate95b95bhat1} \left\| \frac{1}{\Delta t}\int_0^{\Delta t} (e^{A(\Delta t-\tau)}-I_d)B\,d\tau \right\| \le \frac{\|B\|}{\Delta t}\int_0^{\Delta t} (\Delta t-\tau)c_Ae^{c_A(\Delta t-\tau)}\,d\tau \le \frac{1}{2} c_A e^{c_A\Delta t}\|B\|\Delta t. \end{align}\tag{54}\] For the second term on the right-hand side of 53 , since \[\|M_{\tilde{A},\tau}\| = \left\| \int_0^\tau e^{\tilde{A}(\tau-r)}\,dr \right\| \le \int_0^\tau e^{c_{\tilde{A}}(\tau-r)}\,dr \le \tau e^{c_{\tilde{A}}\tau},\] we have \[\begin{align} \label{eq:estimate95b95bhat2} \frac{1}{\Delta t}\int_0^{\Delta t} \|e^{A(\Delta t-\tau)}\bar A M_{\tilde{A},\tau}B\|\,d\tau \le \frac{c_{\bar A}\|B\|}{\Delta t} \int_0^{\Delta t} e^{c_A(\Delta t-\tau)}\tau e^{c_{\tilde{A}}\tau}\,d\tau \le \frac{1}{2} c_{\bar A}e^{(c_A+c_{\tilde{A}})\Delta t}\|B\|\Delta t. \end{align}\tag{55}\] Combining both estimates 54 and 55 , and using \(e^{c_A\Delta t}\le e^{(c_A+c_{\tilde{A}})\Delta t}\), gives \[\|\hat{B}-B\| \le \frac{1}{2} \left( c_A e^{c_A\Delta t} + c_{\bar A}e^{(c_A+c_{\tilde{A}})\Delta t} \right) \|B\|\Delta t\leq \frac{1}{2} (c_A+c_{\bar A})e^{(c_A+c_{\tilde{A}})\Delta t} \|B\|\Delta t.\] ◻

Proof. We verify the two parts of 22. Since \(I=0\), condition (H1) reduces to \(N\in\mathbb{S}_{++}^m,\) \(Q\in\mathbb{S}_{++}^d,\) \(Q+\bar Q\in\mathbb{S}_{++}^d.\) This condition does not depend on \(A,\bar A\), or \(B\), and therefore is unchanged when \((A,\bar A,B)\) is replaced by \((\hat{A},\hat{\bar A},\hat{B})\).

We now prove (H2). Since \(\sigma(x,\mu,a)=\gamma\), we have \(D=F=0\). Set \(A_\beta:=A-\frac{\beta}{2}I_d\), \(\tilde{A}_\beta:=A+\bar A-\frac{\beta}{2}I_d\), \(\hat{A}_\beta:=\hat{A}-\frac{\beta}{2}I_d\), and \(\widehat{\tilde{A}}_\beta:=\hat{A}+\hat{\bar A}-\frac{\beta}{2}I_d\). By 6, \(\|\hat{A}-A\|\), \(\|\hat{\bar A}-\bar A\|\), and \(\|\hat{B}-B\|\) converge to zero as \(\Delta t\to0\). Hence \[\|\widehat{\tilde{A}}_\beta-\tilde{A}_\beta\| \le \|\hat{A}-A\|+\|\hat{\bar A}-\bar A\| \to0.\]

We first prove that \((\widehat{\tilde{A}}_\beta,\hat{B})\) is stabilizable. Since 22 holds for the original coefficients, \((\tilde{A}_\beta,B)\) is stabilizable. Thus there exists \(\Theta_m\in\mathbb{R}^{m\times d}\) such that \(L_m:=\tilde{A}_\beta+B\Theta_m\) is Hurwitz. Hence there exists \(P_m\in\mathbb{S}_{++}^d\) such that \[L_m^\top P_m+P_mL_m\preceq -I_d.\] Using the same feedback \(\Theta_m\), define \(\hat{L}_m:=\widehat{\tilde{A}}_\beta+\hat{B}\Theta_m\). Then \(\hat{L}_m=L_m+E_m\), where \[E_m:=(\widehat{\tilde{A}}_\beta-\tilde{A}_\beta)+(\hat{B}-B)\Theta_m.\] Since \(E_m\to0\) as \(\Delta t\to0\), there exists \(\Delta t_1>0\) such that \(2\|P_m\|\,\|E_m\|\le\frac{1}{2}\) for every \(\Delta t\in(0,\Delta t_1]\). Therefore \[\hat{L}_m^\top P_m+P_m\hat{L}_m = L_m^\top P_m+P_mL_m+E_m^\top P_m+P_mE_m \preceq -I_d+2\|P_m\|\,\|E_m\|I_d \preceq -\frac{1}{2} I_d.\] Thus \(\hat{L}_m\) is Hurwitz, and \((\widehat{\tilde{A}}_\beta,\hat{B})\) is stabilizable.

It remains to verify the Lyapunov inequality in (H2). Since 22 holds for the original coefficients and \(D=F=0\), there exist \(\Theta\in\mathbb{R}^{m\times d}\) and \(P_0\in\mathbb{S}_{++}^d\) such that \[(A_\beta+B\Theta)^\top P_0+P_0(A_\beta+B\Theta)\preceq -I_d.\] Set \(L:=A_\beta+B\Theta\) and \(\hat{L}:=\hat{A}_\beta+\hat{B}\Theta\). Then \(\hat{L}=L+E\), with \(E:=(\hat{A}-A)+(\hat{B}-B)\Theta\). Since \(E\to0\) as \(\Delta t\to0\), there exists \(\Delta t_2>0\) such that \(2\|P_0\|\,\|E\|\le\frac{1}{2}\) for every \(\Delta t\in(0,\Delta t_2]\). Hence \[\hat{L}^\top P_0+P_0\hat{L} = L^\top P_0+P_0L+E^\top P_0+P_0E \preceq -I_d+2\|P_0\|\,\|E\|I_d \preceq -\frac{1}{2} I_d.\] Multiplying the matrix \(P_0\) by \(2\), we obtain \[\hat{L}^\top(2P_0)+(2P_0)\hat{L}\preceq -I_d.\] Since \(2P_0\in\mathbb{S}_{++}^d\), the Lyapunov inequality in (H2) holds for the estimated coefficients.

Choosing \(\Delta t_0>0\) small enough so that both previous arguments apply, conditions (H1) and (H2) hold for \((\hat{A},\hat{\bar A},\hat{B})\) for every \(\Delta t\in(0,\Delta t_0]\). ◻

Proof. By 7, for every \(\Delta t\in(0,\Delta t_0]\), 22 holds with \((A,\bar A,B)\) replaced by \((\hat{A},\hat{\bar A},\hat{B})\). Moreover, by 19 , the auxiliary dynamics ?? are an entropy-regularized McKean-Vlasov LQR system with drift matrices \((\hat{A},\hat{\bar A},\hat{B})\) and zero diffusion. Hence 24, applied to this auxiliary LQR problem, yields the Gaussian feedback ?? and the optimal value \(\widehat{\mathcal{V}}^{\mathrm{aux}}\).

It remains to identify this auxiliary value with the optimal LQR-MF-PhiBE value in ?? . For a fixed policy \(\pi\in\Pi_{\mathrm{add}}\), let \[\widehat{\mathcal{J}}^{\pi}(\mu) := \mathbb{E}\left[ \int_0^\infty e^{-\beta t}r_\lambda^\pi(\hat{\mu}_t)\,dt \,\bigg|\, \hat{\mu}_0=\mu \right],\] where \((\hat{\mu}_t)_{t\ge0}\) is generated by ?? . Since this law flow has generator \(\mathcal{L}_{\hat{b},0}^{\pi}\), the policy-evaluation equation associated with \(\widehat{\mathcal{J}}^{\pi}\) is \[(\mathcal{L}_{\hat{b},0}^{\pi}-\beta)\widehat{\mathcal{J}}^{\pi}(\mu) + r_\lambda^\pi(\mu) = 0, \qquad \mu\in\mathcal{P}_2(\mathbb{R}^d).\] This is precisely the LQR-MF-PhiBE evaluation equation ?? . Therefore, by the notation introduced in 6, we have \(\widehat{\mathcal{J}}^{\pi}(\mu)=\hat{\mathcal{V}}^{\pi}(\mu)\) for every \(\mu\in\mathcal{P}_2(\mathbb{R}^d).\) Taking the supremum over \(\pi\in\Pi_{\mathrm{add}}\) gives \[\widehat{\mathcal{V}}^{\mathrm{aux}}(\mu) = \sup_{\pi\in\Pi_{\mathrm{add}}}\widehat{\mathcal{J}}^{\pi}(\mu) = \sup_{\pi\in\Pi_{\mathrm{add}}}\hat{\mathcal{V}}^{\pi}(\mu) = \hat{\mathcal{V}}^*(\mu).\] Since ?? attains the first supremum, it also attains the supremum in ?? . Hence ?? is an optimal LQR-MF-PhiBE policy. ◻

Proof. Since \(X\) solves ?? , we have \[A_\beta^\top X+XA_\beta-XGX=-Q_0.\] Therefore, with \(\Phi_X=A_\beta-GX\), \[\label{eq:closed95loop95lyapunov95identity95X} \Phi_X^\top X+X\Phi_X = -Q_0-XGX \preceq -Q_0.\tag{56}\] Let \(z(t):=e^{\Phi_Xt}z_0\). From 56 , \[\frac{d}{dt}\,z(t)^\top Xz(t) = z(t)^\top \bigl(\Phi_X^\top X+X\Phi_X\bigr) z(t) \le -\lambda_{\min}(Q_0)\|z(t)\|^2.\] Since \(z^\top Xz\le \|X\|\|z\|^2\), it follows that \[\frac{d}{dt}\,z(t)^\top Xz(t) \le - \frac{\lambda_{\min}(Q_0)}{\|X\|} z(t)^\top Xz(t).\] Thus \[z(t)^\top Xz(t) \le e^{-\lambda_{\min}(Q_0)t/\|X\|} z_0^\top Xz_0.\] Using \(z(t)^\top Xz(t)\ge\lambda_{\min}(X)\|z(t)\|^2\) and \(z_0^\top Xz_0\le \|X\|\|z_0\|^2\), we obtain \[\label{eq:semigroup95bound95X95metric} \|e^{\Phi_Xt}\|^2 \le \frac{\|X\|}{\lambda_{\min}(X)} e^{-\lambda_{\min}(Q_0)t/\|X\|}, \qquad t\ge0.\tag{57}\]

For every symmetric matrix \(Y\), the solution of \(\mathcal{L}_XW=Y\) is \[W = -\int_0^\infty e^{\Phi_X^\top t} Y e^{\Phi_Xt}\,dt.\] Therefore, by 57 , \[\|\mathcal{L}_X^{-1}Y\| \le \int_0^\infty \|e^{\Phi_Xt}\|^2\,dt\,\|Y\| \le \frac{\|X\|^2}{\lambda_{\min}(X)\lambda_{\min}(Q_0)} \|Y\|.\] Hence \[\|\mathcal{L}_X^{-1}\| \le \frac{\|X\|^2}{\lambda_{\min}(X)\lambda_{\min}(Q_0)},\] which gives the first term in ?? .

Next, for any \(Z\in\mathbb{R}^{d\times d}\), \[\|\mathcal{L}_X^{-1}(XZ+Z^\top X)\| \le \|\mathcal{L}_X^{-1}\| \bigl(\|XZ\|+\|Z^\top X\|\bigr) \le \frac{2\|X\|^3}{\lambda_{\min}(X)\lambda_{\min}(Q_0)} \|Z\|.\] This proves the estimate for \(p_X\). Similarly, for any \(H\in\mathbb{S}^d\), \[\|\mathcal{L}_X^{-1}(XHX)\| \le \|\mathcal{L}_X^{-1}\|\,\|X\|^2\,\|H\| \le \frac{\|X\|^4}{\lambda_{\min}(X)\lambda_{\min}(Q_0)} \|H\|.\] This proves the estimate for \(q_X\).

It remains to prove ?? . Let \(v\) be a unit eigenvector of \(X\) associated with \(\lambda_{\min}(X)\). Evaluating ?? on \(v\) gives \[0 = v^\top Q_0v + 2v^\top XA_\beta v - v^\top XGXv.\] Since \(Xv=\lambda_{\min}(X)v\) and \(A_\beta=A_0-\frac{\beta}{2}I_d\), we get \[0 = v^\top Q_0v + 2\lambda_{\min}(X)v^\top A_0v - \beta\lambda_{\min}(X) - \lambda_{\min}(X)^2v^\top Gv.\] Hence \[v^\top Q_0v = \beta\lambda_{\min}(X) + \lambda_{\min}(X)^2v^\top Gv - 2\lambda_{\min}(X)v^\top A_0v.\] Using \(v^\top Q_0v\ge\lambda_{\min}(Q_0)\), \(v^\top Gv\le\|G\|\), and \(-v^\top A_0v\le\|A_0\|\), we obtain \[\lambda_{\min}(Q_0) \le \bigl(\beta+2\|A_0\|\bigr)\lambda_{\min}(X) + \|G\|\lambda_{\min}(X)^2.\] Since \(\lambda_{\min}(X)\le\|X\|\), this implies \[\lambda_{\min}(Q_0) \le \bigl(\beta+2\|A_0\|+\|G\|\|X\|\bigr) \lambda_{\min}(X),\] which proves ?? . Finally, ?? follows by combining the first term in ?? with ?? , and ?? ?? follow by combining ?? with ?? . ◻

Proof of 27. We use the spectral norm. We first choose \(\Delta t_*\in(0,1]\) sufficiently small so that the perturbation estimates below and the smallness conditions in [42] and 7 hold for every \(\Delta t\in(0,\Delta t_*]\). Then, 30 characterizes the optimal policy \(\hat{\pi}\) as ?? . Define \(C_K:=2\max\{1,\|K\|\}\), \(C_\Lambda:=2\max\{1,\|\Lambda\|\}\), and \[\begin{align} \label{eq:constant95C95A95b} C_{\tilde{A}} := \frac{C}{2}(\|A\|+\|A+\bar A\|), \qquad C_B := \frac{C}{2}\|B\|. \end{align}\tag{58}\] We also set \[\label{eq:constant95C95K952} C_{K,2} := 2C_K^3 C_{\tilde{A}} + C_K^4 \Bigl[ 2C_B\|N^{-1}\|\|B\| + C_B^2\|N^{-1}\|\Delta t_* \Bigr],\tag{59}\] \[\label{eq:constant95C95lambd952} C_{\Lambda,2} := 4C_\Lambda^3 C_{\tilde{A}} + C_\Lambda^4 \Bigl[ 2C_B\|N^{-1}\|\|B\| + C_B^2\|N^{-1}\|\Delta t_* \Bigr],\tag{60}\] and \[\label{eq:constant95C95K95195C95Lambd952} C_{K,1} := \bigl(2\|A\|+\|BN^{-1}B^\top\|C_K\bigr)C_{K,2}, \qquad C_{\Lambda,1} := \bigl(2\|A+\bar A\|+\|BN^{-1}B^\top\|C_\Lambda\bigr) C_{\Lambda,2},\tag{61}\] together with \(C_{K,3}:=C_BC_K\) and \(C_{\Lambda,3}:=C_BC_\Lambda\). The constants above are fixed throughout the proof.

Step 1: coefficient estimates. By 6, and by the definition of \(C_{\tilde{A}}\), we have \[\label{eq:coeff95errors95beta95explicit95proof} \|\hat{A}-A\|\le C_{\tilde{A}}\Delta t, \qquad \|\hat{\bar A}-\bar A\|\le C_{\tilde{A}}\Delta t, \qquad \|\hat{B}-B\|\le C_B\Delta t.\tag{62}\] For the proof only, set \[G:=BN^{-1}B^\top, \qquad \hat{G}:=\hat{B} N^{-1}\hat{B}^\top.\] Then \[\hat{G}-G = (\hat{B}-B)N^{-1}B^\top + BN^{-1}(\hat{B}-B)^\top + (\hat{B}-B)N^{-1}(\hat{B}-B)^\top.\] Using 62 , we obtain \[\|\hat{G}-G\| \le 2C_B\|N^{-1}\|\|B\|\Delta t + C_B^2\|N^{-1}\|\Delta t^2.\] Since \(\Delta t\le \Delta t_*\), this yields \[\label{eq:G95error95beta95explicit95proof} \|\hat{G}-G\| \le \Bigl[ 2C_B\|N^{-1}\|\|B\| + C_B^2\|N^{-1}\|\Delta t_* \Bigr]\Delta t.\tag{63}\]

Step 2: perturbation of the centered Riccati equation. Since \(D=F=I=0\), [eq:ARE95SIMPLIFY95RICCATI95K] reduces to \[\beta K = Q+KA+A^\top K-KBN^{-1}B^\top K.\] Equivalently, \[\label{eq:K95exact95CARE95beta95explicit95proof} Q + \left(A-\frac{\beta}{2}I_d\right)^\top K + K\left(A-\frac{\beta}{2}I_d\right) - KGK = 0.\tag{64}\] Similarly, [eq:LQR-MF-PhiBE-ARE-K] is equivalent to \[\label{eq:K95hat95CARE95beta95explicit95proof} Q + \left(\hat{A}-\frac{\beta}{2}I_d\right)^\top\hat{K} + \hat{K}\left(\hat{A}-\frac{\beta}{2}I_d\right) - \hat{K}\hat{G}\hat{K} = 0.\tag{65}\] We now apply [42]. This theorem compares two algebraic Riccati equations of the form \[Q_0+A_0^\top X+XA_0-XG_0X=0\] and \[(Q_0+\Delta Q)+(A_0+\Delta A)^\top \hat{X}+\hat{X}(A_0+\Delta A)-\hat{X}(G_0+\Delta G)\hat{X}=0.\] For 64 , we take \(Q_0=Q\), \(A_0=A-\frac{\beta}{2}I_d\), and \(G_0=G\). Then 65 corresponds to the coefficient differences \[\Delta Q_K:=0, \qquad \Delta A_K:=\hat{A}-A, \qquad \Delta G_K:=\hat{G}-G.\] Let \(l_K,p_K,q_K\) be the corresponding perturbation constants in [42]. Set \[\varepsilon_K:=p_K\|\Delta A_K\|+q_K\|\Delta G_K\|, \qquad \delta_K:=\|\Delta A_K\|+\|\Delta G_K\|\,\|K\|, \qquad \hat{g}_K:=\|G+\Delta G_K\|.\] By 62 and 63 , we have \(\varepsilon_K\to0\) and \(\delta_K\to0\) as \(\Delta t\to0\). Hence, after reducing \(\Delta t_*\) if necessary, the smallness condition in [42] holds. Moreover, the denominator in the bound of [42] is \[l_K-2\delta_K+ \sqrt{(l_K-2\delta_K)^2-4l_K\hat{g}_K\varepsilon_K},\] which converges to \(2l_K\) as \(\Delta t\to0\). Reducing \(\Delta t_*\) again if necessary, this denominator is bounded from below by \(l_K\) for every \(\Delta t\in(0,\Delta t_*]\). Therefore [42] gives \[\label{eq:K95error95sun95prelemma} \|\hat{K}-K\| \le 2\left(p_K\|\Delta A_K\|+q_K\|\Delta G_K\|\right).\tag{66}\] By 8, \[p_K \le \frac{ 2\|K\|^3 \bigl(\beta+2\|A\|+\|G\|\|K\|\bigr) }{ \lambda_{\min}(Q)^2 },\] and \[q_K \le \frac{ \|K\|^4 \bigl(\beta+2\|A\|+\|G\|\|K\|\bigr) }{ \lambda_{\min}(Q)^2 }.\] Since \(\|K\|\le C_K\), \(\|G\|=\|BN^{-1}B^\top\|\), and 62 63 hold, 66 gives \[\label{eq:K95error95beta95explicit95proof} \|\hat{K}-K\| \le \frac{2}{\lambda_{\min}(Q)^2} \bigl(C_{K,1}+\beta C_{K,2}\bigr)\Delta t.\tag{67}\]

Step 3: perturbation of the mean-field Riccati equation. Since \(D=\bar D=F=I=0\), [eq:ARE95SIMPLIFY95RICCATI95LAMBDA] reduces to \[\beta\Lambda = Q+\bar Q + \Lambda\tilde{A}+\tilde{A}^\top\Lambda - \Lambda BN^{-1}B^\top\Lambda, \qquad \tilde{A}:=A+\bar A.\] Equivalently, \[\label{eq:Lambda95exact95CARE95beta95explicit95proof} Q+\bar Q + \left(\tilde{A}-\frac{\beta}{2}I_d\right)^\top\Lambda + \Lambda\left(\tilde{A}-\frac{\beta}{2}I_d\right) - \Lambda G\Lambda = 0.\tag{68}\] The LQR-MF-PhiBE Riccati equation [eq:LQR-MF-PhiBE-ARE-Lambda] is equivalent to \[\label{eq:Lambda95hat95CARE95beta95explicit95proof} Q+\bar Q + \left(\hat{A}+\hat{\bar A}-\frac{\beta}{2}I_d\right)^\top \hat{\Lambda} + \hat{\Lambda} \left(\hat{A}+\hat{\bar A}-\frac{\beta}{2}I_d\right) - \hat{\Lambda}\hat{G}\hat{\Lambda} = 0.\tag{69}\] We apply [42] again. For 68 , we take \(Q_0=Q+\bar Q\), \(A_0=A+\bar A-\frac{\beta}{2}I_d\), and \(G_0=G\). Then 69 corresponds to the coefficient differences \[\Delta Q_\Lambda:=0, \qquad \Delta A_\Lambda:=(\hat{A}-A)+(\hat{\bar A}-\bar A), \qquad \Delta G_\Lambda:=\hat{G}-G.\] By 62 and 63 , we have \(\|\Delta A_\Lambda\|\le 2C_{\tilde{A}}\Delta t\) and \[\|\Delta G_\Lambda\| \le \Bigl[ 2C_B\|N^{-1}\|\|B\| + C_B^2\|N^{-1}\|\Delta t_* \Bigr]\Delta t.\] Let \(l_\Lambda,p_\Lambda,q_\Lambda\) be the corresponding perturbation constants in [42]. Set \[\varepsilon_\Lambda:=p_\Lambda\|\Delta A_\Lambda\|+q_\Lambda\|\Delta G_\Lambda\|, \qquad \delta_\Lambda:=\|\Delta A_\Lambda\|+\|\Delta G_\Lambda\|\,\|\Lambda\|, \qquad \hat{g}_\Lambda:=\|G+\Delta G_\Lambda\|.\] The previous estimates give \(\varepsilon_\Lambda\to0\) and \(\delta_\Lambda\to0\) as \(\Delta t\to0\). Hence, after reducing \(\Delta t_*\) if necessary, the smallness condition in [42] holds. Moreover, the denominator in the bound of [42] is \[l_\Lambda-2\delta_\Lambda+ \sqrt{(l_\Lambda-2\delta_\Lambda)^2-4l_\Lambda\hat{g}_\Lambda\varepsilon_\Lambda},\] which converges to \(2l_\Lambda\) as \(\Delta t\to0\). Reducing \(\Delta t_*\) again if necessary, this denominator is bounded from below by \(l_\Lambda\) for every \(\Delta t\in(0,\Delta t_*]\). Therefore [42] gives \[\label{eq:Lambda95error95sun95prelemma} \|\hat{\Lambda}-\Lambda\| \le 2\left(p_\Lambda\|\Delta A_\Lambda\|+q_\Lambda\|\Delta G_\Lambda\|\right).\tag{70}\] By 8, \[p_\Lambda \le \frac{ 2\|\Lambda\|^3 \bigl(\beta+2\|A+\bar A\|+\|G\|\|\Lambda\|\bigr) }{ \lambda_{\min}(Q+\bar Q)^2 },\] and \[q_\Lambda \le \frac{ \|\Lambda\|^4 \bigl(\beta+2\|A+\bar A\|+\|G\|\|\Lambda\|\bigr) }{ \lambda_{\min}(Q+\bar Q)^2 }.\] Since \(\|\Lambda\|\le C_\Lambda\) and \(\|G\|=\|BN^{-1}B^\top\|\), using 62 and 63 in 70 gives \[\label{eq:Lambda95error95beta95explicit95proof} \|\hat{\Lambda}-\Lambda\| \le \frac{2}{\lambda_{\min}(Q+\bar Q)^2} \bigl(C_{\Lambda,1}+\beta C_{\Lambda,2}\bigr)\Delta t.\tag{71}\]

Step 4: estimate of the policy means. Define \[\label{eq:def95policy95mean95exact95lqr95phibe} m^*(s,\mu) := -N^{-1}B^\top \bigl(K(s-\bar\mu)+\Lambda\bar\mu\bigr), \qquad \hat{m}^*(s,\mu) := -N^{-1}\hat{B}^\top \bigl(\hat{K}(s-\bar\mu)+\hat{\Lambda}\bar\mu\bigr),\tag{72}\] then we have \[\begin{align} \label{eq:policy95difference95decomposition95beta95explicit95proof} \hat{m}^*(s,\mu)-m^*(s,\mu) &= N^{-1}(B^\top K-\hat{B}^\top\hat{K})(s-\bar\mu) + N^{-1}(B^\top\Lambda-\hat{B}^\top\hat{\Lambda})\bar\mu . \end{align}\tag{73}\] For the first coefficient, \[B^\top K-\hat{B}^\top\hat{K} = B^\top(K-\hat{K})+(B-\hat{B})^\top\hat{K}.\] Using 62 and 67 , \[\begin{align} \|B^\top K-\hat{B}^\top\hat{K}\| &\le \|B\|\,\|\hat{K}-K\| + \|\hat{B}-B\|\,\|\hat{K}\| \\ &\le 2\frac{\|B\|}{\lambda_{\min}(Q)^2} \bigl(C_{K,1}+\beta C_{K,2}\bigr)\Delta t + C_B\Delta t \bigl(\|K\|+\|\hat{K}-K\|\bigr) \\ &\le \left[ 2\frac{\|B\|}{\lambda_{\min}(Q)^2} \bigl(C_{K,1}+\beta C_{K,2}\bigr) + C_{K,3} \right]\Delta t + 2\frac{C_B}{\lambda_{\min}(Q)^2} \bigl(C_{K,1}+\beta C_{K,2}\bigr)\Delta t^2. \end{align}\] Similarly, \[B^\top\Lambda-\hat{B}^\top\hat{\Lambda} = B^\top(\Lambda-\hat{\Lambda})+(B-\hat{B})^\top\hat{\Lambda},\] and 71 gives \[\begin{align} \|B^\top\Lambda-\hat{B}^\top\hat{\Lambda}\| &\le \left[ 2\frac{\|B\|}{\lambda_{\min}(Q+\bar Q)^2} \bigl(C_{\Lambda,1}+\beta C_{\Lambda,2}\bigr) + C_{\Lambda,3} \right]\Delta t + 2\frac{C_B}{\lambda_{\min}(Q+\bar Q)^2} \bigl(C_{\Lambda,1}+\beta C_{\Lambda,2}\bigr)\Delta t^2. \end{align}\] Substituting the last two estimates into 73 and taking norms gives exactly ?? .

Step 5: cancellation of the averaged feedback for \(\beta=0\). Assume now that \(\beta=0\), \(d=m\), and \(B\) is invertible. After reducing \(\Delta t_*\) if necessary, \(\hat{B}\) is invertible for every \(\Delta t\in(0,\Delta t_*]\). We prove that the coefficient multiplying \(\bar\mu\) in the feedback mean is exactly preserved.

Let \(\widehat{\tilde{A}}:=\hat{A}+\hat{\bar A}\). From the definitions of \(\hat{A}\), \(\hat{\bar A}\), and \(\hat{B}\) in 18 , the variation-of-constants formula gives \[\label{eq:tildeA-hatB-identities-beta-zero} \widehat{\tilde{A}} = \frac{e^{\tilde{A}\Delta t}-I_d}{\Delta t}, \qquad \hat{B} = \frac{1}{\Delta t}M_{\tilde{A},\Delta t}B.\tag{74}\] Set \(\Phi_{\Delta t}:=\Delta t^{-1}M_{\tilde{A},\Delta t}\). Since \(\Phi_{\Delta t}\to I_d\) as \(\Delta t\to0\), we may reduce \(\Delta t_*\) so that \(\Phi_{\Delta t}\) is invertible. Then 74 gives \(\hat{B}=\Phi_{\Delta t}B\) and \(\widehat{\tilde{A}}=\tilde{A}\Phi_{\Delta t}\). Since \(\Phi_{\Delta t}\) is a matrix function of \(\tilde{A}\), it commutes with \(\tilde{A}\), and therefore \[\label{eq:structural-identity-Lambda-beta-zero} \hat{B}^{-1}\widehat{\tilde{A}} = B^{-1}\Phi_{\Delta t}^{-1}\tilde{A}\Phi_{\Delta t} = B^{-1}\tilde{A}.\tag{75}\]

Define \(M:=\Lambda B\) and \(\hat{M}:=\hat{\Lambda}\hat{B}\). Since \(\beta=0\), multiplying neither equation by an additional factor and using 75 , the exact equation 68 and the estimated equation 69 can be rewritten as the same algebraic equation \[\label{eq:M-Lambda-common-equation-beta-zero} Q+\bar Q + M D + D^\top M^\top - MN^{-1}M^\top = 0, \qquad D:=B^{-1}\tilde{A},\tag{76}\] with \(M\) replaced by \(M\) in the exact case and by \(\hat{M}\) in the estimated case.

It remains only to identify the stabilizing solution. The exact closed-loop matrix associated with \(M\) is \[\tilde{A}-BN^{-1}M^\top = B(D-N^{-1}M^\top),\] and it is Hurwitz because \(\Lambda\) is the stabilizing solution of [eq:ARE95SIMPLIFY95RICCATI95LAMBDA]. The estimated closed-loop matrix associated with the same \(M\) is \[\widehat{\tilde{A}}-\hat{B} N^{-1}M^\top = \hat{B}(D-N^{-1}M^\top),\] where we used 75 . Since \(\hat{B}\to B\) as \(\Delta t\to0\), this matrix is also Hurwitz for \(\Delta t_*\) sufficiently small, by openness of the set of Hurwitz matrices. Hence \(M\) is the stabilizing solution of 76 for the estimated problem. By uniqueness of the stabilizing solution, \(\hat{M}=M\), that is, \[\label{eq:exact-preservation-B-Lambda-beta-zero} \hat{\Lambda}\hat{B}=\Lambda B.\tag{77}\]

Finally, using 72 and the identity \(\int_{\mathbb{R}^d}(s-\bar\mu)\,\mu(ds)=0\), we obtain \[\mathbb{E}_{s\sim\mu}\bigl[m^*(s,\mu)\bigr] = -N^{-1}B^\top\Lambda\bar\mu = -N^{-1}M^\top\bar\mu,\quad\text{and}\quad \mathbb{E}_{s\sim\mu}\bigl[\hat{m}^*(s,\mu)\bigr] = -N^{-1}\hat{B}^\top\hat{\Lambda}\bar\mu = -N^{-1}\hat{M}^\top\bar\mu.\] By 77 , the two quantities coincide. This proves ?? when \(\beta=0\).

Step 6: Value function error. We now estimate the value gap between the optimal LQR policy and the LQR-MF-PhiBE policy, both evaluated in the original dynamics. Let \(L:=N^{-1}B^\top K\), \(\hat{L}:=N^{-1}\hat{B}^\top\hat{K}\), \(U:=N^{-1}B^\top\Lambda\), and \(\hat{U}:=N^{-1}\hat{B}^\top\hat{\Lambda}\). By 72 , the feedback means of \(\pi^*\) and \(\hat{\pi}^*\) are \[\begin{align} \label{eq:value-step-feedback-means} m^*(x,\nu)=-L(x-\bar\nu)-U\bar\nu, \qquad \hat{m}^*(x,\nu)=-\hat{L}(x-\bar\nu)-\hat{U}\bar\nu. \end{align}\tag{78}\] The two Gaussian policies have covariance \(\frac{\lambda}{2}N^{-1}\). Hence, using 14 and 3 , there exists a constant \(c_{\lambda,N}\), independent of \(\nu\), such that \[\begin{align} \label{eq:value-step-reward-pistar} r_\lambda^{\pi^*}(\nu) = -\operatorname{Tr}\bigl((Q+L^\top NL)\operatorname{Cov}(\nu)\bigr) -\bar\nu^\top(Q+\bar Q+U^\top NU)\bar\nu +c_{\lambda,N}, \end{align}\tag{79}\] and \[\begin{align} \label{eq:value-step-reward-pihat} r_\lambda^{\hat{\pi}^*}(\nu) = -\operatorname{Tr}\bigl((Q+\hat{L}^\top N\hat{L})\operatorname{Cov}(\nu)\bigr) -\bar\nu^\top(Q+\bar Q+\hat{U}^\top N\hat{U})\bar\nu +c_{\lambda,N}. \end{align}\tag{80}\]

Let \((\hat{\mu}_t)_{t\ge0}\) be the law flow generated by the original dynamics under \(\hat{\pi}^*\), and set \(\hat{C}_t:=\operatorname{Cov}(\hat{\mu}_t)\). From 78 , the mean and covariance satisfy \[\begin{align} \label{eq:value-step-flow-pihat} \dot{\bar{\hat{\mu}}}_t=(A+\bar A-B\hat{U})\bar{\hat{\mu}}_t, \qquad \dot{\hat{C}}_t=(A-B\hat{L})\hat{C}_t+\hat{C}_t(A-B\hat{L})^\top+\gamma\gamma^\top . \end{align}\tag{81}\] Moreover, by the explicit LQR value formula, \[\begin{align} \label{eq:value-step-optimal-value} \mathcal{V}^{\pi^*}(\nu) = -\operatorname{Tr}\bigl(K\operatorname{Cov}(\nu)\bigr) -\bar\nu^\top\Lambda\bar\nu -R. \end{align}\tag{82}\] Using 82 along the flow \((\hat{\mu}_t)_{t\ge0}\) and integrating by parts, we obtain \[\begin{align} \label{eq:value-step-integration-by-parts} \mathcal{V}^{\pi^*}(\mu) = \int_0^\infty e^{-\beta t} \left[ -\frac{d}{dt}\mathcal{V}^{\pi^*}(\hat{\mu}_t) +\beta\mathcal{V}^{\pi^*}(\hat{\mu}_t) \right]dt . \end{align}\tag{83}\] Therefore, by the definition of \(\mathcal{V}^{\hat{\pi}^*}\), \[\begin{align} \label{eq:value-step-gap-before-integrand} \mathcal{V}^{\pi^*}(\mu)-\mathcal{V}^{\hat{\pi}^*}(\mu) = \int_0^\infty e^{-\beta t} \left[ -\frac{d}{dt}\mathcal{V}^{\pi^*}(\hat{\mu}_t) +\beta\mathcal{V}^{\pi^*}(\hat{\mu}_t) -r_\lambda^{\hat{\pi}^*}(\hat{\mu}_t) \right]dt . \end{align}\tag{84}\]

We compute the integrand in 84 . Using 81 , 80 , and 82 , we get \[\begin{align} \label{eq:value-step-integrand-expanded} &-\frac{d}{dt}\mathcal{V}^{\pi^*}(\hat{\mu}_t) +\beta\mathcal{V}^{\pi^*}(\hat{\mu}_t) -r_\lambda^{\hat{\pi}^*}(\hat{\mu}_t)= \operatorname{Tr}\Bigl( \bigl[ (A-B\hat{L})^\top K+K(A-B\hat{L})-\beta K+Q+\hat{L}^\top N\hat{L} \bigr]\hat{C}_t \Bigr) \notag\\ &\qquad\qquad + \bar{\hat{\mu}}_t^\top \bigl[ (A+\bar A-B\hat{U})^\top\Lambda+\Lambda(A+\bar A-B\hat{U})-\beta\Lambda+Q+\bar Q+\hat{U}^\top N\hat{U} \bigr] \bar{\hat{\mu}}_t \notag\\ &\qquad\qquad -\beta R+\gamma^\top K\gamma-c_{\lambda,N}. \end{align}\tag{85}\] The scalar term in 85 vanishes by [eq:ARE95SIMPLIFY95RICCATI95R]. We next simplify the two matrix residuals in 85 . Since [eq:ARE95SIMPLIFY95RICCATI95K] is equivalent to \[\begin{align} \label{eq:value-step-K-residual-zero} (A-BL)^\top K+K(A-BL)-\beta K+Q+L^\top NL=0, \end{align}\tag{86}\] inserting \(K\) into the closed-loop equation associated with \(\hat{L}\) gives \[\begin{align} \label{eq:value-step-K-residual-expansion} (A-B\hat{L})^\top K+K(A-B\hat{L})-\beta K+Q&+\hat{L}^\top N\hat{L}= (A-BL)^\top K+K(A-BL)-\beta K+Q+L^\top NL \notag\\ & -(\hat{L}-L)^\top B^\top K-KB(\hat{L}-L)+\hat{L}^\top N\hat{L}-L^\top NL . \end{align}\tag{87}\] The first line on the right-hand side of 87 vanishes by 86 . Since \(L=N^{-1}B^\top K\), we have \(B^\top K=NL\) and \(KB=L^\top N\). Hence \[\begin{align} \label{eq:value-step-K-residual-quadratic} (A-B\hat{L})^\top K+K(A-B\hat{L})-\beta K+Q+\hat{L}^\top N\hat{L} = (\hat{L}-L)^\top N(\hat{L}-L). \end{align}\tag{88}\] Similarly, [eq:ARE95SIMPLIFY95RICCATI95LAMBDA] is equivalent to \[\begin{align} \label{eq:value-step-Lambda-residual-zero} (A+\bar A-BU)^\top\Lambda+\Lambda(A+\bar A-BU)-\beta\Lambda+Q+\bar Q+U^\top NU=0. \end{align}\tag{89}\] Inserting \(\Lambda\) into the closed-loop equation associated with \(\hat{U}\) and using \(U=N^{-1}B^\top\Lambda\), we obtain \[\begin{align} \label{eq:value-step-Lambda-residual-quadratic} (A+\bar A-B\hat{U})^\top\Lambda+\Lambda(A+\bar A-B\hat{U})-\beta\Lambda+Q+\bar Q+\hat{U}^\top N\hat{U} = (\hat{U}-U)^\top N(\hat{U}-U). \end{align}\tag{90}\] Substituting 88 and 90 into 85 , and then into 84 , yields \[\begin{align} \label{eq:value-step-quadratic-representation} \mathcal{V}^{\pi^*}(\mu)-\mathcal{V}^{\hat{\pi}^*}(\mu) &= \int_0^\infty e^{-\beta t} \operatorname{Tr}\bigl((\hat{L}-L)^\top N(\hat{L}-L)\hat{C}_t\bigr)dt + \int_0^\infty e^{-\beta t} \bar{\hat{\mu}}_t^\top(\hat{U}-U)^\top N(\hat{U}-U)\bar{\hat{\mu}}_t\,dt . \end{align}\tag{91}\] In particular, \(\mathcal{V}^{\pi^*}(\mu)-\mathcal{V}^{\hat{\pi}^*}(\mu)\ge0.\) It remains to estimate the right-hand side of 91 . From Step 4, we have already proved \[\begin{align} \label{eq:value-step-product-errors} \|B^\top K-\hat{B}^\top\hat{K}\|\le C\Delta t, \qquad \|B^\top\Lambda-\hat{B}^\top\hat{\Lambda}\|\le C\Delta t . \end{align}\tag{92}\] Therefore, by the definitions of \(L,\hat{L},U,\hat{U}\), \[\begin{align} \label{eq:value-step-L-U-errors} \|\hat{L}-L\|+\|\hat{U}-U\| &\le \|N^{-1}\| \left( \|\hat{B}^\top\hat{K}-B^\top K\| + \|\hat{B}^\top\hat{\Lambda}-B^\top\Lambda\| \right) \le C\Delta t . \end{align}\tag{93}\] Moreover, by the initial choice of \(\Delta t_*\) and 7, the matrices \(A-B\hat{L}\) and \(A+\bar A-B\hat{U}\) are uniformly Hurwitz for \(\Delta t\in(0,\Delta t_*]\). Hence the solutions of 81 satisfy, for some constants \(C,\kappa>0\) independent of \(\Delta t\), \[|\bar{\hat{\mu}}_t|^2\le Ce^{-2\kappa t}|\bar\mu|^2, \qquad \operatorname{Tr}(\hat{C}_t) \le Ce^{-2\kappa t}\operatorname{Tr}(\operatorname{Cov}(\mu))+C .\] The constant term in the second estimate comes from the diffusion term \(\gamma\gamma^\top\). Multiplying by \(e^{-\beta t}\) and integrating over \((0,\infty)\) gives \[\begin{align} \label{eq:value-step-second-moment-bound} \int_0^\infty e^{-\beta t} \left( \operatorname{Tr}(\hat{C}_t)+|\bar{\hat{\mu}}_t|^2 \right)dt \le C\bigl(1+m_2(\mu)\bigr). \end{align}\tag{94}\] Combining 91 , 93 , and 94 , we get \[\begin{align} 0 \le \mathcal{V}^{\pi^*}(\mu)-\mathcal{V}^{\hat{\pi}^*}(\mu) &\le C\left(\|\hat{L}-L\|^2+\|\hat{U}-U\|^2\right) \int_0^\infty e^{-\beta t} \left( \operatorname{Tr}(\hat{C}_t)+|\bar{\hat{\mu}}_t|^2 \right)dt\le C\Delta t^2\bigl(1+m_2(\mu)\bigr). \end{align}\] This proves ?? . ◻

Proof. Let \(\bar\mu_t:=\int_{\mathbb{R}^d}x\,\mu_t(dx)\) and set \(\xi_t:=s_t-\bar\mu_t\). Under ?? , \(\int_{\mathbb{R}^m}a\pi^*(da\mid x,\mu) = -N^{-1}B^\top\bigl(K(x-\bar\mu)+\Lambda\bar\mu\bigr).\) Since \(\sigma(x,\mu,a)=\gamma\), the closed-loop MKV dynamics is \[d s_t = \bigl[A_K(s_t-\bar\mu_t)+A_\Lambda\bar\mu_t\bigr]dt + \gamma\,dW_t,\] where \(A_K:=A-BN^{-1}B^\top K\) and \(A_\Lambda:=A+\bar A-BN^{-1}B^\top\Lambda.\) Taking expectation gives \(\dot{\bar\mu}_t=A_\Lambda\bar\mu_t\). Since \(\Lambda\) is the stabilizing solution of [eq:ARE95SIMPLIFY95RICCATI95LAMBDA], \(A_\Lambda\) is Hurwitz. Hence there exist \(C_\Lambda\ge1\) and \(\kappa_\Lambda>0\) such that \[|\bar\mu_t|^2 \le C_\Lambda e^{-2\kappa_\Lambda t}|\bar\mu_0|^2.\]

The estimate for \(\xi_t\) follows from the same quadratic energy argument used in [37]. Indeed, \(\xi_t\) solves \[d\xi_t=A_K\xi_t\,dt+\gamma\,dW_t,\] and the stabilizing property of \(K\) gives a Lyapunov function for \(A_K\). Applying Itô’s formula to this quadratic Lyapunov function and using Gronwall’s inequality yields constants \(C_K\ge1\) and \(\kappa_K>0\) such that \[\mathbb{E}|\xi_t|^2 \le C_K e^{-2\kappa_K t}\mathbb{E}|\xi_0|^2 + C_K|\gamma|^2.\] Finally, since \(m_2(\mu_t)=\mathbb{E}|s_t|^2=\mathbb{E}|\xi_t|^2+|\bar\mu_t|^2\) and \(m_2(\mu_0)=\mathbb{E}|\xi_0|^2+|\bar\mu_0|^2,\) taking \(\kappa:=\min\{\kappa_K,\kappa_\Lambda\}\) and increasing \(C\) if necessary gives \[m_2(\mu_t) \le C e^{-2\kappa t}m_2(\mu_0)+C|\gamma|^2.\] ◻