ACTUAL SELF-SINKHORN STOPPING UNDER UNIFORM IID SAMPLING
Research proof, 18 October 2026. Independent review records are in reviews/.
This file extends the separate proof_candidate.txt; it does not change a paper.

1. Statement and precise scope

Fix d>=1 and an integer m>d/2+1. On the unit torus T^d use the integral-one
periodic Gaussian k_epsilon with Fourier multipliers
  lambda_k=exp(-2 pi^2 epsilon |k|^2).
Let X_1,...,X_n be iid uniform, mu_n=n^-1 sum delta_Xi, and epsilon=epsilon_n
decrease to zero with
  n epsilon_n^(6m) -> infinity.                                (H)
This deliberately conservative condition is sufficient, not asserted optimal.

Let a_n be the exact positive empirical self-scale, extended off sample by
  a_n(x) [K_mu a_n](x)=1.
Starting from b_0(x)=1 use the actual update
  b_(j+1)(x)=sqrt(b_j(x)/[K_mu b_j](x)).
Its restrictions to the observations are exactly the usual finite-vector
damped self-Sinkhorn iterates. Define
  v_j=K_mu b_j, s_j=2 grad log v_j,
  m_j=b_j K_mu b_j, S_j=int sqrt(m_j) dmu_n,
  r_j^2=int (sqrt(m_j)/S_j-1)^2 dmu_n.
The exact empirical score is s_*=2 grad log(K_mu a_n)=-2 grad log a_n.
All error assertions below compare the same data and same kernel.

Put alpha=d/4+1/2, A_n=sqrt(n) epsilon^alpha, and
  J_n=floor[(d/(16 log 2)) log(1/epsilon)],
  t_n=n^-1/2 epsilon^(-d/16), L_n=ceil(4 log n).
For n>=2 let tau_n be the first j in {0,...,L_n} with r_j<=t_n; if the
set is empty, set tau_n=L_n. This budget makes the rule defined on every
sample. The budget is reached without a threshold hit with probability o(1).

THEOREM. For every fixed query x in T^d:
(i) A_n |s_(J_n)(x)-s_*(x)| ->0 in probability, whereas
    sqrt(n) r_(J_n) -> infinity in probability.
(ii) With probability tending to one the first threshold hit exists before
     the budget, tau_n>J_n, and
       t_n/8 <= r_(tau_n) <= t_n.                              (1)
     Moreover A_n |s_(tau_n)(x)-s_*(x)| ->0 in probability.
     In particular sqrt(n) r_(tau_n) -> infinity in probability.
(iii) Both the exact and stopped scores have the same matched-model limit:
       A_n s_*(x) => N(0,V_d), A_n s_(tau_n)(x) => N(0,V_d),
     where, with Lambda(xi)=exp(-2 pi^2 |xi|^2),
       V_d=4 int_Rd (2 pi xi)(2 pi xi)^T
                         [Lambda(xi)/(1+Lambda(xi))]^2 dxi.    (2)
     The matrix is finite, positive definite and a scalar multiple of I_d.

Thus r=o_p(n^-1/2) is not necessary along these actual standard-start
trajectories, although it remains a valid worst-case sufficient condition
over arbitrary candidates. The explicit first residual threshold is looser
by the diverging factor epsilon^(-d/16) and preserves the pointwise limit.
No nonuniform density, Euclidean boundary model, optimal bandwidth window,
universal stopping tolerance, runtime gain or minimax claim is made.
Part (iii) is a supporting kernel CLT in this uniform periodic model; no
general novelty is claimed for a kernel CLT or for Slutsky's theorem.

2. A quantitative deterministic perturbation lemma

Write H=H_epsilon, delta=mu-dx and zeta=||delta||_(H^-m). The Fourier
Sobolev weight is (1+4 pi^2 |k|^2)^m. Constants below depend only on d,m
and the fixed small neighborhood, not epsilon in (0,1] or j>=0.
H^m is a Banach algebra and embeds continuously into C^1. By duality,
  ||h delta||_(H^-m) <= C ||h||_(H^m) ||delta||_(H^-m).
Also ||H||_(H^m -> H^m)<=1, and
  ||H||_(H^-m -> H^m)
    =sup_k (1+4 pi^2 |k|^2)^m lambda_k <= C epsilon^-m.
Consequently E_delta h:=H(h delta) obeys
  ||E_delta||_(H^m -> H^m) <= C epsilon^-m zeta.
For brevity let eta=C_0 epsilon^-m zeta, choosing C_0 sufficiently large.
Assume eta<=eta_0 for a fixed sufficiently small eta_0, and that mu is a
positive probability measure. The empirical measures satisfy this assumption
with probability tending to one; no signed-measure global existence is used.

The scaling equation for a=1+u expands as
  0=E_delta 1+(I+H)u+E_delta u+u E_delta 1+u Hu+u E_delta u.
The inverse R=(I+H)^-1 has H^m norm at most one. Applying -R to the
nonlinear terms gives a contraction on ||u||_(H^m)<=C eta: its derivative
there is bounded by C eta, and its value at zero has norm at most eta.
After fixing eta_0 small, the ball is invariant and the derivative norm is
less than 1/2. This proves existence, uniqueness in the ball, positivity by
Sobolev embedding, and the uniform estimates
  a=1-R Hdelta+O_(H^m)(eta^2), ||a-1||_(H^m)<=C eta.           (3)
On finite positive empirical measures, the restrictions solve the unique
finite self-scaling system. The equation gives its exact off-sample extension.
Thus this local function solution is the exact empirical solution on the
event being analyzed, not a different approximate solution.

The associated Markov operator is
  P_delta h=a (H+E_delta)(a h).
It satisfies P_delta 1=1 and
  ||P_delta-H||_(H^m -> H^m)<=C eta.                           (4)
Set B=(I-H)/2, whose H^m operator norm is at most 1/2. If
e_j=log(b_j/a), the exact relative-log update is
  e_(j+1)=Psi_delta(e_j)
           :=[e_j-log(P_delta exp(e_j))]/2.                  (5)
Analytic composition in a fixed small H^m ball, (4), and P_delta 1=1 give
  ||Psi_delta(e)-B e||_(H^m)
      <=C(eta ||e||_(H^m)+||e||_(H^m)^2).
The derivative is bounded by 1/2+C eta+C||e||_(H^m).
From (3), e_0=R Hdelta+O(eta^2). Choosing eta_0 smaller if necessary
makes the whole trajectory stay in the ball and gives q=3/4 with
  ||e_j||_(H^m)<=C eta q^j, all j>=0.                         (6)

Let ell_j=B^j R Hdelta. If w_j=e_j-ell_j, (5)--(6) imply
  ||w_(j+1)||<= (1/2)||w_j||+C eta^2 q^j, ||w_0||<=C eta^2.
Summing this geometric recursion gives
  ||e_j-ell_j||_(H^m)<=C eta^2 q^j, all j>=0.                (7)
All constants are uniform in j and epsilon. This is the quantitative step
not supplied by a fixed-epsilon implicit-function theorem alone.

3. Uniform score and residual remainders

The exact score difference at e is
  2 grad log[P_delta exp(e)].
Expanding it in H^m and using the C^1 embedding, (4), and (7) gives
  s_j-s_*=v_j^lin+R_j,
  v_j^lin=2 grad H^2 (I+H)^-1 B^j delta,
  sup_(j>=0) ||R_j||_infinity <= C eta^2.                    (8)
The stronger bound C eta^2 q^j also follows, but is not needed.

For the square-root mass field, exactly
  F_j:=sqrt(m_j)=[exp(e_j) P_delta exp(e_j)]^(1/2).
The same expansion yields
  F_j=1+z_j+h_j,
  z_j=(1/2)B^j Hdelta,
  ||z_j||_(H^m)<=C eta 2^-j, ||h_j||_(H^m)<=C eta^2 q^j.
The Lebesgue mean of z_j is zero because delta has zero mass. Therefore
  S_j=1+O(eta^2)
uniformly j: int z_j dmu=pair(delta,z_j)=O(zeta eta), and h_j is
O(eta^2) in C^0. After centering and dividing by S_j, the field differs
from z_j by O(eta^2) in H^m. Expanding its square under the positive
probability measure mu gives an O(eta^3) cross-term bound. Finally
  |int z_j^2 dmu-int z_j^2 dx|
      <=C zeta ||z_j||_(H^m)^2 <=C eta^3.
It follows, uniformly in every j>=0, that
  r_j^2=Q_j+O(eta^3), Q_j:=||z_j||_(L^2(dx))^2.              (9)
This bound uses a squared residual; it does not incorrectly take a Taylor
derivative of a norm at zero.

For uniform iid data,
  E zeta_n^2=n^-1 sum_(k!=0)(1+4 pi^2|k|^2)^-m=O(n^-1).
Thus eta_n=O_p(n^-1/2 epsilon^-m), and condition (H) gives
  eta_n=o_p(1),
  sup_j A_n ||R_j||_infinity=o_p(1),
  D_n:=sup_j n|r_j^2-Q_j|=o_p(1).                           (10)
Indeed n eta_n^3=O_p(n^-1/2 epsilon^-3m)=o_p(1), and
A_n eta_n^2=O_p(n^-1/2 epsilon^(alpha-2m))=o_p(1).
All deterministic estimates hold on events of probability tending to one;
this suffices for the probability conclusions even though the complement
is not claimed to lie in the local contraction basin.

4. Exact proxy variance, and a simultaneous score tail bound

Put beta_k=(1-lambda_k)/2 for k!=0 and beta_0=0. For j=0 use B^0=I;
delta has no zero mode, so the convention at k=0 is immaterial. Parseval gives
  Q_j=(1/4) sum_(k!=0) lambda_k^2 beta_k^(2j) |delta_hat(k)|^2.
Write a_k=lambda_k^2 beta_k^(2j), T_j=sum_(k!=0) a_k.
Uniform iid Fourier coefficients have covariance delta_kl/n, hence
  E Q_j=T_j/(4n),
  Var Q_j=(n-1)/(8n^3) sum_(k!=0) a_k^2.                    (11)
For completeness, the diagonal terms of the corresponding quadratic
U-statistic are constant. The off-diagonal kernel has Fourier coefficients
a_k, zero mean and squared L^2 norm sum a_k^2. Only equal and reversed
ordered index pairs contribute to its variance; their count is 2n(n-1).
The prefactor 1/(4n^2) gives exactly (11). The same identity holds for any
frequency set closed under k -> -k.

At a fixed query x,
  E|v_j^lin(x)|^2=(4/n) sum_(k!=0) |2 pi k|^2
                           lambda_k^4 beta_k^(2j)/(1+lambda_k)^2.
Since beta_k<=1/2 and
  epsilon^(d/2+1) sum |k|^2 exp(-8 pi^2 epsilon |k|^2)<=C,
  E[A_n^2 |v_j^lin(x)|^2]<=C 4^-j.
Consequently, for any deterministic J>=0,
  E[sup_(j>=J) A_n^2 |v_j^lin(x)|^2]
      <=sum_(j>=J) E[A_n^2 |v_j^lin(x)|^2]<=C 4^-J.          (12)
The first inequality is the elementary supremum <= sum inequality; it does
not require independence between iterations. Combining (10)--(12) proves
  sup_(j>=J_n) A_n |s_j(x)-s_*(x)| ->0 in probability.         (13)

5. A diverging root-n residual at a deterministic growing iteration

On the annulus where 1/4<=lambda_k<=1/2, beta_k>=1/4. Counting lattice
points after rescaling by sqrt(epsilon) gives at least c epsilon^-d/2
points for small epsilon. This uses only a fixed annulus with positive
volume and boundary of measure zero. Thus
  T_J>=c epsilon^-d/2 16^-J, max_k a_k<=4^-J.
For J=J_n, floor rounding changes constants only, giving
  T_J>=c epsilon^-d/4 ->infinity,
  sum a_k^2/T_J^2<=max a_k/T_J<=C epsilon^(3d/8)->0.
Equation (11) then gives
  4n Q_J/T_J ->1 in probability.                             (14)
By (10), n r_J^2 also diverges and is at least c epsilon^-d/4
with probability tending to one. This and (13) prove part (i).

6. The first residual hit and control of overshoot

Each Q_j is a nonnegative spectral sum with factors beta_k^(2j).
Consequently Q_j is decreasing in j. From (10), for every j<=J,
  n r_j^2 >=n Q_J-D_n.
By (14), the right side is at least c epsilon^-d/4 with high probability,
which exceeds n t_n^2=epsilon^-d/8. Hence the first hit is later than J
with probability tending to one.

The local contraction (6) implies r_j<=C eta_n q^j. Assumption (H) gives
eta_n=o_p(n^-1/3). Thus at L=ceil(4 log n), r_L=o_p(n^-1/2), since
q^L<=n^(4 log(3/4)) and 1/3-4 log(3/4)>1/2. As t_n>=n^-1/2,
the threshold is reached by L with probability tending to one. The identical
estimate at L-1 proves that the hit is strictly before L with probability
tending to one. Also J<L
eventually: (H) implies log(1/epsilon)<(log n)/(6m) eventually.
This proves the first-hit and budget assertions and, using (13), score
negligibility at the data-dependent tau_n. No optional-stopping theorem or
independence of tau_n is being assumed.

To show that the attained residual itself stays above the root-n scale,
split Q_J into low modes beta_k<1/4 and the other modes. Their low trace
is at most C epsilon^-d/2 16^-J. In the fixed shell
  1/8<=lambda_k<=1/4,
beta_k lies in [3/8,7/16]. Its trace is at least
  c epsilon^-d/2 (9/64)^J.
By (11), its relative variance is at most
  C epsilon^(d/2) (7/6)^(2J) ->0.
The latter limit follows directly from J<=d log(1/epsilon)/(16 log2).
The low expectation divided by this high expectation is at most
  C(4/9)^J ->0.
Chebyshev for the high shell and Markov for the low part show that low
modes contribute o_p(Q_J). Hence, with probability tending to one,
  Q_(J+1)>=Q_J/32.

For each fixed sample, positive spectral sums are log-convex:
  Q_(j+1)^2<=Q_j Q_(j+2)
by Cauchy--Schwarz. When Q_J>0, the ratios Q_(j+1)/Q_j are therefore
nondecreasing, so on the same event
  Q_(j+1)>=Q_j/32 for all j>=J.
At the first hit tau>J, r_(tau-1)>t_n. Equations (9)--(10) imply
  n r_tau^2 >=(1/32) epsilon^-d/8-(33/32)D_n.
The right side exceeds (1/64)epsilon^-d/8 with high probability, whereas
the threshold gives n r_tau^2<=epsilon^-d/8. Taking square roots yields
(1), and proves the remaining assertion of part (ii).

7. The matched-model pointwise central limit law

Equation (3), smooth logarithm and the C^1 embedding give
  s_*(x)=2 grad R Hdelta_n(x)+O_p(n^-1 epsilon^-2m).
The normalized remainder is o_p(1) under (H). The leading field is the
average of iid, mean-zero vectors g_epsilon(x-X_i), where
  g_epsilon(y)=2 sum_(k!=0) (2 pi i k)
                  [lambda_k/(1+lambda_k)] exp(2 pi i k.y).
It is real by conjugate symmetry. Its covariance after scaling is
  epsilon^(d/2+1) E[g_epsilon g_epsilon^T] ->V_d
by a Gaussian-dominated lattice Riemann sum. The limit is (2).
The deterministic bound
  sup_y |g_epsilon(y)|<=C epsilon^(-(d+1)/2)
implies that each normalized triangular-array summand has norm at most
  C/[sqrt(n) epsilon^(d/4)] ->0.
Here (H) implies n epsilon^(d/2)->infinity. Thus the Lindeberg condition
holds, and the multivariate triangular-array CLT proves the exact-score
limit. Equation (13) at tau_n and Slutsky give its stopped counterpart.
This proves (iii) within the explicitly stated uniform periodic model.

8. Mathematical increment and limits to retain in review

The known half-step filter B, geometric contraction, Hellinger residual,
delta method, kernel CLT and Slutsky are tools, not originality claims.
The candidate increment is their quantitative combination proving a
shrinking-bandwidth, uniform-in-iteration nonlinear expansion, an actual
first residual hitting-time guarantee, and a diverging attained root-n
residual despite the same pointwise score law. It answers the necessity
question in one nonempty model/window; it does not decide it for every
density, every bandwidth, or the Euclidean assumptions in Mordant's theorem.

For example d=1 permits m=2, condition n epsilon^12->infinity, and any
epsilon=n^-kappa with 0<kappa<1/12. The threshold exponent d/16 and all
other constants were selected for a short conservative proof. They are
not claimed sharp, practically calibrated, or evidence of a major advance.
The existing epsilon=.64/sqrt(n) calculations lie OUTSIDE (H); they only
diagnose finite trajectories and do not validate this asymptotic theorem.
