1 Introduction↩︎

Sampling from a target probability density \(\pi\), known only up to normalization, is a basic problem in Bayesian inference and computational statistics. At the population level, overdamped Langevin dynamics realizes the Wasserstein gradient flow of relative entropy through independent score-driven diffusions. The corresponding deterministic continuity equation has velocity \(-\nabla\log(\rho/\pi)\), which depends on the unknown evolving density \(\rho\) and therefore cannot be evaluated directly on an empirical measure.

Stein variational gradient descent (SVGD), introduced by [1], takes a different route. It projects this velocity onto a kernel space and uses integration by parts to remove \(\nabla\log\rho\). Write \(\pi=Z^{-1}e^{-V}\) on a \(d\)-dimensional state space, where \(Z\) is the normalizing constant, and let \(b_\pi:=\nabla\log\pi=-\nabla V\). For a smooth symmetric positive-definite kernel \(k\), the standard continuous-time particle system is \[\label{eq:standard-svgd-intro} \dot{X}_i(t) =\frac{1}{N}\sum_{j=1}^N \left\{\nabla_1 k(X_j(t),X_i(t)) +k(X_j(t),X_i(t))b_\pi(X_j(t))\right\}.\tag{1}\] Here \(\nabla_1\) differentiates the first kernel argument, and the sum includes self-interaction. The projected velocity depends on the target only through its score, so the normalizing constant is not needed. The score-weighted term transports the ensemble toward the target, while the kernel derivative supplies a repulsive interaction. Thus SVGD moves the empirical measure \(\mu_t^N:=N^{-1}\sum_{i=1}^N\delta_{X_i(t)}\) without estimating its density or choosing a parametric approximating family. It trades independent stochastic trajectories for a coordinated deterministic approximation of the target law.

Our concern is the long-time behavior of this coordinated particle transport. Unlike Langevin dynamics, the finite system has no diffusive mechanism enforcing ergodicity, and the particles cannot be analyzed as independent samples. The resulting question is how to obtain a many-particle, long-time sampling statement from deterministic interacting dynamics whose fixed-\(N\) trajectories need not converge to the target. The quantity connecting the dynamics to convergence is the kernel Stein discrepancy. If \(\mathcal{H}_k\) is the reproducing-kernel Hilbert space of \(k\), then \[\label{eq:standard-ksd-intro} \begin{align} \operatorname{KSD}_k^2(\mu\mid\pi) &:= \left\|\int \bigl\{\nabla_x k(x,\cdot)+k(x,\cdot)b_\pi(x)\bigr\} \,\,\mathrm d\mu(x)\right\|_{\mathcal{H}_k^d}^{2}. \end{align}\tag{2}\] It vanishes at the target and measures the size of the kernel-projected transport velocity. Along the population SVGD flow, it is exactly the rate at which relative entropy decreases. For \(P\ll Q\), write \[\label{eq:relative-entropy} \operatorname{Ent}(P\mid Q):=\int\log(\,\mathrm dP/\,\mathrm dQ)\,\,\mathrm dP.\tag{3}\] If \(\rho_t^k\) is the smooth-kernel population SVGD flow, then \[\label{eq:standard-population-dissipation-intro} \frac{\,\mathrm d}{\,\mathrm dt}\operatorname{Ent}(\rho_t^k\mid\pi) =-\operatorname{KSD}_k^2(\rho_t^k\mid\pi).\tag{4}\] The population gradient-flow theory was developed by [2] and [3], while [4] established the mean-field limit. For smooth kernels, [5] subsequently obtained longer-time stability estimates for particle systems initialized near the target. The KSD is also used by [6] to test whether a finite sample is consistent with a specified target distribution and by [7] to quantify how well samples approximate the target.

This perspective has led both to new analyses and to modifications of the algorithm. [8] introduced matrix-valued kernels that incorporate problem geometry through preconditioning, while [9] used random batches to reduce the cost of the all-to-all interaction. Kernel design can also change the underlying dissipation: [10] recast SVGD as a kernelized Wasserstein gradient flow of the chi-squared divergence and introduced Laplacian Adjusted Wasserstein Gradient Descent (LAWGD), whose target-adapted kernel represents the inverse Langevin generator and yields strong population convergence under a Poincaré inequality. The required spectral decomposition, however, makes LAWGD a different computational regime from the translation-invariant scalar kernels considered here.

For smooth kernels, the KSD can be evaluated directly on empirical measures, including the self-interaction terms. [11] established non-asymptotic descent estimates for population SVGD and quantitative bounds on its time-averaged KSD, [12] developed a variational description of the many-particle and long-time regimes, and [13] obtained the first explicit finite-particle rate for driving the KSD toward zero. Most closely related to our argument, [14] introduced the normalized joint-law entropy \[H_N^{\mathrm{sm}}(t) :=\frac{1}{N} \operatorname{Ent}\!\left(P_t^{N,\mathrm{sm}}\mid\pi^{\otimes N}\right),\] where \(P_t^{N,\mathrm{sm}}\) is the joint law of the smooth-kernel particle system. Their entropy calculation yields \[\frac{1}{T}\int_0^T \mathbb{E}\!\left[\operatorname{KSD}_k^2(\mu_t^N\mid\pi)\right]\,\,\mathrm dt \le \frac{H_N^{\mathrm{sm}}(0)}{T}+\frac{C_k}{N}.\] Consequently, for uniformly bounded initial normalized entropy, \[\lim_{T\to\infty}\limsup_{N\to\infty} \frac{1}{T}\int_0^T \mathbb{E}\!\left[\operatorname{KSD}_k^2(\mu_t^N\mid\pi)\right]\,\,\mathrm dt =0.\] Under their additional assumptions, this also yields convergence of time-averaged particle marginals. More recently, [15] proved propagation-of-chaos bounds, uniform over the averaging horizon, for time-averaged smooth-kernel SVGD; their uniform-in-physical-time parametric rates concern special finite-rank systems.

Smooth-kernel results show that time-averaged many-particle convergence can be obtained despite the absence of fixed-\(N\) ergodicity. At the population level, general results for standard smooth scalar kernels give qualitative last-iterate convergence and quantitative time-averaged or best-iterate KSD bounds, but no general quantitative rate is known for the last iterate in a strong topology [11], [16]. This contrast motivates kernels with finite-order Fourier decay.

On \(\mathbb{T}^d\), the mean-zero periodic Riesz kernel \(\kappa_a\) is characterized by \[\label{eq:kappa-fourier} \widehat\kappa_a(0)=0, \qquad \widehat\kappa_a(\ell)=(2\pi|\ell|)^{-2a} \quad(\ell\in\mathbb{Z}^d\setminus\{0\}),\tag{5}\] or equivalently by \((-\Delta)^a\kappa_a=\delta_0-1\). Its Fourier multiplier makes the associated Stein discrepancy coercive in a negative Sobolev norm, while the Riesz force provides short-range repulsion. [16] show that this reduced smoothing can yield quantitative population last-iterate convergence. Under their Sobolev regularity and small-initial-entropy assumptions, for every \(a>1\) and a regularity index \(\gamma>\max\{d/2,a-1\}\), \[\label{eq:population-polynomial-rate} \operatorname{Ent}(\rho_t\mid\pi) \le \left( \operatorname{Ent}(\rho_0\mid\pi)^{-(a-1)/\gamma}+\frac{t}{C} \right)^{-\gamma/(a-1)},\tag{6}\] together with corresponding decay in strong Sobolev norms; in particular, the representative entropy and squared \(L^2\) rates are \(t^{-\gamma/(a-1)}\). At the inverse-Laplacian endpoint \(a=1\), [17] establish global exponential convergence for a class of kernels with quadratic Fourier decay, while [16] recover the corresponding Coulomb result on the torus. In a different singular limit, [18] show that concentrating a regular kernel leads to a local quadratic-mobility flow.

The same finite-order Fourier structure creates the particle-level obstruction. Applying the Stein operator in both variables produces a scalar energy kernel with an infinite diagonal, whereas the singular particle dynamics excludes self-interaction. Thus the nonnegative population dissipation cannot simply be evaluated on an empirical measure. The issue is not singularity by itself: reduced smoothing strengthens population coercivity while simultaneously exposing the diagonal at finite \(N\).

1.1 The particle system and its Stein energy↩︎

We now take the state space to be the flat torus \(\mathbb{T}^d\) and specialize to \(k(y,x)=\kappa_a(y-x)\). The first-order Stein operator is \[\label{eq:A-intro} \mathcal{A}_{\pi,y}k(y,x) :=\nabla_y k(y,x)+b_\pi(y)k(y,x).\tag{7}\] The Riesz interaction is not smooth on the diagonal, and its self-interaction is not consistently defined across the range considered below. We therefore study the self-interaction-free system \[\label{eq:particle-system} \dot{X}_i = \frac{1}{N}\sum_{j\ne i}\mathcal{A}_{\pi,X_j}k(X_j,X_i) = \frac{1}{N}\sum_{j\ne i} \left\{\nabla\kappa_a(X_j-X_i) +\kappa_a(X_j-X_i)b_\pi(X_j)\right\}.\tag{8}\] Write \(\boldsymbol{X}(t):=(X_1(t),\ldots,X_N(t))\). For a configuration \(\boldsymbol{x}=(x_1,\ldots,x_N)\), let \[\mu_{\boldsymbol{x}}^N:=\frac{1}{N}\sum_{i=1}^N\delta_{x_i}\] be its empirical measure, and let \[\mathcal{D}_N:=\{\boldsymbol{x}\in(\mathbb{T}^d)^N:x_i\ne x_j\text{ for }i\ne j\}\] be the collision-free configuration space. The first term is the repulsive kernel force; the second transports that interaction through the target score. With \(\kappa_a\) defined by 5 , we work in the finite-energy regime \[\label{eq:parameter-range-intro} 1<a<1+\frac{d}{2}, \qquad \sigma:=d+2-2a\in(0,d).\tag{9}\]

To connect the particle velocity with entropy dissipation, we use divergence relative to the target measure. For a vector field \(u\), this is \[\mathcal{S}_{\pi,x}u :=\frac{1}{\pi(x)}\nabla_x\cdot\bigl(\pi(x)u\bigr) =\nabla_x\cdot u+b_\pi(x)\cdot u.\] Applying this operator to the velocity generated at \(x\) by a source point \(y\) defines the scalar Stein kernel \[\label{eq:G-intro} G_\pi(x,y):=\mathcal{S}_{\pi,x}\mathcal{A}_{\pi,y}k(y,x).\tag{10}\] Thus \(G_\pi(x,y)\) measures the contribution of the interaction between \(x\) and \(y\) to compression of the flow relative to the target. For \(x\ne y\), a direct calculation gives \[\label{eq:G-expansion} \begin{align} G_\pi(x,y) ={}&-\Delta\kappa_a(y-x) +\bigl(b_\pi(x)-b_\pi(y)\bigr)\cdot\nabla\kappa_a(y-x)\\ &+b_\pi(x)\cdot b_\pi(y)\,\kappa_a(y-x). \end{align}\tag{11}\]

At the population level, the quadratic form generated by \(G_\pi\) is nonnegative. Indeed, for a smooth probability measure \(\mu\), \[\label{eq:F-cont} \begin{align} \operatorname{KSD}_{\kappa_a}^2(\mu\mid\pi) &:= \left\| \int_{\mathbb{T}^d}\mathcal{A}_{\pi,y}k(y,\cdot)\,\,\mathrm d\mu(y) \right\|_{\mathcal{H}_{\kappa_a}^d}^{2}\\ &=\iint_{\mathbb{T}^d\times\mathbb{T}^d}G_\pi(x,y) \,\,\mathrm d(\mu-\pi)(x)\,\,\mathrm d(\mu-\pi)(y) =:2\mathcal{F}^\pi(\mu). \end{align}\tag{12}\] We call \(\mathcal{F}^\pi\) the Stein energy. The space \(\mathcal{H}_{\kappa_a}\) is the native Fourier space of \(\kappa_a\), defined in 22 . Along the population SVGD flow, this energy is exactly the entropy dissipation: \[\label{eq:population-dissipation-intro} \frac{\,\mathrm d}{\,\mathrm dt}\operatorname{Ent}(\rho_t\mid\pi) =-\operatorname{KSD}_{\kappa_a}^2(\rho_t\mid\pi) =-2\mathcal{F}^\pi(\rho_t).\tag{13}\] Identity 12 is proved for smooth measures in 5, and the population entropy-dissipation identity 13 is proved in 11. Both arguments use Fourier duality and therefore do not require pointwise reproduction for the singular kernel.

The role of \(G_\pi\) at the particle level is analogous. When the Liouville equation for the joint particle law is tested against relative entropy with respect to \(\pi^{\otimes N}\), each ordered pair contributes \(G_\pi(x_i,x_j)\). The leading term \(-\Delta\kappa_a(y-x)\) in 11 diverges on the diagonal, so the full empirical Stein energy is undefined. Since the particle system omits self-interaction, its entropy production instead involves the off-diagonal empirical Stein energy \[\label{eq:F-discrete} \mathcal{F}_N^\pi(\boldsymbol{x}) :=\frac{1}{2N^2}\sum_{i\ne j}G_\pi(x_i,x_j).\tag{14}\]

1.2 The singular entropy argument↩︎

The smooth-kernel entropy method suggests applying the same normalized joint-law entropy to the off-diagonal Riesz system. Let \(P_t^N\) be its joint particle law and write \[H_N(t):=\frac{1}{N}\operatorname{Ent}(P_t^N\mid\pi^{\otimes N}).\] Removing self-interaction gives the off-diagonal identity \[\label{eq:entropy-intro} \begin{align} \frac{\,\mathrm d}{\,\mathrm dt}H_N(t) &=-\frac{1}{N^2}\mathbb{E}_{P_t^N}\!\left[\sum_{i\ne j} G_\pi(X_i(t),X_j(t))\right]\\ &=-2\mathbb{E}_{P_t^N}\!\left[\mathcal{F}_N^\pi(\boldsymbol{X}(t))\right]. \end{align}\tag{15}\] The problem is to recover the nonnegativity lost by deleting an infinite diagonal.

For comparison, modulated-energy methods compare an empirical measure with a continuum density without evaluating the singular energy on the particle diagonal. The framework was developed for Coulomb-type flows by [19] and extended to translation-invariant Riesz flows by [20], using renormalization and smearing ideas related to [21]. The resulting corrected modulated energy and commutator estimate give finite-time mean-field convergence for Riesz flows. In SVGD, however, the velocity decomposes as \[\label{eq:velocity-intro} v_\mu=-\nabla\kappa_a*\mu+\kappa_a*(b_\pi\mu).\tag{16}\] The first term has the commutator structure treated by that theory within its admissible range; the target-dependent second term creates a discrete empirical trace for which the same estimate is unavailable. This prevents a direct application of existing Riesz-flow theory to a last-iterate particle limit. The time-averaged argument below avoids such a particle-to-population comparison.

Our argument combines the joint-law entropy identity from smooth-kernel SVGD with a direct capped-diagonal passage from particle energies to the continuum KSD. Making this work for the target-dependent, self-interaction-free Riesz system requires three steps:

1. Singular entropy identity. We first prove global collision avoidance and then justify the joint entropy calculation for the singular flow.

2. Capped-diagonal empirical liminf. Capping the full Stein kernel below its singular diagonal gives a direct lower bound from the off-diagonal particle energy to the continuum KSD. This proves that the negative part is bounded by a correction \(c_N\to0\) and permits simultaneous averaging horizons.

3. Quantitative correction below the logarithmic threshold. When \(0<\sigma<2\), an operator factorization and a positive Bessel minorant give the explicit rate \(c_N\lesssim N^{-1+\sigma/d}\).

The resulting finite-\(N\) renormalized entropy method is the central contribution: it converts the singular off-diagonal entropy production into a nonnegative continuum discrepancy up to a vanishing \(N\)-dependent correction. Because the argument works directly with the joint particle law, it does not require a propagation-of-chaos comparison with the population flow.

1.3 Main result and proof mechanism↩︎

The exact off-diagonal identity 15 does not yet provide dissipation because \(\mathcal{F}_N^\pi\) need not be nonnegative. Set \[\label{eq:onsager-intro} c_N:=\left(-\inf_{\boldsymbol{x}\in\mathcal{D}_N}\mathcal{F}_N^\pi(\boldsymbol{x})\right)_+.\tag{17}\] The capped-diagonal liminf proves, throughout \(0<\sigma<d\), that \(c_N\to0\). Therefore the renormalized energy \(\mathcal{E}_N^\pi:=\mathcal{F}_N^\pi+c_N\) is nonnegative. Integrating 15 and using \(H_N(T)\ge0\) gives \[\label{eq:rd-intro} \frac{1}{T}\int_0^T \mathbb{E}_{P_t^N}\!\left[\mathcal{E}_N^\pi(\boldsymbol{X}(t))\right]\,\,\mathrm dt \le \frac{H_N(0)}{2T}+c_N.\tag{18}\] For \(0<\sigma<2\), the quantitative refinement gives \[\label{eq:quantitative-onsager-intro} c_N\le C_{\rm ons}N^{-1+\sigma/d}.\tag{19}\] Thus 18 is the singular replacement for the finite-diagonal correction in the smooth-kernel entropy method, with an explicit rate for \(0<\sigma<2\).

The entropy identity is therefore naturally suited to time-averaged convergence. Let \(P_{t,i}^N\) be the \(i\)-th marginal of \(P_t^N\) and set \[\label{eq:tagged-average-intro} \overline{P}_{\mathrm{av},T}^{N} :=\frac{1}{T}\int_0^T\frac{1}{N}\sum_{i=1}^N P_{t,i}^N\,\,\mathrm dt.\tag{20}\] Writing \(d_{\mathrm{BL}}\) for the bounded-Lipschitz distance, we prove that whenever \(T_N\to\infty\) and \(H_N(0)/T_N\to0\), the time-averaged empirical-measure law converges to \(\delta_\pi\), the Dirac mass at the point \(\pi\in\mathcal{P}(\mathbb{T}^d)\): \[\label{eq:ordered-result-intro} \frac{1}{T_N}\int_0^{T_N} (\boldsymbol{x}\mapsto\mu_{\boldsymbol{x}}^N)_\#P_t^N\,\,\mathrm dt \rightharpoonup\delta_\pi, \qquad d_{\mathrm{BL}}\!\left(\overline{P}_{\mathrm{av},T_N}^{N},\pi\right)\to0.\tag{21}\] The first convergence therefore takes place in \(\mathcal{P}(\mathcal{P}(\mathbb{T}^d))\). In particular, every sequence \(T_N\to\infty\) is admissible under a uniform initial normalized-entropy bound. No exchangeability is required. We also prove that invariant laws of finite relative entropy have empirical-measure laws converging to \(\delta_\pi\), without any uniform entropy bound. Exchangeability is needed only to conclude that their first marginals converge to \(\pi\).

The proof of 21 has four conceptual steps. First, the Riesz repulsion prevents collisions, so the singular flow and its Liouville equation are globally defined. Second, the joint-entropy identity gives an upper bound on the averaged off-diagonal energy. Third, a bounded continuous cap of \(G_\pi\) gives a lower bound for every sequence of empirical-measure laws, including laws averaged over \(T_N\). Finally, compactness and the coercive identity \(\mathcal{F}^\pi(\mu)=0\) if and only if \(\mu=\pi\) identify the limit. The Bessel-minorant argument is an additional quantitative refinement and is not needed for 21 .

2 Setting and theorem statements↩︎

2.1 The torus, the target, and the functional setting↩︎

We identify \(\mathbb{T}^d\) with \([-\tfrac12,\tfrac12)^d\) and use normalized Lebesgue measure. Fourier coefficients are \[\widehat f(\ell)=\int_{\mathbb{T}^d}f(x)e^{-2\pi i\ell\cdot x}\,\,\mathrm dx.\] We use the bounded-Lipschitz metric \[d_{\mathrm{BL}}(\mu,\nu) := \sup_{\|f\|_\infty+\operatorname{Lip}(f)\le1} \left|\int_{\mathbb{T}^d}f\,\,\mathrm d(\mu-\nu)\right|,\] which metrizes weak convergence on \(\mathcal{P}(\mathbb{T}^d)\).

The Sobolev assumptions in 1 are used in the continuum identification of the Riesz KSD with a negative Sobolev norm and in the operator estimates for \(0<\sigma<2\) in 2. Collision avoidance and the entropy calculation use only the bounded derivatives of the target score displayed in their proofs.

The Riesz kernel \(\kappa_a\) was defined in 5 . Its native space is the mean-zero Fourier space \[\label{eq:native-space} \mathcal{H}_{\kappa_a} :=\left\{f:\widehat f(0)=0,\; \|f\|_{\mathcal{H}_{\kappa_a}}^2 :=\sum_{\ell\ne0}\widehat\kappa_a(\ell)^{-1} |\widehat f(\ell)|^2<\infty\right\} =\dot{H}^a_0(\mathbb{T}^d).\tag{22}\] If \(a>d/2\), point evaluation is continuous and this is the usual mean-zero reproducing-kernel Hilbert space generated by \(\kappa_a\). The range considered here may also contain \(a\le d/2\), where point evaluation on \(\dot{H}^a\) is not continuous and \(\kappa_a\) is singular on the diagonal. The literal pointwise reproducing property is then unavailable. We therefore use 22 through Fourier duality; no point evaluation of the singular kernel is needed.

The main finite-particle theorems assume 9 . In this range the kernel is smooth away from the origin, its force is less singular than the Coulomb force \((a=1,\sigma=d)\), and the scalar kernel \(G_\pi\) has locally integrable singularity \(\operatorname{dist}(x,y)^{-\sigma}\).

The identity 12 initially defines \(\mathcal{F}^\pi\) for smooth measures. Later we apply it to weak limits of empirical measures, which need not have densities. We therefore use \(\mathcal{F}^\pi\) for the lower-semicontinuous relaxation of the smooth Gram energy to \(\mathcal{P}(\mathbb{T}^d)\), allowing the value \(+\infty\). The diagonal-cap construction in 6 gives a concrete representation of this relaxation.

2.2 Entropy dissipation and long-time sampling↩︎

Recall from the introduction that \(\mathcal{D}_N\) denotes the collision-free configuration space on which the particle vector field is defined.

The first theorem establishes the analytic foundation for the entropy method. Its three conclusions have distinct roles: collision avoidance makes the singular flow global; the Liouville calculation identifies its exact entropy production; and a qualitative correction turns that sign-indefinite off-diagonal production into a nonnegative energy up to a vanishing finite-\(N\) error.

For a smooth kernel, the analogue of the third conclusion follows by adding the finite diagonal of the Gram matrix. Here that diagonal is infinite. The correction in 24 is shown to vanish by the direct capped-diagonal liminf of 9; no smearing estimate is needed for this qualitative conclusion.

Below the logarithmic threshold \(\sigma=2\), the correction has an explicit natural-scale bound.

The next theorem converts the averaged energy bound into concentration of the time-averaged empirical-measure law at the target, and hence convergence of its label-averaged marginal.

The theorem concerns the full law of the empirical measure and therefore also the label-averaged marginal. It does not require the joint law to be exchangeable. Condition 29 shows that a uniform entropy bound is sufficient but not necessary.

At stationarity, the entropy identity forces the mean off-diagonal dissipation to vanish. Passing this fact to the continuum yields the following consequence.

3 Stein–Riesz structure↩︎

This section establishes the two properties of \(G_\pi\) used throughout the proof. Locally, it has a positive Riesz singularity; this drives collision avoidance and permits a monotone cap of the diagonal. Globally, it is a centered Gram kernel; this makes the continuum energy nonnegative and identifies it as a KSD. The finite-particle difficulty is precisely that these global properties do not make sense on the infinite diagonal.

3.1 Local singularity↩︎

The proposition shows that the target dependence does not alter the leading singularity. The positivity in ?? will be used twice: first to repel collapsing clusters, and later to cap the singular diagonal from below without losing an empirical lower bound.

3.2 Stein centering and positivity↩︎

Local positivity controls the singularity but does not yet explain why the continuum energy measures discrepancy from \(\pi\). That information comes from the global Stein identities.

4 Collision-free particle dynamics and Liouville structure↩︎

Before differentiating the joint entropy, we must know that the singular ODE is globally well-posed. The only possible finite-time obstruction is a collision. Below the logarithmic threshold, the radius of a collapsing cluster has a strictly positive derivative at each sufficiently small running minimum. At and above the threshold, the Riesz pair energy diverges at collision and is decreasing near the collision set. Once collisions are excluded, the Liouville calculation can be performed on the collision-free configuration space.

4.1 Cluster repulsion below the logarithmic threshold↩︎

4.2 A pair-energy barrier at and above the logarithmic threshold↩︎

Define \[\label{eq:pair-energy} \mathcal{W}_N(\boldsymbol{x}):=\sum_{1\le i<j\le N}\kappa_a(x_j-x_i)\tag{36}\] and let \(\boldsymbol{C}^N=(C_1^N,\ldots,C_N^N)\), where \[C_i^N(\boldsymbol{x}) :=\frac{1}{N}\sum_{j\ne i}\kappa_a(x_j-x_i)b_\pi(x_j).\] Then the particle system can be written as \[\label{eq:pair-energy-flow} \dot{\boldsymbol{X}}=-\frac{1}{N}\nabla\mathcal{W}_N(\boldsymbol{X})+\boldsymbol{C}^N(\boldsymbol{X}).\tag{37}\]

Having ruled out collisions, we pass from individual trajectories to the law of the full particle configuration. The appropriate equation is the Liouville equation, the continuity equation for a probability density transported by a deterministic vector field. This is useful here because differentiating relative entropy along a continuity equation produces the divergence of that vector field relative to the reference measure \(\pi^{\otimes N}\). The next identity recognizes this relative divergence as the off-diagonal Stein energy.

4.3 Liouville equation↩︎

From this point on, denote the right-hand side of 8 by \[B_i^N(\boldsymbol{X}) :=\frac{1}{N}\sum_{j\ne i}\mathcal{A}_{\pi,X_j}k(X_j,X_i), \qquad \dot{X}_i=B_i^N(\boldsymbol{X}).\]

Thus \(G_\pi\) is exactly the pairwise contribution to the compressibility of the joint flow relative to \(\pi^{\otimes N}\). Identity 39 is what turns the Liouville transport equation into the entropy-dissipation identity in 8.

If \(P_t^N\) has density \(p_t^N\), it solves \[\label{eq:liouville} \partial_t p_t^N+\operatorname{div}_{\boldsymbol{x}}(p_t^N\boldsymbol{B}^N)=0\tag{40}\] on \(\mathcal{D}_N\). No boundary condition is imposed at collisions: by 7, the flow never reaches that boundary.

5 Joint entropy and Onsager renormalization↩︎

We now turn the deterministic flow into a dissipative estimate for its joint law. The calculation has two parts. The Liouville equation gives an exact identity with \(\mathcal{F}_N^\pi\) as entropy production. Since this off-diagonal quantity can be negative, the direct empirical liminf in 9 supplies a qualitative vanishing correction. For \(0<\sigma<2\), the Bessel-minorant argument of 8 gives its explicit rate.

5.1 Exact entropy production↩︎

The entropy identity below is the structural ingredient inherited from the smooth-kernel joint-law method. The work here is to justify it despite the singular collision set; the no-collision theorem and the approximation in 12 provide that justification.

5.2 Renormalized dissipation↩︎

The entropy identity alone is insufficient because its production term is the off-diagonal energy, not the nonnegative continuum KSD. The empirical liminf in 6 supplies the vanishing correction throughout \(0<\sigma<d\); 8 later gives its explicit rate for \(0<\sigma<2\).

6 Passage to continuum KSD↩︎

Estimate 25 controls a scalar particle energy. To obtain convergence of measures, we must show that this energy cannot disappear in the many-particle limit. We consider laws of the empirical measure \(\mu_{\boldsymbol{x}}^N\) defined in the introduction. Compactness of \(\mathcal{P}(\mathbb{T}^d)\) gives subsequential limits automatically. The work is to pass the singular off-diagonal energy to those limits. We do this by capping \(G_\pi\) at a finite height, passing the resulting bounded continuous kernel to the limit, and then raising the cap. The finite diagonal introduced by the cap costs only \(L/(2N)\). Since this comparison holds for every configuration law, it also applies when the averaging horizon depends on \(N\).

6.1 A capped-diagonal representation↩︎

By 4, \(G_\pi(x,y)\to+\infty\) uniformly as \(\operatorname{dist}(x,y)\to0\). Since it is continuous away from the diagonal, it is also bounded below on \(\mathbb{T}^d\times\mathbb{T}^d\setminus\operatorname{diag}\), where \(\operatorname{diag}:=\{(x,x):x\in\mathbb{T}^d\}\); fix \(C_0\) such that \(G_\pi\ge-C_0\). For \(L\ge1\), define the capped kernel \[G_\pi^{(L)}(x,y) := \begin{cases} \min\{G_\pi(x,y),L\},&x\ne y,\\ L,&x=y, \end{cases}\] and \[\label{eq:F-cap} \mathcal{F}_L^\pi(\mu) :=\frac{1}{2}\iint G_\pi^{(L)}(x,y)\,\,\mathrm d\mu(x)\,\,\mathrm d\mu(y).\tag{41}\] The value on the diagonal makes \(G_\pi^{(L)}\) bounded and continuous on the full product space. Thus \(\mathcal{F}_L^\pi\) is continuous under weak convergence. The family increases with \(L\), and \(\mathcal{F}_L^\pi\ge-C_0/2\).

6.2 Microscopic-to-continuum lower bound↩︎

The next proposition is the bridge from particle dissipation to continuum dissipation. Its lower-bound direction is especially important: it prevents microscopic concentration near the singular diagonal from creating an artificially small limiting energy.

6.3 Time-averaged empirical laws↩︎

For \(T>0\), define \[\label{eq:Lambda-NT} \Lambda_{N,T} :=\frac{1}{T}\int_0^T (\boldsymbol{x}\mapsto\mu_{\boldsymbol{x}}^N)_\#P_t^N\,\,\mathrm dt.\tag{44}\] Its barycenter is exactly the measure in 20 : \[\label{eq:barycenter} \int_{\mathcal{P}(\mathbb{T}^d)}\mu\,\,\mathrm d\Lambda_{N,T}(\mu) =\overline{P}_{\mathrm{av},T}^{N}.\tag{45}\]

7 Target identification and long-time sampling↩︎

The preceding section passes the off-diagonal particle energy to every weak limit of empirical-measure laws. It remains to identify the zero set of the continuum energy. The Riesz Fourier multiplier identifies \(\mathcal{F}^\pi\) with a negative Sobolev norm of \(\mu-\pi\), so \(\pi\) is its unique zero.

7.1 The zero set of the Riesz KSD↩︎

The following coercive identification is the point at which the Riesz Fourier structure matters globally.

7.2 Proofs of the long-time conclusions↩︎

8 Quantitative renormalization below the logarithmic threshold↩︎

The qualitative correction \(c_N\to0\) follows from compactness and applies throughout \(0<\sigma<d\). For \(0<\sigma<2\), the target-dependent part of the Stein kernel can also be separated from a positive translation-invariant minorant. This gives the explicit rate in 2. Throughout this section, assume \[0<\sigma<2, \qquad s:=d-\sigma=2a-2\in(0,d).\] The operator factorization and Bessel minorant below apply throughout this range.

8.1 Normalized Stein operator↩︎

Let \(\mathsf L=(-\Delta)^{1/2}\) on mean-zero periodic distributions, with multiplier \(\lambda_\ell:=2\pi|\ell|\) for \(\ell\ne0\), and let \(\Pi\) denote the mean-zero projection, componentwise for vector fields. Set \[D_{b_\pi}:=\nabla-b_\pi, \qquad T_\pi:=D_{b_\pi}^*\mathsf L^{-2a}D_{b_\pi}.\] The inverse power is the mean-zero pseudoinverse. The Gram calculation in 5 can equivalently be written, for smooth signed \(f\), as \[\label{eq:operator-form} \iint G_\pi(x,y)f(x)f(y)\,\,\mathrm dx\,\,\mathrm dy =\left\langle D_{b_\pi}f,\mathsf L^{-2a}D_{b_\pi}f\right\rangle.\tag{47}\] For \(b_\pi=0\), the corresponding operator is \[T_0=\mathsf L^{2-2a}=\mathsf L^{-s}.\]

We normalize the principal part to the identity. On the mean-zero subspace \(L_0^2(\mathbb{T}^d)\), define \[\mathcal{R}:=\nabla\mathsf L^{-1}, \qquad \mathcal{C}:=\mathsf L^{-a}\Pi M_{b_\pi}\mathsf L^{a-1}.\] Here \(M_{b_\pi}\) denotes multiplication by the vector field \(b_\pi\). Since \(\mathcal{R}^*\mathcal{R}=I\) on \(L_0^2\), \[\label{eq:normalized-B} \mathcal{B} :=\mathsf L^{s/2}T_\pi\mathsf L^{s/2} =(\mathcal{R}-\mathcal{C})^*(\mathcal{R}-\mathcal{C}) =I+\mathcal{K}.\tag{48}\] The next lemma shows that the normalized target perturbation gains two derivatives as a quadratic form.

8.2 A positive Bessel minorant↩︎

For \(M\ge1\), let \(J_M\) be the periodic Bessel kernel with Fourier coefficients \[\label{eq:Bessel-multiplier} \widehat J_M(\ell) =\bigl(M^2+\lambda_\ell^2\bigr)^{-s/2}, \qquad \ell\in\mathbb{Z}^d,\tag{54}\] where \(\lambda_0=0\). Its heat representation is \[\label{eq:Bessel-heat} J_M(z) =\frac{1}{\Gamma(s/2)} \int_0^\infty t^{s/2-1}e^{-M^2t}p_t(z)\,\,\mathrm dt.\tag{55}\] The periodic heat kernel is nonnegative, so \(J_M\ge0\). Moreover, \[J_M\in L^1(\mathbb{T}^d), \qquad \int_{\mathbb{T}^d}J_M\,\,\mathrm dx=M^{-s}.\]

Center this kernel with respect to \(\pi\): \[\label{eq:centered-Bessel} \overline{J}_M^\pi(x,y) :=J_M(x-y)-u_M(x)-u_M(y)+\gamma_M,\tag{56}\] where \[u_M:=J_M*\pi, \qquad \gamma_M:=\int_{\mathbb{T}^d}u_M\,\,\mathrm d\pi.\] Then \(\overline{J}_M^\pi\) is centered in each variable, \[0\le u_M(x)\le\|\pi\|_\infty M^{-s}, \qquad \gamma_M\ge0,\] and \(u_M\) is continuous because \(J_M\in L^1\) and \(\pi\) is continuous.

8.3 A continuous positive-semidefinite remainder↩︎

In this subsection and the next, assume \(0<\sigma<2\).

For \(x\ne y\), define \[\label{eq:KM-definition} K_M^\pi(x,y):=G_\pi(x,y)-\overline{J}_M^\pi(x,y).\tag{62}\] Both terms are singular on the diagonal; the next lemma shows that their difference has a continuous extension.

8.4 Completion of the quantitative estimate↩︎

9 Open directions↩︎

Three extensions remain particularly natural.

1. The qualitative correction \(c_N\to0\) holds for the full range \(0<\sigma<d\), but an explicit rate in the range \(2\le\sigma<d\) remains open. At and above the endpoint, the post-cancellation remainder no longer has a finite diagonal. The ball-growth and smearing methods of [20], [21] control the translation-invariant principal Riesz interaction, but the target-dependent remainder is variable-coefficient and sign-indefinite. An explicit rate would require a corresponding commutator or pseudodifferential truncation estimate.

2. A last-iterate result on logarithmic time scales would follow from a singular propagation-of-chaos estimate with controlled time dependence. Existing Riesz commutator estimates address the translation-invariant term \(-\nabla\kappa_a*(\mu^N-\rho)\) in their admissible range, but not the discrete trace generated by the target-weighted term \(\kappa_a*(b_\pi(\mu^N-\rho))\). Controlling that trace would connect the population decay of [16] to a particle last-iterate limit.

3. It would be useful to admit deterministic atomic initialization. Such initial laws have infinite entropy relative to \(\pi^{\otimes N}\), so the joint-entropy argument would need an initial-layer estimate or a replacement that measures a finite defect.

An extension to \(\mathbb{R}^d\) is also natural, but it brings separate questions of tightness, confinement, and long-range control.

Acknowledgments↩︎

The authors used OpenAI’s ChatGPT during the preparation of this manuscript to assist with mathematical exploration, proof checking, and improvements to the language and readability. The authors reviewed and verified all resulting content and take full responsibility for the manuscript.

10 Periodic Riesz asymptotics↩︎

10.1 The pair-energy gradient near collisions↩︎

The next estimate rules out cancellation of the singular pair forces inside a multiscale cluster. It is the input used in 2.

11 Population entropy dissipation for the singular kernel↩︎

We justify 13 without using pointwise reproduction or evaluating the singular kernel on the diagonal. Let \((\rho_t)_{t\in[0,T]}\) be a smooth, strictly positive solution of \[\label{eq:population-flow} \partial_t\rho_t+\nabla\cdot(\rho_t v_t)=0, \qquad v_t(x):= \int_{\mathbb{T}^d}\mathcal{A}_{\pi,y}k(y,x)\rho_t(y)\,\,\mathrm dy.\tag{72}\] Set \[q_t:=b_\pi\rho_t-\nabla\rho_t =-\rho_t\nabla\log\frac{\rho_t}{\pi}.\] Because \(\rho_t\) is smooth, the convolution of the distribution \(\kappa_a\) with \(q_t\) is smooth. Distributional integration by parts in the source variable gives \[v_t =\int_{\mathbb{T}^d}\mathcal{A}_{\pi,y}k(y,\cdot)\rho_t(y)\,\,\mathrm dy =\kappa_a*q_t.\] Differentiating the relative entropy along 72 and integrating by parts on the torus therefore yields \[\begin{align} \frac{\,\mathrm d}{\,\mathrm dt}\operatorname{Ent}(\rho_t\mid\pi) &=\int_{\mathbb{T}^d}\log\frac{\rho_t}{\pi}\,\partial_t\rho_t\,\,\mathrm dx\\ &=\int_{\mathbb{T}^d}\rho_t\nabla\log\frac{\rho_t}{\pi}\cdot v_t\,\,\mathrm dx\\ &=-\int_{\mathbb{T}^d}q_t\cdot(\kappa_a*q_t)\,\,\mathrm dx. \end{align}\] The Fourier identity established in the proof of 5 gives \[\int_{\mathbb{T}^d}q_t\cdot(\kappa_a*q_t)\,\,\mathrm dx =\left\| \int_{\mathbb{T}^d}\mathcal{A}_{\pi,y}k(y,\cdot)\,\,\mathrm d\rho_t(y) \right\|_{\mathcal{H}_{\kappa_a}^d}^{2} =\operatorname{KSD}_{\kappa_a}^2(\rho_t\mid\pi) =2\mathcal{F}^\pi(\rho_t).\] This proves 13 .

12 Entropy approximation↩︎

We give the localization argument behind 8. Fix \(N\) and a terminal time \(t>0\), let \(\Phi_s^N\) denote the particle flow, and write \[h_0:=\frac{\,\mathrm dP_0^N}{\,\mathrm d\pi^{\otimes N}}.\] For \(m\ge1\), set \[A_m:=\left\{\boldsymbol{x}: m^{-1}\le h_0(\boldsymbol{x})\le m, \quad \min_{0\le s\le t}\min_{i\ne j} \operatorname{dist}\bigl(\Phi_s^N(\boldsymbol{x})_i,\Phi_s^N(\boldsymbol{x})_j\bigr) \ge m^{-1} \right\}.\] The sets \(A_m\) increase. Moreover, \(P_0^N(\bigcup_mA_m)=1\): the density ratio is finite and positive \(P_0^N\)-almost surely, and 7 implies that the minimum pair distance along each trajectory is positive on the compact interval \([0,t]\). Put \(\alpha_m:=P_0^N(A_m)\) and, after discarding finitely many empty sets, define \[P_0^{N,m}:=\frac{P_0^N(\,\cdot\cap A_m)}{\alpha_m}, \qquad P_s^{N,m}:=(\Phi_s^N)_\#P_0^{N,m}.\]

On the trajectory tube generated by \(A_m\), the vector field and its first derivatives are bounded. The flow is a \(C^1\) one-to-one map there, so the area formula applies to the (not necessarily smooth) bounded density of \(P_0^{N,m}\). Equivalently, one may first smooth inside this tube and then remove the smoothing. The Jacobian satisfies \[\frac{\,\mathrm d}{\,\mathrm ds}\log\det D\Phi_s^N(\boldsymbol{x}) =\operatorname{div}_{\boldsymbol{x}}\boldsymbol{B}^N(\Phi_s^N(\boldsymbol{x})),\] whereas \[\frac{\,\mathrm d}{\,\mathrm ds}\log\pi^{\otimes N}(\Phi_s^N(\boldsymbol{x})) =\sum_{i=1}^Nb_\pi(\Phi_s^N(\boldsymbol{x})_i) \cdot B_i^N(\Phi_s^N(\boldsymbol{x})).\] Combining these two formulas with 39 gives \[\label{eq:localized-entropy} \operatorname{Ent}(P_t^{N,m}\mid\pi^{\otimes N}) +2N\int_0^t\mathbb{E}_{P_s^{N,m}}\mathcal{F}_N^\pi\,\,\mathrm ds =\operatorname{Ent}(P_0^{N,m}\mid\pi^{\otimes N}).\tag{73}\] This is the same calculation as in the smooth proof, now with every integration performed on a region uniformly separated from the singular set.

It remains to pass \(m\to\infty\). The local positivity of \(G_\pi\) and compactness away from the diagonal give a constant \(C_G\) such that \(G_\pi(x,y)\ge-C_G\) for \(x\ne y\). Hence the elementary shift \[\widehat{\mathcal{E}}_N^\pi:=\mathcal{F}_N^\pi+\frac{C_G}{2}\] is nonnegative for every configuration; the sharper Onsager estimate is not needed for this fixed-\(N\) approximation. Equation 73 becomes \[\label{eq:localized-renormalized-entropy} \operatorname{Ent}(P_t^{N,m}\mid\pi^{\otimes N}) +2N\int_0^t\mathbb{E}_{P_s^{N,m}}\widehat{\mathcal{E}}_N^\pi\,\,\mathrm ds =\operatorname{Ent}(P_0^{N,m}\mid\pi^{\otimes N}) +NC_Gt.\tag{74}\] Because \(\alpha_m\uparrow1\), conditioning on \(A_m\) gives \[\operatorname{Ent}(P_0^{N,m}\mid\pi^{\otimes N}) =\frac{1}{\alpha_m}\int_{A_m}h_0\log h_0\,\,\mathrm d\pi^{\otimes N} -\log\alpha_m \longrightarrow \operatorname{Ent}(P_0^N\mid\pi^{\otimes N}).\] The right-hand side of 74 is therefore bounded. Lower semicontinuity first shows that \(P_t^N=(\Phi_t^N)_\#P_0^N\) has finite entropy. Since the flow is injective, \(P_t^{N,m}\) is precisely \(P_t^N\) conditioned on \(\Phi_t^N(A_m)\). The same conditional-entropy formula at time \(t\) now gives \[\operatorname{Ent}(P_t^{N,m}\mid\pi^{\otimes N}) \longrightarrow\operatorname{Ent}(P_t^N\mid\pi^{\otimes N}).\] Finally, Tonelli’s theorem and monotone convergence apply to the energy: \[\int_0^t\mathbb{E}_{P_s^{N,m}}\widehat{\mathcal{E}}_N^\pi\,\,\mathrm ds =\frac{1}{\alpha_m}\int_{A_m}\int_0^t \widehat{\mathcal{E}}_N^\pi(\Phi_s^N(\boldsymbol{x}))\,\,\mathrm ds\,\,\mathrm dP_0^N(\boldsymbol{x}) \longrightarrow \int_0^t\mathbb{E}_{P_s^N}\widehat{\mathcal{E}}_N^\pi\,\,\mathrm ds.\] Passing to the limit in 74 and subtracting the deterministic correction \(NC_Gt\) proves the exact identity 23 .

References↩︎