Proof. The full curvature-dimension condition implies the reduced condition used in some of the references, by (Bacher and Sturm 2010, Proposition 2.5(i)). The interior Poisson regularity and gradient estimates apply to a constant right hand side; see (Kell 2013; Jiang 2014). More precisely, if \(8B_R(z)\) is compactly contained in the domain, the quantitative estimate in (Jiang 2014, Theorem 3.1) has the form \[\|\nabla w\|_{L^\infty(B_R(z))} \le C\bigl(R^{-1}+\sqrt{c_\kappa(R^2)}\bigr) \bigl(\|w\|_{L^\infty(B_{8R}(z))} +R^2\|\Delta w\|_{L^\infty(B_{8R}(z))}\bigr).\] Here \(c_\kappa\) is the heat-semigroup variance constant in that estimate; the RCD curvature bound controls it, and the doubling and Poincare constants in \(C\) are uniform on the rescaled balls in use. With bounded rescaled curvature this gives the usual \(R^{-1}\|w\|_\infty+R|a|\) bound. The interior local-boundedness estimate first supplies the \(L^\infty\) norm from an \(L^2\) bound. Our later uses always leave the required eightfold interior room. The local test class has the characterization \[\operatorname{Test}_{\rm loc}(Y) =\{f\in D_{\rm loc}(\Delta)\cap L^\infty_{\rm loc}: |\nabla f|\in L^\infty_{\rm loc},\ \Delta f\in W^{1,2}_{\rm loc}\},\] with the analogous definition on \(\Omega\); see (Connell et al. 2021, Lemmas 2.6–2.7). Thus constant-Laplacian solutions satisfy the stated local test condition. Equivalently, one may first multiply by a cutoff to obtain Laplacian-domain and Hessian regularity, and then use a test cutoff on a smaller ball to obtain the local product calculus. The relevant Hessian and product rules are those of (Gigli 2018).

The noncollapsed hypothesis implies that the tangent module has dimension \(n\) almost everywhere, and that \(\operatorname{tr}\operatorname{Hess}w=\Delta w\) there. The dimensional Bochner inequality therefore implies \(M(w)\ge|\operatorname{Hess}w|^2\mathfrak m\); see (De Philippis and Gigli 2018; Han 2018). Expanding the square \[|\operatorname{Hess}w-bg|^2=|\operatorname{Hess}w|^2-2b\Delta w+nb^2\] proves (2). Local finiteness and the measure formulation follow from the same local test calculus. Polarization of a nonnegative quadratic form gives the last assertion about \((w,b)\).

The cutoff, Poincare, and Sobolev statements are local RCD estimates. For the heat-regularized cutoff construction, see (Mondino and Naber 2019, Lemma 3.1) and (Connell et al. 2021, Lemma 2.6). To obtain the zero-boundary coercivity in (3), extend \(h\) by zero to a fixed larger ball. Poincare on that ball and the positive lower bound for the measure of its annulus give the \(L^2\) estimate. The local Sobolev inequality followed by this estimate gives the other term. One may use an effective doubling dimension larger than both \(n\) and \(2\), so that \(d_*>1\) also when \(n=2\); see (Hajłasz and Koskela 2000, Theorem 5.1 and the following Remark 3). The annulus lower bound and all constants are uniform for the almost Euclidean balls used below. ◻

Proof. Take rescalings tending to zero. Compactness and the regular-point hypothesis imply pointed convergence to \(\mathbb R^n\), after passage to a subsequence, and noncollapsed volume convergence identifies the limit measure with \(\mathcal H^n\); see (De Philippis and Gigli 2018). The rescaled curvature lower bound can be taken with \(\kappa_\nu=\max\{-K,0\}s_\nu^2\), enlarged by an arbitrarily small positive number if necessary.

Fix for the moment a large finite radius. At every bounded center, volumes of balls at a fixed positive radius converge to their Euclidean values. This convergence is uniform in the center: otherwise a subsequence of centers converges in the Euclidean limit and contradicts volume convergence at that limit center. For a noncollapsed RCD space the volume density is at most one at every point, not merely almost everywhere. Bishop–Gromov therefore traps every smaller model-normalized volume ratio between its almost maximal value at the fixed radius and one. The model-volume correction tends uniformly to zero because the rescaled curvature does. This proves the first assertion in (6), uniformly over all smaller radii. The almost-volume-rigidity theorem in (De Philippis and Gigli 2018), applied on a fixed enlargement of the ball, proves the second assertion. Its interior-radius conclusion is why an enlargement is taken first. No openness of the set of exactly regular points is used.

Interior harmonic approximation to Euclidean coordinate functions is available by (Ambrosio and Honda 2018, Corollary 4.12). For each fixed large radius construct the harmonic functions on a strictly larger ball. Interior energy convergence, local boundedness, and the interior gradient estimate give convergence to coordinates and a dimensional gradient bound on the smaller ball. The latter bound is independent of the fixed large radius: the coordinate functions have size \(O(r)\) on a ball of radius \(r\), and the interior gradient estimate includes the factor \(r^{-1}\).

Finally use a diagonal subsequence, allowing the finite radii just chosen to increase sufficiently slowly. It can be chosen so that all the estimates hold up to the multiples of \(T_\nu\) displayed above, with \(\delta_\nu\to0\) and \(\kappa_\nu T_\nu^2\to0\). This gives (7) and the expanding domains of harmonicity simultaneously. ◻

Proof. We give the compactness induction, including the order of the smallness choices. Say that a harmonic column is \(\tau\)-normalizable at \((y,r)\) if there is an invertible matrix \(T\) for which \(T(u-u(y))/r\), on \(B_{2r}(y)\) in units \(r\), differs by at most \(\tau\) from a Euclidean coordinate map under a Gromov–Hausdorff approximation. Replacing the uniform error in this definition by an interior normalized \(L^2\) error gives an equivalent induction, after changing its tolerance by a fixed factor; we use harmonic interior estimates on a slightly larger ball whenever needed. Rotations of the Euclidean coordinates are allowed.

Choose a small fixed tolerance \(\tau\). In an exactly Euclidean ball, a harmonic column \(V\) that is \(\tau\)-close to the coordinate functions has \[\|DV(0)-\mathrm{Id}\|\le C_n\tau, \qquad \sup_{B_{2b}}|DV(0)^{-1}(V-V(0))-x| \le C_n\tau b^2.\] These follow by applying interior derivative estimates to \(V-x\). Choose a dyadic number \(b=2^{-m}\) so small that \(2C_n b<1/4\), and then choose \(\tau\) small enough that the inverse derivative is bounded. Thus affine renormalization improves the error in units \(b\) by a strict factor.

The same propagation holds on all sufficiently flat balls in our sequence. Indeed, if it failed, rescale a failing ball to unit size and apply its assumed normalization. The spaces converge to a Euclidean ball by (6). The normalized harmonic columns have uniform local \(L^2\), energy, and gradient bounds and therefore converge locally uniformly and strongly in interior \(W^{1,2}\) to a Euclidean harmonic column. These are interior convergence statements from (Ambrosio and Honda 2018, sec. 4); no convergence of a full boundary-value problem is required. The limit has the \(\tau\)-normalizability property. Its inverse derivative supplies a fixed matrix on the sequence which gives error strictly less than \(\tau\) at radius \(b\), a contradiction. This additional correction to the already normalized target, and its inverse, have norm at most \(1+C_n\tau\), after allowing the orthogonal target changes in the definition.

At each bounded center the initial normalizability is supplied by Lemma 4, uniformly in that center. Run the preceding step down each of the finitely many dyadic residue classes modulo \(m\). This proves normalizability at every dyadic scale. The mean Gram matrix of a normalized column on its unit ball is \(\mathrm{Id}+O(\tau)\), uniformly once the space is sufficiently flat: a contrary sequence again converges strongly in energy to a Euclidean harmonic column, whose derivative is \(O(\tau)\)-close to \(\mathrm{Id}\) on that interior ball. Consequently these Gram matrices are nonsingular and replacing an inductively chosen \(T\) by the canonical normalization \(R_i(y)^{-1/2}\) changes the normalized map only by an orthogonal map and an \(O(\tau)\) error.

To see explicitly why fixed enlargements are permitted, fix their largest factor first. Start the preceding induction on a larger fixed radius, allowing finitely many negative dyadic indices. For any scale \(r_i\), an ancestor radius contains the desired enlargement, and at most a fixed number of the contraction steps separate that ancestor from \(r_i\). Their normalizing matrices and inverse matrices are uniformly comparable. Apply the interior estimates on the ancestor ball and identify the mean Gram normalization on the descendant ball by strong energy convergence. All factors and the number of steps are fixed before \(\tau\) and the geometric flatness are chosen. This proves normalizability and bounded normalized gradients on every fixed enlargement. As \(\tau\) is arbitrary, the errors tend uniformly to zero along the original sequence, including when the chosen center and scale depend on the sequence index.

The same argument on an enlargement containing both balls in (10) proves that their Gram averages in the first normalization both tend to \(\mathrm{Id}\). This gives (10). Iterating its same-center adjacent-scale case, with error chosen so that each factor is at most \(2^{2a}\), gives (11). The upper bound in (12) also follows directly from (7).

For completeness, the strong convergence of the Gram functions can be seen without assuming convergence of pointwise gradient representatives. A Bochner cutoff estimate first bounds the local normalized Hessian energy. The Hessian product rule and bounded normalized gradients then bound the Gram functions in local \(W^{1,2}\). Sobolev compactness and the convergence of integrals of gradient scalar products identify their limit with \(\mathrm{Id}\); see (Ambrosio and Honda 2018, Corollary 4.3 and Theorem 4.6). Their uniform boundedness upgrades the resulting strong \(L^2\) convergence to every finite \(L^p\). ◻

Proof. In units \(r_i\), apply Lemma 3 to a scalar combination of \(U=U_{i,y}\). The rescaled Bochner matrix is \[ M_{r_i}(U)=r_i^{2-n}R_i(y)^{-1/2}M(u)R_i(y)^{-1/2}, \tag{14}\] where the curvature term in its definition is \(\kappa r_i^2\). Choose a nonnegative cutoff equal to one on \(B_C\) and supported on a larger fixed ball. Its gradient and Laplacian are uniformly bounded. For a fixed vector \(b\), \[\int\phi\,\mathrm dM_{r_i}(b\cdot U) =\frac12\int\bigl(\Gamma(b\cdot U)-|b|^2\bigr) \Delta\phi\,\mathrm d\mathfrak m_{r_i} +\kappa r_i^2\int\phi\Gamma(b\cdot U)\,\mathrm d\mathfrak m_{r_i}.\] Here \(\int\Delta\phi\,\mathrm d\mathfrak m_{r_i}=0\). Both terms tend uniformly to zero by Lemma 5. Positivity of the Bochner measure bounds its mass on \(B_C\). Polarization in the fixed finite dimension, or the same argument for unit vectors, gives the matrix bound. Conjugating back proves (13). ◻

Proof. The approximation in Lemma 5 gives (16). Test a pair of distinct points at a dyadic scale comparable with its distance. Once the error is less than one half of that distance, the coordinate values are distinct. Distances larger than the initial unit scale in a fixed ball are separated directly by the initial uniform Euclidean approximation. Thus \(u\) is injective on \(B_A(p)\).

We prove (15) with its interior margin. Fix a target \(P\) in the indicated ellipsoid and put \(x_0=x\), \(\rho_0=r_i\). Write \(R(\rho,y)\) for the Gram average at a dyadic radius. Suppose inductively that \[|P-u(x_j)|_{R(\rho_j,x_j)^{-1}}<c\rho_j, \qquad \rho_j=4^{-j}\rho_0.\] Approximate surjectivity of the normalized Euclidean coordinate map provides \(x_{j+1}\) with \[d(x_j,x_{j+1})\le(c+\eta)\rho_j, \qquad |P-u(x_{j+1})|_{R(\rho_j,x_j)^{-1}}\le\eta\rho_j.\] Here \(\eta\to0\) is uniform in the scale. For example, choose a point under the Gromov–Hausdorff correspondence whose Euclidean coordinate is the normalized target. Its distance from the center and the error in its coordinate value have the displayed bounds. By (10), the second bound is at most \(2\eta\rho_j\) when measured in \(R(\rho_j/4,x_{j+1})^{-1}\). Choose \(c<1/8\), and then choose \(\eta<c/8\). It is at most \(c\rho_{j+1}\), so the induction continues.

The total displacement is at most \((c+\eta)\sum_j\rho_j<\rho_0/2\). All centers therefore stay in the region where the estimates are valid. Completeness gives a limit \(x_\infty\in B_{\rho_0/2}(x)\), with a strict margin because of the preceding choices. The upper bound \(R(\rho,y) \le C\mathrm{Id}\) in the original coordinate units shows that the Euclidean residual \(|P-u(x_j)|\) tends to zero. Continuity of \(u\) now gives \(u(x_\infty)=P\).

Apply the ellipsoid statement at each point with a sufficiently small radius contained in \(B_A(p)\). Its image is open. Finally, a continuous injection is a homeomorphism locally on each compact interior ball, since these spaces are proper. Alternatively the same separation-scale estimate gives continuity of the inverse. This proves the assertion on the whole open domain. ◻

Proof. When \(D\le 2(C_0+2)r\), symmetry and total heat mass bound the integral of \(p_t(a,z)\) over the ball by one. Dividing by its volume gives \(Cr^{-n}\), which is the first bound with adjusted constants. If \(D>2(C_0+2)r\), then \(d(a,z)\ge D/2\) for all \(a\in B_r(y)\). The Gaussian upper bound gives \(Ct^{-n/2}e^{-cD^2/t}\). For \(t\le r^2\) and \(D/r\) bounded below by this fixed constant, its small-time power is absorbed in the exponential, giving (19).

For the second bound, \(d(a,x)\le(C_0+1)r\). Integrate the first spatial gradient estimate along a minimizing curve from \(x\) to \(a\); that curve stays in a fixed testing enlargement. The inequality \(d(b,z)^2\ge D^2/2-d(b,x)^2\), with \(t\ge r^2\), absorbs the variation of the Gaussian along the curve. This gives (20), initially for almost every fixed source pole \(z\). For each such pole, local Sobolev-to-Lipschitz gives the difference inequality for the continuous representative in the first variable. By symmetry, continuity of the common heat kernel, and full support, a sequence of these full-measure poles converges to any prescribed \(z\). Passing to the limit, with relaxed Gaussian constants if necessary, proves the bounds at every pole. No gradient representative has been paired with a singular measure.

For (21), first test the distributional Laplacian against positive-time heat regularizations of compactly supported test functions. Symmetry and integration by parts give the asserted identity in distributions. Equivalently, differentiate \(P_\tau f\) for positive \(\tau\), pair the Laplacian with the continuous Sobolev function \(p_\tau(x,\cdot)\), and integrate. A cutoff equal to one on the compact supports justifies this pairing; the kernel belongs locally to \(W^{1,2}\) at positive time. The measure and its Sobolev-dual realization agree on this continuous representative. The kernel bounds give continuity of both sides in \(x\), so the identity holds pointwise.

In our application the short-time averaged integral of the total variation of the source is finite: (19), integrated up to \(r^2\), bounds it by a finite Gaussian-weighted multiple of \(r^{2-n}|\mu|\). Also \(P_s f\to f\) in \(L^1\) on the averaging ball. Thus the averaged identity passes to \(s=0\) by absolute convergence. These statements also follow by first using positive lower time endpoints throughout and then applying the bounds just proved. ◻

Proof. The original gradient bound and a Bochner cutoff on a ball of radius \(s\ge1\) give, for centers in a fixed bounded set and up to the support region, \[ M(B_s(x))[v]\le C(s^{n-2}+\kappa s^n)|v|^2 \le Cs^{n-2}|v|^2. \tag{28}\] For radii larger than the support radius the same bound holds with \(\chi M\) in place of \(M\), by monotonicity and \(n\ge2\). The cutoff for this estimate lies inside the expanding domain where (7) holds. The same estimate for the Hessian, followed by Cauchy–Schwarz in the product rule, gives \[\int|E_{\rm cut}|\,\mathrm d\mathfrak m\le CT^{n-2}.\] Since its support is a distance at least \(T/2\) from the fixed testing centers for large \(\nu\), (17) yields \[ \sup_{0<t\le1}\|P_t E_{\rm cut}(x)\| \le C T^{n-2}\sup_{0<t\le1}t^{-n/2}e^{-cT^2/t} =o(1). \tag{29}\] Choose \(\lambda\) at least this supremum plus \(2\kappa\|G\|_\infty\), and add a positive number tending to zero. Equations (21) and (22) give \[\partial_t(P_t(\chi G)+\lambda t\mathrm{Id}) =2P_t(\chi M)-2\kappa P_t(\chi G) +P_tE_{\rm cut}+\lambda\mathrm{Id}.\] The last three terms have nonnegative sum and bounded norm. This proves (23), in fact with lower bound \(2P_t(\chi M)\).

We verify relative, rather than merely absolute, comparison with the ball averages. Fix \(i,x\) and conjugate by \(R_i(x)^{-1/2}\). On every fixed ball \(B_{Lr_i}(x)\), the conjugated Gram field converges in normalized \(L^1\) to \(\mathrm{Id}\), by Lemma 5. The heat kernel at time \(r_i^2\) is at most \(Cr_i^{-n}\), so the integral of its local error tends to zero. The tails are uniform in \(i\): on the shell with radius \(2^m r_i\le1\), slow growth bounds the conjugated Gram average by \(C2^{2am}\), while the Gaussian bounds its contribution by \[C2^{(n+2a)m}e^{-c4^m}.\] When \(2^m r_i>1\), use the original gradient bound and (12); its conjugation cost \(Cr_i^{-2a}\) is at most \(C2^{2am}\). Volume comparison and the same Gaussian again give the displayed summable bound. This treats the entire support of \(\chi\). The missing mass from \(1-\chi\) is bounded in the same way, using total heat mass one. First choose \(L\) large enough to make the tail small and then choose the rescaling late enough to make the local error small. Beyond the expanding almost Euclidean region, the global volume bound adds at most \(e^{C2^m r_i}\) to a shell estimate; since \(r_i\le1\), this factor is absorbed by decreasing the Gaussian constant in \(e^{-c4^m}\). Thus the assertion about the missing heat mass is uniform there as well. Finally, \[\lambda r_i^2R_i(x)^{-1}\le C_a\lambda r_i^{2-2a}\mathrm{Id}=o(1)\mathrm{Id}\] for \(a<1\). This proves (25); its initial scale assertion follows from Euclidean convergence.

The derivative positivity makes \(Q_i\ge Q_{i+1}\), hence \(H_i\le H_{i+1}\). Relative comparison with \(R_i\) and (10) give the adjacent-scale upper bound and all fixed nearby comparisons. Iteration gives (27), choosing the adjacent-scale tolerance after the prescribed exponent \(a\). ◻

Proof. The semigroup identity. Write \(r=r_i\), suppress \(x\) in \(N_j,Q_j\), and let \(\mathcal A f=\mathop{\mathchoice {\displaystyle\int\mkern-18mu-} {\textstyle\int\mkern-15mu-} {\scriptstyle\int\mkern-12mu-} {\scriptscriptstyle\int\mkern-10mu-}}\nolimits_{B_r(y)}f\,\mathrm d\mathfrak m\). We first prove the more explicit estimate \[ \begin{split} |v^T(R_i(y)-Q_i(x))w| &\le Cr|v||w|\\ &\quad+C\sum_{j=0}^i2^{-(i-j)} \sqrt{N_j[v]N_j[w]}. \end{split} \tag{35}\] This also specifies all the spatial and time tails entering (34).

Put \(f=\chi G\), \(S=\Delta f\), and use the kernel functions from Lemma 8. Since \(\chi=1\) on the averaging ball, the exact semigroup identity is \[ \begin{split} R_i(y)-Q_i(x) ={}&\mathcal A P_1f-P_1f(x)-\lambda r^2\mathrm{Id}\\ &-\int_0^{r^2}\mathcal A P_tS\,\mathrm dt -\int_{r^2}^1\bigl(\mathcal A P_tS-P_tS(x)\bigr)\,\mathrm dt. \end{split} \tag{36}\] The averaged endpoint at zero is justified in Lemma 8; it does not require pointwise convergence of the Gram field.

The local source. Consider the contribution of \(2\chi M\) on \(B_1(x)\). The integrated short-time kernel obeys \[ \int_0^{r^2}K_t(z)\,\mathrm dt \le Cr^{2-n}e^{-c d(x,z)^2/r^2}. \tag{37}\] For \(0\le q<i\), integration of (20) on the time block \(I_q=[r_{q+1}^2,r_q^2]\) gives \[ \int_{I_q}|K_t(z)-p_t(x,z)|\,\mathrm dt \le Cr r_q^{1-n}e^{-c d(x,z)^2/r_q^2}. \tag{38}\] To integrate a Gaussian at scale \(r_q\) against \(M(v,w)\), partition \(B_1(x)\) into \(B_{r_q}(x)\) and the annuli \[B_{r_j}(x)\setminus B_{r_{j+1}}(x),\qquad 0\le j<q.\] Define positive numerical coefficients \[ k_0=1,\qquad k_h=2^{(n-2)h}e^{-c4^{h-1}}\quad(h\ge1). \qquad \sum_{h=0}^\infty2^h k_h<\infty. \tag{39}\] On the annulus indexed by \(j<q\), the distance from \(x\) is at least \(2^{q-j-1}r_q\). Equation (5) and the identity \(M(B_{r_j}(x))=r_j^{n-2}N_j\) show that the contribution of the time block \(I_q\) is at most \[C2^{-(i-q)}\sum_{j=0}^q k_{q-j} \sqrt{N_j[v]N_j[w]}.\] Reversing the finite sums, the sum of all these contributions is at most \[\begin{align*} &C\sum_{j=0}^{i-1}2^{-(i-j)} \left(\sum_{h=0}^{i-1-j}2^h k_h\right) \sqrt{N_j[v]N_j[w]}\\ &\hspace{25mm}\le C\sum_{j=0}^{i-1}2^{-(i-j)}\sqrt{N_j[v]N_j[w]}. \end{align*}\] The short-time term uses the same partition with \(q=i\) and (37). It is bounded by \[C\sum_{j=0}^i k_{i-j}\sqrt{N_j[v]N_j[w]} \le C\sum_{j=0}^i2^{-(i-j)}\sqrt{N_j[v]N_j[w]}.\] This proves the required contribution on \(B_1(x)\), with the original unit-ball \(N_j\) and no unmentioned enlarged mass.

The remaining source terms. We treat all remaining parts of (36). If \(D=d(x,z)\ge1\), the integrated kernel bounds imply \[ \int_0^{r^2}K_t(z)\,\mathrm dt+ \int_{r^2}^1|K_t(z)-p_t(x,z)|\,\mathrm dt \le Cr e^{-c'D^2}. \tag{40}\] For the second term, split the Gaussian into two factors, use one to retain \(e^{-c'D^2}\), and use the other to integrate \(t^{-(n+1)/2}e^{-c''/t}\) on \((0,1)\). For the first term, (37) reduces the assertion to absorbing \(r^{1-n}\) in \(e^{-c''/r^2}\). Both bounds are uniform for \(0<r\le1\).

Equation (28), summed over the annuli \(B_{2^{h+1}}(x)\setminus B_{2^h}(x)\), gives \[\int e^{-c'd(x,z)^2}\,\mathrm d(\chi M[v])(z)\le C|v|^2.\] Therefore (40) and polarized Cauchy–Schwarz bound the exterior positive-source contribution by \(Cr|v||w|\). The total variation of the cutoff source is \(O(T^{n-2})\), is supported a distance comparable with \(T\) from \(x\), and has the same Gaussian moment bound, in fact a bound tending to zero. Its contribution is also at most \(Cr|v||w|\).

The curvature source in (22) has bounded density, of norm at most \(C\kappa\). For its short-time part use heat mass one. For its long-time part, integration of (20) in the source variable gives \(Cr/\sqrt t\). This last Gaussian integral follows from local volume comparison for nearby shells and the global exponential volume bound for remote shells. Consequently the curvature contribution is at most \[C\kappa\left(r^2+\int_{r^2}^1 r t^{-1/2}\,\mathrm dt\right)|v||w| \le Cr|v||w|.\] The time-one heat regularization of the bounded matrix field \(f\) has uniformly bounded spatial gradient in the testing region, by the integrated first kernel gradient estimate. Hence \(|v^T(\mathcal A P_1f-P_1f(x))w|\le Cr|v||w|\). Also \(\lambda r^2|v||w|\le Cr|v||w|\). These facts prove (35).

Normalization and decrements. Choose any sufficiently small fixed growth exponent \(a>0\) in (27). By (31), for \(j\le i\), \[\sqrt{N_j[v]}\le C2^{a(i-j)}|v|_{Q_i}.\] Weighted Cauchy–Schwarz now gives, for instance when \(a<1/8\), \[\begin{align*} \sum_{j=0}^i2^{-(i-j)}\sqrt{N_j[v]N_j[w]} &\le C|v|_{Q_i}\sum_{j=0}^i2^{-(1-a)(i-j)}\sqrt{N_j[w]}\\ &\le C|v|_{Q_i} \left(\sum_{j=0}^i2^{-(i-j)/2}N_j[w]\right)^{1/2}. \end{align*}\] Indeed the unused scalar series has terms \(2^{-(2(1-a)-1/2)(i-j)}\), whose sum is bounded. Since \(Q_i\ge c_a r_i^{2a}\mathrm{Id}\), the raw term in (35) is at most \[Cr_i^{1-a}|v|_{Q_i}|w| \le C|v|_{Q_i}\bigl(2^{-i/2}|w|^2\bigr)^{1/2}.\] This proves (34) with \(\gamma=1/2\).

It remains to prove the decrement estimates. Integrate (23) on \([r_i^2/4,r_i^2]\). On a fixed ball \(B_{Ar_i}(x)\), where \(\chi=1\), the heat lower bound (18) gives the left inequality in (33). The upper bound is \[Q_i-Q_{i+1} \le Cr_i^{2-n}\int e^{-c d(x,z)^2/r_i^2}\,\mathrm d(\chi M)(z) +Cr_i^2\mathrm{Id}.\] On \(B_1(x)\), the shell computation with \(q=i\) bounds the positive term by \(C\sum_{j=0}^i k_{i-j}N_j\), which is at most \(C\sum_{j=0}^i2^{-(i-j)/2}N_j\). Outside \(B_1(x)\), (28) and Gaussian decay give \(Cr_i^2\mathrm{Id}\), after absorbing every fixed power of \(r_i^{-1}\) in \(e^{-c/r_i^2}\). Since \(r_i^2\le2^{-i/2}\), this proves the upper bound with the same \(\gamma\). ◻

Proof. Let \(E_i=Q_i^{-1/2}\Delta_iQ_i^{-1/2}\). By (26), \(0\le E_i\le\varepsilon_\nu\mathrm{Id}\) after increasing the sequence \(\varepsilon_\nu\to0\). The exact inverse identity is \[B_i=Q_i^{-1/2}\bigl((\mathrm{Id}-E_i)^{-1}-\mathrm{Id}\bigr)Q_i^{-1/2}.\] The spectral inequality \(E_i\le(\mathrm{Id}-E_i)^{-1}-\mathrm{Id}\le(1+C\varepsilon_\nu)E_i\) proves (41). This use of the functional calculus is for one conjugated positive matrix at a time; matrices from different scales need not commute.

We spell out the convolution estimate. Since \(0\le B_l\le C\varepsilon_\nu H_l\), conjugation by \(Q_l^{1/2}\) gives \(B_lQ_lB_l\le C B_l\). Consequently, by (27), \[ |B_l\xi|_{Q_j}^2\le C2^{2a(l-j)}B_l[\xi],\qquad j\le l. \tag{44}\] For \(j<i\), set \(m=i-j\). The identity \(a_i-a_j=\sum_{l=j}^{i-1}B_l\xi\) and Cauchy–Schwarz give \[|a_i-a_j|_{Q_j}^2 \le C m2^{2am}\sum_{l=j}^{i-1}B_l[\xi].\] Using \(N_j\le\varepsilon_\nu Q_j\), it follows that \[\begin{align*} &\sum_{i<k}\sum_{j\le i}2^{-\gamma(i-j)}N_j[a_i] \tag{45}\\[-2mm] &\qquad\le C\sum_{j<k}N_j[a_j] +C\varepsilon_\nu \sum_{i<k}\sum_{j<i}(i-j)2^{-(\gamma-2a)(i-j)} \sum_{l=j}^{i-1}B_l[\xi]\\ &\qquad\le C\sum_{j<k}N_j[a_j] +C\varepsilon_\nu H_k[\xi]. \end{align*}\] For the last inequality, fix \(l\). For a specified \(m=i-j\), there are at most \(m\) intervals \([j,i)\) of length \(m\) containing \(l\). Hence its coefficient is bounded by \(C\sum_{m\ge1}m^2 2^{-(\gamma-2a)m}<\infty\). Also \(\sum_{l<k}B_l[\xi]=(H_k-H_0)[\xi]\le H_k[\xi]\). The identity contribution is bounded separately using (27): \[ \sum_{i<k}2^{-\gamma i}|a_i|^2 \le C|\xi|^2\sum_{i\ge0}2^{-(\gamma-4a)i} \le C|\xi|^2. \tag{46}\]

The lower heat decrement with \(A=1\), together with (41), gives \[\sum_{i<k}N_i[a_i]\le C\sum_{i<k}B_i[\xi] \le C H_k[\xi].\] Conversely, the upper heat decrement yields \[H_k[\xi]\le C|\xi|^2+C\sum_{i<k}\mathcal F_i[a_i] \le C|\xi|^2+C\sum_{i<k}N_i[a_i] +C\varepsilon_\nu H_k[\xi].\] Choose the rescaling late enough to absorb the last term. This proves (43), since \(H_0\) is comparable with \(\mathrm{Id}\). Equations (45) and (46) now imply (42). Finally the lower heat decrement for a general fixed \(A\) bounds each enlarged-mass quadratic form on \(a_i\) by a fixed multiple of \(B_i[\xi]\). Summing proves the last assertion. ◻

Proof. Taking the supremum over \(v\) with \(|v|_{Q_i}=1\) in (34) gives \[ A_{i,y}H_iA_{i,y}\le C\mathcal F_i. \tag{50}\] This follows directly from \(\sup_{|v|_{Q_i}=1}|v^Tz|^2=z^TH_i z\). Write \(a_i=H_i\xi\). The inverse identity is \[e_{i,y}=-R_i(y)^{-1}A_{i,y}a_i.\] Relative comparison \(R_i(y)\asymp Q_i(x)\) implies \(R_i(y)^{-1}Q_i(x)R_i(y)^{-1}\le C H_i(x)\). Therefore (50) gives \[|e_{i,y}|_{Q_i}^2 \le C a_i^T A_{i,y}H_iA_{i,y}a_i \le C\mathcal F_i[a_i].\] This proves (47). Sum it, use the bound \(J\), and apply (42) to obtain (48). The last assertion follows by \(C_{i,y}[e_{i,y}]\le\varepsilon_\nu|e_{i,y}|_{Q_i}^2\). No commutation of \(R_i(y)\) with \(Q_i(x)\) has been used. ◻

Proof. Fix \(x\) and \(\xi\), and suppress \(x\) from the notation. By (13) and nearby Gram comparison, \(C_{i,y}\le\varepsilon_\nu Q_i\), uniformly in all the choices. Conjugation by \(R_i(y)^{-1}\), followed by relative comparison, gives \(D_{i,y}\le C\varepsilon_\nu H_i\). Bounded overlap proves the first assertion of (53).

For the upper sum, use \(R_i(y)^{-1}\xi=a_i+e_{i,y}\). Positivity gives \[D_{i,y}[\xi]=C_{i,y}[a_i+e_{i,y}] \le2C_{i,y}[a_i]+2C_{i,y}[e_{i,y}].\] The ball defining \(C_{i,y}\) is contained in \(B_{(L+C_0)r_i}(x)\). The enlarged upper sum assertion of Lemma 12 therefore bounds the sum of its first terms by \(C_J H_k[\xi]\). Equation (49) bounds the sum of its second terms by \(C_J\varepsilon_\nu H_k[\xi]\). Since \(\mathrm{Id}\le C H_0\le C H_k\), this proves the middle assertion in (53).

For the reverse direction, at each scale select one of the weight-one centers whose ball covers \(B_{r_i}(x)\). Then \(N_i\le C_{i,y}\). Write \(a_i=R_i(y)^{-1}\xi-e_{i,y}\) to obtain \[N_i[a_i]\le2C_{i,y}[R_i(y)^{-1}\xi] +2C_{i,y}[e_{i,y}] \le2S_i[\xi]+2C_{i,y}[e_{i,y}].\] The central comparison (43) and the normalization-error bound give \[H_k[\xi]\le C\left(|\xi|^2+\sum_{i<k}S_i[\xi]\right) +C\varepsilon_\nu H_k[\xi].\] Choose the rescaling sufficiently late after all fixed constants and absorb the last term. The result is uniform in \(k,x,\xi\) and the center choices. Taking all vectors \(\xi\) proves the last matrix inequality and completes the lemma. ◻

Proof. We first explain the interpretation of tests against \(\mu\). Since \(\mu=-\Delta h\), its action extends continuously to \(W^{1,2}_0(\Omega)\). This positive functional does not charge a set of zero Sobolev capacity: for a compact such set, test it with nonnegative functions which are at least one near that set and whose \(W^{1,2}_0\) norms tend to zero, and use continuity of the functional. Regularity of the measure gives the same conclusion for a Borel set of zero capacity. Consequently bounded Sobolev functions are integrated using their quasi-continuous representatives. Approximation in \(W^{1,2}_0\), followed by bounded truncation and a subsequence converging outside a set of zero capacity, shows that the functional and the measure integral agree on these tests. In particular all bounded Lipschitz compositions of \(h\) which vanish at zero are admissible. This justifies the truncations used next without making a choice of arbitrary \(\mathfrak m\)-almost-everywhere representatives on the support of \(\mu\).

If \(P=0\), testing with \(h\) proves \(h=0\). Suppose \(P>0\) and put \(H=h/P\). For \(t>0\), the test \(\min\{H,t\}\) gives \[\int_{\{H<t\}}|\nabla H|^2\,\mathrm d\mathfrak m\le t.\] Applying (56) to the same truncation yields \[t^2\mathfrak m(H\ge t)^{1/d}\le Ct, \qquad \mathfrak m(H\ge t)\le Ct^{-d}.\] Choose \(0<\alpha<d-1\). The layer-cake formula, using the volume bound for \(t\le1\) and the last estimate for \(t\ge1\), implies \[ \int_\Omega(1+H)^{1+\alpha}\,\mathrm d\mathfrak m\le C. \tag{59}\] The bounded test \(1-(1+H)^{-\alpha}\) gives \[\alpha\int_\Omega \frac{|\nabla H|^2}{(1+H)^{1+\alpha}}\,\mathrm d\mathfrak m\le1.\] Combining this inequality with (59) by weighted Cauchy–Schwarz proves \(\int|\nabla H|\le C\), which is the first assertion. Finally \(h\) itself is an admissible bounded test, so \(\int|\nabla h|^2=\int h\,\mathrm d\mu\le\|h\|_\infty P\). ◻

Proof. Choose \(C_0\) sufficiently large that Laplacian comparison gives \[\Delta q^a\le(n+C_0\kappa r_i^2)\mathfrak m \quad\hbox{on }B_{16Sr_i}(a).\] This is possible because \(d_a\le16Sr_i\) there. Take the sequence late enough that \(C_0\kappa\le1\). Define the positive measure \[\mu_i^a=a_i\mathfrak m-\Delta q^a \quad\hbox{on }B_{16Sr_i}(a), \qquad P_i(a)=r_i^{-n}\mu_i^a(B_{8Sr_i}(a))+\kappa r_i^2.\] Here the measure comparison for the squared distance is on the whole ball, including its center; see (Gigli 2015, Corollary 5.15, equation (5.29)). The source restricted to \(B_{8Sr_i}(a)\) is a continuous functional on \(W^{1,2}_0\) since it is the difference of a bounded density and the distributional Laplacian of the Sobolev function \(q^a\). The Dirichlet variational problem therefore gives \(h_i^a\in W^{1,2}_0\) with \(-\Delta h_i^a=\mu_i^a\). Testing its negative part gives \(h_i^a\ge0\), and \(f_i^a=q^a-h_i^a\) satisfies (64).

Here is the mass estimate, including its normalization. Set \(V_a(t)=\mathfrak m(B_t(a))\) and \(A_a(t)=t^{-n}V_a(t)\). Radial testing of the measure-valued Laplacian of \(q^a\), with \(|\nabla d_a|=1\) almost everywhere off \(a\), gives \[\Delta q^a(B_t(a))=tV_a'(t) \quad\hbox{for almost every }t.\] Equivalently, this follows by approximating radial indicators in \(-\int t\psi'(t)V_a'(t)\,\mathrm dt\); the boundary term at zero vanishes by the local volume upper bound. Thus \[\mu_i^a(B_t(a))=(a_i-n)V_a(t)-t^{n+1}A_a'(t)\] at those radii. Positivity and integration for \(8Sr_i<t<16Sr_i\) imply \[\begin{align*} r_i^{-n}\mu_i^a(B_{8Sr_i}(a)) &\le C\int_{8Sr_i}^{16Sr_i} \frac{\mu_i^a(B_t(a))}{t^{n+1}}\,\mathrm dt\\ &\le C\{A_a(8Sr_i)-A_a(16Sr_i)+\kappa r_i^2\}. \tag{68}\end{align*}\] The last step uses \(V_a(t)\le Ct^n\) on these fixed rescaled balls. For a finite dyadic sum, the volume terms on the right telescope exactly to \(A_a(8Sr_N)-A_a(16S)\). The noncollapsed Bishop inequality gives \(\lim_{t\downarrow0}A_a(t)\le\omega_n\) at every center, where \(\omega_n\) is the Euclidean unit-ball volume. The fixed-radius volume convergence in the preceding section gives \(A_a(16S)\ge\omega_n-o(1)\) uniformly in \(a\in B_4(p)\). Since \(\sum_i\kappa r_i^2\le4\kappa/3\), this proves (63). In particular this argument does not assume that the neighboring centers are regular.

For completeness, the bound for the Dirichlet solution can be obtained with uniform constants directly from (56). Work in units \(r_i=1\). Let \(\tau\) solve \(-\Delta\tau=1\) with zero boundary values on \(B_{8S}(a)\). If \(Y(k)=\mathfrak m(\tau>k)\), testing with \((\tau-k)_+\) and applying Sobolev and Hölder gives, for \(l>k\), \[Y(l)\le C(l-k)^{-2d}Y(k)^{2d-1}.\] Put \(k_j=K(1-2^{-j})\). Since \(2d-1>1\), choose \(K\) using only \(C,d,Y(0)\) so that this recursion implies \(Y(k_j)\le Y(0)2^{-jd/(d-1)}\) by induction. Then \(\tau\le K\). The harmonic function \(f_i^a+a_i\tau\) has the nonnegative boundary values \(q^a\), so the weak maximum principle gives \(f_i^a\ge-a_i\tau\). Together with \(h_i^a\ge0\) this gives \(0\le h_i^a\le q^a+a_i\tau\le C\). Rescaling proves the first bound in (65).

The local Poisson gradient estimate gives \(|\nabla f_i^a|\le C\) throughout \(B_{4S}(a)\) in these units, since the values and the constant right side are bounded. More explicitly, for each \(q\in B_{4S}(a)\) apply the estimate on \(B_{S/4}(q)\); its eightfold enlargement \(B_{2S}(q)\) is contained in \(B_{6S}(a)\), strictly inside the Dirichlet domain \(B_{8S}(a)\). The estimate and its uniform local dependence follow from (Kell 2013; Jiang 2014); the local test-function interpretation is supplied by Lemma 3 (see also (Connell et al. 2021)). Apply Lemma 15 to \(h_i^a\) on this same rescaled Dirichlet domain. Its source mass is at most \(P_i(a)\). This proves both estimates in (66), including their displayed powers of \(r_i\).

It remains to bound the full shifted measure. Choose a nonnegative test cutoff \(\chi\), equal to one on \(B_{2Sr_i}(a)\) and supported in \(B_{4Sr_i}(a)\), with \(|\Delta\chi|\le Cr_i^{-2}\). Since \(h=q^a-f_i^a\) and \(|\nabla q^a|\le8Sr_i\), we have \[\int\left|\tfrac12|\nabla f_i^a|^2-q^a\right|\,\mathrm d\mathfrak m \le\int\bigl(|\nabla q^a||\nabla h| +\tfrac12|\nabla h|^2\bigr)\,\mathrm d\mathfrak m \le Cr_i^{n+2}P_i(a).\] The integral can be restricted to the cutoff support or taken over the Dirichlet domain. Testing the definition of \(M\) and using \(\Delta q^a=a_i\mathfrak m-\mu_i^a\) now gives \[\begin{align*} \int\chi\,\mathrm d\mathcal B_\kappa(f_i^a,1) &=\frac12\int|\nabla f_i^a|^2\Delta\chi\,\mathrm d\mathfrak m +\kappa\int\chi|\nabla f_i^a|^2\,\mathrm d\mathfrak m -(2a_i-n)\int\chi\,\mathrm d\mathfrak m\\ &=(n-a_i)\int\chi\,\mathrm d\mathfrak m-\int\chi\,\mathrm d\mu_i^a +\kappa\int\chi|\nabla f_i^a|^2\,\mathrm d\mathfrak m +O(r_i^nP_i(a))\\ &\le Cr_i^nP_i(a). \end{align*}\] The first two terms in the penultimate line are nonpositive, and the curvature term is bounded by \(C\kappa r_i^{n+2}\), already included in \(Cr_i^nP_i(a)\). Positivity of \(\mathcal B_\kappa(f_i^a,1)\) proves (67). ◻

Proof. Local test regularity and the Hessian product rule give \[\nabla g_a=(\operatorname{Hess}f)(\nabla a,\cdot) +(\operatorname{Hess}a)(\nabla f,\cdot).\] For a test function \(\phi\), express each of the two integrated Hessian terms by differentiating \(\langle\nabla a,\nabla\phi\rangle\) or \(\langle\nabla f,\nabla\phi\rangle\) and integrating by parts. The terms with \(\operatorname{Hess}\phi\) cancel, leaving \[\begin{align*} &\int\operatorname{Hess}a(\nabla f,\nabla\phi)\,\mathrm d\mathfrak m -\int\operatorname{Hess}f(\nabla a,\nabla\phi)\,\mathrm d\mathfrak m\\ &\hspace{1cm} =-\int\Delta f\,\langle\nabla a,\nabla\phi\rangle\,\mathrm d\mathfrak m +\int\Delta a\,\langle\nabla f,\nabla\phi\rangle\,\mathrm d\mathfrak m =0. \end{align*}\] The constancy of \(\Delta f\) and harmonicity of \(a\) are both used in the last equality. This proves the identity with \(\operatorname{Hess}f\) on the right. Subtracting the metric tensor changes its right side by \(2\int\langle\nabla a,\nabla\phi\rangle=0\). Both sides are continuous in the \(W^{1,2}_0\) norm, by the bounded gradient factors and the local \(L^2\) Hessian bounds. Density therefore extends the identity to the stated tests. ◻

Proof. In a Euclidean ball, interior derivative estimates and Taylor’s formula at the center give \[\|w-w(0)-\nabla w(0)\cdot X\|_{2,B_{S\theta}}^{\mathrm{av}} \le C_n\theta^2\|w\|_{2,B_S}^{\mathrm{av}}, \qquad \theta<1/4.\] Choose \(m\) so that \(C_n\theta^2<\theta^{3/2}/4\). If the asserted conclusion failed along increasingly late spaces, rescale by \(r_i\), normalize the nonzero norm of \(w\) to one, and apply interior harmonic compactness (Ambrosio and Honda 2018). The limit is Euclidean harmonic on \(B_S\), with norm at most one there. The normalized coordinate columns converge locally uniformly to orthonormal affine Euclidean coordinates, by Lemma 5. An affine combination of these columns therefore approximates the affine Taylor polynomial of the limit uniformly on \(B_{S\theta}\). Strong convergence on this smaller ball contradicts the strict Euclidean estimate. This contradiction gives a single late-space threshold for every center, scale, and harmonic \(w\). If the normalizing norm is zero, the assertion is immediate. ◻

Proof. Harmonic correction and affine fits. We first work in the unrescaled units of the current member of the sequence and fix a center \(a\) and a covector \(b\). Define \[g_b=\langle\nabla q^a,\nabla(b\cdot u)\rangle, \qquad g_{j,b}=\langle\nabla f_j^a,\nabla(b\cdot u)\rangle.\] The Poisson error and (57) imply \[ \|g_b-g_{j,b}\|_{2,B_{Sr_j}(a)}^{\mathrm{av}} \le Cr_j\sigma_j^a(b)\sqrt{P_j(a)}. \tag{73}\] By Lemma 18, the equation for \(g_{j,b}\) on this ball has divergence source \(2(\operatorname{Hess}f_j^a-g)(\nabla(b\cdot u),\cdot)\). Its \(L^2\) norm is at most \(Cr_j^{n/2}\sigma_j^a(b)\sqrt{P_j(a)}\). Solve this divergence equation with zero boundary values to obtain \(\xi_{j,b}\in W^{1,2}_0(B_{Sr_j}(a))\). Testing with \(\xi_{j,b}\) and then applying Poincaré gives \[ \|\xi_{j,b}\|_{2,B_{Sr_j}(a)}^{\mathrm{av}} +r_j\|\nabla\xi_{j,b}\|_{2,B_{Sr_j}(a)}^{\mathrm{av}} \le Cr_j\sigma_j^a(b)\sqrt{P_j(a)}. \tag{74}\] In particular \(g_{j,b}-\xi_{j,b}\) is harmonic on that ball.

Let \(F_{j,b}\) be the \(L^2(B_{Sr_j}(a))\) orthogonal projection of \(g_b\) onto the finite-dimensional space of affine combinations of \(u\), and set \[e_j(b)=r_j^{-1} \|g_b-F_{j,b}\|_{2,B_{Sr_j}(a)}^{\mathrm{av}}.\] The affine space is nondegenerate by harmonic normalization, so these projections exist uniquely as functions and depend linearly on \(b\). Apply Lemma 19 to \(g_{j,b}-\xi_{j,b}-F_{j,b}\). The two errors in (73) and (74), restricted to the smaller ball, acquire only the fixed volume factor depending on \(m,S\). Dividing by \(r_{j+m}=\theta r_j\) gives \[ e_{j+m}(b)\le\theta^{1/2}e_j(b) +C\sigma_j^a(b)\sqrt{P_j(a)}. \tag{75}\] The factor \(1/2\) in the harmonic improvement leaves more than the needed margin in the coefficient of \(e_j(b)\).

Qualitative convergence. We will also use the uniform estimate \[ e_j(b)\le\varepsilon\sigma_j^a(b),\qquad\varepsilon=o(1). \tag{76}\] Here are details concerning the distance gradients in this assertion. In units \(r_j=1\), put \(q=q^a/r_j^2\), normalize the coordinate column at \(a\), and write \(V=R_j(a)^{-1/2}(u-u(a))/r_j\). On every fixed interior ball, \(V\) converges to orthonormal Euclidean coordinates and \(q^a/r_j^2-|V|^2/2\to0\) uniformly. The norms of the gradients of these two functions have the same limiting squared integrals: for the distance function use \(|\nabla(d_a^2/2)|=d_a\), and for \(|V|^2/2\) use the strong Gram convergence. Their cross term has the same limit. Indeed, for a compactly supported Lipschitz cutoff \(\chi\), \[\int\chi\langle\nabla(q-|V|^2/2),\nabla(|V|^2/2)\rangle\,\mathrm d\mathfrak m =-\int(q-|V|^2/2) \operatorname{div}(\chi\nabla(|V|^2/2))\,\mathrm d\mathfrak m=o(1),\] because \(\Delta(|V|^2/2)=\sum_k|\nabla V_k|^2\) and the relevant gradients are uniformly bounded. Thus the two gradients are close in local \(L^2\). Taking scalar products with \(\nabla V\) and using strong Gram convergence shows \[\langle\nabla q,\nabla V\rangle-V\longrightarrow0 \quad\hbox{in local }L^2.\] This supplies an affine candidate in the defining infimum for \(e_j(b)\). A contradiction subsequence over the center, scale, and unit normalized covector makes the estimate uniform, proving (76).

Summable squared errors. Choose a positive number \(\eta=o(1)\) with \(\eta\ge\varepsilon^2\); for example, add the reciprocal of the sequence index to \(\varepsilon^2\). Fix \(\gamma=1/8\). Iterating (75) in each of the finitely many residue classes modulo \(m\), and using (57), yields \[ e_i(b)\le C\sigma_i^a(b)\sqrt{\widetilde P_i(a)},\qquad \widetilde P_i(a)=\eta2^{-\gamma i} +\sum_{j=0}^i2^{-\gamma(i-j)}P_j(a). \tag{77}\] To see the exponent explicitly, the coefficient of a source at \(j<i\) before squaring is bounded by \(C_m2^{-(1/2-a_0)(i-j)}\). The initial error has the same decay times \(\varepsilon\). Weighted Cauchy–Schwarz bounds the square of the source sum by an exponential convolution of \(P_j\). Since \(a_0=1/16\) and \(\gamma=1/8<1/2-a_0\), weakening that convolution gives exactly (77). This calculation also covers the finitely many initial indices.

Removing the constant term. The constants in the affine fits must be removed with a quantitative bound. Put \(c_j=F_{j,b}(a)\). On \(B_{Sr_{j+1}}(a)\), the difference of consecutive fits has averaged \(L^2\) norm at most \(Cr_j(e_j(b)+e_{j+1}(b))\). Uniform affine nondegeneracy, or harmonic interior evaluation on this ball, therefore gives \[|c_j-c_{j+1}|\le Cr_j(e_j(b)+e_{j+1}(b)).\] Moreover \(|c_j|\le Cr_j\sigma_j^a(b)\). Indeed the raw radial function has this averaged size, its projection has no larger norm, and affine evaluation at the center is uniformly bounded. The growth exponent \(a_0<1\) implies \(c_j\to0\) as \(j\to\infty\). Consequently, using (77), \[\begin{align*} \frac{|c_i|}{r_i\sigma_i^a(b)} &\le C\sum_{l\ge i}2^{-(1-a_0)(l-i)} \bigl(\sqrt{\widetilde P_l(a)} +\sqrt{\widetilde P_{l+1}(a)}\bigr). \tag{78}\end{align*}\] Fix \(\beta=1/16\) and define, with a sufficiently large fixed constant \(C_E\), \[ E_i(a)=C_E\sum_{l=0}^\infty 2^{-\beta|i-l|}\widetilde P_l(a). \tag{79}\] Weighted Cauchy–Schwarz in (78) then gives \(|c_i|\le Cr_i\sigma_i^a(b)\sqrt{E_i(a)}\). The two geometric convolution kernels have finite sums, so \[\sup_a\sum_iE_i(a) \le C\left(\eta+\sup_a\sum_iP_i(a)\right)=o(1).\] Also \(E_i(a)\ge P_i(a)\) after increasing \(C_E\), and the term \(l=0\) gives the center-independent positive floor \(E_i(a)\ge C_E\eta2^{-\beta i}\).

Gradient control and normalization. Combining the value estimates already proved gives \[\|g_{i,b}-(F_{i,b}-c_i)\|_{2,B_{Sr_i}(a)}^{\mathrm{av}} \le Cr_i\sigma_i^a(b)\sqrt{E_i(a)}.\] The difference \(g_{i,b}-F_{i,b}\) satisfies the divergence equation of Lemma 18, because \(F_{i,b}\) is harmonic. Test it against its product with a squared Lipschitz cutoff, equal to one on \(B_{Sr_i/2}(a)\). Cauchy–Schwarz and absorption give the interior energy estimate \[\|\nabla(g_{i,b}-F_{i,b})\|_{2,B_{Sr_i/2}(a)}^{\mathrm{av}} \le C\sigma_i^a(b)\sqrt{E_i(a)}.\] Here the value term is divided by \(r_i\), and the divergence source has the bound used in (74).

For \(T=R_i(a)^{-1/2}\), apply these scalar estimates to its rows, rescale, and stack the resulting affine fits. Their constants at \(a\) have been removed, so they give a matrix \(A_i^a(T)\) and (72). The qualitative radial convergence proved above, the Poisson error, and affine nondegeneracy imply \(\|A_i^a(T)-\mathrm{Id}\|\le C(\varepsilon+\sqrt{P_i(a)}+ \sqrt{E_i(a)})=o(1)\) uniformly. Finally, replacing \(R_i(a)^{-1/2}\) by \(R_i(y)^{-1/2}\) conjugates this matrix by \(R_i(y)^{-1/2}R_i(a)^{1/2}\). This matrix and its inverse have uniformly bounded norms by fixed-distance normalization comparability. Thus both the quantitative bound and \(A_i^a(T)-\mathrm{Id}=o(1)\) persist. Translating \(U\) changes neither side when its center value is subtracted. ◻

Proof. The first assertions follow from the common value of the Laplacian, the Poisson gradient bounds, and (62); see Remark 17. In particular \[\|\operatorname{Hess}v\|_{L^2(B_{8R}(z))} +\|\operatorname{Hess}f^z-g\|_{L^2(B_{8R}(z))}\le Cs.\] Let \(A_x,A_z\) be the matrices in Lemma 20 for the common normalization \(R_i(y)^{-1/2}\). Put \[\delta=U(x)-U(z),\qquad B=A_z-A_x,\qquad c=A_x\delta.\] Since \(U(z)=0\), subtraction of the two radial-invariance identities gives \[ h=BU+c+e,\qquad \|e\|_{W^{1,2}(B_{8R}(z))}\le Cs. \tag{83}\] All coordinate gradients are bounded on this ball, and the values of \(U\), \(\delta\), and \(c\) are bounded by constants depending on the fixed radii. Also \(A_x,A_z=\mathrm{Id}+o(1)\).

Set \(L_vw=\langle\nabla v,\nabla w\rangle\). Differentiate the radial identity \(L_zU=A_zU+O_{W^{1,2}}(s)\) along \(\nabla v\). Its error remains \(O_{L^2}(s)\) because \(|\nabla v|\) is bounded, so \[L_vL_zU=A_zh+O_{L^2}(s).\] The Hessian product rule separately gives the exact identity \[L_vL_zU-L_zh =\operatorname{Hess}f^z(\nabla v,\nabla U) -\operatorname{Hess}v(\nabla f^z,\nabla U).\] The last displayed Hessian bounds and the bounded gradient factors turn its right side into \(h+O_{L^2}(s)\). Thus \[ L_zh=(A_z-\mathrm{Id})h+O_{L^2}(s). \tag{84}\] Every derivative here is legitimate in \(L^2\) by local test regularity and the gradient-product rule; no classical differentiability of the metric is being used.

Insert (83) on both sides of (84), and use the radial identity once more to differentiate \(BU\). This gives \[\bigl(BA_z-(A_z-\mathrm{Id})B\bigr)U-(A_z-\mathrm{Id})c =O_{L^2}(s).\] Affine nondegeneracy, for example on \(B_R(z)\), bounds its linear and constant coefficients separately: \[ \|BA_z-(A_z-\mathrm{Id})B\|+|(A_z-\mathrm{Id})c|\le Cs. \tag{85}\] Writing \(C_z=A_z-\mathrm{Id}\), the first operator is \(B+BC_z-C_zB\). Consequently \[(1-2\|C_z\|)\|B\| \le\|B+BC_z-C_zB\|\le Cs.\] For sufficiently late spaces, \(\|C_z\|<1/4\), and hence \(\|B\|\le Cs\). Next \[(A_x-\mathrm{Id})c=(A_z-\mathrm{Id})c-Bc=O(s).\] Since \(c=A_x\delta\) and \(A_x^{-1}\) is bounded, this yields \[ \delta-c=-A_x^{-1}(A_x-\mathrm{Id})c=O(s). \tag{86}\] Finally (83), the bounded coordinate values, and (86) prove (82) on the full ball \(B_{8R}(z)\). ◻

Proof. The Sobolev product rule and \(\Delta f^z=a_i\) give \[ \operatorname{div}(w^2\nabla f^z) =(a_i+2)w^2+wF,\qquad F=2(\ell\cdot U+k+e_w), \tag{89}\] with \(\|F\|_{L^2(B_{3R}(z))}\le C(\alpha+s)\). The identity is a weak identity: boundedness of \(w\), its Sobolev regularity, and the bounded gradient of \(f^z\) justify the product and its testing with compactly supported Lipschitz functions.

Here is a particular cutoff with the required sign. For \(t\ge0\) let \(k_R(t)\) equal one for \(t\le2R\), decrease linearly to zero on \([2R,3R]\), and equal zero for \(t\ge3R\). Choose \(A_0=2(n+4)/(R-1)\) and put \[\phi(t)= \begin{cases} 1,&0\le t\le R,\\ e^{-A_0(t-R)}k_R(t),&R<t<3R,\\ 0,&t\ge3R. \end{cases}\] This is a bounded nonnegative Lipschitz function. It is bounded below by \(e^{-A_0R}>0\) on \([0,2R]\), and on \(R<t<3R\) it satisfies \(-\phi'(t)\ge A_0\phi(t)\). Since \(n\le a_i\le n+1\) and \(R>1\), the choice of \(A_0\) implies \[ -t\phi'(t)-(a_i+2)\phi(t) \ge |\phi'(t)|+\phi(t) \quad\hbox{for almost every }R<t<3R. \tag{90}\] The derivative and the cutoff are bounded by fixed constants.

Write \(t=d_z\) and \[e_f=\langle\nabla d_z,\nabla f^z\rangle-d_z, \qquad \mathcal A=\int_{R<d_z<3R}(|\phi'|+\phi)w^2\,\mathrm d\mathfrak m.\] The Poisson gradient error gives \(\|e_f\|_{L^2(B_{3R}(z))}\le Cs\), because \(\nabla(d_z^2/2)=d_z\nabla d_z\). Test (89) against \(\phi(d_z)\). Using (90), and moving the inner-ball term to the right, gives \[ \mathcal A\le C\alpha^2 +\left|\int\phi wF\,\mathrm d\mathfrak m\right| +\left|\int_{R<d_z<3R}\phi'e_f w^2\,\mathrm d\mathfrak m\right|. \tag{91}\] The decisive coefficient error has the bound \[\begin{align*} \left|\int_{R<d_z<3R}\phi'e_f w^2\,\mathrm d\mathfrak m\right| &\le\|w\|_\infty \left(\int|\phi'|w^2\,\mathrm d\mathfrak m\right)^{1/2} \left(\int|\phi'|e_f^2\,\mathrm d\mathfrak m\right)^{1/2}\\ &\le Cs\sqrt{\mathcal A} \le\tfrac14\mathcal A+Cs^2. \tag{92}\end{align*}\] All integrals in the middle line are over the annulus. For the \(F\) term, the inner ball contributes at most \(C\alpha(\alpha+s)\), while the annulus contributes at most \(C\sqrt{\mathcal A}(\alpha+s)\). Another application of Young’s inequality therefore gives \[\left|\int\phi wF\,\mathrm d\mathfrak m\right| \le\tfrac14\mathcal A+C(\alpha^2+s^2).\] Substitution in (91) proves \(\mathcal A\le C(\alpha^2+s^2)\). The positive lower bound for \(\phi\) up to radius \(2R\) now gives (88).

In particular, the fixed bound for \(\|w\|_\infty\) in (92) is used before \(w\) is divided by \(\alpha\). The conclusion is quadratic in \(s\) at the level of squared norms. ◻

Proof. Subtract a constant from \(v\) so that \(v(z)=0\); this does not change any gradient, Laplacian, or Bochner measure. Its bounded gradient gives a uniform bound for its values on \(B_{8R}(z)\). Let \(c_0+b\cdot U\) be its \(L^2(B_R(z))\) orthogonal projection onto the affine coordinate span, and put \[w=v-c_0-b\cdot U,\qquad \alpha=\|w\|_{L^2(B_R(z))}.\] Affine nondegeneracy and the bounds for \(v,U\) give \(|c_0|+|b|\le C\). The function \(w\) is harmonic on \(B_{8R}(z)\) and has uniformly bounded values and gradients there.

The gradient product rule and the Poisson estimates yield \[ \nabla(L_zv-v) =(\operatorname{Hess}f^z-g)(\nabla v,\cdot) +\operatorname{Hess}v(\nabla f^z,\cdot), \qquad \|\nabla(L_zv-v)\|_{L^2(B_{8R}(z))}\le Cs. \tag{94}\] Poincaré on \(B_{4R}(z)\), using the gradient on \(B_{8R}(z)\), therefore gives a constant \(k_v\) such that \(L_zv-v=k_v+O_{L^2(B_{4R}(z))}(s)\). Subtracting the projection, and using \(L_zU=A_zU+O_{W^{1,2}}(s)\), gives \[ L_zw=w+\ell\cdot U+k+e_w, \qquad \|e_w\|_{L^2(B_{4R}(z))}\le Cs, \tag{95}\] where \(\ell^T=b^T(\mathrm{Id}-A_z)\) and \(k=k_v+c_0\). We need a bound for these coefficients in terms of the small residual norm. Caccioppoli for the harmonic function \(w\) gives \[\|\nabla w\|_{L^2(B_{R/2}(z))}\le C\alpha.\] The bounded gradient of \(f^z\) and (95) now bound the \(L^2\) norm of \(\ell\cdot U+k\) on this smaller ball by \(C(\alpha+s)\). Affine nondegeneracy gives \[ |\ell|+|k|\le C(\alpha+s). \tag{96}\] Lemma 22 consequently applies and yields \[ \|w\|_{L^2(B_{2R}(z))}\le C(\alpha+s). \tag{97}\]

We claim that \(\alpha\le Cs\) with a uniform constant for all sufficiently late spaces. If no such constant and threshold existed, choose a sequence of increasingly late spaces and configurations with \(\alpha/s\to\infty\). The positive floor in (71) makes \(s>0\), so \(\widehat w=w/\alpha\) is defined along this sequence. The norm on \(B_R(z)\) is one, and (97) bounds its norm on \(B_{2R}(z)\) uniformly. Interior harmonic compactness therefore gives local uniform and strong \(W^{1,2}\) convergence on \(B_{3R/2}(z)\) to a Euclidean harmonic function \(w_\infty\); see (Ambrosio and Honda 2018). Use a rotation so that the normalized columns \(U\), centered at \(z\), converge to the Euclidean coordinate column \(X\).

The coefficients \(\ell/\alpha,k/\alpha\) are bounded by (96); pass to a convergent subsequence. The error \(e_w/\alpha\) tends to zero in \(L^2\). The gradients of \(f^z\) converge strongly to the radial vector field \(X\): the Poisson gradient error tends to zero, since \(P_i(z)\le s^2\to0\), and the squared-distance gradient convergence was proved in Lemma 20. The normalized harmonic gradients have uniform interior bounds and converge strongly. Thus scalar products of these gradients pass to the limit in (95), which becomes \[ X\cdot\nabla w_\infty =w_\infty+\ell_\infty\cdot X+k_\infty \quad\hbox{on }B_{3R/2}(0). \tag{98}\] One may first pass this identity against interior tests and then use the smoothness of the Euclidean harmonic limit.

Uniform convergence on a neighborhood of \(\overline B_R(0)\) and convergence of the measures imply \[\|w_\infty\|_{L^2(B_R(0))}=1,\qquad \int_{B_R(0)}w_\infty\,\mathrm dX=0,\qquad \int_{B_R(0)}w_\infty X_j\,\mathrm dX=0\quad(1\le j\le n).\] The boundary sphere has zero Euclidean measure, so restriction to the projection ball causes no loss in these limits. Expand \(w_\infty\) in its harmonic Taylor series at zero. The operator \(X\cdot\nabla-1\) multiplies the homogeneous term of degree \(q\) by \(q-1\). Equation (98) therefore forces all terms of degree at least two to vanish. Analyticity makes \(w_\infty\) affine on the ball. Its affine orthogonality would then give zero norm, a contradiction. This proves the claimed uniform bound \(\alpha\le Cs\).

Finally Caccioppoli on the projection ball gives \(\|\nabla w\|_{L^2(B_{R/2}(z))}\le Cs\). Since \(B_{2L+4}(y)\subset B_{R/2}(z)\) and the constant \(c_0\) does not affect gradients, this is (93). ◻

Proof. Continue first in the rescaled units used in the preceding subsection, and choose \(b\) from Lemma 23. The normalization at \(y\) is exactly \[\mathop{\mathchoice {\displaystyle\int\mkern-18mu-} {\textstyle\int\mkern-15mu-} {\scriptstyle\int\mkern-12mu-} {\scriptscriptstyle\int\mkern-10mu-}}\nolimits_{B_1(y)} (\langle\nabla U_j,\nabla U_k\rangle)_{j,k}\,\mathrm d\mathfrak m=\mathrm{Id}.\] Consequently \[\begin{align*} \mathop{\mathchoice {\displaystyle\int\mkern-18mu-} {\textstyle\int\mkern-15mu-} {\scriptstyle\int\mkern-12mu-} {\scriptscriptstyle\int\mkern-10mu-}}\nolimits_{B_1(y)}h\,\mathrm d\mathfrak m &=b+\mathop{\mathchoice {\displaystyle\int\mkern-18mu-} {\textstyle\int\mkern-15mu-} {\scriptstyle\int\mkern-12mu-} {\scriptscriptstyle\int\mkern-10mu-}}\nolimits_{B_1(y)} \langle\nabla U,\nabla(v-b\cdot U)\rangle\,\mathrm d\mathfrak m =b+O(s). \end{align*}\] Here the gradient bound for \(U\) and (93) control the last error. Equation (82) then gives \[ b=U(x)-U(z)+O(s)=\delta+O(s). \tag{100}\]

The harmonic function \(w=v-c_0-b\cdot U\) has gradient energy at most \(Cs^2\) on \(B_{2L+4}(y)\). Choose a nonnegative test cutoff \(\zeta\), equal to one on \(B_L(y)\) and supported in \(B_{2L+2}(y)\), with bounded Laplacian. Testing the definition of the nonnegative measure \(M(w)\) gives \[\begin{align*} M(w)(B_L(y)) &\le\frac12\int|\nabla w|^2\Delta\zeta\,\mathrm d\mathfrak m +\kappa r_i^2\int\zeta|\nabla w|^2\,\mathrm d\mathfrak m\\ &\le Cs^2. \tag{101}\end{align*}\] Together with (81), positivity and polarization imply \(M(b\cdot U)(B_L(y))\le Cs^2\). The normalized mass \(M(U)(B_L(y))\) has bounded trace by Lemma 6. Therefore (100) implies \[M((\delta-b)\cdot U)(B_L(y)) \le|\delta-b|^2\operatorname{tr}M(U)(B_L(y))\le Cs^2.\] Another positive quadratic-form triangle inequality gives \[ M(\delta\cdot U)(B_L(y))\le Cs^2. \tag{102}\] In particular the \(O(s)\) coefficient error has contributed \(O(s^2)\) mass, including the possibly singular part of the Bochner measure.

We finish by recording the scale conversion explicitly. Return to the unrescaled units, set \(T=R_i(y)^{-1/2}\) and \(\delta u=u(x)-u(z)\). Under \(d'=d/r_i\), \(\mathfrak m'=r_i^{-n}\mathfrak m\) and \(U=Tu/r_i\) up to an additive constant, the matrix-valued measure transforms as \[ M'(U)=r_i^{2-n}T M(u)T^T. \tag{103}\] Indeed the squared gradient matrix of \(U\) in the primed metric is \(TGT^T\), the Laplacian gains a factor \(r_i^2\), and the lower curvature parameter becomes \(\kappa r_i^2\). Meanwhile \(\delta=U(x)-U(z)=T\delta u/r_i\). Hence, on the corresponding balls, \[\begin{align*} M'(\delta\cdot U)(B_L'(y)) &=r_i^{-n}\delta u^TR_i(y)^{-1} M(u)(B_{Lr_i}(y))R_i(y)^{-1}\delta u\\ &=r_i^{-2}D_{i,y}[\delta u]. \end{align*}\] Applying (102) and \(s^2=E_i(x)+E_i(z)\) proves (99). ◻

Proof. All sums in this proof are over \(\mathcal A_i(x)\). Put \(w_y=P_x-c_{i,y}\). Equation (109) and the fixed choice (112) give \[ D_{i,y}[w_y]\le C r^2E. \tag{120}\] The same center-distance bounds and harmonic normalization give \(|w_y|_H\le Cr\). Since \(0\le D_{i,y}\le Z\le\varepsilon H\), conjugation by \(H^{-1/2}\) gives \[ |D_{i,y}w|_Q^2\le\varepsilon D_{i,y}[w] \quad\text{and}\quad |\xi^TD_{i,y}w|^2\le D_{i,y}[\xi]D_{i,y}[w]. \tag{121}\] Here the second inequality is Cauchy–Schwarz for a positive semidefinite form and does not require its invertibility. Applying the first inequality to (120) and using the bounded number of summands proves (114).

Differentiation of (113) gives \[e(\xi)=\sum_y (\mathrm D\phi_{i,y})\xi\ D_{i,y}w_y.\] By normalized cutoff bounds and matrix comparability, \[ |(\mathrm D\phi_{i,y})\xi|\le Cr^{-1}|\xi|_H, \qquad |\mathrm D^2\phi_{i,y}(\xi,\xi)|\le Cr^{-2}|\xi|_H^2. \tag{122}\] These inequalities and (121) prove the bound for \(e\) in (116). Since \(0\le S\le Z\le\varepsilon H\), the same conjugation gives \(SQS\le\varepsilon S\le\varepsilon Z\), proving its bound for \(S\). For the sharper diagonal bound, use the second inequality in (121) before summing: \[|\xi^Te(\xi)| \le C\sqrt E\,|\xi|_H\sum_y\sqrt{D_{i,y}[\xi]} \le C\sqrt{Z[\xi]E}\,|\xi|_H.\] Finally, \[\mathrm D^2v_i(\xi,\xi) =\sum_y\left( \mathrm D^2\phi_{i,y}(\xi,\xi)D_{i,y}w_y +2(\mathrm D\phi_{i,y})\xi\ D_{i,y}\xi\right).\] Equations (121)–(122) bound the two sums by \(Cr^{-1}\sqrt\varepsilon\sqrt E\,|\xi|_H^2\) and \(Cr^{-1}\sqrt\varepsilon\sqrt{Z[\xi]}\,|\xi|_H\), respectively. Dropping \(\sqrt\varepsilon\le1\) proves (118).

For the isotropic estimates, \(D_{i,y}\le\varepsilon H\) also gives \(\|D_{i,y}\|_{H\to Q}\le\varepsilon\) and \(|D_{i,y}w_y|_Q\le C\varepsilon r\). In order \(m\ge1\), the product in (113) has only the term with \(m\) derivatives on the cutoff and the \(m\) terms with one derivative on the affine factor. By (110), the frozen source-norm comparison, and bounded overlap, their total norm is at most \[C\varepsilon r^{1-m} \left(C_1^{m+1}(m!)^s +m C_1^m((m-1)!)^s\right).\] Increase \(C_v\) once to obtain (119). The preceding value bound with \(|w_y|_H\le Cr\) proves its \(m=0\) case. ◻

Proof. For \(1\le l\le m-1\) the summand in the first expression of (127) is \[\frac{(l+1)(m-l+1)}{m+1} \left(\frac{(l+1)!(m-l+1)!}{(m+1)!}\right)^\sigma.\] The two terms with \(l=1,m-1\) are \(O_\sigma(m^{-\sigma})\); the terms with \(l=2,m-2\) are \(O_\sigma(m^{-2\sigma})\). For \(3\le l\le m-3\), the factorial ratio is at most \(24/((m+1)m(m-1))\): its maximum in this interval occurs at an endpoint, as is seen by taking the ratio of consecutive terms. The first factor is at most \(m+1\), and there are at most \(m\) such terms. Their sum is therefore \(O_\sigma(m^{2-3\sigma})\), which is \(O_\sigma(m^{-\sigma})\) because \(\sigma\ge1\). This proves \(\delta_m\to0\); the finitely many small values of \(m\) do not matter. Including the terms \(l=0,m\) proves the second assertion of (127).

Let \(K=J^{-1}\). The multilinear product formula, followed by the operator norm inequality and division by \(m!\), gives \[\frac{\|\mathrm D^m K\|}{m!} \le\sqrt A\sum_{l=1}^m \frac{\|\mathrm D^l J\|}{l!} \frac{\|\mathrm D^{m-l}K\|}{(m-l)!}.\] In particular the binomial coefficients in the product formula have been removed exactly. After dividing out \(B^m\), universal majorants \(a_m\) for the positive-order coefficients satisfy \[ a_m\le A L_m+\sqrt A \sum_{l=1}^{m-1}L_l a_{m-l}. \tag{129}\] Choose \(M=M(A,\sigma)\) so that \(\sqrt A\,\delta_m\le1/2\) for \(m\ge M\). The recursion for the finitely many orders \(m<M\) gives finite universal bounds depending only on \(A\) and \(\sigma\). Choose \(C_{A,\sigma}\) to dominate \(\sqrt A\), \(2A\), and all these bounds after division by \(L_m\). If \(a_q\le C_{A,\sigma}L_q\) for \(q<m\), then (129) gives, for \(m\ge M\), \[\frac{a_m}{L_m} \le A+\sqrt A\,C_{A,\sigma}\delta_m \le A+\frac12 C_{A,\sigma} \le C_{A,\sigma}.\] This proves (128) with a single constant for all orders; no factor \(C_{A,\sigma}^m\) has been introduced. Transposition preserves Euclidean operator norms and the same derivative estimates. ◻

Proof. Freeze \(x\) and abbreviate \(r=r_i\), \(H=H_i(x)\), \(Q=H^{-1}\). The linear source and target normalizations are \[\begin{aligned} P&=P_x+rH^{-1/2}X,\\ \widetilde h(X) &=r^{-1}\bigl(h_i(P_x+rH^{-1/2}X)-h_i(P_x)\bigr),\\ \widetilde v(X)&=r^{-1}H^{-1/2}v_i(P_x+rH^{-1/2}X). \end{aligned}\] At \(X=0\), \(\widetilde J=\mathrm D\widetilde h=J_iH^{-1/2}\) has singular values between \(A^{-1/2}\) and \(A^{1/2}\), and \[\widetilde J^{-T}\widetilde v =r^{-1}J_i^{-T}v_i.\] This identity checks the covector normalization in the update. All matrices in the change of variables are constant matrices.

For every \(m\ge1\), (126) gives \[\frac{\|\mathrm D^m\widetilde J(0)\|}{m!} \le\frac{B^m((m+1)!)^{1+\sigma}}{m!} =B^m L_m.\] Lemma 27 therefore bounds the coefficients of \(\widetilde J^{-T}\) by \(C_{A,\sigma}B^mL_m\). On the other hand, (119) gives \[\frac{\|\mathrm D^q\widetilde v(0)\|}{q!} \le\varepsilon C_v^{q+1}(q!)^{s-1} \le\varepsilon C_v B^qL_q\qquad(q\ge0).\] The product rule divided by \(m!\) and (127) now imply \[ \frac{\|\mathrm D^m(\widetilde J^{-T}\widetilde v)(0)\|}{m!} \le\varepsilon C_A B^m L_m. \tag{130}\]

To pass to next-scale jet units, note that \(H_{i+1}(x)\ge H_i(x)\). Thus the norm of any multilinear map with source norm \(H_{i+1}(x)\) is at most its norm with source norm \(H_i(x)\). The change from \(r_i\) to \(r_{i+1}=r_i/2\) supplies the additional factor \(2^{1-m}\). Relative to the allowance \(B^{m-1}(m!)^{1+\sigma}\) at order \(m\ge2\), the old contribution is at most \(2^{1-m}\le1/2\), while (130) contributes at most \[\varepsilon C_A B\,2^{1-m}(m+1)^{1+\sigma}.\] The supremum of \(2^{1-m}(m+1)^{1+\sigma}\) over \(m\ge2\) is finite. One sufficiently small choice of \(\varepsilon\) therefore makes the new contribution at most one half in every order simultaneously. This proves the lemma. ◻

Proof. We give the estimates under the induction hypothesis through level \(k-1\) and then close that induction. Fix \(x\) and a nonzero source vector \(\xi\). Every expression in this proof is evaluated at \(P_x\). Set \[q_i=|\xi|_{H_i(x)},\qquad b_i=\frac{r_i|\mathrm D^2h_i(\xi,\xi)|}{q_i},\qquad K_i=J_i^{-T},\qquad W_i=K_iv_i.\] The metric bound gives \(|K_iw|\le\sqrt A|w|_{Q_i(x)}\) for covectors \(w\), since \((J_i^TJ_i)^{-1}\le A Q_i(x)\). The jet bounds through order three, the ordinary inverse derivative formula, and the duality of \(H_i\) and \(Q_i\) give \[\begin{align*} |(\mathrm D_\xi K_i)w| &\le C_A\frac{q_i}{r_i}|w|_{Q_i(x)}, \tag{131}\\ |(\mathrm D_\xi^2 K_i)w| &\le C_A\frac{q_i^2}{r_i^2}|w|_{Q_i(x)}. \tag{132}\end{align*}\] For completeness, after the frozen normalization in the preceding lemma, the first formula is \(\mathrm D_\xi K_i=-K_i(\mathrm D_\xi J_i)^TK_i\). The second derivative is a sum of two products containing two first derivatives of \(J_i\) and one product containing \(\mathrm D_\xi^2J_i\). Bounds for \(\mathrm D^2h_i\) and \(\mathrm D^3h_i\), with the bounded normalized inverse, give exactly (131)–(132) after undoing the normalization. Thus these estimates are also valid for arbitrarily anisotropic \(H_i\).

Twice differentiating \(W_i\) gives all three types of terms: \[ \mathrm D^2W_i(\xi,\xi) =K_i\mathrm D^2v_i(\xi,\xi) +2(\mathrm D_\xi K_i)(\mathrm Dv_i\xi) +(\mathrm D_\xi^2 K_i)v_i. \tag{133}\] By Lemma 26 and (131)–(132), multiplication by \(r_i/q_i\) bounds each of these terms by a constant times \[f_i:=\sqrt{Z_i(x)[\xi]}+\sqrt{E_i(x)}\,q_i.\] More explicitly, the first term uses (118); the second uses \(|\mathrm Dv_i\xi|_{Q_i}\le C\sqrt\varepsilon f_i\); and the third uses \(|v_i|_{Q_i}\le C\sqrt\varepsilon r_i\sqrt{E_i(x)}\). Because \(q_{i+1}\ge q_i\) and \(r_{i+1}=r_i/2\), the old second derivative contracts by at least a factor of one half. Hence \[ b_{i+1}\le\tfrac12 b_i+C_A f_i, \qquad b_0=0. \tag{134}\] Convolution with the sequence \(1,1/2,1/4,\ldots\), whose sum is two, is bounded on \(\ell^2\). Thus, for every finite \(k\) covered by the induction hypothesis, \[\begin{align*} \sum_{i<k}b_i^2 &\le C_A\sum_{i<k}\left(Z_i(x)[\xi]+E_i(x)q_i^2\right) \\ &\le C_A H_k(x)[\xi]. \tag{135}\end{align*}\] The last line uses (107), monotonicity \(H_i\le H_k\), and (108). In particular it does not require that the matrices \(H_i\) remain bounded as \(i\to\infty\).

The first derivative of the increment satisfies \[\Delta_iJ\,\xi:=(J_{i+1}-J_i)\xi =K_i(\mathrm Dv_i\xi)+(\mathrm D_\xi K_i)v_i.\] The same estimates give the pointwise square bound \[ |\Delta_iJ\,\xi|^2 \le C_A\varepsilon \left(Z_i(x)[\xi]+E_i(x)q_i^2\right). \tag{136}\] Consequently their sum for \(i<k\) is at most \(C_A\varepsilon H_k(x)[\xi]\).

It remains to control the linear part of the increment of the squared length. The undifferentiated inverse gives \[\langle J_i\xi,K_i\mathrm Dv_i\xi\rangle =\xi^T\mathrm Dv_i\xi =S_i(x)[\xi]+\xi^Te_i(\xi).\] For the differentiated inverse, direct matrix multiplication gives \[\begin{align*} \langle J_i\xi,(\mathrm D_\xi K_i)v_i\rangle &=-\xi^T(\mathrm D_\xi J_i)^TJ_i^{-T}v_i \\ &=-\langle\mathrm D^2h_i(\xi,\xi),K_iv_i\rangle_{\mathbb R^n}. \tag{137}\end{align*}\] The inner product on the right is the ordinary Euclidean target inner product. Its absolute value is bounded by \(C_A b_iq_i\sqrt{E_i(x)}\), using (114) and dropping the additional factor \(\sqrt\varepsilon\le1\).

Now expand \(|J_{i+1}\xi|^2-|J_i\xi|^2\), sum over \(i<k\), and use (117), (136), and (137). The two linear error sums satisfy, respectively, \[\begin{align*} \sum_{i<k}|\xi^Te_i(\xi)| &\le C\left(\sum_{i<k}Z_i(x)[\xi]\right)^{1/2} \left(\sum_{i<k}E_i(x)q_i^2\right)^{1/2} \le C\sqrt\varepsilon H_k(x)[\xi],\\ \sum_{i<k}b_iq_i\sqrt{E_i(x)} &\le\left(\sum_{i<k}b_i^2\right)^{1/2} \left(\sum_{i<k}E_i(x)q_i^2\right)^{1/2} \le C_A\sqrt\varepsilon H_k(x)[\xi]. \end{align*}\] We have proved the quantitative cumulative error estimate \[ \left|\,|J_k\xi|^2-|\xi|^2 -2\sum_{i<k}S_i(x)[\xi]\,\right| \le C_A(\sqrt\varepsilon+\varepsilon)H_k(x)[\xi]. \tag{138}\] Only the sum of \(E_i\) has been used, never the sum of its square roots.

We now fix the constants and complete the simultaneous induction. By (107) and \(0\le S_i\le Z_i\), \[ C_*^{-1}H_k \le\mathrm{Id}+2\sum_{i<k}S_i \le2C_*H_k. \tag{139}\] Choose \(A\ge4C_*\), increasing it if necessary for the initial comparability of \(H_0\) and \(\mathrm{Id}\). Fix \(\sigma,B\) as in Lemma 28. Only now choose the data sufficiently flat that its smallness requirement holds and \[C_A(\sqrt\varepsilon+\varepsilon)<\frac1{2C_*}.\] Then (138) and (139) imply (125) with slack at level \(k\).

Initially \(h_0\) is affine, so all its positive-order derivatives beyond the first vanish, and both bootstrap statements hold. Suppose they hold through level \(i\). The pointwise inverse of \(J_i\) exists everywhere on \(\Omega\) and is smooth there; hence the update defines a smooth \(h_{i+1}\) on this same open set. The proof of (138) at \(k=i+1\) uses inverse and jet estimates only through level \(i\): the error sum uses \(b_j\) only for \(j\le i\), whose recurrence involves still earlier increments. It therefore proves the metric bound and invertibility at level \(i+1\) without assuming them. Lemma 28 proves all jet bounds at level \(i+1\) using the old-level data. This is an induction on the finite scale index with an all-orders assertion at each step, and completes the proof. ◻

Proof. From (114) and (125), \[|h_{i+1}(P_x)-h_i(P_x)| \le C_A\sqrt\varepsilon r_i\sqrt{E_i(x)} \le C_A\varepsilon r_i.\] The geometric series gives uniform convergence on \(B_1(p)\) and a tail estimate \[ \sup_{x\in B_1(p)}|F(x)-h_i(P_x)|\le\tau r_i, \qquad \tau=C_A\varepsilon\longrightarrow0. \tag{140}\] Each approximating map is continuous, so \(F\) is continuous.

Let \(K_2=B(2!)^{1+\sigma}\), the fixed second-derivative constant in (126). By (105), along any value-space segment contained in the central ellipsoid at \(x\), with preimages in \(B_{r_i}(x)\subset B_1(p)\), \[ |\mathrm D^2h_i(P_y)(\zeta,\zeta)| \le\frac{K_2C_{\mathrm{cmp}}}{r_i} |\zeta|_{H_i(x)}^2. \tag{141}\] Put \(K_T=K_2C_{\mathrm{cmp}}/2\). Choose a fixed \(D\ge1\) so large that \[ D>4/a,\qquad 4K_T/D\le1/(8\sqrt A). \tag{142}\] These choices depend only on constants already fixed in the construction.

Take a further sufficiently late rescaling, without changing those constants, so that the normalized harmonic approximation at the fixed scale ratio \(D\) gives \[ \tfrac12 d(x,z)\le|P_z-P_x|_{H_i(x)}\le2d(x,z) \quad\text{whenever }Dr_i^{-1}d(x,z)\in[1/2,1] \tag{143}\] on the interior balls in use. To justify this precise consequence, normalized harmonic approximation controls distances in value space up to \(o(1)r_i\) on the fixed ball; since the displayed distances are at least \(r_i/(2D)\), its relative error tends to zero after \(D\) has been fixed. The relative comparison of \(Q_i\) with the normalization matrices preserves these bounds. Also require \[ 4D\tau\le1/(8\sqrt A). \tag{144}\]

Set, for example, \(\rho=(16D+4)^{-1}\). If \(x,z\in B_\rho(p)\) are distinct and \(t=d(x,z)\), then \(t\le2\rho\). Choose a dyadic scale with \[Dt\le r_i<2Dt.\] It has \(r_i<4D\rho<1\), and \(B_{r_i}(x)\subset B_1(p)\). Thus every pair in the final ball, not merely sufficiently generic pairs, is covered by an available scale. Write \(\Delta=P_z-P_x\). Equation (143) gives \(t/2\le|\Delta|_{H_i(x)}\le2t\). By (142), \(|\Delta|_{H_i(x)}<a r_i\). The entire straight segment \[P_x+s\Delta,\qquad 0\le s\le1,\] therefore lies in the central ellipsoid in (105), and hence in \(u(B_{r_i}(x))\subset\Omega\). In particular every point where a derivative is integrated has an interior preimage and obeys the matrix comparison used in (141).

Taylor’s formula along this segment yields \[h_i(P_z)-h_i(P_x)=J_i(P_x)\Delta+\mathcal R, \qquad |\mathcal R|\le\frac{K_T}{r_i}|\Delta|_{H_i(x)}^2 \le\frac{4K_T}{D}t\le\frac{t}{8\sqrt A}.\] The metric bound gives \[\frac{t}{2\sqrt A}\le|J_i(P_x)\Delta|\le2\sqrt A\,t.\] Finally the two tails in (140) have total size at most \(2\tau r_i\le4D\tau t\le t/(8\sqrt A)\). Combining the last three inequalities proves, for every \(x,z\in B_\rho(p)\), \[ \frac{1}{4\sqrt A}d(x,z) \le|F(x)-F(z)|\le4\sqrt A\,d(x,z). \tag{145}\] The case \(x=z\) is immediate. In particular \(F\) is injective there and has a continuous inverse onto its image.

The set \(u(B_\rho(p))\) is open, by the previously established open coordinates. The continuous injective map \[F\circ u^{-1}:u(B_\rho(p))\longrightarrow\mathbb R^n\] therefore has open image by invariance of domain (Hatcher 2002, Theorem 2B.3). This proves the proposition. All scale choices used \(t=d(x,z)\) in the ambient space; no intrinsic path distance on the chosen neighborhood was substituted. ◻

Proof of Theorem 1. Let \(p\) be the prescribed \(n\)-regular point of the noncollapsed \(\mathrm{RCD}(K,n)\) space. The preceding sections provide the geometric, harmonic, heat, and directional-pinching conclusions on sufficiently late rescalings at \(p\). First fix the finitely many patch enlargement constants using only \(n\). Choose the cutoff with \(s=2\) and fix \(\sigma=2\). The cumulative constant \(C_*\), the bootstrap constants \(A,B\), and then \(D\) in (142) consequently depend only on \(n\). Choose a rescaling sufficiently late that all their flatness requirements, the simultaneous induction, and (144) hold. Proposition 31 gives an open ball about \(p\) and a chart with the required two-sided distance bound on that rescaled space. If its metric is \(d/r_*\), multiply the coordinate values of that chart by \(r_*\) when returning to the original metric. Both sides of (145) then scale by the same factor, while the domain remains an open neighborhood and the image remains an open subset of \(\mathbb R^n\). Thus Theorem 1 holds with \(L_n=4\sqrt{A(n)}\). The required rescaling, and hence the radius of the neighborhood in the original metric, may depend on \(X\) and \(p\). ◻