64c3814
==== TABLE OF CONTENTS ====
Bi-Lipschitz-Coordinates-at-Regular-Points-of-Noncollapsed-RCD-Spaces//coordinates.tex
Bi-Lipschitz-Coordinates-at-Regular-Points-of-Noncollapsed-RCD-Spaces//introduction.tex
Bi-Lipschitz-Coordinates-at-Regular-Points-of-Noncollapsed-RCD-Spaces//paper.tex
Bi-Lipschitz-Coordinates-at-Regular-Points-of-Noncollapsed-RCD-Spaces//pinching.tex
Bi-Lipschitz-Coordinates-at-Regular-Points-of-Noncollapsed-RCD-Spaces//renormalization.tex
Bi-Lipschitz-Coordinates-at-Regular-Points-of-Noncollapsed-RCD-Spaces//references.bib
===========================

==== BEGIN Bi-Lipschitz-Coordinates-at-Regular-Points-of-Noncollapsed-RCD-Spaces//coordinates.tex ====
\section{Uniform harmonic coordinates}
\label{co:section}

We first obtain an open coordinate domain.  The estimates in this
section allow the derivatives of these coordinates to degenerate as
the scale decreases.  Their degeneration will be measured by positive
matrices in the next section.

\subsection{Local calculus and rescaling}

For locally Sobolev functions write
\(\Gamma(f,g)=\langle\nabla f,\nabla g\rangle\) and
\(\Gamma(f)=\Gamma(f,f)\).  Laplacians of gradient squares will be
understood as measures.  Products with such measures always use the
continuous representative of a Lipschitz function, or the capacity
representative of a Sobolev function when that is the appropriate
meaning of the pairing.

We write \(g\) for the canonical metric tensor.

\begin{lemma}[Local analytic facts]\label{co:local-calculus}
Let \((Y,d,\m)\) be a noncollapsed \(\RCD(-\kappa,n)\) space and let
\(\Omega\subset Y\) be open.  A locally bounded solution of
\(\Delta w=a\) in \(\Omega\), with \(a\) constant, is locally Lipschitz
and is a local test function.  In particular it has a local
\(L^2\) Hessian, \(\Gamma(w_1,w_2)\in W^{1,2}_{\rm loc}(\Omega)\)
for two such functions, and the Hessian product rule holds locally.
The measure
\begin{equation}\label{co:M-definition}
 M(w)=\frac12\Delta\Gamma(w)+\kappa\Gamma(w)\m
\end{equation}
is locally finite and nonnegative.  For every constant \(b\),
\begin{equation}\label{co:shifted-bochner}
 M(w)-(2b\Delta w-nb^2)\m
 \ge |\Hess w-bg|^2\m.
\end{equation}
On the real vector space of local constant-Laplacian functions,
the left side is a nonnegative quadratic form in the pair
\((w,b)\), with both entries scaled together.
All these statements hold on interior balls for weak solutions with
locally bounded \(L^2\) norm, after using the interior estimates to
obtain local boundedness.

On fixed interior scales, bounded Lipschitz cutoffs with bounded
Laplacian are available.  Under uniform local volume comparison and
a uniform rescaled curvature bound, their bounds, local Poincare
and Sobolev constants, and interior harmonic and Poisson estimates
are uniform.  In particular, on a fixed-radius ball in the rescaled
spaces below there is a number \(d_*>1\) and a constant \(C\) such
that
\begin{equation}\label{co:dirichlet-sobolev}
 \|h\|_{L^2}+\|h\|_{L^{2d_*}}
 \le C\|\nabla h\|_{L^2},\qquad h\in W^{1,2}_0(B),
\end{equation}
whenever the ball and a fixed larger ball are contained in the
working region.
\end{lemma}

\begin{proof}
The full curvature-dimension condition implies the reduced condition
used in some of the references, by \cite[Proposition~2.5(i)]{BS10}.
The interior Poisson regularity and gradient estimates apply to a
constant right hand side; see \cite{Kell13,Jiang14}.  More precisely,
if \(8B_R(z)\) is compactly contained in the domain, the quantitative
estimate in \cite[Theorem~3.1]{Jiang14} 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
\cite[Lemmas~2.6--2.7]{CDNZPSW}.  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 \cite{Gigli18}.

The noncollapsed hypothesis implies that the tangent module has
dimension \(n\) almost everywhere, and that
\(\tr\Hess w=\Delta w\) there.  The dimensional Bochner inequality
therefore implies \(M(w)\ge|\Hess w|^2\m\); see
\cite{DPG18,Han18}.  Expanding the square
\[
 |\Hess w-bg|^2=|\Hess w|^2-2b\Delta w+nb^2
\]
proves \eqref{co:shifted-bochner}.  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
\cite[Lemma~3.1]{MN19} and \cite[Lemma~2.6]{CDNZPSW}.
To obtain the zero-boundary coercivity in
\eqref{co:dirichlet-sobolev}, 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
\cite[Theorem~5.1 and the following Remark~3]{HK00}. The annulus lower
bound and all constants are uniform for the almost Euclidean balls
used below.
\end{proof}

For a column \(u=(u_1,\ldots,u_n)^T\) of harmonic functions, define
\begin{equation}\label{co:G-M-matrix}
 G=(\Gamma(u_j,u_k))_{j,k=1}^n,
 \qquad M=(M(u_j,u_k))_{j,k=1}^n,
\end{equation}
where \(M(\cdot,\cdot)\) denotes polarization of
\eqref{co:M-definition}.  Thus \(M[b]=M(b\cdot u)\) is a positive
measure for every \(b\in\R^n\).  If \(E\) is Borel, positivity
implies
\begin{equation}\label{co:matrix-measure-cs}
 |M(v,w)|(E)\le \sqrt{M[v](E)M[w](E)}.
\end{equation}
Here \(|M(v,w)|\) is total variation.  To see this, take the trace
measure as a dominating measure and apply matrix Cauchy--Schwarz
to the positive semidefinite matrix density.  The same argument
applies to integrals with a nonnegative weight.

We record precisely the rescaling regime used throughout the proof.
\(\omega_n\) denotes the measure of the unit ball in \(\R^n\).

\begin{lemma}[An expanding almost Euclidean regime]\label{co:flat-setup}
Fix an \(n\)-regular point \(p\) in the space of
Theorem~\ref{thm:main}.  There are rescalings
\[
 (X_\nu,d_\nu,\m_\nu,p_\nu)
 =(X,s_\nu^{-1}d,s_\nu^{-n}\mathcal H_d^n,p),
 \qquad s_\nu\downarrow0,
\]
numbers \(\kappa_\nu\downarrow0\), and radii \(T_\nu\uparrow\infty\)
with \(\kappa_\nu T_\nu^2\to0\), having the following properties.
The space \(X_\nu\) is \(\RCD(-\kappa_\nu,n)\).  For a sequence
\(\delta_\nu\downarrow0\), all balls with center in
\(B_{64T_\nu}(p_\nu)\) and radius at most \(64T_\nu\) satisfy
\begin{equation}\label{co:uniform-flatness}
 \left|\frac{\m_\nu(B_r(y))}{\omega_n r^n}-1\right|
 \le\delta_\nu,
 \qquad
 d_{GH}\bigl(\overline B_r(y),\overline B_r(0^n)\bigr)
 \le\delta_\nu r.
\end{equation}
There are harmonic columns \(u_\nu\) on
\(B_{32T_\nu}(p_\nu)\), converging on every fixed ball locally
uniformly and strongly in \(W^{1,2}\) to Euclidean coordinates, and
\begin{equation}\label{co:original-gradient-bound}
 |\nabla u_\nu|\le C_n
 \quad\text{on }B_{16T_\nu}(p_\nu).
\end{equation}
\end{lemma}

\begin{proof}
Take rescalings tending to zero.  Compactness and the regular-point
hypothesis imply pointed convergence to \(\R^n\), after passage to
a subsequence, and noncollapsed volume convergence identifies the
limit measure with \(\mathcal H^n\); see \cite{DPG18}.  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
\eqref{co:uniform-flatness}, uniformly over all smaller radii.
The almost-volume-rigidity theorem in \cite{DPG18}, 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 \cite[Corollary~4.12]{AH18}.  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
\eqref{co:original-gradient-bound} and the expanding domains of
harmonicity simultaneously.
\end{proof}

From now on the sequence index is suppressed.  Every assertion on
a fixed ball, or with finitely many fixed enlargement factors, is
understood to hold for all sufficiently large \(\nu\).  A quantity
denoted by \(o(1)\) is uniform in the centers in that ball and in
all the smaller scales.  We set
\begin{equation}\label{co:dyadic-R}
 r_i=2^{-i},\qquad
 R_i(y)=\fint_{B_{r_i}(y)}G\dd\m,
 \qquad i=0,1,\ldots.
\end{equation}

\subsection{Transformation at every center and scale}

Harmonic normalization across scales is a central tool in
\cite[Theorem~1.32]{CN15}, \cite[Section~7]{CJN21}, and
\cite[Theorem~4.1]{HH24}. The noncollapsed RCD transformation framework of
Bru\`e--Naber--Semola \cite[Section~3.3, Proposition~3.13 and Corollary~3.16]{BNS20}
is a further methodological antecedent. We give the all-center statement and
its proof needed on the expanding domains used here.

\begin{lemma}[Uniform harmonic normalization]\label{co:harmonic-normalization}
For every pair of fixed constants \(A,C\ge1\), the matrices \(R_i(y)\) are
positive definite for \(y\in B_A(p)\) and \(i\ge0\).  In the
rescaled ball \(B_{Cr_i}(y)\), the map
\begin{equation}\label{co:normalized-column}
 U_{i,y}=\frac{R_i(y)^{-1/2}(u-u(y))}{r_i}
\end{equation}
is, after an orthogonal change of target coordinates, an
\(o(1)\) approximation to Euclidean coordinates, uniformly and
strongly in interior \(W^{1,2}\).  Its gradients are uniformly
bounded on each fixed interior enlargement.  Its Gram matrix
converges strongly in every finite \(L^p\) on such enlargements
to \(\Id\).

For every fixed \(C\ge1\), if \(d(y,z)\le Cr_i\) and
\(C^{-1}r_i\le r_j\le Cr_i\), then
\begin{equation}\label{co:nearby-R}
 (1-o(1))R_i(y)\le R_j(z)\le(1+o(1))R_i(y).
\end{equation}
For every \(a>0\), after choosing a sufficiently late rescaling,
\begin{equation}\label{co:slow-R}
 C_a^{-1}2^{-2a|i-j|}R_j(y)
 \le R_i(y)\le C_a2^{2a|i-j|}R_j(y).
\end{equation}
The constants are uniform in the center and the two indices.  In
particular, since \(R_0(y)=\Id+o(1)\) on fixed balls,
\begin{equation}\label{co:R-lower}
 c_a r_i^{2a}\Id\le R_i(y)\le C\Id.
\end{equation}
\end{lemma}

\begin{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)-\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 \eqref{co:uniform-flatness}.  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 \cite[Section~4]{AH18}; 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~\ref{co:flat-setup}, 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 \(\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 \(\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
\eqref{co:nearby-R} proves that their Gram averages in the first
normalization both tend to \(\Id\).  This gives
\eqref{co:nearby-R}.  Iterating its same-center adjacent-scale
case, with error chosen so that each factor is at most \(2^{2a}\),
gives \eqref{co:slow-R}.  The upper bound in
\eqref{co:R-lower} also follows directly from
\eqref{co:original-gradient-bound}.

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
\(\Id\); see \cite[Corollary~4.3 and Theorem~4.6]{AH18}.  Their
uniform boundedness upgrades the resulting strong \(L^2\)
convergence to every finite \(L^p\).
\end{proof}

\begin{lemma}[Small normalized Bochner mass]\label{co:bochner-small}
For each pair of fixed constants \(A,C\ge1\), uniformly for \(y\in B_A(p)\) and
\(i\ge0\),
\begin{equation}\label{co:small-mass}
 r_i^{2-n}M(B_{Cr_i}(y))\le o(1)R_i(y).
\end{equation}
\end{lemma}

\begin{proof}
In units \(r_i\), apply Lemma~\ref{co:local-calculus} to a scalar
combination of \(U=U_{i,y}\).  The rescaled Bochner matrix is
\begin{equation}\label{co:M-scaling}
 M_{r_i}(U)=r_i^{2-n}R_i(y)^{-1/2}M(u)R_i(y)^{-1/2},
\end{equation}
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\dd M_{r_i}(b\cdot U)
 =\frac12\int\bigl(\Gamma(b\cdot U)-|b|^2\bigr)
                  \Delta\phi\dd\m_{r_i}
  +\kappa r_i^2\int\phi\Gamma(b\cdot U)\dd\m_{r_i}.
\]
Here \(\int\Delta\phi\dd\m_{r_i}=0\).  Both terms tend uniformly
to zero by Lemma~\ref{co:harmonic-normalization}.  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
\eqref{co:small-mass}.
\end{proof}

\subsection{Open coordinate images}

\begin{lemma}[Open harmonic coordinates]\label{co:open-coordinates}
For every fixed \(A\ge1\), the map \(u\) is a homeomorphism
from \(B_A(p)\) onto an open subset of \(\R^n\), once the
rescaling is sufficiently late.  There is a fixed \(c>0\) such
that, whenever the balls involved stay in a fixed interior
working region,
\begin{equation}\label{co:ellipsoid-containment}
 \{P:|P-u(x)|_{R_i(x)^{-1}}<c r_i\}
 \subset u(B_{r_i/2}(x)).
\end{equation}
Moreover, for every fixed \(D\ge1\), if
\(D^{-1}r_i\le d(x,z)\le Dr_i\), then
\begin{equation}\label{co:separation-scale}
 \bigl||u(z)-u(x)|_{R_i(x)^{-1}}-d(x,z)\bigr|
 \le o_D(1)r_i.
\end{equation}
\end{lemma}

\begin{proof}
The approximation in Lemma~\ref{co:harmonic-normalization} gives
\eqref{co:separation-scale}.  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 \eqref{co:ellipsoid-containment} 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 \eqref{co:nearby-R}, 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\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.
\end{proof}

\section{Heat averages and cumulative matrix control}
\label{co:heat-section}

We next replace the ball Gram averages by matrices that are
monotone as the spatial scale decreases.  The local mass forms use
inverse ball averages at nearby centers.  Relative comparability
of these averages with the heat averages does not by itself compare
the twice-conjugated mass forms.
We will bound the accumulated cost of changing the normalization by
a small multiple of the final inverse heat matrix.
Heat Gram normalization and quadratic telescoping are used in
\cite[Sections~8.1--8.2]{CJN21}. The estimates below give the cumulative
comparison for the inverse matrices without a commutativity
assumption.

\subsection{Continuous kernels and compactly supported sources}

Let \(p_t(x,z)\) be the common continuous heat kernel and \(P_t\)
the heat semigroup.  We use the Gaussian and first spatial
gradient estimates of \cite[Theorem~1.2 and Corollary~1.2]{JLZ16}.
On the fixed testing regions, for \(0<t\le1\), volume comparison
allows their form
\begin{equation}\label{co:kernel-bounds}
 p_t(x,z)\le C t^{-n/2}e^{-d(x,z)^2/(Ct)},
 \qquad
 |\nabla_xp_t(x,z)|\le C t^{-(n+1)/2}e^{-d(x,z)^2/(Ct)}.
\end{equation}
The gradient statement initially holds almost everywhere.  The
lower bound, used when \(d(x,z)\le A\sqrt t\), is
\begin{equation}\label{co:kernel-lower}
 p_t(x,z)\ge c_A t^{-n/2}.
\end{equation}
We may use one fixed weaker negative curvature bound in these
estimates.  The constants are then uniform in the sequence.

The following interpretation will be useful because \(M\) can
have a singular part.

\begin{lemma}[Kernel differences and measure sources]\label{co:kernel-measures}
Fix \(C_0\ge1\), \(x,y\) in a fixed testing region, and
\(d(x,y)\le C_0r\), where \(0<r\le1\).  Put
\[
 K_t(z)=\fint_{B_r(y)}p_t(a,z)\dd\m(a),
 \qquad D=d(x,z).
\]
After enlarging fixed constants, the following bounds hold at
every source pole \(z\):
\begin{align}
 K_t(z)&\le C r^{-n}e^{-cD^2/r^2},&&0<t\le r^2,
       \label{co:short-kernel}\\
 |K_t(z)-p_t(x,z)|&\le C r t^{-(n+1)/2}e^{-cD^2/t},
       &&r^2\le t\le1.\label{co:long-kernel}
\end{align}
If \(f\in W^{1,2}\cap L^\infty\) has compact support and its
Laplacian is a finite signed measure \(\mu\) of compact support,
then
\begin{equation}\label{co:semigroup-measure}
 P_t f-P_s f=\int_s^t P_\tau\mu\dd\tau,
 \qquad 0<s<t,
\end{equation}
using the continuous representatives.  In the application
\(f=\chi\Gamma(u_j,u_k)\), this identity can also be averaged over
a ball and passed to \(s\downarrow0\).
\end{lemma}

\begin{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 \eqref{co:short-kernel}.

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
\eqref{co:long-kernel}, 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 \eqref{co:semigroup-measure}, 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: \eqref{co:short-kernel},
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.
\end{proof}

Choose a nonnegative cutoff \(\chi\) equal to one on \(B_T(p)\),
supported in \(B_{2T}(p)\), with
\[
 |\nabla\chi|\le C/T,\qquad |\Delta\chi|\le C/T^2,
\]
where \(T=T_\nu\) comes from Lemma~\ref{co:flat-setup}.  Extend
\(\chi G\) by zero outside its support.  Its distributional
Laplacian is the compactly supported matrix measure
\begin{equation}\label{co:cutoff-source}
 \Delta(\chi G)=2\chi M-2\kappa\chi G\m+E_{\rm cut}\m,
 \qquad E_{\rm cut}=(\Delta\chi)G+2\Gamma(\chi,G).
\end{equation}
The last expression is a matrix of \(L^1\) functions supported in
the cutoff annulus.

\begin{lemma}[Monotone heat Gram matrices]\label{co:heat-matrices}
There is \(\lambda=\lambda_\nu>0\), with \(\lambda_\nu\to0\),
such that, at all fixed testing centers and \(0<t\le1\),
\begin{equation}\label{co:heat-derivative}
 P_t(\chi M)
 \le\partial_t\bigl(P_t(\chi G)+\lambda t\Id\bigr)
 \le 3P_t(\chi M)+C\Id.
\end{equation}
Define
\begin{equation}\label{co:QH-definition}
 Q_i(x)=P_{r_i^2}(\chi G)(x)+\lambda r_i^2\Id,
 \qquad H_i(x)=Q_i(x)^{-1}.
\end{equation}
Then
\begin{align}
 &(1-o(1))R_i(x)\le Q_i(x)\le(1+o(1))R_i(x),
       \qquad Q_0(x)=\Id+o(1),\label{co:QR-comparison}\\
 &H_i(x)\le H_{i+1}(x)\le(1+o(1))H_i(x).
       \label{co:H-monotone}
\end{align}
The nearby-center and slow-growth estimates of
Lemma~\ref{co:harmonic-normalization} hold for \(Q_i\) as well.
In particular, for every prescribed \(a>0\), with a sufficiently
late rescaling,
\begin{equation}\label{co:slow-QH}
 Q_j(x)\le C_a2^{2a(i-j)}Q_i(x)\quad(j\le i),
 \qquad \|H_i(x)\|\le C_a2^{2ai}.
\end{equation}
\end{lemma}

\begin{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,
\begin{equation}\label{co:raw-mass-growth}
 M(B_s(x))[v]\le C(s^{n-2}+\kappa s^n)|v|^2
 \le Cs^{n-2}|v|^2.
\end{equation}
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 \eqref{co:original-gradient-bound} holds.  The same estimate
for the Hessian, followed by Cauchy--Schwarz in the product rule,
gives
\[
 \int|E_{\rm cut}|\dd\m\le CT^{n-2}.
\]
Since its support is a distance at least \(T/2\) from the fixed
testing centers for large \(\nu\), \eqref{co:kernel-bounds} yields
\begin{equation}\label{co:cutoff-small}
 \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).
\end{equation}
Choose \(\lambda\) at least this supremum plus
\(2\kappa\|G\|_\infty\), and add a positive number tending to
zero.  Equations \eqref{co:semigroup-measure} and
\eqref{co:cutoff-source} give
\[
 \partial_t(P_t(\chi G)+\lambda t\Id)
 =2P_t(\chi M)-2\kappa P_t(\chi G)
     +P_tE_{\rm cut}+\lambda\Id.
\]
The last three terms have nonnegative sum and bounded norm.  This
proves \eqref{co:heat-derivative}, 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 \(\Id\), by
Lemma~\ref{co:harmonic-normalization}.  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
\eqref{co:R-lower}; 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}\Id=o(1)\Id
\]
for \(a<1\).  This proves \eqref{co:QR-comparison}; 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
\eqref{co:nearby-R} give the adjacent-scale upper bound and all
fixed nearby comparisons.  Iteration gives
\eqref{co:slow-QH}, choosing the adjacent-scale tolerance after
the prescribed exponent \(a\).
\end{proof}

\subsection{The dyadic heat comparison}

The matrices
\begin{equation}\label{co:N-definition}
 N_j(x)=r_j^{2-n}M(B_{r_j}(x))
\end{equation}
are masses in the original harmonic coordinate components.  They
have not been conjugated by a Gram inverse.  In particular,
Lemma~\ref{co:bochner-small} and \eqref{co:QR-comparison} give
\begin{equation}\label{co:N-small}
 N_j(x)\le\varepsilon_\nu Q_j(x),
 \qquad \varepsilon_\nu\longrightarrow0,
\end{equation}
uniformly in \(j\) and in fixed testing regions.  We allow
\(\varepsilon_\nu\) to increase by fixed multiplicative constants
in what follows.

\begin{lemma}[Dyadic heat comparison]\label{co:heat-comparison}
Fix constants \(A,C_0\ge1\), where \(C_0\) bounds
\(d(x,y)/r_i\), and fix a bounded region for \(x\).
There is a fixed exponent \(\gamma>0\), which
may be taken to be \(1/2\), with the following properties for all
sufficiently late rescalings.  Set
\begin{equation}\label{co:F-definition}
 \mathcal F_i(x)=2^{-\gamma i}\Id+
       \sum_{j=0}^i2^{-\gamma(i-j)}N_j(x).
\end{equation}
For the prescribed enlargement \(A\),
\begin{equation}\label{co:heat-decrement}
 c_A r_i^{2-n}M(B_{Ar_i}(x))
 \le Q_i(x)-Q_{i+1}(x)\le C\mathcal F_i(x).
\end{equation}
For \(d(x,y)\le C_0r_i\) and \(v,w\in\R^n\),
\begin{equation}\label{co:bilinear-heat-comparison}
 |v^T(R_i(y)-Q_i(x))w|
 \le C|v|_{Q_i(x)}\sqrt{\mathcal F_i(x)[w]}.
\end{equation}
All constants are independent of \(i\) and of the expanding
cutoff radius.  They may depend on the stated fixed parameters.
\end{lemma}

\begin{proof}
\emph{The semigroup identity.}
Write \(r=r_i\), suppress \(x\) in \(N_j,Q_j\), and let
\(\mathcal A f=\fint_{B_r(y)}f\dd\m\).  We first prove the more
explicit estimate
\begin{equation}\label{co:raw-bilinear-heat}
 \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}
\end{equation}
This also specifies all the spatial and time tails entering
\eqref{co:bilinear-heat-comparison}.

Put \(f=\chi G\), \(S=\Delta f\), and use the kernel functions
from Lemma~\ref{co:kernel-measures}.  Since \(\chi=1\) on the
averaging ball, the exact semigroup identity is
\begin{equation}\label{co:average-heat-identity}
 \begin{split}
 R_i(y)-Q_i(x)
 ={}&\mathcal A P_1f-P_1f(x)-\lambda r^2\Id\\
 &-\int_0^{r^2}\mathcal A P_tS\dd t
  -\int_{r^2}^1\bigl(\mathcal A P_tS-P_tS(x)\bigr)\dd t.
 \end{split}
\end{equation}
The averaged endpoint at zero is justified in
Lemma~\ref{co:kernel-measures}; it does not require pointwise
convergence of the Gram field.

\smallskip\noindent
\emph{The local source.}
Consider the contribution of \(2\chi M\) on \(B_1(x)\).
The integrated short-time kernel obeys
\begin{equation}\label{co:integrated-short-kernel}
 \int_0^{r^2}K_t(z)\dd t
 \le Cr^{2-n}e^{-c d(x,z)^2/r^2}.
\end{equation}
For \(0\le q<i\), integration of
\eqref{co:long-kernel} on the time block
\(I_q=[r_{q+1}^2,r_q^2]\) gives
\begin{equation}\label{co:integrated-long-kernel}
 \int_{I_q}|K_t(z)-p_t(x,z)|\dd t
 \le Cr r_q^{1-n}e^{-c d(x,z)^2/r_q^2}.
\end{equation}
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
\begin{equation}\label{co:gaussian-shell-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.
\end{equation}
On the annulus indexed by \(j<q\), the distance from \(x\) is at
least \(2^{q-j-1}r_q\).  Equation
\eqref{co:matrix-measure-cs} 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
\eqref{co:integrated-short-kernel}.  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.

\smallskip\noindent
\emph{The remaining source terms.}
We treat all remaining parts of
\eqref{co:average-heat-identity}.  If \(D=d(x,z)\ge1\), the
integrated kernel bounds imply
\begin{equation}\label{co:exterior-integrated-kernel}
 \int_0^{r^2}K_t(z)\dd t+
 \int_{r^2}^1|K_t(z)-p_t(x,z)|\dd t
 \le Cr e^{-c'D^2}.
\end{equation}
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,
\eqref{co:integrated-short-kernel} reduces the assertion to
absorbing \(r^{1-n}\) in \(e^{-c''/r^2}\).  Both bounds are uniform
for \(0<r\le1\).

Equation \eqref{co:raw-mass-growth}, summed over the annuli
\(B_{2^{h+1}}(x)\setminus B_{2^h}(x)\), gives
\[
 \int e^{-c'd(x,z)^2}\dd(\chi M[v])(z)\le C|v|^2.
\]
Therefore \eqref{co:exterior-integrated-kernel} 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 \eqref{co:cutoff-source} 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
\eqref{co:long-kernel} 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}\dd t\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
\eqref{co:raw-bilinear-heat}.

\smallskip\noindent
\emph{Normalization and decrements.}
Choose any sufficiently small fixed growth exponent \(a>0\) in
\eqref{co:slow-QH}.  By \eqref{co:N-small}, 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}\Id\), the raw term in
\eqref{co:raw-bilinear-heat} 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 \eqref{co:bilinear-heat-comparison} with
\(\gamma=1/2\).

It remains to prove the decrement estimates.  Integrate
\eqref{co:heat-derivative} on
\([r_i^2/4,r_i^2]\).  On a fixed ball \(B_{Ar_i}(x)\), where
\(\chi=1\), the heat lower bound
\eqref{co:kernel-lower} gives the left inequality in
\eqref{co:heat-decrement}.  The upper bound is
\[
 Q_i-Q_{i+1}
 \le Cr_i^{2-n}\int e^{-c d(x,z)^2/r_i^2}\dd(\chi M)(z)
        +Cr_i^2\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)\),
\eqref{co:raw-mass-growth} and Gaussian decay give \(Cr_i^2\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\).
\end{proof}

\begin{remark}[Fixed enlargements and indexing]\label{co:enlargement-indexing}
No enlarged mass was needed in the preceding bilinear estimate.
When one is used elsewhere, it can be returned to the same
notation without losing the decay.  Indeed, for fixed \(A\ge1\)
put \(h=\lceil\log_2 A\rceil\).  For \(j\ge h\),
\[
 r_j^{2-n}M(B_{Ar_j}(x))\le2^{h(n-2)}N_{j-h}(x).
\]
The finitely many \(j<h\) are bounded by \(C_A\Id\) in the
original coordinate components, by a fixed-scale cutoff estimate.
Thus an exponentially weighted sum of these enlarged masses is
bounded by a fixed multiple of
\(2^{-\gamma i}\Id+\sum_{j\le i}2^{-\gamma(i-j)}N_j\).
A center shift bounded by \(C_0r_j\) only increases \(A\) by
\(C_0\).  In contrast, the mass matrices \(D_{i,y}\) below include
two Gram inverses; their comparison requires the argument that
follows.
\end{remark}

\subsection{Accumulation without commuting matrices}

We separate the finite-dimensional calculation from its geometric
application.  All inequalities in this subsection are inequalities
of quadratic forms on the fixed space \(\R^n\).

\begin{lemma}[Weighted inverse increments]\label{co:inverse-increments}
Fix a testing center \(x\), and suppress it in the notation.
Let \(a>0\) in \eqref{co:slow-QH} be chosen so that
\(4a<\gamma\).  Set
\[
 B_i=H_{i+1}-H_i,\qquad \Delta_i=Q_i-Q_{i+1},
 \qquad a_i=H_i\xi
\]
for a fixed \(\xi\in\R^n\).  Then, uniformly for every integer
\(k\ge0\),
\begin{align}
 &H_i\Delta_iH_i\le B_i\le(1+C\varepsilon_\nu)H_i\Delta_iH_i,
             \label{co:inverse-increment-order}\\
 &\sum_{i<k}\mathcal F_i[a_i]\le C H_k[\xi],
             \label{co:weighted-total}\\
 &C^{-1}H_k\le \Id+\sum_{i<k}H_iN_iH_i\le C H_k.
             \label{co:central-cumulative}
\end{align}
For every fixed \(A\ge1\), the upper sum bound also holds with
\(r_i^{2-n}M(B_{Ar_i}(x))\) in place of \(N_i\).
\end{lemma}

\begin{proof}
Let \(E_i=Q_i^{-1/2}\Delta_iQ_i^{-1/2}\).
By \eqref{co:H-monotone}, \(0\le E_i\le\varepsilon_\nu\Id\)
after increasing the sequence \(\varepsilon_\nu\to0\).
The exact inverse identity is
\[
 B_i=Q_i^{-1/2}\bigl((\Id-E_i)^{-1}-\Id\bigr)Q_i^{-1/2}.
\]
The spectral inequality
\(E_i\le(\Id-E_i)^{-1}-\Id\le(1+C\varepsilon_\nu)E_i\)
proves \eqref{co:inverse-increment-order}.  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
\eqref{co:slow-QH},
\begin{equation}\label{co:B-vector-bound}
 |B_l\xi|_{Q_j}^2\le C2^{2a(l-j)}B_l[\xi],\qquad j\le l.
\end{equation}
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]
       \label{co:convolution}\\[-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]\notag\\
 &\qquad\le C\sum_{j<k}N_j[a_j]
             +C\varepsilon_\nu H_k[\xi].\notag
\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
\eqref{co:slow-QH}:
\begin{equation}\label{co:background-sum}
 \sum_{i<k}2^{-\gamma i}|a_i|^2
 \le C|\xi|^2\sum_{i\ge0}2^{-(\gamma-4a)i}
 \le C|\xi|^2.
\end{equation}

The lower heat decrement with \(A=1\), together with
\eqref{co:inverse-increment-order}, 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 \eqref{co:central-cumulative}, since \(H_0\) is comparable
with \(\Id\).  Equations \eqref{co:convolution} and
\eqref{co:background-sum} now imply
\eqref{co:weighted-total}.  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.
\end{proof}

\begin{lemma}[Changing the mean normalization]\label{co:normalization-error}
Fix \(C_0\ge1\).  For \(d(x,y)\le C_0r_i\), put
\[
 A_{i,y}=R_i(y)-Q_i(x),\qquad
 e_{i,y}=R_i(y)^{-1}\xi-H_i(x)\xi.
\]
Then
\begin{equation}\label{co:normalization-error-pointwise}
 |e_{i,y}|_{Q_i(x)}^2
 \le C\mathcal F_i(x)[H_i(x)\xi].
\end{equation}
If at each \(i<k\) at most \(J\) such centers are chosen, then
\begin{equation}\label{co:normalization-error-sum}
 \sum_{i<k}\sum_y|e_{i,y}|_{Q_i(x)}^2
 \le C_J H_k(x)[\xi].
\end{equation}
For any positive semidefinite matrices \(C_{i,y}\) satisfying
\(C_{i,y}\le\varepsilon_\nu Q_i(x)\), their error cost is
\begin{equation}\label{co:normalization-mass-error}
 \sum_{i<k}\sum_y C_{i,y}[e_{i,y}]
 \le C_J\varepsilon_\nu H_k(x)[\xi].
\end{equation}
\end{lemma}

\begin{proof}
Taking the supremum over \(v\) with \(|v|_{Q_i}=1\) in
\eqref{co:bilinear-heat-comparison} gives
\begin{equation}\label{co:bilinear-dual-form}
 A_{i,y}H_iA_{i,y}\le C\mathcal F_i.
\end{equation}
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 \eqref{co:bilinear-dual-form} 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 \eqref{co:normalization-error-pointwise}.  Sum it,
use the bound \(J\), and apply
\eqref{co:weighted-total} to obtain
\eqref{co:normalization-error-sum}.  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.
\end{proof}

\begin{lemma}[Cumulative comparison for local mass matrices]
\label{co:cumulative-comparison}
Fix constants \(A_0,C_0,L\ge1\) and a positive integer \(J\).
At each center
\(x\in B_{A_0}(p)\) and each scale \(i\ge0\), choose at most \(J\)
centers \(y\) with \(d(x,y)\le C_0r_i\).  The choice may depend
arbitrarily on \(x,i\).  Set
\begin{equation}\label{co:D-definition}
 C_{i,y}=r_i^{2-n}M(B_{Lr_i}(y)),
 \qquad D_{i,y}=R_i(y)^{-1}C_{i,y}R_i(y)^{-1}.
\end{equation}
Give the chosen centers weights \(\theta_{i,y}\in[0,1]\), and
suppose that for every \(i\) at least one has weight one and
satisfies \(B_{Lr_i}(y)\supset B_{r_i}(x)\).  Write
\begin{equation}\label{co:ZS-definition}
 Z_i=\sum_yD_{i,y},\qquad
 S_i=\sum_y\theta_{i,y}D_{i,y}.
\end{equation}
There is \(C_*<\infty\), depending only on \(n,A_0,C_0,L,J\),
such that for all sufficiently late rescalings, every such choice,
and every integer \(k\ge0\),
\begin{equation}\label{co:cumulative-estimates}
 Z_i\le o(1)H_i(x),\qquad
 \Id+\sum_{i<k}Z_i\le C_*H_k(x),\qquad
 H_k(x)\le C_*\left(\Id+\sum_{i<k}S_i\right).
\end{equation}
The \(o(1)\) is uniform over the stated centers, all scales, and
all choices satisfying these conditions.
\end{lemma}

\begin{proof}
Fix \(x\) and \(\xi\), and suppress \(x\) from the notation.
By \eqref{co:small-mass} 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
\eqref{co:cumulative-estimates}.

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~\ref{co:inverse-increments} therefore bounds the sum of
its first terms by \(C_J H_k[\xi]\).
Equation \eqref{co:normalization-mass-error} bounds the sum of
its second terms by \(C_J\varepsilon_\nu H_k[\xi]\).
Since \(\Id\le C H_0\le C H_k\), this proves the middle
assertion in \eqref{co:cumulative-estimates}.

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 \eqref{co:central-cumulative} 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.
\end{proof}
==== END Bi-Lipschitz-Coordinates-at-Regular-Points-of-Noncollapsed-RCD-Spaces//coordinates.tex ====

==== BEGIN Bi-Lipschitz-Coordinates-at-Regular-Points-of-Noncollapsed-RCD-Spaces//introduction.tex ====
\section{Introduction}

Euclidean tangent spaces describe the geometry of a metric space at a
point, but a coordinate chart must compare all pairs of points in a
neighborhood. Under a lower Ricci curvature bound these are distinct
regularity questions.

Throughout, noncollapsed means that the reference measure is exactly
\(\mathcal H_d^n\). We use the full \(\CD(K,n)\) condition together with
quadratic Cheeger energy. Spaces are complete, separable, geodesic, and
have full-support reference measure finite on bounded sets. A point is
\emph{\(n\)-regular} if every pointed Gromov--Hausdorff tangent there is
pointed-isometric to \((\R^n,|\cdot|,0)\). Here a tangent is a pointed
limit of rescalings by inverse radii tending to infinity.
We prove the following open-neighborhood result.

\begin{theorem}\label{thm:main}
For every integer \(n\ge2\) there is a constant \(L_n\ge1\) with
the following property. Let \(K\in\R\), and let
\((X,d,\mathcal H_d^n)\) be a noncollapsed \(\RCD(K,n)\) space.
Let \(p\in X\).
If every pointed Gromov--Hausdorff tangent of \(X\) at \(p\) is
pointed-isometric to \((\R^n,|\cdot|,0)\), then there exist an open
neighborhood \(U\) of \(p\), an open set \(V\subset\R^n\), and a
homeomorphism \(F:U\to V\) such that
\[
 L_n^{-1}d(x,z)\le |F(x)-F(z)|\le L_nd(x,z)
 \qquad (x,z\in U).
\]
The distance on the left and right is the original distance \(d\)
restricted to \(U\).
\end{theorem}

The neighborhood may depend on \(X\) and \(p\); the distortion
constant depends only on the dimension. The theorem concerns every
prescribed regular point, not only almost every point of the reference
measure. It also yields the following manifold region.

\begin{corollary}\label{cor:full-measure-atlas}
Under the hypotheses of Theorem~\ref{thm:main}, there is an open set
\(\mathcal U\subset X\) containing every \(n\)-regular point, with
\(\mathcal H_d^n(X\setminus\mathcal U)=0\), and a countable atlas
of \(L_n\)-bi-Lipschitz charts from ambient-open subsets of
\(\mathcal U\) onto open subsets of \(\R^n\). Each chart uses the
restricted ambient distance, and its transition maps are
\(L_n^2\)-bi-Lipschitz.
\end{corollary}

\begin{proof}
The \(n\)-regular set has full \(\mathcal H_d^n\) measure by
\cite[Theorem~1.12(ii),(vi)]{DPG18}, applied with reference measure
\(\mathcal H_d^n\). Take \(\mathcal U\) to be the union of the chart
domains in Theorem~\ref{thm:main}. As a separable metric space,
\(\mathcal U\) is second countable, so these domains admit a
countable subcover. On an overlap, composing one chart with the
inverse of the other gives the stated transition bound.
\end{proof}

The open region \(\mathcal U\) may contain points outside the exact
regular locus; the corollary does not assert that this locus is open.

\subsection{History and significance}

The structure theory of noncollapsed limits of Riemannian manifolds
with lower Ricci bounds provides an early form of this question.
Cheeger and Colding obtained manifold neighborhoods with bi-H\"older
parametrizations and asked about a bi-Lipschitz upgrade
\cite[Theorem~5.14 and Remark~5.15]{CC97}. Their question concerns
quantitatively regular Ricci-limit neighborhoods and is a historical
antecedent of the synthetic regular-point problem considered here.

Synthetic lower Ricci bounds were developed through optimal transport
by Lott and Villani and by Sturm \cite{LV09,Sturm06II}. The
Riemannian refinement of Ambrosio, Gigli and Savar\'e distinguishes
quadratic Cheeger energy and linear heat flow \cite{AGS14}; the
finite-dimensional connections with Bochner inequalities were
developed by Erbar, Kuwada and Sturm \cite{EKS15}. In the noncollapsed
setting, De Philippis and Gigli established volume convergence,
volume rigidity, and quantitative Euclidean approximation
\cite{DPG18}.

Mondino and Naber proved countable bi-Lipschitz rectifiability of
RCD spaces: measurable pieces cover almost every point and map to
measurable Euclidean subsets \cite[Theorem~1.1]{MN19}. This does not
give an ambient-open neighborhood of each prescribed regular point.
Kapovitch and Mondino proved that the entire regular set is contained
in an open manifold region with bi-H\"older charts
\cite[Theorem~3.1 and Corollary~3.2]{KM21}. This statement does not
require the exact regular locus itself to be open. An open
bi-Lipschitz chart requires stronger control across infinitely many
scales.

The conjecture that every regular point in the noncollapsed setting
admits a bi-Lipschitz chart is explicitly recalled in the recent work
of Honda and Zhang \cite[Section~1.2]{HZ26}. Theorem~\ref{thm:main}
answers this question under the precise noncollapsed \(\RCD(K,n)\)
hypotheses stated above. Honda and Zhang's geometric regularity
theorems assume a positive lower bound for injectivity radius on
neighborhoods \cite[Theorems~1.1 and 3.35]{HZ26}. Their analytic
criterion also permits Morrey \(L^2\) control of the Laplacians,
together with a neighborhood-uniform positive Gram determinant
\cite[Theorem~3.9 and Proposition~3.10]{HZ26}. Neither additional
hypothesis is assumed here.

Several nearby results clarify the strength of the conclusion.
Colding and Naber proved that bounded subsets of quantitative
Reifenberg loci in noncollapsed Ricci limits, with sufficiently small
flatness and scale parameters, admit bi-Lipschitz embeddings into a
Euclidean space of possibly higher dimension
\cite[Theorem~1.1]{CN13}. They also constructed examples with
Euclidean tangents that admit no coordinates inducing a H\"older
continuous metric tensor \cite[Theorem~1.2]{CN13}. The conclusion here
is an open chart into \(\R^n\); it does not assert H\"older continuity
of the metric tensor in that chart. Likewise, the restrictions on
bi-Lipschitz regular maps in \cite[Proposition~5.1]{HP23} concern
maps whose pullbacks of smooth target test functions have locally
bounded Laplacians \cite[Definition~3.1]{HP23}. No such analytic
regularity is asserted for the limiting coordinates constructed
below; in particular, their components are not asserted to have
bounded Laplacian.

\subsection{The two main estimates}

After rescaling at \(p\), volume comparison makes all sufficiently
small balls at nearby centers quantitatively Euclidean. Harmonic
functions approximating Euclidean coordinates then give an open
coordinate map \(u=(u_1,\ldots,u_n)^T\). The scale-dependent mean Gram matrices of \(u\)
are invertible, but they need not be uniformly nondegenerate as the
scale tends to zero. We compensate for this degeneration in the
\emph{space of coordinate values}.

Two kinds of control make this possible.  Accumulated Bochner mass
records how much the coordinates must be stretched to compensate for
their degeneration; this accumulation may be unbounded.  A separate
estimate evaluates the same mass forms on actual coordinate differences
$u(x)-u(z)$ and controls them by defects whose sum is small.  The
construction uses the cumulative mass to prescribe the stretching and
the directional estimate to keep its metric error summable.

Write \(r_i=2^{-i}\), and let
\(G=(\langle\nabla u_j,\nabla u_k\rangle)_{j,k=1}^n\) be the Gram
matrix field. Its mean on \(B_{r_i}(y)\) is denoted by \(R_i(y)\).
Let \(\mathfrak m\) denote Hausdorff measure in the current length
units. If the rescaled curvature lower bound is \(-\kappa\), the
matrix-valued Bochner measure is
\(M=\tfrac12\Delta G+\kappa G\mathfrak m\), with the Laplacian
understood entrywise. For a fixed enlargement \(L>1\), define the
positive semidefinite matrix
\[
 D_{i,y}=R_i(y)^{-1}
       \bigl[r_i^{2-n}M(B_{Lr_i}(y))\bigr]R_i(y)^{-1},
 \qquad D[\xi]=\xi^TD\xi.
\]
Corrected heat averages of \(G\) yield positive definite forms
\(H_i(x)\), increasing with \(i\), on the space of coordinate
values. Lemma~\ref{co:cumulative-comparison} compares these forms
with accumulated bounded-overlap sums of the local matrices
\(D_{i,y}\), including a reverse comparison for suitable weighted
sums. This is the cumulative input to the construction.

The first main estimate gives directional control of the same local
matrices. For each fixed \(L,C_{\mathrm{sep}}>1\), there are
nonnegative defect functions \(E_i\) on the working region such that,
whenever \(x,z,y\in B_4(p)\) and
\(d(x,y),d(z,y)\le C_{\mathrm{sep}}r_i\),
\[
 D_{i,y}[u(x)-u(z)]
 \le C r_i^2\bigl(E_i(x)+E_i(z)\bigr),
 \qquad \sup_{x\in B_4(p)}\sum_{i\ge0}E_i(x)\longrightarrow0
\]
as the initial rescaling becomes flat
(Proposition~\ref{pi:directional-pinching}). The proof uses Poisson
replacements of squared distance, affine approximation of radial
derivatives, and a sharp projection estimate for the difference of
two replacements. A radial transport argument preserves the
\emph{squared} defect in this last estimate. Poisson replacement,
harmonic transformation, heat Gram normalization, and sharp
squared-Hessian pinching have substantive antecedents in the work
of Cheeger and Colding and of Cheeger, Jiang and Naber
\cite{CC96,CJN21}. The two-center directional bound displayed above
is proved here in the form needed for the value-space construction.

Scale-dependent transformations are also central in
\cite[Theorem~1.32]{CN15} and \cite[Theorems~1.7 and 4.1]{HH24}.
Potentially unbounded, noncommuting normalization matrices already
occur in these methods and in \cite[Theorem~3.8]{HZ26}; those
freedoms alone do not remove the H\"older loss. The issue here is
how to use the directional estimate to control the cumulative
metric error without assuming summability of its square root.

The second main ingredient is the analytic construction in
Proposition~\ref{re:metric-iteration}. On the common open value
domain \(\Omega=u(B_1(p))\), the construction produces smooth maps
\(h_i:\Omega\to\R^n\) whose ordinary Jacobians \(Dh_i\) remain
invertible. From the forms \(D_{i,y}\) and smooth cutoffs we build
covector fields \(v_i\) and, for \(P\in\Omega\), iterate
\[
 h_0(P)=P,\qquad
 h_{i+1}(P)=h_i(P)+Dh_i(P)^{-T}v_i(P).
\]
The leading change in the pullback metric is a positive sum of the
forms \(D_{i,y}\). When the inverse transpose is differentiated, its
linear contribution to a squared length contains
\(D^2h_i(\xi,\xi)\). A recurrence for these diagonal derivatives
allows the errors to be summed using \(\sum_i E_i\), even when
\(\sum_i\sqrt{E_i}\) diverges. A separate estimate for all derivatives,
with Gevrey growth in their order, closes the iteration despite the
loss of one derivative in each update.

This part of the argument uses only positive matrix comparisons,
directional smallness, and smooth cutoffs; it does not differentiate
the scale metrics \(H_i(x)\) or assume they remain bounded. It isolates
an analytic way to turn summable directional information into
coordinates, separate from the geometric origin of that information.
The limit \(F=\lim_i h_i\circ u\) has the asserted metric bounds.
The openness of its image follows from the open coordinates \(u\)
and invariance of domain.

Summability also governs bi-Lipschitz Reifenberg parametrizations.
For extrinsic sets, David and Toro control coherent changes of
approximating planes by square-summable errors
\cite[Proposition~8.34]{DT12}. For metric spaces, Gigli and Violo
obtain bi-Lipschitz parametrizations under summability of ordinary
dyadic Gromov--Hausdorff errors \cite[Theorems~2.2 and 2.7]{GV25}.
Regularity at a point does not itself supply that summability.
Here the summable quantity is instead the directional Bochner
defect, and its conversion into metric control is proved by the
explicit iteration above.

Section~\ref{co:section} constructs the open harmonic domain and
its scale normalizations. Section~\ref{co:heat-section} proves the
cumulative matrix comparison, and Section~\ref{pi:section} proves
the summable directional estimate. Section~\ref{re:section} uses
these two inputs to construct the renormalization and prove the
all-pairs metric bounds and openness of the limiting chart.

\subsection{Conventions}

All balls use the ambient metric. Every enlargement appearing in a
local estimate is fixed before a sufficiently late rescaling is
chosen. We write \(C\) for constants depending on the dimension and
on specified fixed enlargement factors. Uniform quantities tending
to zero along the rescalings will be denoted by \(o(1)\); this
notation includes uniformity over the centers and all dyadic scales
under consideration. A subscript, as in \(o_A(1)\), permits
dependence on constants fixed earlier.

For a positive definite matrix \(A\), set
\[
 A[\xi]=\xi^T A\xi,\qquad |\xi|_A=A[\xi]^{1/2}.
\]
Matrix inequalities are inequalities of quadratic forms. Vectors
and covectors are represented using the fixed Euclidean pairing,
and their length norms are specified explicitly. The norm
\(\|T\|_{H\to Q}\) of a multilinear map uses \(|\cdot|_H\) on each
argument and \(|\cdot|_Q\) on its value. Derivatives of maps on
coordinate value space are ordinary Euclidean derivatives. We use
\(\fint_B\) for the integral divided by the measure of \(B\).
==== END Bi-Lipschitz-Coordinates-at-Regular-Points-of-Noncollapsed-RCD-Spaces//introduction.tex ====

==== BEGIN Bi-Lipschitz-Coordinates-at-Regular-Points-of-Noncollapsed-RCD-Spaces//paper.tex ====
\documentclass[11pt]{article}
\usepackage[margin=1in]{geometry}
\usepackage{amsmath,amssymb,amsthm}
\usepackage{microtype}
\usepackage[colorlinks=true,allcolors=blue]{hyperref}
\hypersetup{pdftitle={Bi-Lipschitz Coordinates at Regular Points of Noncollapsed RCD Spaces},pdfauthor={OpenAI}}
\setlength{\emergencystretch}{2em}
\numberwithin{equation}{section}
\newtheorem{theorem}{Theorem}[section]
\newtheorem{proposition}[theorem]{Proposition}
\newtheorem{lemma}[theorem]{Lemma}
\newtheorem{corollary}[theorem]{Corollary}
\theoremstyle{definition}
\newtheorem{definition}[theorem]{Definition}
\theoremstyle{remark}
\newtheorem{remark}[theorem]{Remark}
\newcommand{\R}{\mathbb R}
\newcommand{\N}{\mathbb N}
\newcommand{\m}{\mathfrak m}
\newcommand{\Id}{\mathrm{Id}}
\newcommand{\CD}{\mathrm{CD}}
\newcommand{\RCD}{\mathrm{RCD}}
\newcommand{\Hess}{\operatorname{Hess}}
\newcommand{\Ric}{\operatorname{Ric}}
\newcommand{\tr}{\operatorname{tr}}
\newcommand{\supp}{\operatorname{supp}}
\newcommand{\Lip}{\operatorname{Lip}}
\newcommand{\vol}{\operatorname{vol}}
\newcommand{\diam}{\operatorname{diam}}
\newcommand{\dd}{\,\mathrm d}
\newcommand{\fint}{\mathop{\mathchoice
 {\displaystyle\int\mkern-18mu-}
 {\textstyle\int\mkern-15mu-}
 {\scriptstyle\int\mkern-12mu-}
 {\scriptscriptstyle\int\mkern-10mu-}}\nolimits}
\DeclareMathOperator{\diver}{div}
\title{Bi-Lipschitz Coordinates at Regular Points\\
of Noncollapsed \(\RCD\) Spaces}
\author{OpenAI}
\date{September 25, 2026}
\begin{document}
\maketitle
\begin{abstract}
We resolve the regular-point bi-Lipschitz conjecture in the noncollapsed
setting. For every integer \(n\ge2\) and every real \(K\), every regular
point of a noncollapsed \(\RCD(K,n)\) space has an open neighborhood
bi-Lipschitz homeomorphic to an open subset of \(\R^n\). The bi-Lipschitz
constant depends only on \(n\), and the neighborhood uses the restricted
ambient distance.
\end{abstract}
\input{introduction}
\input{coordinates}
\input{pinching}
\input{renormalization}
\bibliographystyle{plain}
\bibliography{references}
\end{document}
==== END Bi-Lipschitz-Coordinates-at-Regular-Points-of-Noncollapsed-RCD-Spaces//paper.tex ====

==== BEGIN Bi-Lipschitz-Coordinates-at-Regular-Points-of-Noncollapsed-RCD-Spaces//pinching.tex ====
\section{Summable directional pinching}
\label{pi:section}

We prove the estimate which relates a difference of coordinate values to
the Bochner measure.  The summable quantity in this section is an energy;
its square root will occur in individual estimates, but will never be
summed without an exponentially decreasing weight.

\subsection{Conventions and domains}

We use the sequence of rescaled spaces and harmonic columns constructed
in Section~\ref{co:section}.  We suppress the index of that sequence.
Thus the lower Ricci bound is $-\kappa$, with $\kappa\to0$, and $o(1)$
is uniform over the centers and dyadic scales under consideration.
The measure is always the Hausdorff measure in the current length units.
We write $r_i=2^{-i}$, $i\ge0$, and use the matrices $R_i$ and the
matrix-valued positive measure $M(u)$ from the preceding section.

Fix numbers $L,C_{\mathrm{sep}}>1$.  The final estimate will involve the mass on
$B_{Lr_i}(y)$ and centers satisfying
\begin{equation}
 x,z,y\in B_4(p),\qquad
 d(x,y),d(z,y)\le C_{\mathrm{sep}}r_i.
 \label{pi:center-configuration}
\end{equation}
There is no restriction in replacing an initially prescribed number by
a larger number greater than one.  Set
\begin{equation}
 R=16(L+C_{\mathrm{sep}}+2),\qquad S=32R.
 \label{pi:radii}
\end{equation}
All the balls used below are contained in $B_{16S+4}(p)$ before the
additional dyadic rescaling.  We choose this fixed enlargement first
and then take a sufficiently late member of the sequence.  In
particular, the harmonic functions are defined there and all the
conclusions of Lemma~\ref{co:harmonic-normalization} apply on its
fixed enlargements.  Poisson solutions centered at $a$ will be defined
on $B_{8Sr_i}(a)$.  Their estimates will be used on $B_{2Sr_i}(a)$;
the affine fits will be made on $B_{Sr_i}(a)$; radial invariance will
be used only on $B_{Sr_i/2}(a)$.  In the two-center argument, its
projection ball is $B_{Rr_i}(z)$ and the transport calculation uses
$B_{3Rr_i}(z)$.  These choices leave strict margins between every
testing ball and its domain of definition.

For a measurable ball $B$ of positive measure, put
\[
 \|f\|_{2,B}^{\mathrm{av}}
 =\left(\fint_B |f|^2\dd\m\right)^{1/2}.
\]
Unadorned $L^2$ norms below use integrals, rather than averages.
After a dyadic rescaling, the volumes of all fixed-radius balls in
use are bounded above and below by positive constants.  Thus these
two norm conventions differ by uniformly bounded factors at each of
the fixed radii.

We use two consequences of the preceding section.  First, the
Dirichlet inequalities on the rescaled balls in use include, for
some $d>1$ and a constant depending on the fixed radius,
\begin{equation}
 \|\psi\|_{L^{2d}}^2\le C\int |\nabla\psi|^2\dd\m,
 \qquad \psi\in W^{1,2}_0(B).
 \label{pi:dirichlet-sobolev}
\end{equation}
as provided by Lemma~\ref{co:local-calculus}, specifically
\eqref{co:dirichlet-sobolev}.
The corresponding Poincar\'e inequality is included.  In dimension
two one may choose any one of the available exponents strictly
greater than two.  The domains are proper subsets of larger balls
whose complements have a fixed positive volume, so the Dirichlet
inequality is coercive.  Second, if
\[
 \sigma_i^a(b)=|b|_{R_i(a)},\qquad b\in\R^n,
\]
then, after fixing for example $a_0=1/16$ in the growth estimate,
\begin{equation}
 \frac{\sigma_j^a(b)}{\sigma_i^a(b)}
 \le C2^{a_0|i-j|}\quad(b\ne0),\qquad
 \sup_{B_{2Sr_i}(a)}|\nabla(b\cdot u)|
 \le C\sigma_i^a(b).
 \label{pi:coordinate-input}
\end{equation}
Each constant in this section may depend on $n,L,C_{\mathrm{sep}}$, and consequently
on $R,S$, but not on $i$, the centers, or the late member of the sequence.

\subsection{Positive sources and Poisson pinching}

Poisson replacements of squared distance follow the almost-rigidity
method of Cheeger and Colding \cite{CC96}; see also the quantitative
regularization in \cite[Section~6.1]{CJN21}. We retain the full
positive Bochner measure and prove the uniform estimates needed at
all testing centers.

The following elementary Dirichlet-form observation is useful because
an application of Cauchy--Schwarz to the energy estimate alone would
not give the required dependence on the mass.

\begin{lemma}[A linear gradient bound for a positive source]
\label{pi:positive-source}
Let $\Omega$ be one of the bounded rescaled domains satisfying
\eqref{pi:dirichlet-sobolev}, and suppose that $h\in W^{1,2}_0(\Omega)$
is nonnegative and bounded, with
\[
 \int\langle\nabla h,\nabla\psi\rangle\dd\m
 =\int\psi\dd\mu
\]
for compactly supported Lipschitz tests, where $\mu$ is a finite
positive measure.  If $\mu(\Omega)\le P$, then
\begin{equation}
 \int_\Omega|\nabla h|\dd\m\le CP,
 \qquad
 \int_\Omega|\nabla h|^2\dd\m\le\|h\|_\infty P.
 \label{pi:positive-source-bounds}
\end{equation}
The constant in the first inequality depends only on the Sobolev
constant, its exponent, and an upper bound for $\m(\Omega)$.
\end{lemma}

\begin{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 $\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\dd\m\le t.
\]
Applying \eqref{pi:dirichlet-sobolev} to the same truncation yields
\[
 t^2\m(H\ge t)^{1/d}\le Ct,
 \qquad \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
\begin{equation}
 \int_\Omega(1+H)^{1+\alpha}\dd\m\le C.
 \label{pi:source-moment}
\end{equation}
The bounded test $1-(1+H)^{-\alpha}$ gives
\[
 \alpha\int_\Omega
 \frac{|\nabla H|^2}{(1+H)^{1+\alpha}}\dd\m\le1.
\]
Combining this inequality with \eqref{pi:source-moment} 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\dd\mu\le\|h\|_\infty P$.
\end{proof}

For a function with constant Laplacian and a scalar $b$, write
\begin{equation}
 \mathcal B_\kappa(w,b)
 =M(w)-(2b\Delta w-nb^2)\m.
 \label{pi:shifted-measure}
\end{equation}
By Lemma~\ref{co:local-calculus}, this is a nonnegative measure and
\begin{equation}
 \mathcal B_\kappa(w,b)
 \ge |\Hess w-bg|^2\m.
 \label{pi:shifted-lower-bound}
\end{equation}
It is quadratic in the pair $(w,b)$.  Positivity, applied to the sum
and difference of two such pairs, gives the measure inequality
\begin{equation}
 \mathcal B_\kappa(w_1-w_2,b_1-b_2)
 \le2\mathcal B_\kappa(w_1,b_1)
    +2\mathcal B_\kappa(w_2,b_2).
 \label{pi:shifted-triangle}
\end{equation}
These statements use the noncollapsed dimension and trace identity
and the improved Bochner inequality; see \cite{DPG18,Han18,Gigli18}.

\begin{lemma}[Poisson replacements with summable defects]
\label{pi:poisson-pinching}
There are nonnegative numbers $P_i(a)$, defined for every
$a\in B_4(p)$ and $i\ge0$, with
\begin{equation}
 \sup_{a\in B_4(p)}\sum_{i=0}^{\infty}P_i(a)=o(1),
 \label{pi:poisson-sum}
\end{equation}
and functions $f_i^a$ on $B_{8Sr_i}(a)$ such that, with
$q^a=d_a^2/2$,
\begin{equation}
 \Delta f_i^a=a_i:=n+C_0\kappa r_i^2,
 \qquad f_i^a-q^a\in W^{1,2}_0(B_{8Sr_i}(a)).
 \label{pi:poisson-definition}
\end{equation}
Here $C_0$ depends only on $n,S$ and is the same for every center
and every scale.  The following bounds hold:
\begin{align}
 0\le q^a-f_i^a&\le Cr_i^2,
 &\sup_{B_{4Sr_i}(a)}|\nabla f_i^a|&\le Cr_i,
 \label{pi:poisson-boundedness}\\
 \int_{B_{8Sr_i}(a)}|\nabla(f_i^a-q^a)|^2\dd\m
 &\le Cr_i^{n+2}P_i(a),
 &\int_{B_{8Sr_i}(a)}|\nabla(f_i^a-q^a)|\dd\m
 &\le Cr_i^{n+1}P_i(a),
 \label{pi:poisson-errors}\\
 \mathcal B_\kappa(f_i^a,1)(B_{2Sr_i}(a))
 &\le Cr_i^nP_i(a).
 \label{pi:poisson-shifted-mass}
\end{align}
\end{lemma}

\begin{proof}
Choose $C_0$ sufficiently large that Laplacian comparison gives
\[
 \Delta q^a\le(n+C_0\kappa r_i^2)\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\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
\cite[Corollary~5.15, equation~(5.29)]{Gigli15}.
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 \eqref{pi:poisson-definition}.

Here is the mass estimate, including its normalization.  Set
$V_a(t)=\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)\dd t$; 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}}\dd t\notag\\
 &\le C\{A_a(8Sr_i)-A_a(16Sr_i)+\kappa r_i^2\}.
 \label{pi:volume-drop}
\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
\eqref{pi:poisson-sum}.  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
\eqref{pi:dirichlet-sobolev}.  Work in units $r_i=1$.  Let $\tau$
solve $-\Delta\tau=1$ with zero boundary values on $B_{8S}(a)$.
If $Y(k)=\m(\tau>k)$, testing with $(\tau-k)_+$ and applying
Sobolev and H\"older 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 \eqref{pi:poisson-boundedness}.

 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 \cite{Kell13,Jiang14};
the local test-function interpretation is supplied by
Lemma~\ref{co:local-calculus} (see also \cite{CDNZPSW}).
Apply Lemma~\ref{pi:positive-source} to $h_i^a$ on this same
rescaled Dirichlet domain.  Its source mass is at most $P_i(a)$.
This proves both estimates in \eqref{pi:poisson-errors}, 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|\dd\m
 \le\int\bigl(|\nabla q^a||\nabla h|
                  +\tfrac12|\nabla h|^2\bigr)\dd\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\m-\mu_i^a$ now gives
\begin{align*}
 \int\chi\dd\mathcal B_\kappa(f_i^a,1)
 &=\frac12\int|\nabla f_i^a|^2\Delta\chi\dd\m
   +\kappa\int\chi|\nabla f_i^a|^2\dd\m
   -(2a_i-n)\int\chi\dd\m\\
 &=(n-a_i)\int\chi\dd\m-\int\chi\dd\mu_i^a
   +\kappa\int\chi|\nabla f_i^a|^2\dd\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
\eqref{pi:poisson-shifted-mass}.
\end{proof}

\begin{remark}[The units of the Poisson estimates]
\label{pi:poisson-units}
Under $d'=r_i^{-1}d$, $\m'=r_i^{-n}\m$ and
$\widehat f_i^a=r_i^{-2}f_i^a$, one has
$\Delta'\widehat f_i^a=a_i$ and $\kappa'=\kappa r_i^2$.
The measure $\mathcal B_{\kappa'}(\widehat f_i^a,1)$ is
$r_i^{-n}\mathcal B_\kappa(f_i^a,1)$.  Consequently all the
Poisson bounds above have size $CP_i(a)$ in their natural rescaled
units.  In particular, by \eqref{pi:shifted-triangle}, differences
of two replacements at the same scale are harmonic and have
Bochner mass at most $C(P_i(x)+P_i(z))$ on any of the common
interior balls.  This estimate controls the full measure, including
its possible singular part.
\end{remark}

\subsection{Radial differentiation of the coordinate span}

The Poisson estimates provide a summable scalar budget. We next
fit the radial derivatives by affine combinations of the harmonic
coordinates and distribute the errors over neighboring scales.

We first record the weak cancellation used to replace a radial
derivative by a harmonic function.

\begin{lemma}[Weak radial cancellation]
\label{pi:weak-cancellation}
Suppose that $f$ has constant Laplacian and $a$ is harmonic on a
ball on which both functions have bounded gradients.  Put
$g_a=\langle\nabla f,\nabla a\rangle$.  For interior Sobolev tests
$\phi\in W^{1,2}_0$ one has
\begin{equation}
 \int\langle\nabla g_a,\nabla\phi\rangle\dd\m
 =2\int(\Hess f-g)(\nabla a,\nabla\phi)\dd\m.
 \label{pi:weak-cancellation-equation}
\end{equation}
The identity is local, with the right side interpreted on a
slightly larger interior ball if a cutoff is needed.
\end{lemma}

\begin{proof}
Local test regularity and the Hessian product rule give
\[
 \nabla g_a=(\Hess f)(\nabla a,\cdot)
             +(\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 $\Hess\phi$ cancel, leaving
\begin{align*}
 &\int\Hess a(\nabla f,\nabla\phi)\dd\m
  -\int\Hess f(\nabla a,\nabla\phi)\dd\m\\
 &\hspace{1cm}
 =-\int\Delta f\,\langle\nabla a,\nabla\phi\rangle\dd\m
   +\int\Delta a\,\langle\nabla f,\nabla\phi\rangle\dd\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 $\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.
\end{proof}

\begin{lemma}[Harmonic affine improvement]
\label{pi:affine-improvement}
There is an integer $m\ge3$, depending only on $n,S$, such that,
with $\theta=2^{-m}$, the following holds for all sufficiently
late spaces, all $i$ and all $a\in B_4(p)$.  If $w$ is harmonic
on $B_{Sr_i}(a)$, then there is an affine combination $F$ of the
components of $u$ with
\begin{equation}
 \|w-F\|_{2,B_{S\theta r_i}(a)}^{\mathrm{av}}
 \le\tfrac12\theta^{3/2}
       \|w\|_{2,B_{Sr_i}(a)}^{\mathrm{av}}.
 \label{pi:affine-improvement-equation}
\end{equation}
\end{lemma}

\begin{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 \cite{AH18}.  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~\ref{co:harmonic-normalization}.  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.
\end{proof}

\begin{lemma}[Summable radial invariance]
\label{pi:radial-invariance}
There are functions $E_i(a)>0$, for $a\in B_4(p)$, such that
\begin{equation}
 E_i(a)\ge P_i(a),\qquad
 \sup_a\sum_{i=0}^\infty E_i(a)=o(1),\qquad
 \inf_a E_i(a)>0\quad\hbox{for each fixed }i.
 \label{pi:pinching-summability}
\end{equation}
Let $T=R_i(y)^{-1/2}$ with $d(a,y)\le C_{\mathrm{sep}}r_i$, and put
$U=T(u-u(a))/r_i$.  In the metric $r_i^{-1}d$, with
$\widehat f_i^a=f_i^a/r_i^2$, there is a matrix $A_i^a(T)$
satisfying $\|A_i^a(T)-\Id\|=o(1)$ and
\begin{equation}
 \left\|\langle\nabla\widehat f_i^a,\nabla U\rangle
       -A_i^a(T)(U-U(a))\right\|_{W^{1,2}(B_{S/2}(a))}
 \le C\sqrt{E_i(a)}.
 \label{pi:radial-invariance-equation}
\end{equation}
The ball and all derivatives in this display use the rescaled
metric.  Arbitrary additive constants may be used in $U$.
\end{lemma}

\begin{proof}
\emph{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 \eqref{pi:coordinate-input} imply
\begin{equation}
 \|g_b-g_{j,b}\|_{2,B_{Sr_j}(a)}^{\mathrm{av}}
 \le Cr_j\sigma_j^a(b)\sqrt{P_j(a)}.
 \label{pi:raw-radial-error}
\end{equation}
By Lemma~\ref{pi:weak-cancellation}, the equation for $g_{j,b}$
on this ball has divergence source
$2(\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\'e gives
\begin{equation}
 \|\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)}.
 \label{pi:radial-harmonic-correction}
\end{equation}
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~\ref{pi:affine-improvement} to
$g_{j,b}-\xi_{j,b}-F_{j,b}$.  The two errors in
\eqref{pi:raw-radial-error} and
\eqref{pi:radial-harmonic-correction}, restricted to the smaller
ball, acquire only the fixed volume factor depending on $m,S$.
Dividing by $r_{j+m}=\theta r_j$ gives
\begin{equation}
 e_{j+m}(b)\le\theta^{1/2}e_j(b)
                  +C\sigma_j^a(b)\sqrt{P_j(a)}.
 \label{pi:fit-recurrence}
\end{equation}
The factor $1/2$ in the harmonic improvement leaves more than the
needed margin in the coefficient of $e_j(b)$.

\smallskip\noindent
\emph{Qualitative convergence.}
We will also use the uniform estimate
\begin{equation}
 e_j(b)\le\varepsilon\sigma_j^a(b),\qquad\varepsilon=o(1).
 \label{pi:qualitative-radial-error}
\end{equation}
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\dd\m
 =-\int(q-|V|^2/2)
       \operatorname{div}(\chi\nabla(|V|^2/2))\dd\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
\eqref{pi:qualitative-radial-error}.

\smallskip\noindent
\emph{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
\eqref{pi:fit-recurrence} in each of the finitely many residue
classes modulo $m$, and using \eqref{pi:coordinate-input}, yields
\begin{equation}
 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).
 \label{pi:smoothed-fit-error}
\end{equation}
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 \eqref{pi:smoothed-fit-error}.  This
calculation also covers the finitely many initial indices.

\smallskip\noindent
\emph{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 \eqref{pi:smoothed-fit-error},
\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).
 \label{pi:inward-constant-sum}
\end{align}
Fix $\beta=1/16$ and define, with a sufficiently large fixed
constant $C_E$,
\begin{equation}
 E_i(a)=C_E\sum_{l=0}^\infty
                  2^{-\beta|i-l|}\widetilde P_l(a).
 \label{pi:definition-E}
\end{equation}
Weighted Cauchy--Schwarz in \eqref{pi:inward-constant-sum}
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}$.

\smallskip\noindent
\emph{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~\ref{pi:weak-cancellation}, 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 \eqref{pi:radial-harmonic-correction}.

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
\eqref{pi:radial-invariance-equation}.  The qualitative radial
convergence proved above, the Poisson error, and affine
nondegeneracy imply
$\|A_i^a(T)-\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)-\Id=o(1)$ persist.  Translating $U$ changes neither
side when its center value is subtracted.
\end{proof}

\subsection{Two centers and the sharp coordinate projection}

In Euclidean space, the difference $q^z-q^x$ of the squared-distance
potentials is affine, with gradient $x-z$.  We need a
quantitative version of both facts for the equal-Laplacian Poisson
replacements: their harmonic difference must be close to an affine
coordinate function, and its coefficient must be close to the
coordinate difference of the centers.

Fix a configuration \eqref{pi:center-configuration}.  In this
subsection take $r_i$ as the length unit and replace the measure
by $r_i^{-n}\m$.  Write
\[
 f^x=f_i^x/r_i^2,\qquad f^z=f_i^z/r_i^2,
 \qquad U=R_i(y)^{-1/2}(u-u(z))/r_i,
 \qquad s^2=E_i(x)+E_i(z).
\]
All balls and derivatives until the final scale conversion are in
these units.  We write $M$ for the Bochner measure with lower-bound
parameter $\kappa r_i^2$ in these units, and $g$ for the rescaled
metric tensor.  The equal-source construction gives
$\Delta f^x=\Delta f^z=a_i$.  Put
\begin{equation}
 v=f^z-f^x,\qquad
 L_a w=\langle\nabla f^a,\nabla w\rangle\ (a=x,z),\qquad
 h=\langle\nabla v,\nabla U\rangle.
 \label{pi:pair-notation}
\end{equation}
The domain nesting in \eqref{pi:radii} gives
\[
 B_{8R}(z)\subset B_{S/2}(x)\cap B_{S/2}(z),
 \qquad
 B_{2L+4}(y)\subset B_{R/2}(z).
\]
All these inclusions are strict with a fixed margin, since
$d(x,z)\le2C_{\mathrm{sep}}$ and $S/2=16R$.

The commutator below first identifies $h$ with $U(x)-U(z)$ up to
$O_{L^2}(s)$.  The sharp projection then gives an affine function
$c_0+b\cdot U$ whose gradient differs from that of $v$ by
$O_{L^2}(s)$ on the required interior ball.  Averaging against the
normalized Gram matrix recovers $b=U(x)-U(z)+O(s)$.  These two
conclusions transfer the $Cs^2$ Bochner bound for $v$ to the
coordinate function in the direction $U(x)-U(z)$.

\begin{lemma}[The two-center commutator]
\label{pi:commutator}
On $B_{8R}(z)$ one has
\begin{equation}
 \Delta v=0,\qquad
 |\nabla v|\le C,\qquad
 M(v)(B_{8R}(z))\le Cs^2,
 \label{pi:pair-poisson-bounds}
\end{equation}
and
\begin{equation}
 \|h-(U(x)-U(z))\|_{L^2(B_{8R}(z))}\le Cs.
 \label{pi:constant-direction}
\end{equation}
\end{lemma}

\begin{proof}
The first assertions follow from the common value of the Laplacian,
the Poisson gradient bounds, and
\eqref{pi:shifted-triangle}; see
Remark~\ref{pi:poisson-units}.  In particular
\[
 \|\Hess v\|_{L^2(B_{8R}(z))}
 +\|\Hess f^z-g\|_{L^2(B_{8R}(z))}\le Cs.
\]
Let $A_x,A_z$ be the matrices in
Lemma~\ref{pi:radial-invariance} 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
\begin{equation}
 h=BU+c+e,\qquad
 \|e\|_{W^{1,2}(B_{8R}(z))}\le Cs.
 \label{pi:commutator-affine}
\end{equation}
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=\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
 =\Hess f^z(\nabla v,\nabla U)
       -\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
\begin{equation}
 L_zh=(A_z-\Id)h+O_{L^2}(s).
 \label{pi:commutator-equation}
\end{equation}
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 \eqref{pi:commutator-affine} on both sides of
\eqref{pi:commutator-equation}, and use the radial identity once
more to differentiate $BU$.  This gives
\[
 \bigl(BA_z-(A_z-\Id)B\bigr)U-(A_z-\Id)c
 =O_{L^2}(s).
\]
Affine nondegeneracy, for example on $B_R(z)$, bounds its linear
and constant coefficients separately:
\begin{equation}
 \|BA_z-(A_z-\Id)B\|+|(A_z-\Id)c|\le Cs.
 \label{pi:commutator-coefficients}
\end{equation}
Writing $C_z=A_z-\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-\Id)c=(A_z-\Id)c-Bc=O(s).
\]
Since $c=A_x\delta$ and $A_x^{-1}$ is bounded, this yields
\begin{equation}
 \delta-c=-A_x^{-1}(A_x-\Id)c=O(s).
 \label{pi:delta-from-c}
\end{equation}
Finally \eqref{pi:commutator-affine}, the bounded coordinate
values, and \eqref{pi:delta-from-c} prove
\eqref{pi:constant-direction} on the full ball $B_{8R}(z)$.
\end{proof}

The next lemma supplies the norm bound on a ball larger than the
projection domain.  This is the step which permits normalization
of a potentially very small projected residual.

\begin{lemma}[Radial transport of an unnormalized residual]
\label{pi:radial-transport}
Work in the rescaled units above.  Suppose that $w$ is a bounded
Sobolev function on $B_{4R}(z)$ with $\|w\|_\infty\le C_1$,
and that
\begin{equation}
 L_zw=w+\ell\cdot U+k+e_w,
 \qquad \|e_w\|_{L^2(B_{4R}(z))}\le C_1s.
 \label{pi:transport-hypothesis}
\end{equation}
Put $\alpha=\|w\|_{L^2(B_R(z))}$.  If
$|\ell|+|k|\le C_1(\alpha+s)$, then
\begin{equation}
 \|w\|_{L^2(B_{2R}(z))}\le C(\alpha+s),
 \label{pi:transport-conclusion}
\end{equation}
where $C$ is independent of $\alpha,s$.  The constants may depend
on $C_1$ and the already fixed radii and coordinate bounds.
\end{lemma}

\begin{proof}
The Sobolev product rule and $\Delta f^z=a_i$ give
\begin{equation}
 \operatorname{div}(w^2\nabla f^z)
 =(a_i+2)w^2+wF,\qquad
 F=2(\ell\cdot U+k+e_w),
 \label{pi:transport-divergence}
\end{equation}
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
\begin{equation}
 -t\phi'(t)-(a_i+2)\phi(t)
 \ge |\phi'(t)|+\phi(t)
 \quad\hbox{for almost every }R<t<3R.
 \label{pi:transport-cutoff-sign}
\end{equation}
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\dd\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
\eqref{pi:transport-divergence} against $\phi(d_z)$.
Using \eqref{pi:transport-cutoff-sign}, and moving the inner-ball
term to the right, gives
\begin{equation}
 \mathcal A\le C\alpha^2
   +\left|\int\phi wF\dd\m\right|
   +\left|\int_{R<d_z<3R}\phi'e_f w^2\dd\m\right|.
 \label{pi:transport-energy}
\end{equation}
The decisive coefficient error has the bound
\begin{align}
 \left|\int_{R<d_z<3R}\phi'e_f w^2\dd\m\right|
 &\le\|w\|_\infty
       \left(\int|\phi'|w^2\dd\m\right)^{1/2}
       \left(\int|\phi'|e_f^2\dd\m\right)^{1/2}\notag\\
 &\le Cs\sqrt{\mathcal A}
 \le\tfrac14\mathcal A+Cs^2.
 \label{pi:transport-square-error}
\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\dd\m\right|
 \le\tfrac14\mathcal A+C(\alpha^2+s^2).
\]
Substitution in \eqref{pi:transport-energy} proves
$\mathcal A\le C(\alpha^2+s^2)$.  The positive lower bound for
$\phi$ up to radius $2R$ now gives
\eqref{pi:transport-conclusion}.

In particular, the fixed bound for $\|w\|_\infty$ in
\eqref{pi:transport-square-error} is used before $w$ is divided
by $\alpha$.  The conclusion is quadratic in $s$ at the level
of squared norms.
\end{proof}

\begin{lemma}[Sharp approximation by the harmonic coordinate span]
\label{pi:sharp-projection}
For all sufficiently late spaces, there is a vector $b\in\R^n$
with $|b|\le C$ such that
\begin{equation}
 \|\nabla(v-b\cdot U)\|_{L^2(B_{2L+4}(y))}\le Cs.
 \label{pi:sharp-projection-equation}
\end{equation}
The constants and the late-space threshold are uniform in the
configuration \eqref{pi:center-configuration} and the dyadic scale.
\end{lemma}

\begin{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
\begin{equation}
 \nabla(L_zv-v)
 =(\Hess f^z-g)(\nabla v,\cdot)
    +\Hess v(\nabla f^z,\cdot),
 \qquad
 \|\nabla(L_zv-v)\|_{L^2(B_{8R}(z))}\le Cs.
 \label{pi:euler-gradient}
\end{equation}
Poincar\'e 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
\begin{equation}
 L_zw=w+\ell\cdot U+k+e_w,
 \qquad
 \|e_w\|_{L^2(B_{4R}(z))}\le Cs,
 \label{pi:projected-euler-equation}
\end{equation}
where $\ell^T=b^T(\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
\eqref{pi:projected-euler-equation} now bound the $L^2$ norm
of $\ell\cdot U+k$ on this smaller ball by $C(\alpha+s)$.
Affine nondegeneracy gives
\begin{equation}
 |\ell|+|k|\le C(\alpha+s).
 \label{pi:forcing-coefficient-bound}
\end{equation}
Lemma~\ref{pi:radial-transport} consequently applies and yields
\begin{equation}
 \|w\|_{L^2(B_{2R}(z))}\le C(\alpha+s).
 \label{pi:projection-outer-bound}
\end{equation}

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 \eqref{pi:pinching-summability} makes $s>0$, so
$\widehat w=w/\alpha$ is defined along this sequence.
The norm on $B_R(z)$ is one, and
\eqref{pi:projection-outer-bound} 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 \cite{AH18}.
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
\eqref{pi:forcing-coefficient-bound}; 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~\ref{pi:radial-invariance}.  The normalized harmonic
gradients have uniform interior bounds and converge strongly.
Thus scalar products of these gradients pass to the limit in
\eqref{pi:projected-euler-equation}, which becomes
\begin{equation}
 X\cdot\nabla w_\infty
 =w_\infty+\ell_\infty\cdot X+k_\infty
 \quad\hbox{on }B_{3R/2}(0).
 \label{pi:limiting-euler-equation}
\end{equation}
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\dd X=0,\qquad
 \int_{B_R(0)}w_\infty X_j\dd X=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
\eqref{pi:limiting-euler-equation} 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
\eqref{pi:sharp-projection-equation}.
\end{proof}

\begin{proposition}[Summable directional pinching]
\label{pi:directional-pinching}
For every fixed $L,C_{\mathrm{sep}}>1$, the functions $E_i$ may be chosen
as in \eqref{pi:definition-E}, with the uniform summability and
positive floors in \eqref{pi:pinching-summability}, so that the
following holds for all sufficiently late spaces.  Define
\[
 D_{i,y}=R_i(y)^{-1}
       \bigl[r_i^{2-n}M(u)(B_{Lr_i}(y))\bigr]R_i(y)^{-1}.
\]
For every configuration \eqref{pi:center-configuration},
\begin{equation}
 D_{i,y}[u(x)-u(z)]
 \le Cr_i^2\bigl(E_i(x)+E_i(z)\bigr).
 \label{pi:directional-estimate}
\end{equation}
The same functions work simultaneously for all such centers and
all dyadic scales.  All enlargement constants are fixed before
choosing the late-space threshold.
\end{proposition}

\begin{proof}
Continue first in the rescaled units used in the preceding
subsection, and choose $b$ from
Lemma~\ref{pi:sharp-projection}.  The normalization at $y$ is
exactly
\[
 \fint_{B_1(y)}
   (\langle\nabla U_j,\nabla U_k\rangle)_{j,k}\dd\m=\Id.
\]
Consequently
\begin{align*}
 \fint_{B_1(y)}h\dd\m
 &=b+\fint_{B_1(y)}
       \langle\nabla U,\nabla(v-b\cdot U)\rangle\dd\m
 =b+O(s).
\end{align*}
Here the gradient bound for $U$ and
\eqref{pi:sharp-projection-equation} control the last error.
Equation~\eqref{pi:constant-direction} then gives
\begin{equation}
 b=U(x)-U(z)+O(s)=\delta+O(s).
 \label{pi:coefficient-recovery}
\end{equation}

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\dd\m
    +\kappa r_i^2\int\zeta|\nabla w|^2\dd\m\notag\\
 &\le Cs^2.
 \label{pi:residual-mass}
\end{align}
Together with \eqref{pi:pair-poisson-bounds}, 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~\ref{co:bochner-small}.  Therefore
\eqref{pi:coefficient-recovery} implies
\[
 M((\delta-b)\cdot U)(B_L(y))
 \le|\delta-b|^2\tr M(U)(B_L(y))\le Cs^2.
\]
Another positive quadratic-form triangle inequality gives
\begin{equation}
 M(\delta\cdot U)(B_L(y))\le Cs^2.
 \label{pi:normalized-directional-mass}
\end{equation}
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$, $\m'=r_i^{-n}\m$
and $U=Tu/r_i$ up to an additive constant, the matrix-valued
measure transforms as
\begin{equation}
 M'(U)=r_i^{2-n}T M(u)T^T.
 \label{pi:matrix-measure-scaling}
\end{equation}
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 \eqref{pi:normalized-directional-mass} and
$s^2=E_i(x)+E_i(z)$ proves
\eqref{pi:directional-estimate}.
\end{proof}

\begin{remark}[Use with the heat-normalized matrices]
\label{pi:heat-interface}
The construction of $E_i$ and the radial argument used $R_i$ for
the scalar coordinate norms and did not use the heat matrices.
Lemma~\ref{co:heat-matrices}, specifically
\eqref{co:QR-comparison}, identifies $Q_i$ and $R_i$ up to
uniform relative $o(1)$.  They may therefore be interchanged in
bounded norm comparisons in the next section.  This introduces
only multiplicative constants; it adds no per-scale error term
to the summable quantities in
\eqref{pi:pinching-summability} or
\eqref{pi:directional-estimate}.  In particular the explicit
floor $E_i(a)\ge C_E\eta2^{-\beta i}$ is available when choosing
a single center that nearly minimizes $E_i$ on an active
coordinate patch.
\end{remark}
==== END Bi-Lipschitz-Coordinates-at-Regular-Points-of-Noncollapsed-RCD-Spaces//pinching.tex ====

==== BEGIN Bi-Lipschitz-Coordinates-at-Regular-Points-of-Noncollapsed-RCD-Spaces//renormalization.tex ====
\section{Renormalization in the space of coordinate values}
\label{re:section}

We now turn the harmonic coordinates into a bi-Lipschitz chart.  The
construction takes place in the ordinary Euclidean space of their values.
In particular, the differentiations below do not involve a pointwise inverse
of the differential of the harmonic map.

\subsection{Notation and quantitative inputs}
\label{re:input-data}

Write $r_i=2^{-i}$, $P_x=u(x)$, and
\[
  \Omega=u(B_1(p)).
\]
By Lemma~\ref{co:open-coordinates}, $\Omega$ is open and $u$ is a
homeomorphism from $B_1(p)$ onto $\Omega$.  We use the positive definite
matrices $Q_i(x)$ and $H_i(x)=Q_i(x)^{-1}$ of
Lemma~\ref{co:heat-matrices}.  For a positive definite matrix $H$,
put $|\xi|_H=(\xi^T H\xi)^{1/2}$.  Vectors in the source value space
are measured with $H_i(x)$; covectors in that space are measured with
$Q_i(x)$, its dual matrix.  If $T$ is an $m$-linear map with covector
values, our operator norm is
\[
 \|T\|_{H\to Q}
 =\sup\bigl\{|T(\xi_1,\ldots,\xi_m)|_Q:
                 |\xi_a|_H\le1\ (1\le a\le m)\bigr\}.
\]
The notation $\|T\|_{H\to\R^n}$ has the analogous meaning with the
Euclidean target norm.  For $m=0$ these conventions mean the target norm
of the value.  We denote ordinary Fr\'echet derivatives in value space
by $\mathrm D^m$.  Whenever a derivative norm at $P_x$ uses $H_i(x)$ or
$Q_i(x)$, these matrices are frozen: they are not differentiated.

The fixed constants from the preceding sections, including all required
enlargements of balls, will be chosen before the final rescaling.  We
may then use a small parameter $\varepsilon>0$ tending to zero along
that sequence of rescalings.  Increasing this parameter by a fixed
constant factor if necessary, the following bounds hold:
\begin{equation}
 \label{re:monotone-metrics}
 H_i(x)\le H_{i+1}(x)\le(1+\varepsilon)H_i(x),
 \qquad Q_0(x)=\Id+o(1).
\end{equation}
The matrices $R_i(y)^{-1}$ and $H_i(x)$ are uniformly comparable if
$d(x,y)\le C r_i$, for each of the fixed constants $C$ used below.
Moreover, normalized harmonic approximation gives the separation and
central-ellipsoid conclusions of
Lemmas~\ref{co:harmonic-normalization} and
\ref{co:open-coordinates}.  In particular, there are fixed $a>0$ and
$C_{\mathrm{cmp}}\ge1$ such that, whenever the balls involved are
interior,
\begin{equation}
 \label{re:ellipsoid-input}
 \{P:|P-P_x|_{H_i(x)}<a r_i\}
       \subset u(B_{r_i}(x)),
 \qquad H_i(y)\le C_{\mathrm{cmp}}H_i(x)
       \quad(d(x,y)\le r_i).
\end{equation}
Here and below, the harmonic map is defined on the larger balls supplied
by the preceding sections, even when an image in
\eqref{re:ellipsoid-input} is not contained in $\Omega$.

For a fixed sufficiently large $L$, recall the positive semidefinite
matrices
\begin{equation}
 \label{re:patch-matrix}
 D_{i,y}=R_i(y)^{-1}
       \bigl[r_i^{2-n}M(B_{Lr_i}(y))\bigr]R_i(y)^{-1}.
\end{equation}
Lemma~\ref{co:cumulative-comparison} applies to boundedly many
centers at distance $O(r_i)$ from a testing point, with weights in
$[0,1]$ and with at least one weight equal to one whose mass ball contains
$B_{r_i}(x)$.  With the notation $Z_i$ for the unweighted sum and $S_i$
for the weighted sum, it gives a fixed constant $C_*\ge1$ such that
\begin{equation}
 \label{re:cumulative-input}
 \begin{aligned}
 Z_i(x)&\le\varepsilon H_i(x),\\
 \Id+\sum_{i<k}Z_i(x)&\le C_*H_k(x),\\
 H_k(x)&\le C_*\left(\Id+\sum_{i<k}S_i(x)\right).
 \end{aligned}
\end{equation}
Enlarging $C_*$ if needed, these statements include $k=0$.

By Proposition~\ref{pi:directional-pinching}, for these fixed choices
of $L$ and center separations there are functions $E_i$ on a fixed
domain containing $B_1(p)$ such that
\begin{equation}
 \label{re:pinching-input}
 \sup_{x\in B_1(p)}\sum_{i=0}^\infty E_i(x)\le\varepsilon,
 \qquad E_i(x)\ge e_i^*>0,
\end{equation}
where $e_i^*$ is independent of $x$ at each fixed $i$, and
\begin{equation}
 \label{re:directional-input}
 D_{i,y}[P_x-P_z]
       \le C r_i^2\bigl(E_i(x)+E_i(z)\bigr)
 \quad\text{if }d(x,y),d(z,y)\le C_0r_i.
\end{equation}
We use $D[\xi]=\xi^TD\xi$ for quadratic forms.  The positive floor in
\eqref{re:pinching-input} is supplied by the exponentially decreasing
positive term in the construction of $E_i$.  Such a term can always
be added with arbitrarily small sum.  It is used only to select fixed
centers for the patches.

All later constants depend on dimension and the fixed preceding
constants, but not on the scale, the testing point, or the number of
scales.  Constants denoted by $C_A$ may additionally depend on a
bootstrap constant $A$ and on the jet constants fixed after $A$.

\subsection{Smooth fields with directional control}

For each $i$ choose a maximal $r_i$-separated subset
$\mathcal Y_i\subset B_1(p)$.  Properness and local doubling imply
that this family is finite; maximality implies that every point of
$B_1(p)$ is at distance at most $r_i$ from some member.

Fix a smooth function $\psi:\R^n\to[0,1]$ equal to one on
$\{|X|\le2\}$ and supported in $\{|X|<3\}$, with
\begin{equation}
 \label{re:bump-jets}
 \|\mathrm D^m\psi\|_{\infty}
       \le C_b^{m+1}(m!)^s\qquad(m\ge0)
\end{equation}
for a fixed $s>1$.  One may, for example, take $s=2$.  To recall why
such bounds are available, the function $e^{-1/t}$ for $t>0$, extended
by zero for $t\le0$, has on each fixed bounded interval derivative
bounds $C^{m+1}(m!)^2$.  Cauchy's estimate on a complex disk of radius
a fixed small multiple of $t$ bounds its $m$th derivative by
$m!(C/t)^m e^{-c/t}$, and optimization in $t$ gives this assertion.
Products, integration, and composition with a fixed quadratic polynomial
give a radial cutoff with the specified radii and the same factorial
order, by the finite product and chain formulas.

Define
\[
 \phi_{i,y}(P)=
 \psi\left(r_i^{-1}R_i(y)^{-1/2}(P-P_y)\right).
\]
For $x\in B_1(p)$ let
\[
 \mathcal A_i(x)=\{y\in\mathcal Y_i:P_x\in\supp\phi_{i,y}\}.
\]
The harmonic separation estimates show that every member of
$\mathcal A_i(x)$ satisfies $d(x,y)\le C_0r_i$ for a fixed $C_0$.
Indeed, fix a growth exponent $\alpha=1/16$ in~\eqref{co:slow-R}.
For $r_i\le t=d(x,y)\le1$, choose $j\le i$ with $t\le r_j<2t$.
The separation estimate at scale $r_j$, followed by inverse matrix
comparison in~\eqref{co:slow-R}, gives
\[
 |P_x-P_y|_{R_i(y)^{-1}}
 \ge c(r_i/r_j)^\alpha |P_x-P_y|_{R_j(y)^{-1}}
 \ge c r_i(t/r_i)^{1-\alpha}.
\]
Membership in the cutoff support bounds the left side by $3r_i$, so
$t/r_i$ is bounded. For $t>1$ in $B_1(p)$, the initial harmonic
approximation at scale one gives the same conclusion by comparison
with $R_0(y)$. Packing
therefore gives $|\mathcal A_i(x)|\le N$ for a fixed $N$.
For sufficiently small flatness, a net point with $d(x,y)\le r_i$
satisfies $|P_x-P_y|_{R_i(y)^{-1}}<2r_i$, so that
$\phi_{i,y}(P_x)=1$.  Choose $L>C_0+2$ in
\eqref{re:patch-matrix}.  In particular the mass ball of this
weight-one center contains $B_{r_i}(x)$.

These facts justify the application of
\eqref{re:cumulative-input} to
\begin{equation}
 \label{re:active-sums}
 Z_i(x)=\sum_{y\in\mathcal A_i(x)}D_{i,y},\qquad
 S_i(x)=\sum_{y\in\mathcal A_i(x)}
                     \phi_{i,y}(P_x)D_{i,y}.
\end{equation}
We fix these sums from now on.  They are used only pointwise as
quadratic forms; no spatial differentiability of them is asserted.

For a patch let
\[
 T_{i,y}=\{x\in B_1(p):P_x\in\supp\phi_{i,y}\}.
\]
It is nonempty, since it contains $y$.  All its points lie within
$C_0r_i$ of $y$.  By the positive floor in
\eqref{re:pinching-input}, its infimum of $E_i$ is positive.
Choose one point $z_{i,y}\in T_{i,y}$ satisfying
\begin{equation}
 \label{re:fixed-center}
 E_i(z_{i,y})<2\inf_{x\in T_{i,y}}E_i(x),
 \qquad c_{i,y}=P_{z_{i,y}}.
\end{equation}
Thus a single fixed point satisfies
$E_i(z_{i,y})\le2E_i(x)$ at every active testing point of the patch,
including the support of every derivative of its cutoff.  No continuity
of $E_i$ or measurable selection is needed.  Restricting the minimization
to $T_{i,y}$ also keeps its center in the stated domain at the coarsest
scales.

The following is the covector field used in the iteration:
\begin{equation}
 \label{re:field-definition}
 v_i(P)=\sum_{y\in\mathcal Y_i}
           \phi_{i,y}(P)D_{i,y}(P-c_{i,y}).
\end{equation}
It is a smooth field on all of value space, since the sum is finite.

\begin{lemma}[Field estimates]
 \label{re:field-estimates}
At $P=P_x$, $x\in B_1(p)$, write $r=r_i$, $H=H_i(x)$,
$Q=Q_i(x)$, $Z=Z_i(x)$, $S=S_i(x)$, and $E=E_i(x)$.
There is a linear error $e(\xi)$ such that
\begin{align}
 r^{-1}|v_i|_Q&\le C\sqrt\varepsilon\sqrt E,
       \label{re:field-value}\\
 (\mathrm Dv_i)\xi&=S\xi+e(\xi),
       \label{re:field-linear}\\
 |S\xi|_Q&\le\sqrt\varepsilon\sqrt{Z[\xi]},
       &|e(\xi)|_Q&\le C\sqrt\varepsilon\sqrt E\,|\xi|_H,
       \label{re:field-small-linear}\\
 |\xi^Te(\xi)|&\le C\sqrt{Z[\xi]E}\,|\xi|_H,
       \label{re:field-diagonal}\\
 r|\mathrm D^2v_i(\xi,\xi)|_Q
       &\le C\left(\sqrt{Z[\xi]}\,|\xi|_H
                          +\sqrt E\,|\xi|_H^2\right).
       \label{re:field-second}
\end{align}
There is a fixed $C_v\ge1$ such that, in every order $m\ge0$,
\begin{equation}
 \label{re:field-jets}
 r_i^{m-1}\|\mathrm D^m v_i(P_x)\|_{H_i(x)\to Q_i(x)}
       \le\varepsilon C_v^{m+1}(m!)^s.
\end{equation}
\end{lemma}

\begin{proof}
All sums in this proof are over $\mathcal A_i(x)$.  Put
$w_y=P_x-c_{i,y}$.  Equation~\eqref{re:directional-input} and the
fixed choice~\eqref{re:fixed-center} give
\begin{equation}
 \label{re:center-energy}
 D_{i,y}[w_y]\le C r^2E.
\end{equation}
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
\begin{equation}
 \label{re:positive-matrix-bound}
 |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].
\end{equation}
Here the second inequality is Cauchy--Schwarz for a positive
semidefinite form and does not require its invertibility.  Applying
the first inequality to~\eqref{re:center-energy} and using the bounded
number of summands proves~\eqref{re:field-value}.

Differentiation of~\eqref{re:field-definition} gives
\[
 e(\xi)=\sum_y
          (\mathrm D\phi_{i,y})\xi\ D_{i,y}w_y.
\]
By normalized cutoff bounds and matrix comparability,
\begin{equation}
 \label{re:cutoff-derivative}
 |(\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.
\end{equation}
These inequalities and~\eqref{re:positive-matrix-bound} prove the
bound for $e$ in~\eqref{re:field-small-linear}.  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
\eqref{re:positive-matrix-bound} 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~\eqref{re:positive-matrix-bound}--\eqref{re:cutoff-derivative}
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~\eqref{re:field-second}.

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~\eqref{re:field-definition} has only the term with $m$
derivatives on the cutoff and the $m$ terms with one derivative on
the affine factor.  By~\eqref{re:bump-jets}, 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~\eqref{re:field-jets}.  The preceding
value bound with $|w_y|_H\le Cr$ proves its $m=0$ case.
\end{proof}

\subsection{The iteration and its derivative bounds}

On the common open set $\Omega$ define, as long as the Jacobians are
invertible there,
\begin{equation}
 \label{re:iteration}
 h_0(P)=P,
 \qquad J_i(P)=\mathrm Dh_i(P),
 \qquad h_{i+1}(P)=h_i(P)+J_i(P)^{-T}v_i(P).
\end{equation}
The inverse transpose is chosen to turn the leading derivative of
$v_i$ directly into a change of metric.  To see both the gain and the
error, put $W=J^{-T}v$ for a smooth map $h$ with $J=\mathrm Dh$.
Differentiating $J^TJ^{-T}=\Id$ gives the exact identity
\begin{equation}
 \label{re:metric-step}
 \begin{aligned}
 |(J+\mathrm DW)\xi|^2-|J\xi|^2
   ={}&2\xi^T\mathrm Dv\,\xi
       -2\langle\mathrm D^2h(\xi,\xi),J^{-T}v\rangle\\
     &+|\mathrm DW\,\xi|^2.
 \end{aligned}
\end{equation}
By the field estimates, $\mathrm Dv_i\,\xi=S_i\xi+e(\xi)$.
The first term therefore supplies the desired positive increment
$2S_i[\xi]$, up to the directional error.  The middle term explains
why we will estimate the diagonal derivative
$\mathrm D^2h_i(\xi,\xi)$ by a scale recurrence.  Bounds on all
higher derivatives serve a different purpose: they keep the smooth
update controlled despite its loss of one derivative.

We shall prove simultaneously a metric bound and all higher-derivative
bounds.  Choose $A\ge1$, subsequently fixed in terms of $C_*$, and
fix $\sigma\ge\max\{1,s\}$.  The desired bounds are
\begin{align}
 A^{-1}H_i(x)&\le J_i(P_x)^TJ_i(P_x)\le A H_i(x),
      \label{re:metric-bootstrap}\\
 r_i^{m-1}\|\mathrm D^m h_i(P_x)\|_{H_i(x)\to\R^n}
       &\le B^{m-1}(m!)^{1+\sigma}\qquad(m\ge2),
      \label{re:jet-bootstrap}
\end{align}
where $B$ will depend on the fixed data and $A,\sigma$, but not on
$\varepsilon,i,x,m$.

The next lemma shows that matrix inversion preserves the exponential
base $B$ in these derivative bounds.  A new factor growing exponentially
with the order at each step would prevent the simultaneous induction.

\begin{lemma}[Factorial convolution and inverse coefficients]
 \label{re:inverse-majorant}
Fix $\sigma\ge1$ and put $L_0=1$ and
\[
 L_m=(m+1)((m+1)!)^\sigma\qquad(m\ge1).
\]
Then
\begin{equation}
 \label{re:convolution-decay}
 \delta_m:=\sum_{l=1}^{m-1}\frac{L_lL_{m-l}}{L_m}
       \longrightarrow0,
 \qquad
 \sum_{l=0}^m L_lL_{m-l}\le C_\sigma L_m.
\end{equation}
Suppose a smooth matrix-valued function $J$ is invertible at a point
and, in Euclidean operator norms at that point,
\[
 \|J^{-1}\|\le\sqrt A,\qquad
 \frac{\|\mathrm D^mJ\|}{m!}\le B^mL_m\quad(m\ge1),
 \qquad A,B\ge1.
\]
There is a constant $C_{A,\sigma}$, independent of $m$ and $B$, such that
\begin{equation}
 \label{re:inverse-jets}
 \frac{\|\mathrm D^m(J^{-1})\|}{m!}
       \le C_{A,\sigma}B^mL_m\qquad(m\ge0).
\end{equation}
The assertion also holds for the inverse transpose.
\end{lemma}

\begin{proof}
For $1\le l\le m-1$ the summand in the first expression of
\eqref{re:convolution-decay} 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~\eqref{re:convolution-decay}.

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
\begin{equation}
 \label{re:inverse-recursion}
 a_m\le A L_m+\sqrt A
                    \sum_{l=1}^{m-1}L_l a_{m-l}.
\end{equation}
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~\eqref{re:inverse-recursion} 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~\eqref{re:inverse-jets} 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.
\end{proof}

\begin{lemma}[One step for all jets]
 \label{re:jet-step}
Choose $B\ge C_v$ sufficiently large and then $\varepsilon$ sufficiently
small, depending on $A,\sigma,B$ and the fixed data.  If
\eqref{re:metric-bootstrap} and~\eqref{re:jet-bootstrap} hold at level
$i$, the update~\eqref{re:iteration} satisfies
\eqref{re:jet-bootstrap} at level $i+1$.  This assertion does not
require the new-level metric bound.
\end{lemma}

\begin{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$,~\eqref{re:jet-bootstrap} gives
\[
 \frac{\|\mathrm D^m\widetilde J(0)\|}{m!}
 \le\frac{B^m((m+1)!)^{1+\sigma}}{m!}
 =B^m L_m.
\]
Lemma~\ref{re:inverse-majorant} therefore bounds the coefficients
of $\widetilde J^{-T}$ by $C_{A,\sigma}B^mL_m$.
On the other hand,~\eqref{re:field-jets} 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
\eqref{re:convolution-decay} now imply
\begin{equation}
 \label{re:increment-jets}
 \frac{\|\mathrm D^m(\widetilde J^{-T}\widetilde v)(0)\|}{m!}
       \le\varepsilon C_A B^m L_m.
\end{equation}

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~\eqref{re:increment-jets} 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.
\end{proof}

\subsection{Diagonal second derivatives and the metric bootstrap}

\begin{proposition}[Uniform metric and jet bounds]
 \label{re:metric-iteration}
There are fixed $A,B$ and a sufficiently small flatness threshold
such that the maps~\eqref{re:iteration} are smooth on all of $\Omega$
for every finite $i$, their Jacobians are invertible there, and
\eqref{re:metric-bootstrap}--\eqref{re:jet-bootstrap} hold for all
$x\in B_1(p)$ and every $i$.
\end{proposition}

\begin{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)},
       \label{re:first-inverse-cost}\\
 |(\mathrm D_\xi^2 K_i)w|
     &\le C_A\frac{q_i^2}{r_i^2}|w|_{Q_i(x)}.
       \label{re:second-inverse-cost}
\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
\eqref{re:first-inverse-cost}--\eqref{re:second-inverse-cost} 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:
\begin{equation}
 \label{re:second-increment-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.
\end{equation}
By Lemma~\ref{re:field-estimates} and
\eqref{re:first-inverse-cost}--\eqref{re:second-inverse-cost},
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~\eqref{re:field-second}; 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
\begin{equation}
 \label{re:diagonal-recurrence}
 b_{i+1}\le\tfrac12 b_i+C_A f_i,
 \qquad b_0=0.
\end{equation}
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)
       \notag\\
 &\le C_A H_k(x)[\xi].
       \label{re:diagonal-square-sum}
\end{align}
The last line uses~\eqref{re:cumulative-input}, monotonicity
$H_i\le H_k$, and~\eqref{re:pinching-input}.  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
\begin{equation}
 \label{re:squared-increment}
 |\Delta_iJ\,\xi|^2
 \le C_A\varepsilon
           \left(Z_i(x)[\xi]+E_i(x)q_i^2\right).
\end{equation}
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
       \notag\\
 &=-\langle\mathrm D^2h_i(\xi,\xi),K_iv_i\rangle_{\R^n}.
       \label{re:inverse-transpose-contraction}
\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~\eqref{re:field-value} 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~\eqref{re:field-diagonal},
\eqref{re:squared-increment}, and
\eqref{re:inverse-transpose-contraction}.  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
\begin{equation}
 \label{re:metric-error}
 \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].
\end{equation}
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~\eqref{re:cumulative-input} and $0\le S_i\le Z_i$,
\begin{equation}
 \label{re:leading-metric-comparison}
 C_*^{-1}H_k
       \le\Id+2\sum_{i<k}S_i
       \le2C_*H_k.
\end{equation}
Choose $A\ge4C_*$, increasing it if necessary for the initial
comparability of $H_0$ and $\Id$.  Fix $\sigma,B$ as in
Lemma~\ref{re:jet-step}.  Only now choose the data sufficiently flat
that its smallness requirement holds and
\[
 C_A(\sqrt\varepsilon+\varepsilon)<\frac1{2C_*}.
\]
Then~\eqref{re:metric-error} and
\eqref{re:leading-metric-comparison} imply
\eqref{re:metric-bootstrap} 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
\eqref{re:metric-error} 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~\ref{re:jet-step} 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.
\end{proof}

\begin{remark}
The construction neither differentiates the fields $H_i$ nor assumes
that they have a finite limit.  It also does not require the maps
$h_i$ themselves to be globally injective.  Their ordinary Jacobians
are invertible on the fixed open value-space domain, which is exactly
what the update needs.
\end{remark}

\subsection{Passage to an open bi-Lipschitz chart}

\begin{proposition}[The limiting chart]
 \label{re:limit-chart}
For sufficiently flat data in the preceding construction, there is
$\rho>0$ such that the uniform limit
\[
 F(x)=\lim_{i\to\infty}h_i(u(x)),\qquad x\in B_1(p),
\]
restricts to a bi-Lipschitz homeomorphism from $B_\rho(p)$ onto an
open subset of $\R^n$.  Its distance comparisons use the ambient
distance $d$ restricted to $B_\rho(p)$.
\end{proposition}

\begin{proof}
From~\eqref{re:field-value} and
\eqref{re:metric-bootstrap},
\[
 |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
\begin{equation}
 \label{re:uniform-tail}
 \sup_{x\in B_1(p)}|F(x)-h_i(P_x)|\le\tau r_i,
 \qquad \tau=C_A\varepsilon\longrightarrow0.
\end{equation}
Each approximating map is continuous, so $F$ is continuous.

Let $K_2=B(2!)^{1+\sigma}$, the fixed second-derivative constant in
\eqref{re:jet-bootstrap}.  By~\eqref{re:ellipsoid-input}, along
any value-space segment contained in the central ellipsoid at $x$,
with preimages in $B_{r_i}(x)\subset B_1(p)$,
\begin{equation}
 \label{re:taylor-constant}
 |\mathrm D^2h_i(P_y)(\zeta,\zeta)|
       \le\frac{K_2C_{\mathrm{cmp}}}{r_i}
                         |\zeta|_{H_i(x)}^2.
\end{equation}
Put $K_T=K_2C_{\mathrm{cmp}}/2$.  Choose a fixed $D\ge1$ so large
that
\begin{equation}
 \label{re:choice-D}
 D>4/a,\qquad 4K_T/D\le1/(8\sqrt A).
\end{equation}
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
\begin{equation}
 \label{re:pair-normalization}
 \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]
\end{equation}
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
\begin{equation}
 \label{re:choice-tail}
 4D\tau\le1/(8\sqrt A).
\end{equation}

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~\eqref{re:pair-normalization}
gives $t/2\le|\Delta|_{H_i(x)}\le2t$.
By~\eqref{re:choice-D}, $|\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
\eqref{re:ellipsoid-input}, 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~\eqref{re:taylor-constant}.

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~\eqref{re:uniform-tail} 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)$,
\begin{equation}
 \label{re:final-bilipschitz}
 \frac{1}{4\sqrt A}d(x,z)
       \le|F(x)-F(z)|\le4\sqrt A\,d(x,z).
\end{equation}
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\R^n
\]
therefore has open image by invariance of domain
\cite[Theorem~2B.3]{Hatcher02}. 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.
\end{proof}

\begin{proof}[Proof of Theorem~\ref{thm:main}]
Let $p$ be the prescribed $n$-regular point of the noncollapsed
$\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 \eqref{re:choice-D} consequently depend only
on $n$. Choose a rescaling sufficiently late that all
their flatness requirements, the simultaneous induction, and
\eqref{re:choice-tail} hold.  Proposition~\ref{re:limit-chart} 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
\eqref{re:final-bilipschitz} then scale by the same factor, while
the domain remains an open neighborhood and the image remains an
open subset of $\R^n$. Thus Theorem~\ref{thm:main} 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$.
\end{proof}
==== END Bi-Lipschitz-Coordinates-at-Regular-Points-of-Noncollapsed-RCD-Spaces//renormalization.tex ====

==== BEGIN Bi-Lipschitz-Coordinates-at-Regular-Points-of-Noncollapsed-RCD-Spaces//references.bib ====
@article{DPG18,
  author = {De Philippis, Guido and Gigli, Nicola},
  title = {Non-collapsed spaces with {Ricci} curvature bounded from below},
  journal = {Journal de l'{\'E}cole polytechnique --- Math{\'e}matiques},
  volume = {5},
  pages = {613--650},
  year = {2018},
  doi = {10.5802/jep.80},
  url = {https://www.numdam.org/articles/10.5802/jep.80/},
  eprint = {1708.02060v3},
  archivePrefix = {arXiv},
  note = {Revised version: \href{https://arxiv.org/abs/1708.02060v3}{arXiv:1708.02060v3}, 23 June 2025}
}

@article{Han18,
  author = {Han, Bang-Xian},
  title = {{Ricci} tensor on {$\mathrm{RCD}^{*}(K,N)$} spaces},
  journal = {The Journal of Geometric Analysis},
  volume = {28},
  number = {2},
  pages = {1295--1314},
  year = {2018},
  doi = {10.1007/s12220-017-9863-7},
  url = {https://doi.org/10.1007/s12220-017-9863-7},
  eprint = {1412.0441v4},
  archivePrefix = {arXiv},
  note = {\href{https://arxiv.org/abs/1412.0441v4}{arXiv:1412.0441v4}}
}

@article{Gigli18,
  author = {Gigli, Nicola},
  title = {Nonsmooth differential geometry---An approach tailored for spaces
           with {Ricci} curvature bounded from below},
  journal = {Memoirs of the American Mathematical Society},
  volume = {251},
  number = {1196},
  year = {2018},
  doi = {10.1090/memo/1196},
  url = {https://doi.org/10.1090/memo/1196},
  eprint = {1407.0809v1},
  archivePrefix = {arXiv},
  note = {\href{https://arxiv.org/abs/1407.0809v1}{arXiv:1407.0809v1}}
}

@article{Gigli15,
  author = {Gigli, Nicola},
  title = {On the differential structure of metric measure spaces and applications},
  journal = {Memoirs of the American Mathematical Society},
  volume = {236},
  number = {1113},
  year = {2015},
  doi = {10.1090/memo/1113},
  url = {https://doi.org/10.1090/memo/1113},
  eprint = {1205.6622v2},
  archivePrefix = {arXiv},
  note = {\href{https://arxiv.org/abs/1205.6622v2}{arXiv:1205.6622v2}}
}

@article{CM20,
  author = {Cavalletti, Fabio and Mondino, Andrea},
  title = {New formulas for the {Laplacian} of distance functions and applications},
  journal = {Analysis \& PDE},
  volume = {13},
  number = {7},
  pages = {2091--2147},
  year = {2020},
  doi = {10.2140/apde.2020.13.2091},
  url = {https://doi.org/10.2140/apde.2020.13.2091},
  eprint = {1803.09687},
  archivePrefix = {arXiv}
}

@article{AH18,
  author = {Ambrosio, Luigi and Honda, Shouhei},
  title = {Local spectral convergence in {$\mathrm{RCD}^{*}(K,N)$} spaces},
  journal = {Nonlinear Analysis},
  volume = {177},
  pages = {1--23},
  year = {2018},
  doi = {10.1016/j.na.2017.04.003},
  url = {https://doi.org/10.1016/j.na.2017.04.003},
  eprint = {1703.04939v1},
  archivePrefix = {arXiv},
  note = {\href{https://arxiv.org/abs/1703.04939v1}{arXiv:1703.04939v1}}
}

@article{JLZ16,
  author = {Jiang, Renjin and Li, Huaiqian and Zhang, Huichun},
  title = {Heat kernel bounds on metric measure spaces and some applications},
  journal = {Potential Analysis},
  volume = {44},
  number = {3},
  pages = {601--627},
  year = {2016},
  doi = {10.1007/s11118-015-9521-2},
  url = {https://doi.org/10.1007/s11118-015-9521-2},
  eprint = {1407.5289v2},
  archivePrefix = {arXiv},
  note = {\href{https://arxiv.org/abs/1407.5289v2}{arXiv:1407.5289v2}}
}

@article{CDNZPSW,
  author = {Connell, Chris and Dai, Xianzhe and
            N{\'u}{\~n}ez-Zimbr{\'o}n, Jes{\'u}s and Perales, Raquel and
            Su{\'a}rez-Serrato, Pablo and Wei, Guofang},
  title = {Maximal volume entropy rigidity for
           {$\mathsf{RCD}^{*}(-(N-1),N)$} spaces},
  journal = {Journal of the London Mathematical Society},
  volume = {104},
  number = {4},
  pages = {1615--1681},
  year = {2021},
  doi = {10.1112/jlms.12470},
  url = {https://doi.org/10.1112/jlms.12470},
  eprint = {1809.06909v2},
  archivePrefix = {arXiv},
  note = {\href{https://arxiv.org/abs/1809.06909v2}{arXiv:1809.06909v2}}
}

@misc{Kell13,
  author = {Kell, Martin},
  title = {A note on {Lipschitz} continuity of solutions of {Poisson}
           equations in metric measure spaces},
  year = {2013},
  eprint = {1307.2224v2},
  archivePrefix = {arXiv},
  primaryClass = {math.MG},
  doi = {10.48550/arXiv.1307.2224},
  url = {https://arxiv.org/abs/1307.2224},
  note = {\href{https://arxiv.org/abs/1307.2224v2}{arXiv:1307.2224v2},
          17 July 2013}
}

@article{Jiang14,
  author = {Jiang, Renjin},
  title = {{Cheeger}-harmonic functions in metric measure spaces revisited},
  journal = {Journal of Functional Analysis},
  volume = {266},
  number = {3},
  pages = {1373--1394},
  year = {2014},
  doi = {10.1016/j.jfa.2013.11.022},
  url = {https://doi.org/10.1016/j.jfa.2013.11.022},
  eprint = {1307.1334v2},
  archivePrefix = {arXiv},
  note = {\href{https://arxiv.org/abs/1307.1334v2}{arXiv:1307.1334v2}}
}

@article{BS10,
  author = {Bacher, Kathrin and Sturm, Karl-Theodor},
  title = {Localization and tensorization properties of the
           curvature-dimension condition for metric measure spaces},
  journal = {Journal of Functional Analysis},
  volume = {259},
  number = {1},
  pages = {28--56},
  year = {2010},
  doi = {10.1016/j.jfa.2010.03.024},
  url = {https://doi.org/10.1016/j.jfa.2010.03.024},
  eprint = {1003.2116v1},
  archivePrefix = {arXiv},
  note = {\href{https://arxiv.org/abs/1003.2116v1}{arXiv:1003.2116v1}}
}

@article{CJN21,
  author = {Cheeger, Jeff and Jiang, Wenshuai and Naber, Aaron},
  title = {Rectifiability of singular sets of noncollapsed limit spaces
           with {Ricci} curvature bounded below},
  journal = {Annals of Mathematics},
  volume = {193},
  number = {2},
  pages = {407--538},
  year = {2021},
  doi = {10.4007/annals.2021.193.2.2},
  url = {https://annals.math.princeton.edu/2021/193-2/p02},
  eprint = {1805.07988v1},
  archivePrefix = {arXiv},
  note = {\href{https://arxiv.org/abs/1805.07988v1}{arXiv:1805.07988v1}}
}

@article{CN13,
  author = {Colding, Tobias Holck and Naber, Aaron},
  title = {Lower {Ricci} curvature, branching and the bilipschitz structure
           of uniform {Reifenberg} spaces},
  journal = {Advances in Mathematics},
  volume = {249},
  pages = {348--358},
  year = {2013},
  doi = {10.1016/j.aim.2013.09.005},
  url = {https://doi.org/10.1016/j.aim.2013.09.005},
  eprint = {1111.2184v1},
  archivePrefix = {arXiv},
  note = {\href{https://arxiv.org/abs/1111.2184v1}{arXiv:1111.2184v1}}
}

@article{KM21,
  author = {Kapovitch, Vitali and Mondino, Andrea},
  title = {On the topology and the boundary of {$N$}-dimensional
           {$\mathrm{RCD}(K,N)$} spaces},
  journal = {Geometry \& Topology},
  volume = {25},
  number = {1},
  pages = {445--495},
  year = {2021},
  doi = {10.2140/gt.2021.25.445},
  url = {https://doi.org/10.2140/gt.2021.25.445},
  eprint = {1907.02614v3},
  archivePrefix = {arXiv},
  note = {\href{https://arxiv.org/abs/1907.02614v3}{arXiv:1907.02614v3}}
}

@article{HP23,
  author = {Honda, Shouhei and Peng, Yuanlin},
  title = {A note on the topological stability theorem from
           {$\mathrm{RCD}$} spaces to {Riemannian} manifolds},
  journal = {Manuscripta Mathematica},
  volume = {172},
  number = {3--4},
  pages = {971--1007},
  year = {2023},
  doi = {10.1007/s00229-022-01418-7},
  url = {https://doi.org/10.1007/s00229-022-01418-7},
  eprint = {2202.06500v2},
  archivePrefix = {arXiv},
  note = {\href{https://arxiv.org/abs/2202.06500v2}{arXiv:2202.06500v2}}
}

@misc{HZ26,
  author = {Honda, Shouhei and Zhang, Ruobing},
  title = {Regularity and structure of spaces with synthetic {Ricci}
           bounds and positive injectivity radius},
  year = {2026},
  eprint = {2608.16021v1},
  archivePrefix = {arXiv},
  primaryClass = {math.DG},
  doi = {10.48550/arXiv.2608.16021},
  url = {https://arxiv.org/abs/2608.16021},
  note = {\href{https://arxiv.org/abs/2608.16021v1}{arXiv:2608.16021v1},
          17 August 2026}
}

@article{CC96,
  author = {Cheeger, Jeff and Colding, Tobias H.},
  title = {Lower bounds on {Ricci} curvature and the almost rigidity of warped products},
  journal = {Annals of Mathematics},
  volume = {144},
  number = {1},
  pages = {189--237},
  year = {1996},
  doi = {10.2307/2118589},
  note = {\href{https://doi.org/10.2307/2118589}{doi:10.2307/2118589}}
}

@article{CC97,
  author = {Cheeger, Jeff and Colding, Tobias H.},
  title = {On the structure of spaces with {Ricci} curvature bounded below. {I}},
  journal = {Journal of Differential Geometry},
  volume = {46},
  number = {3},
  pages = {406--480},
  year = {1997},
  doi = {10.4310/jdg/1214459974},
  note = {\href{https://doi.org/10.4310/jdg/1214459974}{doi:10.4310/jdg/1214459974}}
}

@article{LV09,
  author = {Lott, John and Villani, Cedric},
  title = {{Ricci} curvature for metric-measure spaces via optimal transport},
  journal = {Annals of Mathematics},
  volume = {169},
  number = {3},
  pages = {903--991},
  year = {2009},
  doi = {10.4007/annals.2009.169.903},
  note = {\href{https://doi.org/10.4007/annals.2009.169.903}{doi:10.4007/annals.2009.169.903}}
}

@article{Sturm06II,
  author = {Sturm, Karl-Theodor},
  title = {On the geometry of metric measure spaces. {II}},
  journal = {Acta Mathematica},
  volume = {196},
  pages = {133--177},
  year = {2006},
  doi = {10.1007/s11511-006-0003-7},
  note = {\href{https://doi.org/10.1007/s11511-006-0003-7}{doi:10.1007/s11511-006-0003-7}}
}

@article{AGS14,
  author = {Ambrosio, Luigi and Gigli, Nicola and Savar{\'e}, Giuseppe},
  title = {Metric measure spaces with {Riemannian Ricci} curvature bounded from below},
  journal = {Duke Mathematical Journal},
  volume = {163},
  number = {7},
  pages = {1405--1490},
  year = {2014},
  doi = {10.1215/00127094-2681605},
  note = {\href{https://doi.org/10.1215/00127094-2681605}{doi:10.1215/00127094-2681605}}
}

@article{EKS15,
  author = {Erbar, Matthias and Kuwada, Kazumasa and Sturm, Karl-Theodor},
  title = {On the equivalence of the entropic curvature-dimension condition and
           {Bochner}'s inequality on metric measure spaces},
  journal = {Inventiones Mathematicae},
  volume = {201},
  number = {3},
  pages = {993--1071},
  year = {2015},
  doi = {10.1007/s00222-014-0563-7},
  note = {\href{https://doi.org/10.1007/s00222-014-0563-7}{doi:10.1007/s00222-014-0563-7}}
}

@article{MN19,
  author = {Mondino, Andrea and Naber, Aaron},
  title = {Structure theory of metric measure spaces with lower {Ricci} curvature bounds},
  journal = {Journal of the European Mathematical Society},
  volume = {21},
  number = {6},
  pages = {1809--1854},
  year = {2019},
  doi = {10.4171/JEMS/874},
  eprint = {1405.2222v3},
  archivePrefix = {arXiv},
  note = {\href{https://arxiv.org/abs/1405.2222v3}{arXiv:1405.2222v3}}
}

@article{CN15,
  author = {Cheeger, Jeff and Naber, Aaron},
  title = {Regularity of {Einstein} manifolds and the codimension 4 conjecture},
  journal = {Annals of Mathematics},
  volume = {182},
  number = {3},
  pages = {1093--1165},
  year = {2015},
  doi = {10.4007/annals.2015.182.3.5},
  note = {\href{https://doi.org/10.4007/annals.2015.182.3.5}{doi:10.4007/annals.2015.182.3.5}}
}

@article{HH24,
  author = {Huang, Hongzhi and Huang, Xian-Tao},
  title = {Almost splitting maps, transformation theorems and smooth fibration theorems},
  journal = {Advances in Mathematics},
  volume = {457},
  pages = {109914},
  year = {2024},
  doi = {10.1016/j.aim.2024.109914},
  eprint = {2207.10029v3},
  archivePrefix = {arXiv},
  note = {\href{https://arxiv.org/abs/2207.10029v3}{arXiv:2207.10029v3}}
}

@article{GV25,
  author = {Gigli, Nicola and Violo, Ivan Yuri},
  title = {Notes on the {Cheeger} and {Colding} version of the {Reifenberg}
           theorem for metric spaces},
  journal = {Bollettino dell'Unione Matematica Italiana},
  volume = {18},
  number = {1},
  pages = {217--271},
  year = {2025},
  doi = {10.1007/s40574-025-00458-7},
  eprint = {2406.13771v1},
  archivePrefix = {arXiv},
  note = {\href{https://doi.org/10.1007/s40574-025-00458-7}{doi:10.1007/s40574-025-00458-7}}
}

@article{DT12,
  author = {David, Guy and Toro, Tatiana},
  title = {{Reifenberg} parameterizations for sets with holes},
  journal = {Memoirs of the American Mathematical Society},
  volume = {215},
  number = {1012},
  year = {2012},
  eprint = {0910.4869v1},
  archivePrefix = {arXiv},
  note = {\href{https://arxiv.org/abs/0910.4869v1}{arXiv:0910.4869v1}}
}

@book{Hatcher02,
  author = {Hatcher, Allen},
  title = {Algebraic Topology},
  publisher = {Cambridge University Press},
  year = {2002},
  url = {https://pi.math.cornell.edu/~hatcher/AT/AT.pdf},
  note = {\href{https://pi.math.cornell.edu/~hatcher/AT/AT.pdf}{Author's electronic edition}}
}

@article{HK00,
  author = {Haj{\l}asz, Piotr and Koskela, Pekka},
  title = {Sobolev met {Poincar\'e}},
  journal = {Memoirs of the American Mathematical Society},
  volume = {145},
  number = {688},
  year = {2000},
  doi = {10.1090/memo/0688},
  note = {\href{https://doi.org/10.1090/memo/0688}{doi:10.1090/memo/0688}}
}

@misc{BNS20,
  author = {Bru{\`e}, Elia and Naber, Aaron and Semola, Daniele},
  title = {Boundary Regularity and Stability for Spaces with {Ricci} Bounded Below},
  year = {2020},
  eprint = {2011.08383v1},
  archivePrefix = {arXiv},
  url = {https://arxiv.org/abs/2011.08383v1},
  note = {\href{https://arxiv.org/abs/2011.08383v1}{arXiv:2011.08383v1}, 17 November 2020; locators refer to this version}
}
==== END Bi-Lipschitz-Coordinates-at-Regular-Points-of-Noncollapsed-RCD-Spaces//references.bib ====