Lemma B.2. Recall the eigenmatrix $V\in\mathbb{R}^{p\times p}$ defined in
(2.1). Let

$\mathbf{V}:=\begin{pmatrix}V&0\\\
0&I_{n}\end{pmatrix}\in\mathbb{R}^{(n+p)\times(n+p)}.$

We then rotate $\mathsf{H}$ and $\mathsf{G}$ in (B.2) as

$H(z):=\mathbf{V}^{\top}\mathsf{H}(z)\mathbf{V},\quad G(z):=H(z)^{-1}.$ (B.13)

Under this definition, with (B.3), we have that

$H=\begin{pmatrix}-\Lambda_{0}^{-1}&V^{\top}X\\\
X^{\top}V&-z\end{pmatrix},\quad
G=\begin{pmatrix}z\Lambda_{0}^{1/2}\mathcal{G}_{1}\Lambda_{0}^{1/2}&\Lambda_{0}V^{\top}X\mathcal{G}_{2}\\\
\mathcal{G}_{2}X^{\top}V\Lambda_{0}&\mathcal{G}_{2}\end{pmatrix}.$ (B.14)

where the resolvents $\mathcal{G}_{1}$ and $\mathcal{G}_{2}$ are defined as

$\mathcal{G}_{1}(z):=(\Lambda_{0}^{1/2}V^{\top}XX^{\top}V\Lambda_{0}^{1/2}-z)^{-1},\quad\mathcal{G}_{2}(z):=\mathsf{G}_{2}(z).$
(B.15)

In this case, it is easy to see that Lemmas B.2 and B.3 hold for $G(z)$ if we
replace $\Pi_{\mathsf{G}}(z)$ by

$\Pi(z)=\begin{pmatrix}-\Lambda_{0}(1+m(z)\Lambda_{0})^{-1}&0\\\
0&m(z)\end{pmatrix}.$ (B.16)

The main advantage of $\Pi$ over $\Pi_{\mathsf{G}}$ is that $\Pi$ is a
diagonal matrix. Similarly, we can define $\widetilde{\mathcal{G}}_{1}$,
$\widetilde{\mathcal{G}}_{2}$, $\widetilde{G}$, and $\widetilde{\Pi}$ by
replacing $\Lambda_{0}$ with $\Lambda$ in the above definitions.

The following lemma outlines some basic estimates regarding $m(z)$ and
$\operatorname{Im}m(z)$. Notably, it shows that the entries of $\Pi$ in (B.16)
are of order 1. These estimates will be used tacitly in the subsequent proofs.

######  Lemma B.5.

Suppose items (i) and (iii) of Assumption 2.3 hold. Then, for any small
constant $c>0$, the following estimates hold uniformly for $z$ satisfying
$c\leq|z|\leq c^{-1}$:

$|m(z)|\asymp
1,\quad\operatorname{Im}m(z)\asymp\begin{cases}\sqrt{\kappa+\eta},\
&\text{if}\ E\in\operatorname{supp}(\varrho)\\\
\frac{\eta}{\sqrt{\kappa+\eta}},\ &\text{if}\
E\notin\operatorname{supp}(\varrho)\end{cases},$ (B.17)
$\min_{j=1}^{p}|1+m(z)\sigma_{j}|\gtrsim 1.$ (B.18)

For sufficiently small constant $c>0$, the following estimate holds uniformly
for $z$ satisfying $|z|\leq c$:

$|zm(z)|\asymp\begin{cases}1,&\ \text{if}\ c_{n}\leq 1-\tau\\\ |z|,&\
\text{if}\ c_{n}\geq 1+\tau\end{cases}.$ (B.19)

###### Proof.

The estimates in (B.17) are proved in Lemma 4.10 and equation (A.7) of [45].
The estimate (B.18) is proved in Lemma A.6 of [28]. Finally, the estimate
(B.19) follows from the fact that the deformed MP law has support within
$[\lambda_{-},\lambda_{+}]$ with a delta mass $(1-c_{n})\delta_{0}$ when
$c_{n}<1$. ∎

As an important consequence of the averaged local laws in Lemma B.2, we have
the following rigidity estimate of the eigenvalues
$\lambda_{j}\equiv\lambda_{j}(\mathcal{Q}_{2})$.

######  Lemma B.6 (Theorem 3.12 of [45]).

Under the assumptions of Lemma B.2, for $\gamma_{j}\in[a_{2k},a_{2k-1}]$, we
have

$|\lambda_{j}-\gamma_{j}|\prec[(n_{k}+1-j)\wedge(j-n_{k-1})]^{-1/3}n^{-{2}/{3}},$
(B.20)

where we recall the notations in (2.11) and Definition 2.2.

###  B.2 Eigenvalues and eigenvectors of the spiked covariance model

In this subsection, we present a collection of useful results regarding the
eigenvalues and eigenvectors of the spiked model $\widetilde{\mathcal{Q}}_{1}$
in (2.4). Recall that $\widetilde{\lambda}_{i}$ and $\lambda_{i}$ represent
the eigenvalues of $\widetilde{\mathcal{Q}}_{1}$ and $\mathcal{Q}_{1}$,
respectively, while $\mathfrak{a}_{i}$ and $\mathfrak{b}_{i}$ are defined in
(2.20). Furthermore, we emphasize that the ${\bm{u}_{i}}$ mentioned in Lemmas
B.7 and B.8 below refers to the eigenvectors of $\widetilde{\mathcal{Q}}_{1}$,
rather than $\mathcal{Q}_{1}$. First, Lemma B.7 gives the first-order limits
of the outlier eigenvalues and eigenvectors, along with nearly optimal
convergence rates.

######  Lemma B.7 (Theorems 3.2 and 3.3 of [27]).

Suppose Assumption 2.3 holds. For any $i\in[\\![{r}]\\!]$, we have that

$\left|\widetilde{\lambda}_{i}-\mathfrak{a}_{i}\right|\prec n^{-1/2}.$ (B.21)

Moreover, for the corresponding outlier eigenvectors, we have that for any
$i,j\in[\\![{r}]\\!]$,

$\left|\langle\bm{u}_{i},\bm{v}_{j}\rangle\right|^{2}-\frac{\mathfrak{b}_{i}}{\widetilde{\sigma}_{i}}\delta_{ij}\prec
n^{-1/2}.$ (B.22)

Next, we state some estimates for the non-outlier eigenvalues and
eigenvectors.

######  Lemma B.8.

Suppose Assumption 2.3 holds. For $i\in[\\![{\mathsf{K}-r}]\\!]$, we have the
eigenvalue sticking estimate:

$|\widetilde{\lambda}_{i+r}-\lambda_{i}|\prec n^{-1},$ (B.23)

and the eigenvector delocalization estimate:

$\left|\langle\bm{u}_{i+r},{\bf{v}}\rangle\right|^{2}\prec n^{-1}$ (B.24)

for any deterministic unit vector ${\bf{v}}\in\mathbb{R}^{p}$. When $c_{n}>1$,
we have that

$\|U_{0}^{\top}{\bf{v}}\|_{2}^{2}\prec 1.$ (B.25)

Finally, for $i\in[\\![{r}]\\!]$ and $j\in[\\![{r+1,p}]\\!],$ we have that

$\left|\langle\bm{u}_{i},\bm{v}_{j}\rangle\right|^{2}-\frac{\mathfrak{c}_{i}}{\sigma_{j}}G_{ij}(\mathfrak{a}_{i})G_{ji}(\mathfrak{a}_{i})\prec
n^{-3/2},\quad\text{with}\quad\mathfrak{c}_{i}:=\frac{d_{i}^{2}\mathfrak{b}_{i}}{\widetilde{\sigma}_{i}^{2}}.$
(B.26)

###### Proof.

The results (B.23) and (B.24) have been proved in Theorems 3.7 and 3.14 of
[29] for $i\leqslant\tau p$ with $\tau\in(0,1)$ being a small constant. The
reason for considering a subset of indices is that only the edge regularity
condition around $\lambda_{+}=a_{1}$ was assumed in [29]. However, with the
stronger regularity conditions in (2.16) and (2.17), we can extend (B.23) and
(B.24) to all $i\in[\\![{\mathsf{K}-r}]\\!]$ following the proof in [29]. We
omit the details. As for (B.25), we have the trivial bound

${\|U_{0}^{\top}{\bf{v}}\|_{2}^{2}}\leq\eta\operatorname{Im}[{\bf{v}}^{\top}\widetilde{\mathsf{G}}_{1}(\mathrm{i}\eta){\bf{v}}],$

which follows from the spectral decomposition of $\widetilde{\mathsf{G}}_{1}$.
The RHS can be written as

${\bf{v}}^{\top}\widetilde{\mathsf{G}}_{1}(\mathrm{i}\eta){\bf{v}}=\left({\bf{v}}^{\top}\Sigma^{-1/2},\mathbf{0}_{n}^{\top}\right)\frac{\widetilde{\mathsf{G}}(\mathrm{i}\eta)}{\mathrm{i}\eta}\begin{pmatrix}\Sigma^{-1/2}{\bf{v}}\\\
\mathbf{0}_{n}\end{pmatrix}.$

Then, (B.25) follows immediately from the local law on
$\widetilde{\mathsf{G}}$ in (B.12) by taking $\eta=(\log n)^{-1}$.

It remains to prove (B.26). The proof relies on a representation of
$\left|\langle\bm{u}_{i},\bm{v}_{j}\rangle\right|^{2}$ in terms of contour
integrals of the resolvent $G$ (see e.g., [13, 18, 27, 29]). We only sketch
the proof since many arguments below are similar to those in [18, 27]. Let
$\\{{\bf{u}}_{i}\\}_{i=1}^{p}=\\{V^{\top}\bm{u}_{i}\\}_{i=1}^{p}$ be the left
singular eigenvectors of $\Lambda^{1/2}V^{\top}X$. Due to the relation

$\langle\bm{v}_{j},\bm{u}_{i}\rangle=(\bm{v}_{j}^{\top}V){\bf{u}}_{i}={\bf{e}}_{j}^{\top}{\bf{u}}_{i},$
(B.27)

it suffices to focus on the entries ${\bf{u}}_{i}(j)$ of ${\bf{u}}_{i}$ in the
following proof. Let $\Gamma_{i}$ denote the boundary of
$B_{\rho}(\mathfrak{a}_{i})$, which represents a disk of radius $\rho$
centered at $\mathfrak{a}_{i}$. Here, we choose $\rho$ to be sufficiently
small so that $B_{2\rho}(\mathfrak{a}_{i})$ does not contain $\lambda_{+}$ or
any other outlier $\mathfrak{a}_{k}$, $k\in[\\![{r}]\\!]\setminus\\{i\\}$.
Then, by (B.21) and item (iv) of Assumption 2.3, we know that with high
probability, $\Gamma_{i}$ only encloses $\widetilde{\lambda}_{i}$ and no other
eigenvalues of $G(z)$. By applying Cauchy’s integral formula, we obtain that

$|{\bf{u}}_{i}(j)|^{2}=-\frac{1}{2\pi\mathrm{i}}\oint_{\Gamma_{i}}(\widetilde{\mathcal{G}}_{1})_{jj}(z)\mathrm{d}z=-\lim_{\delta\downarrow
0}\frac{1}{2\pi\mathrm{i}}\oint_{\Gamma_{i}}(\widetilde{\mathcal{G}}_{1}^{\delta})_{jj}(z)\mathrm{d}z,$
(B.28)

where $\widetilde{\mathcal{G}}_{1}^{\delta}$ is defined by introducing a
perturbation to the $j$-th eigenvalue of $\Lambda$:

$\widetilde{\mathcal{G}}_{1}^{\delta}(z):=\big{(}\Lambda_{\delta}^{1/2}V^{\top}XX^{\top}V\Lambda_{\delta}^{1/2}-z\big{)}^{-1},\quad\text{with}\quad\Lambda_{\delta}:=\Lambda+\delta\widetilde{\sigma}_{j}\mathbf{e}_{j}\mathbf{e}_{j}^{\top}.$

To compute the right-hand side (RHS) of (B.28), we utilize the Woodbury matrix
identity to express $\widetilde{\mathcal{G}}^{\delta}_{1}$ in terms of
$\mathcal{G}_{1}$. Using the fact that $\mathcal{G}_{1}$ contains no pole
inside $\Gamma_{i}$, we then obtain the following expression:

$\langle{\bf{u}}_{i},{\bf{e}}_{j}\rangle^{2}=\lim_{\delta\downarrow
0}\frac{1+\delta}{\delta^{2}}\frac{1}{2\pi\mathrm{i}}\oint_{\Gamma_{i}}\left[\mathbf{D}_{r,j}(\delta)^{-1}+1+z\mathbf{E}_{r,j}^{\top}\mathcal{G}_{1}(z)\mathbf{E}_{r,j}\right]_{r+1,r+1}^{-1}\frac{\mathrm{d}z}{z}.$
(B.29)

where we introduce the matrices

$\mathbf{E}_{r,j}:=({\bf{e}}_{1},\cdots,{\bf{e}}_{r},{\bf{e}}_{j})\in\mathbb{R}^{p\times(r+1)},\quad\mathbf{D}_{r,j}(\delta):=\operatorname{diag}\left\\{d_{1},d_{2},\cdots,d_{r},\delta\right\\}\in\mathbb{R}^{(r+1)\times(r+1)}.$

For a detailed derivation of (B.29), we refer readers to Lemma 5.7 and
equations (7.5) and (7.6) of [27].

Denote
$\Lambda_{0,j}:=\operatorname{diag}\\{\sigma_{1},\cdots,\sigma_{r},\sigma_{j}\\}$,
and

$\Delta_{j}(z):=-\Upsilon_{j}(z)-z\mathbf{E}_{r,j}^{\top}\mathcal{G}_{1}(z)\mathbf{E}_{r,j},\quad\Upsilon_{j}(z):=\left(1+m(z)\Lambda_{0,j}\right)^{-1}.$
(B.30)

By the local law (B.10) and the estimate on $\operatorname{Im}m(z)$ in (B.17),
we have that uniformly in $z\in\Gamma_{i}$,

$\left\|\Delta_{j}(z)\right\|\prec n^{-1/2}.$ (B.31)

Denote $\Xi_{j}(z):=\mathbf{D}_{r,j}^{-1}+1-\Upsilon_{j}(z).$ By item (iv) of
2.3 and the choice of $\rho$, we have that $\|\Xi_{j}(z)^{-1}\|\lesssim 1$
uniformly in $z\in\Gamma_{i}$. Then, applying Taylor expansion to

$\left(\mathbf{D}_{r,j}^{-1}+1+z\mathbf{E}_{r,j}^{\top}\mathcal{G}_{1}(z)\mathbf{E}_{r,j}\right)^{-1}=\left(\Xi_{j}(z)-\Delta_{j}(z)\right)^{-1}$

till order three, we can expand (B.29) as

$\langle{\bf{u}}_{i},{\bf{e}}_{j}\rangle^{2}=s_{1}+s_{2}+s_{3}+s_{4},$

where $s_{i}$, $i=1,2,3,4,$ are defined as

$\displaystyle s_{1}:=\lim_{\delta\downarrow
0}\frac{1+\delta}{\delta^{2}}\frac{1}{2\pi\mathrm{i}}\oint_{\Gamma_{i}}\frac{1}{1+\delta^{-1}-(1+m(z)\sigma_{j})^{-1}}\frac{\mathrm{d}z}{z},$
$\displaystyle s_{2}:=\lim_{\delta\downarrow
0}\frac{1+\delta}{\delta^{2}}\frac{1}{2\pi\mathrm{i}}\oint_{\Gamma_{i}}\frac{\left(\Delta(z)\right)_{r+1,r+1}}{[1+\delta^{-1}-(1+m(z)\sigma_{j})^{-1}]^{2}}\frac{\mathrm{d}z}{z},$
$\displaystyle s_{3}:=\lim_{\delta\downarrow
0}\frac{1+\delta}{\delta^{2}}\frac{1}{2\pi\mathrm{i}}\oint_{\Gamma_{i}}\left[\Xi_{j}^{-1}(z)\Delta(z)\Xi_{j}^{-1}(z)\Delta(z)\Xi_{j}^{-1}(z)\right]_{r+1,r+1}\frac{\mathrm{d}z}{z},$
$\displaystyle s_{4}:=\lim_{\delta\downarrow
0}\frac{1+\delta}{\delta^{2}}\frac{1}{2\pi\mathrm{i}}\oint_{\Gamma_{i}}\left[\left(\Xi_{j}^{-1}(z)\Delta(z)\right)^{3}\left(\Xi_{j}(z)-\Delta(z)\right)^{-1}\right]_{r+1,r+1}\frac{\mathrm{d}z}{z}.$

Using (B.31), it is easy to see that

$|s_{4}|\prec n^{-3/2}.$

For $s_{1}$ and $s_{2},$ by following a similar argument as in [27,
Proposition 7.2] or [18, Lemma 4.31], we can show that for sufficiently small
$\delta,$ the functions $z^{-1}$,
$[1+\delta^{-1}-(1+m(z)\sigma_{j})^{-1}]^{-1}$, and $(\Delta(z))_{r+1,r+1}$
are all holomorphic on and inside $\Gamma_{i}$ with high probability. Hence,
applying Cauchy’s integral theorem, we get that

$s_{1}=s_{2}=0\quad\text{with high probability}.$ (B.32)

To estimate the leading term $s_{3}$, we employ a modified proof strategy
similar to that of [18, Lemma 4.32]. By definition, we can write that

$\displaystyle\left[\Xi_{j}^{-1}(z)\Delta(z)\Xi_{j}^{-1}(z)\Delta(z)\Xi_{j}^{-1}(z)\right]_{r+1,r+1}=\frac{1}{\left[1+\delta^{-1}-(1+m(z)\sigma_{j})^{-1}\right]^{2}}\sum_{k=1}^{r+1}\Delta_{r+1,k}\Delta_{k,r+1}\Xi_{j}(k,k)^{-1},$

where we have introduced the abbreviation $\Delta\equiv\Delta(z)$ and used
$\Xi_{j}(k,k)$ to denote the $k$-th diagonal entry of $\Xi_{j}(z)$. We
decompose the sum into a diagonal part with $k=r+1$ and an off-diagonal part
with $k\in[\\![{r}]\\!]$:

$\sum_{k=1}^{r+1}\Delta_{r+1,k}\Delta_{k,r+1}\Xi_{j}(k,k)^{-1}=\mathsf{E}_{d}+\mathsf{E}_{o},$

where $\mathsf{E}_{d}$ and $\mathsf{E}_{o}$ are defined as

$\displaystyle\mathsf{E}_{d}:=\Delta^{2}_{r+1,r+1}\Xi_{j}(r+1,r+1),\quad\mathsf{E}_{o}:=z^{2}\sum_{k=1}^{r}\frac{({\bf{e}}_{j}^{\top}\mathcal{G}_{1}(z){\bf{e}}_{k})({\bf{e}}_{k}^{\top}\mathcal{G}_{1}(z){\bf{e}}_{j})}{1+d_{k}^{-1}-(1+m(z)\sigma_{k})^{-1}}.$

By the discussion above (B.32), we have that for sufficiently small $\delta$,
$\mathsf{E}_{d}$ is holomorphic on and inside $\Gamma_{i}$ with high
probability. This implies that

$\frac{1}{2\pi\mathrm{i}}\oint_{\Gamma_{i}}\frac{\mathsf{E}_{d}}{\left[1+\delta^{-1}-(1+m(z)\sigma_{j})^{-1}\right]^{2}}\frac{\mathrm{d}z}{z}=0\quad\text{with
high probability}.$

For the off-diagonal part $\mathsf{E}_{o}$, using (B.14), we can rewrite it as

$\mathsf{E}_{o}=\frac{1}{\sigma_{j}}\sum_{k=1}^{r}\frac{1}{\sigma_{k}}\frac{G_{jk}(z)G_{kj}(z)}{1+d_{k}^{-1}-(1+m(z)\sigma_{k})^{-1}}.$

Under item (iv) of Assumption 2.3, using the definition of $m(z)$ in (2.8) and
the definition of $\mathfrak{a}_{i}$ in (2.20), we find that for
$k\in[\\![{r}]\\!]\setminus\\{i\\}$,

$1+d_{k}^{-1}-(1+m(\mathfrak{a}_{i})\sigma_{k})^{-1}=1+d_{k}^{-1}-\left(1-{\sigma_{k}}/\widetilde{\sigma}_{i}\right)^{-1}\gtrsim
1.$

Consequently, due to the continuity of $m(z)$ near $\mathfrak{a}_{i}$, we can
choose $\rho$ sufficiently small such that
$1+d_{k}^{-1}-(1+m(z)\sigma_{k})^{-1}$ has no zero on and inside $\Gamma_{i}$.
Then, by Cauchy’s integral theorem, we have that for
$k\in[\\![{r}]\\!]\setminus\\{i\\}$,

$\oint_{\Gamma_{i}}\frac{G_{jk}(z)G_{kj}(z)}{\left[1+\delta^{-1}-(1+m(z)\sigma_{j})^{-1}\right]^{2}\left[1+d_{k}^{-1}-(1+m(z)\sigma_{k})^{-1}\right]}\frac{\mathrm{d}z}{z}=0$

with high probability. Finally, we only need to consider the $k=i$ term in
$\mathsf{E}_{o}$. Applying the parametrization $z=h(\zeta),$ we obtain that
with high probability,

$\displaystyle\frac{1+\delta}{\delta^{2}}\frac{1}{2\pi\mathrm{i}}\oint_{\Gamma_{i}}\frac{\mathsf{E}_{o}}{\left[1+\delta^{-1}-(1+m(z)\sigma_{j})^{-1}\right]^{2}}\frac{\mathrm{d}z}{z}$
$\displaystyle=\frac{(1+\delta)d_{i}}{2\pi\mathrm{i}\widetilde{\sigma}_{i}\sigma_{j}\sigma_{i}}\oint_{\gamma_{i}}\frac{h^{\prime}(\zeta)G_{ji}(h(\zeta))G_{ij}(h(\zeta))}{h(\zeta)\left[1+\delta-\delta(1+\zeta\sigma_{j})^{-1}\right]^{2}}\frac{1+\zeta\sigma_{i}}{\zeta+\widetilde{\sigma}_{i}^{-1}}\mathrm{d}\zeta,$
$\displaystyle=(1+\delta)\frac{d_{i}^{2}\mathfrak{b}_{i}}{\widetilde{\sigma}_{i}^{2}\sigma_{j}}\frac{G_{ji}(\mathfrak{a}_{i})G_{ij}(\mathfrak{a}_{i})}{\left[1+\delta-\delta(1-\sigma_{j}/\widetilde{\sigma}_{i})^{-1}\right]^{2}},$

where we used the relation $m(h(\zeta))=\zeta$ as in (2.8) in the first step,
and the residual theorem at $\zeta=-\widetilde{\sigma}_{i}^{-1}$ in the second
step. Finally, taking $\delta\to 0$ and using the dominated convergence
theorem, we conclude that

$s_{3}=\frac{d_{i}^{2}\mathfrak{b}_{i}}{\widetilde{\sigma}_{i}^{2}\sigma_{j}}G_{ji}(\mathfrak{a}_{i})G_{ij}(\mathfrak{a}_{i}),$

which completes the proof of (B.26). ∎

###  B.3 Local laws for the rectangular matrix Dyson Brownian motion

In this subsection, we consider the rectangular matrix DBM of the non-spiked
model:

$W_{0}(t):=W_{0}+\sqrt{t}X^{G},\quad\text{with}\quad
W_{0}:=\Lambda_{0}^{1/2}V^{\top}X,$

where $X^{G}$ is defined above (4.3). In this subsection, we present the local
laws for its resolvent, which will be utilized in our proof. From these local
laws, we can also derive the local laws for the resolvent of the (spiked)
rectangular matrix DBM defined in (4.3) and (4.4).

Let us begin by introducing some notations. We define
$\Lambda_{t}:=\Lambda_{0}+t$ and $m_{t}(z):=m(z,\Lambda_{t})$ as given by
(2.8). In other words, $m_{t}(z)$ is the unique solution to the equation

$z=-\frac{1}{m_{t}}+\frac{1}{n}\operatorname{Tr}\frac{\Lambda_{t}}{1+m_{t}\Lambda_{t}}$
(B.33)

subject to the condition that $\operatorname{Im}m_{t}(z)\geq 0$ whenever
$\operatorname{Im}z\geq 0$. Subsequently, we denote the density function
associated with $m_{t}$ as $\varrho_{t}$, and refer to the spectral edges and
quantiles of $\varrho_{t}$ as $a_{k}(t)$ and $\gamma_{k}(t)$, following the
definitions in Lemma 2.1 and in Definition 2.2. Using the stability of the
self-consistent equation (B.33) established in [45], it is not hard to check
that

$\sup_{x\in\mathbb{R}}|\varrho_{t}(x)-\varrho(x)|=\operatorname{o}(1)$ (B.34)

when $t=\operatorname{o}(1)$. Then, we define $\mathcal{G}_{1,t}$,
$\mathcal{G}_{2,t}$, $G_{t}$, and $\Pi_{t}$ as described in Remark B.4, with
the replacements of $W_{0}$, $\Lambda_{0}$, and $m(z)$ by $W_{0}(t)$,
$\Lambda_{t}$, and $m_{t}(z)$, respectively. Let
$\\{\lambda_{k}(t)\\}_{k=1}^{p}$ and $\\{{\bf{u}}_{k}(t)\\}_{k=1}^{p}$ denote
the eigenvalues and corresponding eigenvectors of $W_{0}(t)W_{0}(t)^{\top}$.
Finally, we introduce the error control parameter

$\Psi_{t}(z):=\sqrt{\frac{\operatorname{Im}m_{t}(z)}{n\eta}}+\frac{1}{n\eta},$

Now, we are ready to state the main result of this subsection.

######  Lemma B.9.

Suppose (i)–(iii) of Assumption 2.3 hold. Fix a time $T\asymp
n^{-1/3+\mathsf{c}}$ for some constant $\mathsf{c}\in(0,1/6).$ The following
estimates hold uniformly in $z\in\mathbf{D}$ and $t\in[0,T]$.

  1. (i)

The support of $\varrho_{t}$ has the same number of connected components as
$\varrho_{t}$, and the spectral edges $a_{k}(t)$ of $\varrho_{t}$ are
perturbations of those of $\varrho$:
$\max_{k=1}^{2q}|a_{k}(t)-a_{k}|=\operatorname{o}(1).$ In particular, it
implies that $|\lambda_{\pm,t}-\lambda_{\pm}|=\operatorname{o}(1)$, where
$\lambda_{+,t}:=a_{1}(t)$ and $\lambda_{-,t}:=a_{2q}(t)$ denote the rightmost
and leftmost spectral edges of $\varrho_{t}$, respectively. The number of
quantiles within each component remains unchanged, that is, $n_{k}(t)$ defined
as in (2.11) is constant in $t$.

  2. (ii)

For any deterministic unit vectors ${\bf{u}},{\bf{v}}\in\mathbb{R}^{p+n}$, we
have the anisotropic local law (recall Lemma B.2)

$\left|{\bf{u}}^{\top}{G}_{t}(z){\bf{v}}-{\bf{u}}^{\top}\Pi_{t}(z){\bf{v}}\right|\prec\Psi_{t}(z).$
(B.35)

For $\mathsf{m}_{n,t}:=n^{-1}\operatorname{Tr}\mathcal{G}_{2,t}$, the
following averaged local law holds:

$|\mathsf{m}_{n,t}(z)-m_{t}(z)|\prec\frac{1}{n\eta}.$ (B.36)

  3. (iii)

For $\gamma_{j}(t)\in[a_{2k}(t),a_{2k-1}(t)]$, we have that (recall Lemma B.6)

$\left|\lambda_{j}(t)-\gamma_{j}(t)\right|\prec[(n_{k}(t)+1-j)\wedge(j-n_{k-1}(t))]^{-1/3}n^{-{2}/{3}}.$
(B.37)

  4. (iv)

For $k\in[\\![{\mathsf{K}}]\\!]$ and any deterministic unit vector
${\bf{v}}\in\mathbb{R}^{p}$, we have that (recall Lemma B.8)

$\left|\langle{\bf{u}}_{k}(t),{\bf{v}}\rangle\right|^{2}\prec n^{-1}.$ (B.38)

###### Proof.

Under the edge regularity condition (2.16), it is easy to see that $a_{k}(t)$
and $n_{k}(t)$ are continuous in $t$ when $t=\operatorname{o}(1)$. This
establishes part (i) of the lemma, taking into account the fact that
$n_{k}(t)$ takes integer values. For the other parts, we notice that
$W_{0}(t)$ can be expressed as

$W_{0}(t)=\Lambda_{t}^{1/2}\mathbf{U}\begin{pmatrix}X\\\ X^{G}\end{pmatrix},$
(B.39)

where $\mathbf{U}\in\mathbb{R}^{p\times 2p}$ is a semi-orthogonal matrix
defined as

$\mathbf{U}:=\begin{pmatrix}(\Lambda_{0}/\Lambda_{t})^{1/2}V^{\top},(t/\Lambda_{t})^{1/2}I_{p}\end{pmatrix}.$

This is the random matrix model studied in [45]. Furthermore, (B.34) and part
(i) together validate the regularity conditions in (2.16) and (2.17) for
$\Lambda_{t}$. Then, for each fixed $t$, (ii) is proved in Theorems 3.6 and
3.7 of [45], (iii) is proved in Theorem 3.12 of [45], and (iv) follows easily
from (ii) and (iii) as explained, for example, in [12, Theorem 2.8]. Finally,
by utilizing a standard $\varepsilon$-net argument, we can extend these
results uniformly to all $t\in[0,T]$. ∎

######  Remark B.10.

The estimate (B.36) provides an averaged local law for $\mathsf{m}_{n,t}.$ In
the following proof, we will also need an averaged local law for
$\underline{\mathsf{m}}_{n,t}(z):=p^{-1}\operatorname{Tr}\mathcal{G}_{1,t}$,
whose classical limit is given by

$m_{1,t}(z):=-\frac{1}{p}\operatorname{Tr}\frac{1}{z(1+m_{t}(z)\Lambda_{t})}.$

From (B.33), we can observe that

$m_{1,t}(z)=\frac{c_{n}^{-1}-1}{z}+c_{n}^{-1}m_{t}(z).$ (B.40)

On the other hand, since $W_{0}(t)W_{0}(t)^{\top}$ share the non-zero same
eigenvalues as $W_{0}(t)^{\top}W_{0}(t)$, we have a similar relation between
$\underline{\mathsf{m}}_{n,t}(z)$ and $\mathsf{m}_{n,t}(z)$:

$\underline{\mathsf{m}}_{n,t}(z)=\frac{c_{n}^{-1}-1}{z}+c_{n}^{-1}\mathsf{m}_{n,t}(z).$
(B.41)

Under these two equations, we can directly deduce from (B.36) that

$|\underline{\mathsf{m}}_{n,t}(z)-m_{1,t}(z)|\prec(n\eta)^{-1}.$ (B.42)

Next, we consider the rectangular matrix DBM of the spiked model:

$W(t):=W+\sqrt{t}X^{G},\quad\text{with}\quad W:=\Lambda^{1/2}V^{\top}X,$

where $\Lambda$ contains a non-spiked component $\Lambda_{0}$ and $r$ spikes
as in (2.3). Let $\widetilde{\Lambda}_{t}:=\Lambda+t$, and we again denote
$\Lambda_{t}=\Lambda_{0}+t$ and $m_{t}(z)=m(z,\Lambda_{t})$. Then, we define
$\widetilde{G}_{t}$ and $\widetilde{\Pi}_{t}$ as described in Remark B.4, with
the replacements of $W_{0}$, $\Lambda_{0}$, and $m(z)$ by $W(t)$,
$\widetilde{\Lambda}_{t}$, and $m_{t}(z)$, respectively. Let
$\\{\widetilde{\lambda}_{k}(t)\\}_{k=1}^{p}$ denote the eigenvalues of
$W(t)W(t)^{\top}$. We define a smaller spectral parameter domain than (B.7)
for arbitrarily small constants $c,\tau\in(0,1)$:

$\mathbf{D}^{\prime}\equiv\mathbf{D}^{\prime}(c,\tau,t,n):=\left\\{z=E+\mathrm{i}\eta\in\mathbb{C}_{+}:\lambda_{-}-c\leqslant
E\leqslant\lambda_{+}+c,n^{-1+\tau}\leqslant\eta\leqslant\tau^{-1}\right\\},$
(B.43)

We show that a similar local law as in (B.35) holds for
$\widetilde{G}_{t}-\widetilde{\Pi}_{t}$.

######  Lemma B.11.

Suppose Assumption 2.3 holds. Fix a time $T\asymp n^{-1/3+\mathsf{c}}$ for
some constant $\mathsf{c}\in(0,1/6).$ Then, there exists a constant $c>0$ such
that the following local law holds uniformly for $z\in\mathbf{D}^{\prime}$ and
$t\in[0,T]$, considering any deterministic unit vectors
${\bf{u}},{\bf{v}}\in\mathbb{R}^{p+n}$:

$\big{|}{\bf{u}}^{\top}\widetilde{G}_{t}(z){\bf{v}}-{\bf{u}}^{\top}\widetilde{\Pi}_{t}(z){\bf{v}}\big{|}\prec\Psi_{t}(z).$
(B.44)

Moreover, for $i\in[\\![{\mathsf{K}-r}]\\!]$, we have the eigenvalue sticking
estimate:

$|\widetilde{\lambda}_{i+r}(t)-\lambda_{i}(t)|\prec n^{-1}.$ (B.45)

###### Proof.

For the proof of (B.44), the only concern arises for $z$ around the $r$
outlier locations. However, by confining ourselves to a narrower spectral
domain $\mathbf{D}^{\prime}$ around $\operatorname{supp}(\varrho_{t})$, we
ensure that $z$ stays at a distance of order $\Omega(1)$ from the outliers. As
a result, (B.44) is essentially a consequence of the local law (B.35). The
proof of (B.45) is the same as that of (B.23). We omit the details. ∎

##  Appendix C Proof of the main results

###  C.1 Proof of the results in Section 2.2

In this subsection, we provide the proofs of Theorem 2.4, Corollary 2.5, and
Corollary 2.6 using the local laws presented in Appendix B and Theorem 2.12.

###### Proof of Theorem 2.4.

Applying Theorem 2.12 with $\omega_{j}=\ell(\sigma_{j})$ and noticing
$p^{-1}\sum_{j=1}^{r}\ell(\sigma_{j})\phi(\bm{v}_{j},\bm{v}_{j},\gamma_{i-r})=\operatorname{O}(p^{-1}),$
we can deduce that

$\mathbb{P}\left(\left|\bm{u}_{i}^{\top}\ell(\Sigma)\bm{u}_{i}-\vartheta(\ell,\gamma_{i-r})\right|>n^{-\mathfrak{d}/4}\right)\leq
n^{-\mathfrak{d}/2},\quad i\in[\\![{r+1,\mathsf{K}}]\\!],$

for a constant $\mathfrak{d}>0$. By combining this estimate with the trivial
bound
$\left|\bm{u}_{i}^{\top}\ell(\Sigma)\bm{u}_{i}-\vartheta(\ell,\gamma_{i-r})\right|\lesssim
1$, we conclude (2.30). To show (2.31), let $\Gamma_{0}\equiv\Gamma_{0}(n)$ be
the contour centered at $0$ with radius $\rho_{n}=(\log n)^{-1}/2$. By the
rigidity estimate (B.20), the eigenvalue sticking (B.23), and the lower bound
(2.19), we know that with high probability, $\Gamma_{0}$ only encloses the
$p-n$ zero eigenvalues of $\mathcal{Q}_{1}$. Then, with the spectral
decomposition of $\widetilde{\mathsf{G}}_{1}$ (recall (B.1)), we obtain that
with high probability,

$\displaystyle\operatorname{Tr}\left[U_{0}^{\top}\ell(\Sigma)U_{0}\right]=-\frac{1}{2\pi\mathrm{i}}\oint_{\Gamma_{0}}\operatorname{Tr}\big{[}\widetilde{\mathsf{G}}_{1}(z)\ell(\Sigma)\big{]}\mathrm{d}z.$

Applying the anisotropic law in Lemma B.3 and using Cauchy’s integral theorem
at $z=0$, we obtain that

$\displaystyle\operatorname{Tr}\left[U_{0}^{\top}\ell(\Sigma)U_{0}\right]$
$\displaystyle=\frac{1}{2\pi\mathrm{i}}\oint_{\Gamma_{0}}\operatorname{Tr}\left[(1+m(z)\Sigma)^{-1}\ell(\Sigma)\right]\frac{\mathrm{d}z}{z}+\operatorname{O}_{\prec}(n^{-1/2})$
$\displaystyle=\sum_{j=1}^{p}\frac{\ell(\widetilde{\sigma}_{j})}{1+m(0)\widetilde{\sigma}_{j}}+\operatorname{O}_{\prec}(n^{-1/2}).$

This yields (2.31) by noticing that

$\displaystyle\vartheta(\ell,0)=\frac{1}{p-n}\sum_{j=r+1}^{p}\frac{\ell(\widetilde{\sigma}_{j})}{1+m(0)\widetilde{\sigma}_{j}}=\frac{1}{p-n}\sum_{j=1}^{p}\frac{\ell(\widetilde{\sigma}_{j})}{1+m(0)\widetilde{\sigma}_{j}}+\operatorname{O}(n^{-1}).$

Now, we focus on the proof of (2.29). Recall the decomposition (4.2):

$\displaystyle\bm{u}_{i}^{\top}\ell(\Sigma)\bm{u}_{i}$
$\displaystyle=\sum_{j=1}^{r}\ell(\widetilde{\sigma}_{j})|\langle\bm{u}_{i},\bm{v}_{j}\rangle|^{2}+\sum_{j=r+1}^{p}\ell(\widetilde{\sigma}_{j})|\langle\bm{u}_{i},\bm{v}_{j}\rangle|^{2}:=\mathsf{w}_{1i}+\mathsf{w}_{2i},\quad
i\in[\\![{r}]\\!].$ (C.1)

Using (B.22), we obtain that

$\mathsf{w}_{1i}=\frac{\ell(\widetilde{\sigma}_{i})}{\widetilde{\sigma}_{i}}\mathfrak{b}_{i}+\operatorname{O}_{\prec}(n^{-1/2}).$
(C.2)

The remaining part of the proof is devoted to estimating $\mathsf{w}_{2i}.$
Since $\widetilde{\sigma}_{j}=\sigma_{j}$ for $j\geqslant r$, by (B.26), we
have

$\displaystyle\mathsf{w}_{2i}$
$\displaystyle=\mathfrak{c}_{i}\sum_{j=r+1}^{p}\frac{\ell(\sigma_{j})}{\sigma_{j}}G_{ij}(\mathfrak{a}_{i})G_{ji}(\mathfrak{a}_{i})+\operatorname{O}_{\prec}(n^{-1/2})$
$\displaystyle=\mathfrak{c}_{i}\sum_{j=1}^{p}\frac{\ell(\sigma_{j})}{\sigma_{j}}G_{ij}(\mathfrak{a}_{i})G_{ji}(\mathfrak{a}_{i})-\mathfrak{c}_{i}\sum_{j=1}^{r}\frac{\ell(\sigma_{j})}{\sigma_{j}}G_{ij}(\mathfrak{a}_{i})G_{ji}(\mathfrak{a}_{i})+\operatorname{O}_{\prec}(n^{-1/2}).$

Applying (B.10) to estimate the second sum on the RHS and recalling the
notations in Remark B.4, we can rewrite the above equation as

$\mathsf{w}_{2i}=\mathfrak{c}_{i}\mathfrak{a}_{i}^{2}\sigma_{i}\sum_{j=1}^{p}\ell(\sigma_{j})(\mathcal{G}_{1})_{ij}(\mathfrak{a}_{i})(\mathcal{G}_{1})_{ji}(\mathfrak{a}_{i})-\frac{\mathfrak{c}_{i}\sigma_{i}\ell(\sigma_{i})}{[1+m(\mathfrak{a}_{i})\sigma_{i}]^{2}}+\operatorname{O}_{\prec}(n^{-1/2}).$
(C.3)

It remains to estimate the quantity

$\mathcal{L}_{i}=\sum_{j=1}^{p}\ell(\sigma_{j})(\mathcal{G}_{1})_{ij}(\mathfrak{a}_{i})(\mathcal{G}_{1})_{ji}(\mathfrak{a}_{i}).$
(C.4)

For a small $\mathsf{t}>0,$ we define a new resolvent:

$\mathcal{G}(\mathsf{t},z):=\left(\Lambda_{0}^{1/2}V^{\top}XX^{\top}V\Lambda_{0}^{1/2}-z-\mathsf{t}L\right)^{-1},$
(C.5)

where $L$ is a diagonal matrix defined as
$L=\operatorname{diag}\\{\ell(\sigma_{1}),\cdots,\ell(\sigma_{p})\\}.$ Under
this definition, we notice that

$\mathcal{L}_{i}=\left.\frac{\partial\mathcal{G}_{ii}(\mathsf{t},\mathfrak{a}_{i})}{\partial\mathsf{t}}\right|_{\mathsf{t}=0}.$
(C.6)

To estimate (C.6), it suffices to establish a local law on
$\mathcal{G}_{ii}(\mathsf{t},\mathfrak{a}_{i})$ for any small $\mathsf{t}>0$.
We rewrite $\mathcal{G}(\mathsf{t},\mathfrak{a}_{i})$ as

$\displaystyle\mathcal{G}(\mathsf{t},\mathfrak{a}_{i})$
$\displaystyle=A(\mathsf{t})^{-1/2}\mathfrak{G}(\mathsf{t},\mathfrak{a}_{i})A(\mathsf{t})^{-1/2},$

where $A(\mathsf{t}):=I+\mathsf{t}L/\mathfrak{a}_{i}$ and

$\mathfrak{G}(\mathsf{t},z):=\left(\Lambda(\mathsf{t})^{1/2}V^{\top}XX^{\top}V\Lambda(\mathsf{t})^{1/2}-z\right)^{-1}\quad\text{with}\quad\Lambda(\mathsf{t}):=A(\mathsf{t})^{-1}\Lambda_{0}.$

Now, define $\widetilde{m}_{\mathsf{t}}(z,\mathfrak{a}_{i})\equiv
m(z,\Lambda(\mathsf{t}))$ as in (2.22) with $x=\mathfrak{a}_{i}$. As
$\mathsf{t}\to 0$, we can utilize the local law (B.10) to
$\mathfrak{G}(\mathsf{t},z)$ at $z=\mathfrak{a}_{i}$, and substitute
$\Lambda_{0}$ and $m(z)$ with $\Lambda(\sf{t})$ and
$\widetilde{m}_{\mathsf{t}}(z,\mathfrak{a}_{i})$, respectively. This leads to
that

$\displaystyle\mathcal{G}_{ii}(\mathsf{t},\mathfrak{a}_{i})$
$\displaystyle=-\frac{A_{ii}(\mathsf{t})^{-1}}{\mathfrak{a}_{i}[1+\widetilde{m}_{\mathsf{t}}(\mathfrak{a}_{i},\mathfrak{a}_{i})\Lambda_{ii}(\mathsf{t})]}+\operatorname{O}_{\prec}(n^{-1/2})$
(C.7)

for all $\mathsf{t}=\operatorname{o}(1)$. We now take $\mathsf{t}=n^{-1/4}$.
Using the eigenvalue rigidity (B.20) for $\lambda_{1}$ and the fact that
$\mathfrak{a}_{i}$ is an outlier of $\varrho$, we get from the definition
(C.5) that $\|\mathcal{G}(\mathsf{t},\mathfrak{a}_{i})\|\lesssim 1$ with high
probability. With this bound, we readily derive that with high probability,

$\mathcal{L}_{i}-\frac{\mathcal{G}_{ii}(\mathsf{t},\mathfrak{a}_{i})-\mathcal{G}_{ii}(0,\mathfrak{a}_{i})}{\mathsf{t}}=\operatorname{O}(\mathsf{t}).$
(C.8)

Combining (C.6), (C.7) and (C.8), we obtain that

$\displaystyle\mathcal{L}_{i}$
$\displaystyle=\left.-\frac{\partial}{\partial\mathsf{t}}\frac{A_{ii}(\mathsf{t})^{-1}}{\mathfrak{a}_{i}[1+\widetilde{m}_{\mathsf{t}}(\mathfrak{a}_{i},\mathfrak{a}_{i})\Lambda_{ii}(\mathsf{t})]}\right|_{\mathsf{t}=0}+\operatorname{O}_{\prec}(\mathsf{t}+n^{-1/2}/\mathsf{t})$
$\displaystyle=\frac{1}{[1+m(\mathfrak{a}_{i})\sigma_{i}]^{2}}\left(\frac{\dot{m}_{0}(\mathfrak{a}_{i})\sigma_{i}}{\mathfrak{a}_{i}}+\frac{\ell(\sigma_{i})}{\mathfrak{a}_{i}^{2}}\right)+\operatorname{O}_{\prec}(n^{-1/4}).$

Plugging it back into (C.3) gives an estimate for $\mathsf{w}_{2i}$. Combining
this estimate with (C.2), we obtain that

$\displaystyle\bm{u}_{i}^{\top}\ell(\Sigma)\bm{u}_{i}=\mathfrak{b}_{i}\left(\frac{\ell(\widetilde{\sigma}_{i})}{\widetilde{\sigma}_{i}}+\mathfrak{a}_{i}\dot{m}_{0}(\mathfrak{a}_{i})\right)+\operatorname{O}_{\prec}(n^{-1/4}),$

where we used that
$\mathfrak{c}_{i}={d_{i}^{2}\mathfrak{b}_{i}}/{\widetilde{\sigma}_{i}^{2}}$
and $m(\mathfrak{a}_{i})=-\widetilde{\sigma}_{i}^{-1}$. This concludes (2.29).
∎

###### Proof of Corollary 2.5.

By Theorem 2.4, it suffices to simplify the formulas of
$\vartheta(x)\equiv\vartheta(\ell,x)$ and $\psi_{i}\equiv\psi_{i}(\ell)$ when
$\ell(x)=x$. When $x=0$, we notice that

$\vartheta(0)=\frac{1}{p(1-c_{n}^{-1})}\sum_{j=r+1}^{p}\frac{\widetilde{\sigma}_{j}}{1+m(0)\widetilde{\sigma}_{j}}=\frac{1}{c_{n}-1}\frac{1}{n}\sum_{j=1}^{p}\frac{\sigma_{j}}{1+m(0)\sigma_{j}}+\operatorname{O}(n^{-1}).$

By (2.8), we have the identity

$\frac{1}{m(0)}=\frac{1}{n}\sum_{j=1}^{p}\frac{\sigma_{j}}{1+m(0)\sigma_{j}}.$

This proves (2.34). Next, when $i\in[\\![{r+1,\mathsf{K}}]\\!]$, we have
$\gamma_{i-r}>0$ and

$\vartheta(\gamma_{i-r})=\frac{c_{n}}{p\gamma_{i-r}}\sum_{j=r+1}^{p}\frac{\widetilde{\sigma}^{2}_{j}}{|1+m(\gamma_{i-r})\widetilde{\sigma}_{j}|^{2}}=\frac{1}{n\gamma_{i-r}}\sum_{j=1}^{p}\frac{\sigma^{2}_{j}}{|1+m(\gamma_{i-r})\sigma_{j}|^{2}}+\operatorname{O}(n^{-1}).$
(C.9)

By taking the imaginary parts of both sides of (2.8), we obtain that

$\frac{1}{n}\sum_{j=1}^{p}\frac{\sigma^{2}_{j}}{|1+m(\gamma_{i-r})\sigma_{j}|^{2}}=\frac{1}{|m(\gamma_{i-r})|^{2}}.$

Plugging it into (C.9) concludes the proof of (2.33). It remains to prove
(2.32). We write $\psi_{i}$ as

$\psi_{i}=\mathfrak{b}_{i}\big{(}1+\mathfrak{a}_{i}\dot{m}_{0}(\mathfrak{a}_{i})\big{)}.$
(C.10)

By setting $\ell(\sigma_{i})=\sigma_{i}$ and $z=x=\mathfrak{a}_{i}$, we can
simplify (2.25) using the first equation in (2.24), which yields that

$\dot{m}_{0}(\mathfrak{a}_{i})=\frac{1}{\mathfrak{a}_{i}}\left(\frac{m^{\prime}(\mathfrak{a}_{i})}{|m(\mathfrak{a}_{i})|^{2}}-1\right).$

Substituting this expression into (C.10) concludes (2.32). ∎

###### Proof of Corollary 2.6.

This corollary follows directly from Lemma A.1, Lemma A.2, and Theorem 2.4. ∎

###  C.2 Proof of the results in Section 2.3

First, we provide the proof of Theorem 2.8 by using Lemma 4.1 and Lemma 4.2.

###### Proof of Theorem 2.8.

Recall that for $D_{0}=\Lambda-\mathfrak{t}$ and $D_{t}=D_{0}+t$, we choose
their non-spiked parts as $D_{00}=\Lambda_{0}-\mathfrak{t}$ and
$D_{t0}=D_{00}+t$, respectively. Then, we define $m_{t}=m(z,D_{t0})$ as in
(B.33) and denote the density function associated with $m_{t}$ as
$\varrho_{t}$. The quantiles $\gamma_{k}(t)$ of $\varrho_{t}$ are defined as
in Definition 2.2. Given that $\mathfrak{t}=\operatorname{o}(1)$, it is
straightforward to show that the eigenvalues of $D_{00}$ and $D_{0}$ satisfy
(iii) and (iv) in 2.3. (In particular, the regularity conditions in (2.16) and
(2.17) can be verified using (B.34) and part (i) of Lemma B.9.) Then,
combining Lemma 4.1 and Lemma 4.2, we obtain that

$\sqrt{p}\begin{pmatrix}\xi_{1}\langle\bm{v},\bm{u}_{i_{1}+r}\rangle\\\
\vdots\\\
\xi_{L}\langle\bm{v},\bm{u}_{i_{L}+r}\rangle\end{pmatrix}\simeq\mathcal{N}(\bm{0},\Xi_{L}(\mathfrak{t})),$
(C.11)

where $\Xi_{L}(\mathfrak{t})$ is defined in (4.6) with
${\bf{v}}=V^{\top}\bm{v}$ (recall that ${\bf{u}}_{i}=V^{\top}\bm{u}_{i}$).
Since $D_{\mathfrak{t}}=\Lambda$, by comparing equations (B.33) and (2.8), we
see that $\varrho_{\mathfrak{t}}=\varrho$. Hence, we have
$\gamma_{i_{k}}(\mathfrak{t})=\gamma_{i_{k}}$ for $k\in[\\![{L}]\\!]$.
Moreover, using the stability of the self-consistent equation (2.8) (see
Appendix A of [45]), it is easy to check that

$m(x,D_{00})=m(x)+\operatorname{O}(\sqrt{\mathfrak{t}}).$ (C.12)

Then, using the definition (4.5), we conclude that

$\|\Xi_{L}(\mathfrak{t})-\Xi_{L}\|\lesssim\sqrt{\mathfrak{t}}.$

This establishes (2.37) together with (C.11). ∎

Recall that the proof of Lemma 4.1 is based on Theorem 4.5 and Lemma 4.6.
Using these two results, we can show that the LHS of (4.7) converge to the RHS
also in the sense of convergence of moments. Then, by applying Lemma 4.2 and
following the above proof of Theorem 2.8, we can deduce that (2.37) holds in
terms of matching moments, that is, for any fixed
$\mathbf{k}=(k_{1},\ldots,k_{L})\in\mathbb{N}^{L}$, there exists a constant
$c>0$ such that

$\mathbb{E}\prod_{j=1}^{L}\left(p|\langle\bm{v},\bm{u}_{i_{j}+r}\rangle|^{2}\right)^{k_{j}}=\prod_{j=1}^{L}\mathbb{E}\left(|\mathcal{N}_{j}|^{2k_{j}}\right)+\operatorname{O}(n^{-c}).$
(C.13)

Now, we proceed to establish Theorem 2.12 by utilizing this estimate. Our
approach follows a similar argument employed for Wigner matrices in [10, 15].

###### Proof of Theorem 2.12.

Due to (B.27), it suffices to prove the following estimate for ${\bf{u}}_{i}$:

$\mathbb{P}\bigg{(}\Big{|}\sum_{j=1}^{p}w_{j}\left|{\bf{u}}_{i}(j)\right|^{2}-p^{-1}\sum_{j=1}^{p}w_{j}\phi_{j}(\gamma_{i-r})\Big{|}>{\varepsilon}\bigg{)}\leqslant
n^{-\mathfrak{d}}/{\varepsilon}^{2},$ (C.14)

where the function $\phi_{j}(x)$ is defined as

$\phi_{j}(x):=\phi(\bm{v}_{j},\bm{v}_{j},x)=\frac{c_{n}\widetilde{\sigma}_{j}}{x|1+m(x)\widetilde{\sigma}_{j}|^{2}},\quad
x>0.$

Applying Markov’s inequality, we can derive (C.14) from the following
estimate:

$p^{-2}\mathbb{E}\Big{|}\sum_{j=1}^{p}w_{j}\left(p\left|{\bf{u}}_{i}(j)\right|^{2}-\phi_{j}(\gamma_{i-r})\right)\Big{|}^{2}\leq
n^{-\mathfrak{d}}.$ (C.15)

We observe that the LHS of (C.15) can be bounded by

$\max_{j\neq
k\in[\\![{p}]\\!]}\mathbb{E}\left(p\left|{\bf{u}}_{i}(j)\right|^{2}-\phi_{j}(\gamma_{i-r})\right)\left(p\left|{\bf{u}}_{i}(k)\right|^{2}-\phi_{k}(\gamma_{i-r})\right)+p^{-1}\max_{j=1}^{p}\mathbb{E}\left(p\left|{\bf{u}}_{i}(j)\right|^{2}-\phi_{j}(\gamma_{i-r})\right)^{2}.$
(C.16)

Using (C.13) with $\bm{v}=\bm{v}_{j}$ and $k_{j}\in\\{1,2\\}$, we can conclude
that

$p\mathbb{E}\left|{\bf{u}}_{i}(j)\right|^{2}=\phi_{j}(\gamma_{i-r})+\operatorname{O}(n^{-c}),\quad
p^{2}\mathbb{E}\left|{\bf{u}}_{i}(j)\right|^{4}=3\phi_{j}(\gamma_{i-r})^{2}+\operatorname{O}(n^{-c}).$
(C.17)

With these two estimates, we can bound the second term in (C.16) as

$p^{-1}\mathbb{E}\left(p\left|{\bf{u}}_{i}(j)\right|^{2}-\phi_{j}(\gamma_{i-r})\right)^{2}=\operatorname{O}(p^{-1}),\quad
j\in[\\![{p}]\\!].$ (C.18)

It remains to control the first term in (C.16). We now apply (C.13) with
$L=1$, $i_{1}=i-r$, $k_{1}=2$, and
$\bm{v}=(\xi_{1}\bm{v}_{j}+\xi_{2}\bm{v}_{k})/\sqrt{2}$, where $\xi_{1}$ and
$\xi_{2}$ are two uniformly random signs independent of $X$. This gives that

$\displaystyle
p^{2}\mathbb{E}\left|{\bf{u}}_{i}(j)\right|^{4}+p^{2}\mathbb{E}\left|{\bf{u}}_{i}(k)\right|^{4}+6p^{2}\mathbb{E}\left(\left|{\bf{u}}_{i}(j)\right|^{2}\left|{\bf{u}}_{i}(k)\right|^{2}\right)=3\left(\phi_{j}(\gamma_{i-r})^{2}+\phi_{k}(\gamma_{i-r})^{2}+2\phi_{j}(\gamma_{i-r})\phi_{k}(\gamma_{i-r})\right).$

Together with (C.17), this equation implies that

$p^{2}\mathbb{E}\left(\left|{\bf{u}}_{i}(j)\right|^{2}\left|{\bf{u}}_{i}(k)\right|^{2}\right)=\phi_{j}(\gamma_{i-r})\phi_{k}(\gamma_{i-r})+\operatorname{O}(n^{-c}).$
(C.19)

Combining (C.17) and (C.19), we conclude that

$\mathbb{E}\left(p\left|{\bf{u}}_{i}(j)\right|^{2}-\phi_{j}(\gamma_{i-r})\right)\left(p\left|{\bf{u}}_{i}(k)\right|^{2}-\phi_{k}(\gamma_{i-r})\right)\lesssim
n^{-c}.$ (C.20)

Plugging (C.18) and (C.20) into (C.16), we obtain (C.15), which further
concludes (C.14). ∎

###  C.3 Proof of the results in Section 3

###### Proof of Lemma 3.1.

The estimate (3.5) can be derived from [42, Section 5] or [46, Section 3],
along with the eigenvalue sticking estimate (B.23). The first estimate in
(3.6) follows directly from (B.21).

To prove the second estimate in (3.6), we apply the inverse function theorem
to (2.8) and use that $m(\mathfrak{a}_{i})=-\widetilde{\sigma}_{i}^{-1}$,
which yields the relation

$h^{\prime}(-\widetilde{\sigma}_{i}^{-1})=\left[m^{\prime}(\mathfrak{a}_{i})\right]^{-1}.$
(C.21)

Next, considering the contour $\Gamma_{i}$ defined below (B.27), we apply
Cauchy’s integral formula to obtain that

$m^{\prime}(\mathfrak{a}_{i})=\frac{1}{2\pi\mathrm{i}}\oint_{\Gamma_{i}}\frac{m(z)}{(z-\mathfrak{a}_{i})^{2}}\mathrm{d}z=\frac{1}{2\pi\mathrm{i}}\oint_{\Gamma_{i}}\frac{\mathsf{m}_{n}(z)}{(z-\mathfrak{a}_{i})^{2}}\mathrm{d}z+\operatorname{O}_{\prec}(n^{-1/2})=\mathsf{m}_{n}^{\prime}(\mathfrak{a}_{i})+\operatorname{O}_{\prec}(n^{-1/2}).$
(C.22)

In the second step, we utilized the local law (B.10), and in the third step,
we used the fact that, with high probability, $\mathsf{m}_{n}$ has no pole
inside $\Gamma_{i}$ due to the rigidity estimate (B.20). Since $r$ is fixed,
using (B.21) and (B.23), we get that

$\mathsf{m}_{n}^{\prime}(\mathfrak{a}_{i})=\frac{1}{n}\sum_{j=1}^{n}\frac{1}{(\lambda_{j}-{\mathfrak{a}}_{i})^{2}}=\frac{1}{n}\sum_{j=r+1}^{n}\frac{1}{(\widetilde{\lambda}_{j}-{\mathfrak{a}}_{i})^{2}}+\operatorname{O}_{\prec}(n^{-1})=\widehat{m}_{i}^{\prime}+\operatorname{O}_{\prec}(n^{-1/2}).$
(C.23)

Plugging (C.23) and (C.22) into (C.21) concludes the second estimate in (3.6).
With a similar argument, employing (B.10) and (B.21), we can deduce that

$-\widetilde{\sigma}_{i}^{-1}=m(\mathfrak{a}_{i})=\mathsf{m}_{n}(\mathfrak{a}_{i})+\operatorname{O}_{\prec}(n^{-1/2})=\widehat{m}_{i}+\operatorname{O}_{\prec}(n^{-1/2}).$
(C.24)

This verifies the third estimate in (3.6).

For the last estimate in (3.6), we first notice that by (2.15) and (3.5),

$K_{\varepsilon}^{+}=p-\operatorname{o}(p)$ (C.25)

for small enough ${\varepsilon}>0$. Since $m$ is a real increasing function
outside $\operatorname{supp}(\varrho)$, we have
$b_{1}=m(\lambda_{+})<m(\mathfrak{a}_{i})<0.$ Combining this with condition
(2.18), we obtain $\min_{j=r+1}^{p}|1+m(\mathfrak{a}_{i})\sigma_{j}|\gtrsim
1,$ which, together with (C.24), implies
$\min_{j=r+1}^{p}|1+\widehat{m}_{i}\sigma_{j}|\gtrsim 1$ with high
probability. Then, from (3.5), we see that

$K_{{\varepsilon}}^{-}=\operatorname{o}(p).$ (C.26)

On the other hand, by employing the weak convergence in (3.5), along with
(B.18) and (C.24)–(C.26), we obtain that for small $\delta>0$,

$\displaystyle\frac{1}{n}\sum_{j=1}^{p}\frac{\ell(\sigma_{j})\sigma_{j}}{[1+m(\mathfrak{a}_{i})\sigma_{j}]^{2}}$
$\displaystyle=\frac{1}{n}\sum_{j=r+1}^{\lceil(1-\delta)p\rceil}\frac{\ell(\sigma_{j})\sigma_{j}}{(1+m(\mathfrak{a}_{i})\sigma_{j})^{2}+\delta}+\operatorname{O}(\delta)$
$\displaystyle=\frac{1}{n}\sum_{j=r+1}^{\lceil(1-\delta)p\rceil}\frac{\ell(\widehat{\sigma}_{j})\widehat{\sigma}_{j}}{(1+\widehat{m}_{i}\widehat{\sigma}_{j})^{2}+\delta}+\operatorname{o}_{\mathbb{P}}(1)+\operatorname{O}(\delta)$
$\displaystyle=\frac{1}{n}\sum_{j=(r+1)\vee K_{\varepsilon}^{-}}^{p\wedge
K_{\varepsilon}^{+}}\frac{\ell(\widehat{\sigma}_{j})\widehat{\sigma}_{j}}{(1+\widehat{m}_{i}\widehat{\sigma}_{j})^{2}+\delta}+\operatorname{o}_{\mathbb{P}}(1)+\operatorname{O}(\delta)$
$\displaystyle=\frac{1}{n}\sum_{j=(r+1)\vee K_{\varepsilon}^{-}}^{p\wedge
K_{\varepsilon}^{+}}\frac{\ell(\widehat{\sigma}_{j})\widehat{\sigma}_{j}}{(1+\widehat{m}_{i}\widehat{\sigma}_{j})^{2}}+\operatorname{o}_{\mathbb{P}}(1)+\operatorname{O}(\delta).$

As $\delta$ is arbitrary, we conclude from this estimate that

$\frac{1}{n}\sum_{j=1}^{p}\frac{\ell(\sigma_{j})\sigma_{j}}{[1+m(\mathfrak{a}_{i})\sigma_{j}]^{2}}=\frac{1}{n}\sum_{j=(r+1)\vee
K_{\varepsilon}^{-}}^{p\wedge
K_{\varepsilon}^{+}}\frac{\ell(\widehat{\sigma}_{j})\widehat{\sigma}_{j}}{(1+\widehat{m}_{i}\widehat{\sigma}_{j})^{2}}+\operatorname{o}_{\mathbb{P}}(1).$
(C.27)

Now, plugging (B.21), (C.22)–(C.24), and (C.27) into (2.25), we obtain that

$\dot{m}_{0}(\mathfrak{a}_{i})=\frac{m^{\prime}(\mathfrak{a}_{i})}{n\mathfrak{a}_{i}}\sum_{j=1}^{p}\frac{\ell(\sigma_{j})\sigma_{j}}{[1+m(\mathfrak{a}_{i})\sigma_{j}]^{2}}=\frac{\widehat{m}_{i}^{\prime}}{n\widehat{\mathfrak{a}}_{i}}\sum_{j=(r+1)\vee
K_{\varepsilon}^{-}}^{p\wedge
K_{\varepsilon}^{+}}\frac{\ell(\widehat{\sigma}_{j})\widehat{\sigma}_{j}}{(1+\widehat{m}_{i}\widehat{\sigma}_{j})^{2}}+\operatorname{o}_{\mathbb{P}}(1).$
(C.28)

This concludes the last estimate in (3.6). ∎

###### Proof of Theorem 3.2.

The estimate (3.11) follows immediately from Lemma 3.1. For the estimate
(3.12), utilizing the eigenvalue rigidity estimate (B.20) and the eigenvalue
sticking estimate (B.23), we obtain that

$|\widetilde{\lambda}_{i}-\gamma_{i-r}|\prec n^{-2/3},\quad
i\in[\\![{r+1,\mathsf{K}}]\\!].$ (C.29)

Using (B.23) again, we can derive the following expression for
$\widehat{m}(\widetilde{\lambda}_{i})$ when $\eta=n^{-1/2}$:

$\displaystyle\widehat{m}(\widetilde{\lambda}_{i})$
$\displaystyle=\frac{1}{n}\sum_{j=r+1:|j-i|\geq
n^{1/4}}^{n}\frac{1}{\widetilde{\lambda}_{j}-\widetilde{\lambda}_{i}-\mathrm{i}\eta}+\frac{1}{n}\sum_{j=r+1:|j-i|<n^{1/4}}^{n}\frac{1}{\widetilde{\lambda}_{j}-\widetilde{\lambda}_{i}-\mathrm{i}\eta}$
$\displaystyle=\frac{1}{n}\sum_{j=1:|j+r-i|\geq
n^{1/4}}^{n-r}\frac{1}{{\lambda}_{j}-\widetilde{\lambda}_{i}-\mathrm{i}\eta}+\frac{1}{n}\sum_{j=1:|j+r-i|<n^{1/4}}^{n-r}\frac{1}{{\lambda}_{j}-\widetilde{\lambda}_{i}-\mathrm{i}\eta}+\operatorname{O}_{\prec}(n^{-1/4})$
$\displaystyle=\frac{1}{n}\sum_{j=1}^{n-r}\frac{1}{{\lambda}_{j}-\widetilde{\lambda}_{i}-\mathrm{i}\eta}+\operatorname{O}_{\prec}(n^{-1/4})=\mathsf{m}_{n}(\widetilde{\lambda}_{i}+\mathrm{i}\eta)+\operatorname{O}_{\prec}(n^{-1/4}).$
(C.30)

In the second step, we utilize the trivial bound
$|{\lambda}_{j}-\widetilde{\lambda}_{i}-\mathrm{i}\eta|^{-1}\leq\eta^{-1}$
when $|j+r-i|<n^{1/4}$, and when $|j+r-i|\geq n^{1/4}$, we have
$|{\lambda}_{j}-\widetilde{\lambda}_{i}|\gtrsim n^{-1/4}$ with high
probability due to the eigenvalue rigidity (B.20). When combined with the
averaged local law (B.9), the two estimates (C.29) and (C.30) imply that

$\widehat{m}(\widetilde{\lambda}_{i})=m(\widetilde{\lambda}_{i}+\mathrm{i}\eta)+\operatorname{O}_{\prec}(n^{-1/4})=m(\gamma_{i-r})+\operatorname{O}_{\prec}(n^{-1/4}),$
(C.31)

where in the second step we used the following bound for any constant $C>0$:

$|m(z_{1})-m(z_{2})|\lesssim|z_{1}-z_{2}|^{1/2}\quad\forall\ \
z_{1},z_{2}\in\mathbb{C}_{+},\ C^{-1}\leq|z_{1}|,|z_{2}|\leq C.$ (C.32)

(This estimate can be proven by employing the Stieltjes transform form of
$m(z)$ and the square root behaviors of $\varrho$ around the spectral edges.)
Similarly, using Lemma B.3, we can derive that

$\widehat{m}_{0}=\operatorname{Re}m(\mathrm{i}\eta)+\operatorname{O}_{\prec}(n^{-1/2})=m(0)+\operatorname{O}_{\prec}(n^{-1/2}).$
(C.33)

With the estimates (C.29), (C.31), and (C.33) at hand, the proof of (3.12) is
then similar to that of the last estimate in (3.6), and we omit the details. ∎

###### Proof of Lemma 3.3.

The estimate (3.13) follows from (3.6) and (C.22)–(C.24). The estimate (3.14)
can be derived from (C.29), (C.31), and (C.33). ∎

##  Appendix D Eigenvector moment flow and proof of Lemma 4.10

In this section, we study the dynamics of the eigenvalues and eigenvectors of
the rectangular DBM and establish Lemma 4.10. We first introduce some new
notations. We denote the processes of eigenvalues and eigenvectors as
$\bm{\lambda}(\cdot)\mathrel{\vbox{\hbox{.}\hbox{.}}}=(\bm{\lambda}(t))_{t\geqslant
0}$ and $U(\cdot)=(U(t))_{t\geqslant 0}$, respectively, where
$\bm{\lambda}(t)=(\lambda_{1}(t),\ldots,\lambda_{p}(t))$ and
$U(t)=({\bf{u}}_{1}(t),\ldots,{\bf{u}}_{p}(t))$. We further denote by
$\mathbb{E}^{U(0)}_{\bm{\lambda}(\cdot)}(\cdot)$ the expectation with respect
to the process $U(\cdot)$ conditioned on the eigenvalue process
$\bm{\lambda}(\cdot)$ and the initial state $U(0)$ of the eigenvector process.
Let $\mathsf{L}(t)$ be the generator associated with the process $U(t)$ (or
equivalently the dynamics given by (4.21)) conditioned on
$\bm{\lambda}(\cdot)$ and $U(0)$. In other words, for any smooth functions
$f:\mathbb{R}^{p^{2}}\rightarrow\mathbb{R}$, we have the equation

$\frac{\mathrm{d}}{\mathrm{d}t}\mathbb{E}_{\bm{\lambda}(\cdot)}^{U(0)}f\bigl{(}{U(t)}\bigr{)}\;=\;\mathbb{E}_{\bm{\lambda}(\cdot)}^{U(0)}(\mathsf{L}(t)f)\bigl{(}{U(t)}\bigr{)}\,,$
(D.1)

where $(\mathsf{L}(t)f)(U(t))$ represents the action of the generator on the
function $f$ evaluated at $U(t)$. Finally, we define the
$\bm{\lambda}(\cdot)$-measurable $p\times p$ matrices
$\Upsilon(t)=(\Upsilon_{kl}(t))$ with entries introduced in (4.25).

######  Lemma D.1.

The generator $\mathsf{L}(t)$ associated with the eigenvector flow (4.21)
conditioned on $\bm{\lambda}(\cdot)$ and $U(0)$ is

$\mathsf{L}(t)=\sum_{k,l=1}^{p}\frac{1}{2}\Upsilon_{kl}(t)(\mathsf{X}_{kl})^{2}\,,$
(D.2)

where $\mathsf{X}_{kl}$ is defined as

$\mathsf{X}_{kl}\;\mathrel{\vbox{\hbox{.}\hbox{.}}}=\;\sum_{j=1}^{p}\biggl{(}{{\bf{u}}_{k}(j)\frac{\partial}{\partial{\bf{u}}_{l}(j)}-{\bf{u}}_{l}(j)\frac{\partial}{\partial{\bf{u}}_{k}(j)}}\biggr{)}\,.$
(D.3)

###### Proof.

This result can be proved in the same way as Lemma 2.4 in [15], employing
(4.21) and Itô’s formula. The details are omitted. ∎

Now, we will follow the argument presented in [15, Section 3.1] to show that
$\mathsf{L}(t)$ acting on $f_{t}(\xi)$ in (4.22) can be interpreted as a
multi-particle random walk in random environments described by the DBM of
eigenvalues. This establishes the desired _eigenvector moment flow_ in Lemma
4.10. Now, we will derive a more general form of EMF that encompasses Lemma
4.10 as a special case. Let $\bm{\mathrm{v}}_{1},\dots,\bm{\mathrm{v}}_{M}$ be
a deterministic sequence of orthonormal vectors in $\mathbb{R}^{p}$, where $M$
is a fixed integer. For $(k,i)\in[\\![{p}]\\!]\times[\\![{M}]\\!],$ we define
the variables
$z_{k}^{i}\;\mathrel{\vbox{\hbox{.}\hbox{.}}}=\;\sqrt{p}\langle{\bm{\mathrm{v}}_{i}}\mspace{2.0mu},{\bm{\mathrm{u}}_{k}}\rangle$
and abbreviate
$Z=(z_{k}^{i}\mathrel{\vbox{\hbox{.}\hbox{.}}}(k,i)\in[\\![{p}]\\!]\times[\\![{M}]\\!])$.
We further define the collection of “particle configurations”

$\mathcal{S}_{n}\;\mathrel{\vbox{\hbox{.}\hbox{.}}}=\;\biggl{\\{}{\xi\in\mathbb{Z}^{[\\![{p}]\\!]\times[\\![{M}]\\!]}\mathrel{\vbox{\hbox{.}\hbox{.}}}\sum_{i=1}^{M}\xi_{k}^{i}\in
2\mathbb{Z}\text{\rm{ for all $k\in[\\![{p}]\\!]$}}}\biggr{\\}}\,.$ (D.4)

The special basis elements $e_{k}^{i}$,
$(k,i)\in[\\![{p}]\\!]\times[\\![{M}]\\!],$ for the configurations in
$\mathcal{S}_{n}$ are defined as
$(e_{k}^{i})_{l}^{j}\mathrel{\vbox{\hbox{.}\hbox{.}}}=\delta_{kl}\delta_{ij}$.
For each $\xi\in\mathcal{S}_{n}$, we assign a polynomial $P_{\xi}$ of the
variables in $Z$, defined as follows:

$P_{\xi}\;\mathrel{\vbox{\hbox{.}\hbox{.}}}=\;\prod_{k=1}^{p}\prod_{i=1}^{M}(z_{k}^{i})^{\xi_{k}^{i}}$
(D.5)

if $\xi_{k}^{i}\geqslant 0$ for all
$(k,i)\in[\\![{p}]\\!]\times[\\![{M}]\\!]$, and $P_{\xi}=0$ otherwise. With
the above notations, we can derive the following key result through direct
calculations.

######  Lemma D.2.

For $k,l\in[\\![{p}]\\!]$ and $i\in[\\![{M}]\\!]$, define the operator
$A_{kl}^{i}\mathrel{\vbox{\hbox{.}\hbox{.}}}\mathcal{S}_{n}\to\mathcal{S}_{n}$
through

$A_{kl}^{i}\xi\;\mathrel{\vbox{\hbox{.}\hbox{.}}}=\;\xi-
e_{k}^{i}+e_{l}^{i}\,.$ (D.6)

Then, we have that

$\displaystyle\mathsf{L}(t)P_{\xi}\;=\;$
$\displaystyle\sum_{k,l}\Upsilon_{kl}(t)\sum_{i}\Bigl{[}{\xi_{k}^{i}(\xi_{k}^{i}-1)P_{A_{kl}^{i}A_{kl}^{i}\xi}-\xi_{k}^{i}(\xi_{l}^{i}+1)P_{\xi}}\Bigr{]}$
$\displaystyle+$ $\displaystyle\sum_{k,l}\Upsilon_{kl}(t)\sum_{i\neq
j}\Bigl{[}{\xi_{k}^{i}\xi_{k}^{j}P_{A_{kl}^{i}A_{kl}^{j}\xi}-\xi_{k}^{i}\xi_{l}^{j}P_{A_{kl}^{i}A_{lk}^{j}\xi}}\Bigr{]}\,.$
(D.7)

###### Proof.

For any differentiable function $f(Z)$ of $Z$, by applying the chain rule to
$z_{k}^{i}=\sqrt{p}\sum_{j=1}^{p}\bm{\mathrm{v}}_{i}(j)\bm{\mathrm{u}}_{k}(j)$,
we obtain that

$\sum_{j=1}^{p}{\bf{u}}_{k}(j)\frac{\partial}{\partial{\bf{u}}_{l}(j)}f\;=\;\sum_{i=1}^{M}z_{k}^{i}\frac{\partial}{\partial
z_{l}^{i}}f\,.$ (D.8)

On the other hand, using the definition in (D.5), we obtain that

$z_{k}^{i}\frac{\partial}{\partial
z_{l}^{i}}P_{\xi}\;=\;\xi_{l}^{i}P_{\xi+e_{k}^{i}-e_{l}^{i}}\,.$ (D.9)

Now, using the above two identities, we can derive the following equation for
$k\not\sim l$:

$\displaystyle(\mathsf{X}_{kl})^{2}P_{\xi}\;=\;$
$\displaystyle\sum_{i}\Bigl{[}\xi_{k}^{i}(\xi_{k}^{i}-1)P_{\xi+2e_{l}^{i}-2e_{k}^{i}}+\xi_{l}^{i}(\xi_{l}^{i}-1)P_{\xi+2e_{k}^{i}-2e_{l}^{i}}-\bigl{(}{\xi_{k}^{i}(\xi_{l}^{i}+1)+\xi_{l}^{i}(\xi_{k}^{i}+1)}\bigr{)}P_{\xi}\Bigr{]}$
$\displaystyle+$ $\displaystyle\sum_{i\neq
j}\Bigl{[}\xi_{l}^{i}\xi_{l}^{j}P_{\xi+e_{k}^{i}+e_{k}^{j}-e_{l}^{i}-e_{l}^{j}}+\xi_{k}^{i}\xi_{k}^{j}P_{\xi+e_{l}^{i}+e_{l}^{j}-e_{k}^{i}-e_{k}^{j}}-\xi_{l}^{i}\xi_{k}^{j}P_{\xi+e_{k}^{i}+e_{l}^{j}-e_{l}^{i}-e_{k}^{j}}-\xi_{k}^{i}\xi_{l}^{j}P_{\xi+e_{l}^{i}+e_{k}^{j}-e_{k}^{i}-e_{l}^{j}}\Bigr{]}\,.$

Finally, utilizing the symmetry $\Upsilon_{kl}(t)=\Upsilon_{lk}(t)$, this
equation yields that

$\displaystyle\mathsf{L}(t)P_{\xi}\;=\;$
$\displaystyle\sum_{k,l}\Upsilon_{kl}(t)\sum_{i}\Bigl{[}{\xi_{k}^{i}(\xi_{k}^{i}-1)P_{\xi+2e_{l}^{i}-2e_{k}^{i}}-\xi_{k}^{i}(\xi_{l}^{i}+1)P_{\xi}}\Bigr{]}$
$\displaystyle+$ $\displaystyle\sum_{k,l}\Upsilon_{kl}(t)\sum_{i\neq
j}\Bigl{[}{\xi_{k}^{i}\xi_{k}^{j}P_{\xi+e_{l}^{i}+e_{l}^{j}-e_{k}^{i}-e_{k}^{j}}-\xi_{k}^{i}\xi_{l}^{j}P_{\xi+e_{l}^{i}+e_{k}^{j}-e_{k}^{i}-e_{l}^{j}}}\Bigr{]}\,.$

We then complete the proof by substituting the definition (D.6) into this
equation. ∎

######  Remark D.3.

In light of Lemma D.2, the generator $\mathsf{L}(t)$ acting on $P_{\xi}$ can
interpreted as a colored multi-particle random walk in random environments. In
this process, particles move on the lattice $[\\![{p}]\\!]$, with each
particle carrying a color from the set $[\\![{M}]\\!]$. The particles at each
site in $[\\![{p}]\\!]$ are unordered. The particle configuration is encoded
by $\xi\in\mathcal{S}_{n}$, where $\xi_{k}^{i}$ represents the number of
particles of color $i$ at site $k$. The operator $A_{kl}^{i}$ describes the
jump of a particle with color $i$ from site $k$ to site $l$. Thus, on the RHS
of (D.2), the first term corresponds to moving two particles of the same color
from site $k$ to site $l$, the third term represents moving two particles of
different colors from site $k$ to site $l$, and the last term involves
exchanging two particles of different colors between site $k$ and site $l$.
The second term on the right-hand side of (D.2) ensures that $\mathsf{L}$ acts
as the generator of a continuous-time Markov process. Note that under the
dynamics described by (D.2), every site $k$ always contains an even number of
particles, and the total number of particles of a given color remains
conserved. By applying the techniques developed in Appendix E below, we can
extend our analysis to the more intricate case of a colored EMF, leading to
the derivation of the more general results presented in (2.38). However, in
the interest of simplicity, we do not delve into these details within the
scope of this paper and defer their exploration to future work.

###### Proof of Lemma 4.10.

Note that $f(\xi)$ can be written as
$f_{t}(\xi)\;\mathrel{\vbox{\hbox{.}\hbox{.}}}=\;\mathbb{E}_{\bm{\lambda}(\cdot)}^{U(0)}Q_{\xi}(U(t))\,,$
where $Q_{\xi}$ is defined as

$Q_{\xi}\;\mathrel{\vbox{\hbox{.}\hbox{.}}}=\;P_{\xi}\prod_{k=i_{0}}^{\mathsf{K}}\frac{1}{(2\xi(k)-1)!!},\quad\text{with}\quad
P_{\xi}\;=\;\prod_{k=i_{0}}^{\mathsf{K}}z_{k}(t)^{2\xi(k)}\,.$

Then, Lemma 4.10 is a simple consequence of Lemmas D.1 and D.2 by taking $M=1$
and $\xi_{k}^{1}=2\xi(k)$. ∎

##  Appendix E Proof of Theorem 4.5

This section is dedicated to proving a key technical result of this
paper—Theorem 4.5. As demonstrated in Lemma 4.10 and Appendix D, the EMF
$f_{t}(\xi)$ can be viewed as a multi-particle random walk in a random
environment with generator (D.2). Through a careful analysis of the EMF, we
will show that $f_{t}(\xi)$ relaxes to the “equilibrium state”
$g_{t}(\xi,{\varepsilon}_{1})$ on the time scale $t\gg n^{-1/3}$. Towards the
end of this section, we will also briefly discuss how to extend our arguments
to the scenario where the eigenvectors are projected onto multiple distinct
directions.

The following proposition is an analog of Proposition 3.2 in [15]. Recall that
a measure $\pi$ on the configuration space is said to be _reversible_ with
respect to a generator $\mathsf{L}$ if, for any functions $\mathrm{f}$ and
$\mathrm{g}$, the following equality holds:

$\sum_{\xi}\pi(\xi)\mathrm{g}(\xi)\mathsf{L}\mathrm{f}(\xi)=\sum_{\xi}\pi(\xi)\mathrm{f}(\xi)\mathsf{L}\mathrm{g}(\xi).$

We then define the Dirichlet form with respect to the reversible measure $\pi$
as

$\mathsf{D}^{\pi}(\mathrm{f})=-\sum_{\xi}\pi(\xi)\mathrm{f}(\xi)\mathsf{L}\mathrm{f}(\xi).$

######  Proposition E.1 (Proposition 3.2 of [15]).

Define the measure $\pi$ on the configuration space by assigning the following
weight to each $\xi=(\xi_{1},\cdots,\xi_{p})$:

$\pi(\xi)=\prod_{i=1}^{p}\mathfrak{L}(\xi_{i}),\quad\mathfrak{L}(k):=\prod_{i=1}^{k}(1-(2i)^{-1}).$
(E.1)

Then, $\pi$ is a reversible measure for $\mathscr{B}$ defined in (4.24),
yielding that

$-\sum_{\xi}\pi(\xi)\mathrm{g}(\xi)\mathscr{B}\mathrm{f}(\xi)=\frac{1}{2}\sum_{\xi}\pi(\xi)\sum_{l\not\sim
k}\Upsilon_{kl}\xi(k)(1+2\xi(l))\bigl{(}{\mathrm{f}(\xi^{k\to
l})-\mathrm{f}(\xi)}\bigr{)}\bigl{(}{\mathrm{g}(\xi^{k\to
l})-\mathrm{g}(\xi)}\bigr{)}.$

In addition, the reversibility of $\pi$ remains valid for the generator
$\mathscr{B}^{\prime}$ obtained by replacing $\Upsilon$ in $\mathscr{B}$ with
any other symmetric matrix $\Upsilon^{\prime}$, i.e.,
$\Upsilon_{kl}^{\prime}=\Upsilon_{lk}^{\prime}$.

Without loss of generality, in the following proof, we will focus on the case
where

$i_{0}=1\quad\text{and}\quad c_{n}<1.$ (E.2)

In simpler terms, we assume that our model is free from outliers or trivial
eigenvalues at 0. In fact, under 2.3, the outlier eigenvalues
$\widetilde{\lambda}_{i}$, $i\in[\\![{r}]\\!]$, and the trivial eigenvalues at
0 are away from $\\{\widetilde{\lambda}_{i}\\}_{i=i_{0}}^{\mathsf{K}}$ with a
distance of order 1. Leveraging the finite speed of propagation estimate
established in Lemma C.14 of [31], we can show that the influence of the
eigenvalues/eigenvectors corresponding to the outliers and the trivial zero
eigenvalues on $f_{t}(\xi)$ is exponentially small. Consequently, the proof of
the general case without assuming (1.6) can be concluded by following the
arguments outlined in Remarks B.15 and C.2 of [31].

Given a small constant $\delta>0,$ we introduce non-negative cutoff functions
$\theta_{-,\delta}$ supported on $[\lambda_{-}/3,\infty)$ and
${\theta}_{+,\delta}$ supported on $(-\infty,3\lambda_{+}]$, satisfying the
following properties:

$\theta_{-,\delta}(x)=\begin{cases}0,&x<\lambda_{-}/3\\\ 1,&x\geqslant
2\lambda_{-}/3\end{cases},\qquad\theta_{+,\delta}=\begin{cases}0,&x>3\lambda_{+}\\\
1,&x\leqslant 2\lambda_{+}\end{cases};$ (E.3)

the derivatives of $\theta_{-}$ and $\theta_{+}$ satisfy that for all
$x\in\mathbb{R}$,

$\displaystyle 0\leq\theta_{-,\delta}^{\prime}(x)\lesssim
n^{\delta}\theta_{-,\delta}(x)+\exp(-n^{\delta/2}),\quad|\theta_{-,\delta}^{\prime\prime}(x)|\lesssim
n^{2\delta}\theta_{-,\delta}(x)+\exp(-n^{\delta/2}),$ (E.4) $\displaystyle
0\leq-\theta_{+,\delta}^{\prime}(x)\lesssim
n^{\delta}\theta_{+,\delta}(x)+\exp(-n^{\delta/2}),\quad|\theta_{+,\delta}^{\prime\prime}(x)|\lesssim
n^{2\delta}\theta_{+,\delta}(x)+\exp(-n^{\delta/2}).$

As an example, we can consider the function $\theta_{-,\delta}$ defined as

$\theta_{-,\delta}(x)=\frac{1}{2}\exp\left(n^{\delta}(x-\lambda_{-}/2)\right)\quad\text{for}\quad\lambda_{-}/2-n^{-\delta/3}\leq
x\leq\lambda_{-}/2.$

Then, we properly choose the value of $\theta_{-,\delta}$ on
$[\lambda_{-}/3,\lambda_{-}/2-n^{-\delta/3}]$ and
$[\lambda_{-}/2,2\lambda_{-}/3]$ such that (E.3) and (E.4) hold. We can choose
the function $\theta_{+,\delta}$ in a similar way. We are now prepared to
prove Theorem 4.5.

###### Proof of Theorem 4.5.

Let ${\varepsilon}>0$ be a sufficiently small constant. We choose
$g_{t}(\xi,{\varepsilon})$ as defined in (4.17), along with the cutoff
functions $\theta_{\pm,{\varepsilon}}$. Next, we introduce the function

$F_{t}\equiv
F_{t}(L,\mathsf{q},{\varepsilon}):=\theta_{-,{\varepsilon}}(\lambda_{\mathsf{K}}(t))\theta_{+,{\varepsilon}}(\lambda_{1}(t))\sum_{|\xi|=L}\pi(\xi)\widetilde{f}_{t}(\xi,{\varepsilon})^{\mathsf{q}},$
(E.5)

where $\mathsf{q}\in 2\mathbb{N}$ is an even integer, $\pi$ is defined in
(E.1), and

$\widetilde{f}_{t}(\xi,{\varepsilon}):=f_{t}(\xi)-g_{t}(\xi,{\varepsilon}).$

To conclude the proof, it suffices to show that for any fixed $K\in\mathbb{N}$
and $\mathsf{p}\in 2\mathbb{N}$, there exist constants
$0<{\varepsilon}<\tilde{\varepsilon}$ such that the following estimate holds
with high probability:

$\sup_{t\in[T_{1}/2,T_{2}]}F_{t}(K,\mathsf{p},\tilde{\varepsilon})\lesssim
n^{-{\varepsilon}\mathsf{p}+K}.$ (E.6)

Given constants $0<{\varepsilon}^{\prime}<{\varepsilon}$, $L\in\mathbb{N}$,
and $\mathsf{q}\in 2\mathbb{N}$, we define the function
$\mathcal{A}_{t}(L,\mathsf{q},{\varepsilon},{\varepsilon}^{\prime}):=\log(F_{t}(L,\mathsf{q},{\varepsilon})+n^{-{\varepsilon}^{\prime}\mathsf{q}+L})$.
Then, the estimate (E.6) can be derived from the following technical lemma on
$\mathcal{A}_{t}$, the proof of which is deferred to Section E.1.

######  Lemma E.2.

Fix any $L\in\mathbb{N}$. Under the assumptions of Theorem 4.5, suppose there
exist constants $\tilde{\varepsilon}>0$,
${\varepsilon}\in(0,\min\\{\tilde{\varepsilon},1/12\\})$, and a fixed integer
$\mathsf{p}_{0}\in 2\mathbb{N}$, such that (E.6) holds with high probability,
uniformly in $t\in[T_{1}/2,T_{2}]$, for $K=L-1$ and any fixed even integer
$\mathsf{p}\geq\mathsf{p}_{0}$. Then, there exists a constant
${\varepsilon}^{\prime\prime}\in(0,\min\\{{\varepsilon},3\mathsf{c}/8\\})$ and
a fixed $\mathsf{q}_{0}\in 2\mathbb{N}$ such that the function
$\mathcal{A}_{t}(L,\mathsf{q},{\varepsilon},{\varepsilon}^{\prime})$ satisfies
the following SDE for any fixed
$0<{\varepsilon}^{\prime}\leq{\varepsilon}^{\prime\prime}$ and
$\mathsf{q}\geq\mathsf{q}_{0}$:

$\mathrm{d}\mathcal{A}_{t}(L,\mathsf{q},{\varepsilon},{\varepsilon}^{\prime})=S(t)\mathrm{d}t+\sigma(t)\mathrm{d}B(t),$
(E.7)

where $B(t)$ denotes a standard Brownian motion, $S(t)$ can be expressed as

$S(t)=-C(t)\frac{F_{t}}{F_{t}+n^{-{\varepsilon}^{\prime}\mathsf{q}+L}}-s_{+}(t)+\mathcal{E}(t),$

and the following estimates hold with high probability, uniformly in
$t\in[T_{1}/2,T_{2}]$: $\sigma(t)=\operatorname{O}(n^{-1/2+2{\varepsilon}})$,
$C(t)\gtrsim n^{1/3-\mathsf{c}/4}$, $s_{+}(t)\geq 0$, and
$\mathcal{E}(t)=\operatorname{O}(n^{2{\varepsilon}})$.

We now proceed to prove (E.6) by utilizing Lemma E.2. First, we note that
(E.6) trivially holds for $K=0$ since
$F_{t}(0,\mathsf{p},{\varepsilon}_{0})=0$. Now, suppose (E.6) holds for
$K=L-1$ and any fixed even integer $\mathsf{p}\geq\mathsf{p}_{0}$. Then, by
applying Lemma E.2, we can deduce that the equation (E.7) holds for some
constant ${\varepsilon}^{\prime}>0$ and $\mathsf{q}\geq\mathsf{q}_{0}$. By
applying the Burkholder-Davis-Gundy inequality to the diffusion term in (E.7),
we obtain that the following event holds with high probability:

$\left\\{\mathcal{A}_{t}\leq{\mathcal{A}}_{s}-\alpha\int_{s}^{t}\frac{F_{t^{\prime}}}{F_{t^{\prime}}+n^{-{\varepsilon}^{\prime}\mathsf{q}+L}}\mathrm{d}t^{\prime}+\operatorname{O}\left(n^{-{1}/{3}+2{\varepsilon}+\mathsf{c}}\right)\
\forall t,s\in[T_{1}/2,T_{2}]\right\\},$ (E.8)

where $\alpha=cn^{1/3-\mathsf{c}/4}$ for a small constant $c>0$. Note that,
according to the definition of $\mathcal{A}_{t}$, we have the deterministic
rough bound $\mathcal{A}_{s}\leq C\log n$ for a constant $C>0$. Therefore, the
above event implies that

$\left\\{{\mathcal{A}}_{t}\leq\min({\mathcal{A}}_{T_{1}/2},C\log
n)-\alpha\int_{T_{1}/2}^{t}\frac{F_{s}}{F_{s}+n^{-{\varepsilon}^{\prime}\mathsf{q}+L}}\mathrm{d}s+\operatorname{o}(1)\
\forall t\in[T_{1}/2,T_{2}]\right\\}$ (E.9)

holds with high probability. Now, we define the following stopping time (with
the convention $\inf\emptyset=T_{2}$):

$\tilde{t}:=(T_{1}/2)\vee\inf\\{t\in[T_{1}/2,T_{2}]:\mathcal{A}_{t}\leq\log(2n^{-{\varepsilon}^{\prime}\mathsf{q}+L})\\}.$

Noticing that ${F_{t}}/(F_{t}+n^{-{\varepsilon}^{\prime}\mathsf{q}+L})\geq
1/2$ for $t\in[T_{1}/2,\tilde{t}]$, on the event (E.9), we have

$\mathcal{A}_{\tilde{t}}\leq\min(\mathcal{A}_{T_{1}/2},C\log
n)-\frac{\alpha}{2}(\tilde{t}-T_{1}/2)+\operatorname{o}(1)\quad\text{w.h.p.}$

Since $\alpha T_{1}/2\gtrsim n^{3\mathsf{c}/4}\gg\log n$, this estimate
implies that

$\tilde{t}\leq T_{1}\quad\text{w.h.p.}$ (E.10)

Now, we apply this fact to (E.8) with $s=\tilde{t}$ and conclude that

$\sup_{t\in[T_{1},T_{2}]}{\mathcal{A}}_{t}\leq{\mathcal{A}}_{\tilde{t}}+1\quad\text{w.h.p.}$
(E.11)

Together with the definition of $F_{t}$, this implies that with high
probability,
$\sup_{t\in[T_{1},T_{2}]}F_{t}\leq(2e-1)n^{-{\varepsilon}^{\prime}\mathsf{q}+L},$
which completes the induction argument and establishes (E.6).

Now, given (E.6) with $K=L$ and $0<{\varepsilon}<\tilde{\varepsilon}$, we
choose ${\varepsilon}_{1}=\tilde{\varepsilon}$ and
${\varepsilon}_{0}={\varepsilon}/2$ and let $\mathsf{p}$ be sufficiently large
so that $-{\varepsilon}\mathsf{p}+L<-{\varepsilon}_{0}\mathsf{p}$. This allows
us to obtain the following inequality:

$F_{t}(L,\mathsf{p},{\varepsilon}_{1})=\theta_{-,{\varepsilon}_{1}}(\lambda_{\mathsf{K}}(t))\theta_{+,{\varepsilon}_{1}}(\lambda_{1}(t))\sum_{|\xi|=L}\pi(\xi)\widetilde{f}_{t}(\xi,{\varepsilon}_{1})^{\mathsf{p}}\ll
n^{-{\varepsilon}_{0}\mathsf{p}}\quad\text{w.h.p.}$

By utilizing the rigidity estimate (B.37) for $\lambda_{1}(t)$ and
$\lambda_{\mathsf{K}}(t)$, and considering the definitions of
$\theta_{\pm,{\varepsilon}_{1}}$, we can deduce that
$\theta_{-,{\varepsilon}_{1}}(\lambda_{\mathsf{K}}(t))=\theta_{+,{\varepsilon}_{1}}(\lambda_{1}(t))=1$
with high probability. Consequently, taking into account the fact that
$\pi(\xi)$ are all of order 1, the above estimate implies that

$\widetilde{f}_{t}(\xi,{\varepsilon}_{1})^{\mathsf{p}}\leq
n^{-{\varepsilon}_{0}\mathsf{p}}\quad\text{w.h.p.}$

for any configuration $\xi$ with $|\xi|=L$. This concludes the proof of
Theorem 4.5.

∎

###  E.1 Proof of Lemma E.2

In this subsection, we present the proof of Lemma E.2. By Itô’s formula, we
have

$d{\mathcal{A}}_{t}=\frac{\mathrm{d}F_{t}}{F_{t}+n^{-{\varepsilon}^{\prime}\mathsf{q}+L}}-\frac{\mathrm{d}[F]_{t}}{2(F_{t}+n^{-{\varepsilon}^{\prime}\mathsf{q}+L})^{2}},$

where the second term on the RHS is negative. Hence, it suffices to prove that
the following equation holds for $t\in[T_{1}/2,T_{2}]$:

$\frac{\mathrm{d}F_{t}}{F_{t}+n^{-{\varepsilon}^{\prime}\mathsf{q}+L}}=-C(t)\frac{F_{t}}{F_{t}+n^{-{\varepsilon}^{\prime}\mathsf{q}+L}}\mathrm{d}t-s_{+}(t)\mathrm{d}t+\mathcal{E}(t)\mathrm{d}t+\sigma(t)\mathrm{d}B(t),$
(E.12)

where $s_{+}(t)\geq 0$ and $C(t)$, $\mathcal{E}(t)$, $\sigma(t)$ satisfy the
desired estimates stated in Lemma E.2.

For simplicity of notation, in the following proof, we will abbreviate
$F\equiv F_{t}$, $f\equiv f_{t}$, $g\equiv g_{t}$,
$\widetilde{f}\equiv\widetilde{f}_{t}$, $\lambda_{k}\equiv\lambda_{k}(t)$,
${\bf{u}}_{k}\equiv{\bf{u}}_{k}(t)$,
$\phi(\lambda_{k}(t),{\varepsilon})\equiv\phi(\lambda_{k})$,
$\theta_{-}\equiv\theta_{-,{\varepsilon}}(\lambda_{\mathsf{K}})$,
$\theta_{+}\equiv\theta_{+,{\varepsilon}}(\lambda_{1})$, and
$\mathsf{m}_{n}\equiv\mathsf{m}_{n,t}$,
$\underline{\mathsf{m}}_{n}\equiv\underline{\mathsf{m}}_{n,t}$ (where we adopt
the notations in Section B.3 with $\Lambda_{0}=D_{00}$ and
$\Lambda_{t}=D_{00}+t$). Using (4.20) and Itô’s formula, we get that

$\displaystyle\mathrm{d}\phi(\lambda_{k})=\phi^{\prime}(\lambda_{k})\mathrm{d}\lambda_{k}+\phi^{\prime\prime}(\lambda_{k})\frac{2\lambda_{k}}{n}\mathrm{d}t,$
(E.13)
$\displaystyle\mathrm{d}g(\xi)=\sum_{k}\Big{(}\partial_{\lambda_{k}}g(\xi)\Big{)}\mathrm{d}\lambda_{k}+\sum_{k}\Big{(}\partial^{2}_{\lambda_{k}}g(\xi)\Big{)}\frac{2\lambda_{k}}{n}\mathrm{d}t,$
(E.14)
$\displaystyle\mathrm{d}\theta_{-}=\theta_{-}^{\prime}\mathrm{d}\lambda_{\mathsf{K}}+\theta_{-}^{\prime\prime}\frac{2\lambda_{\mathsf{K}}}{n}\mathrm{d}t,\quad\mathrm{d}\theta_{+}=\theta_{+}^{\prime}\mathrm{d}\lambda_{1}+\theta_{+}^{\prime\prime}\frac{2\lambda_{1}}{n}\mathrm{d}t.$
(E.15)

Note that by the definition of $g$ in (4.17) and the estimate (4.18), we have

$\theta_{-}\theta_{+}\partial_{\lambda_{k}}g(\xi)=\operatorname{O}(n^{{\varepsilon}}),\quad\theta_{-}\theta_{+}\partial^{2}_{\lambda_{k}}g(\xi)=\operatorname{O}(n^{2{\varepsilon}}).$
(E.16)

Now, applying Itô’s formula to $F$ defined in (E.5) and using (E.13)–(E.15),
we can derive that

$\displaystyle\mathrm{d}F=$
$\displaystyle\sum_{\xi}\pi(\xi)\widetilde{f}(\xi)^{\mathsf{q}}\left(\theta_{-}\theta_{+}^{\prime}\mathrm{d}\lambda_{1}+\theta_{-}\theta_{+}^{\prime\prime}\frac{2\lambda_{1}}{n}\mathrm{d}t+\theta_{-}^{\prime}\theta_{+}\mathrm{d}\lambda_{\mathsf{K}}+\theta_{-}^{\prime\prime}\theta_{+}\frac{2\lambda_{\mathsf{K}}}{n}\mathrm{d}t\right)$
(E.17)
$\displaystyle+\theta_{-}\theta_{+}\sum_{\xi}\pi(\xi)\left(\mathsf{q}\widetilde{f}(\xi)^{\mathsf{q}-1}(\mathrm{d}f(\xi)-\mathrm{d}g(\xi))+\mathsf{q}(\mathsf{q}-1)\widetilde{f}(\xi)^{\mathsf{q}-2}\sum_{k}\big{(}\partial_{\lambda_{k}}g(\xi)\big{)}^{2}\frac{2\lambda_{k}}{n}\mathrm{d}t\right)$
(E.18)
$\displaystyle-\sum_{\xi}\pi(\xi)\left(\mathsf{q}\widetilde{f}(\xi)^{\mathsf{q}-1}\right)\left(\partial_{\lambda_{\mathsf{K}}}g(\xi)\theta_{-}^{\prime}\theta_{+}\frac{4\lambda_{\mathsf{K}}}{n}+\partial_{\lambda_{1}}g(\xi)\theta_{-}\theta_{+}^{\prime}\frac{4\lambda_{1}}{n}\right)\mathrm{d}t.$
(E.19)

Here, we have used the convention that $\sum_{\xi}\equiv\sum_{\xi:|\xi|=L}.$
Next, we will control these terms one by one.

First, considering the term (E.17), we utilize (E.3) and (E.4) with
$\delta={\varepsilon}$, along with the SDE (4.20), to obtain that

$\displaystyle\theta_{-}\theta_{+}^{\prime}\mathrm{d}\lambda_{1}+\theta_{-}\theta_{+}^{\prime\prime}\frac{2\lambda_{1}}{n}\mathrm{d}t=$
$\displaystyle~{}\theta_{-}\theta_{+}^{\prime}2\sqrt{\lambda_{1}}\frac{\mathrm{d}B_{11}}{\sqrt{n}}+\theta_{-}\theta_{+}^{\prime}\bigg{(}1+\frac{1}{n}\sum_{k\neq
1}\frac{\lambda_{1}+\lambda_{k}}{\lambda_{1}-\lambda_{k}}\bigg{)}\mathrm{d}t$
(E.20)
$\displaystyle~{}+\operatorname{O}\left(n^{-1+2{\varepsilon}}\theta_{-}\theta_{+}+\exp(-n^{{\varepsilon}/2})\right)\mathrm{d}t.$

It then follows from (E.20) and the bound (E.4) that

$\sum_{\xi}\frac{\pi(\xi)\widetilde{f}(\xi)^{\mathsf{q}}}{F+n^{-{\varepsilon}^{\prime}\mathsf{q}+L}}\left(\theta_{-}\theta_{+}^{\prime}\mathrm{d}\lambda_{1}+\theta_{-}\theta_{+}^{\prime\prime}\frac{2\lambda_{1}}{n}\mathrm{d}t\right)=\operatorname{O}(n^{-1+2{\varepsilon}})\mathrm{d}t+\operatorname{O}(n^{-1/2+{\varepsilon}})\mathrm{d}B_{11}-s_{+}^{(1)}(t)\mathrm{d}t,$
(E.21)

where $s_{+}^{(1)}(t)$ represents a positive term arising from the second term
on the RHS of (E.20). Using a similar argument, we can also show that

$\displaystyle\sum_{\xi}\frac{\pi(\xi)\widetilde{f}(\xi)^{\mathsf{q}}}{F+n^{-{\varepsilon}^{\prime}\mathsf{q}+L}}\left(\theta_{-}^{\prime}\theta_{+}\mathrm{d}\lambda_{\mathsf{K}}+\theta_{-}^{\prime\prime}\theta_{+}\frac{2\lambda_{\mathsf{K}}}{n}\mathrm{d}t\right)$
$\displaystyle=\operatorname{O}(n^{-1+2{\varepsilon}})\mathrm{d}t+\operatorname{O}(n^{-1/2+{\varepsilon}})\mathrm{d}B_{\mathsf{K}\mathsf{K}}+\frac{\theta_{-}^{\prime}\theta_{+}\sum_{\xi}\pi(\xi)\widetilde{f}(\xi)^{\mathsf{q}}}{F+n^{-{\varepsilon}^{\prime}\mathsf{q}+L}}\bigg{(}1+\frac{1}{n}\sum_{k\neq\mathsf{K}}\frac{\lambda_{\mathsf{K}}+\lambda_{k}}{\lambda_{\mathsf{K}}-\lambda_{k}}\bigg{)}\mathrm{d}t$
$\displaystyle=\operatorname{O}(n^{-1+2{\varepsilon}})\mathrm{d}t+\operatorname{O}(n^{-1/2+{\varepsilon}})\mathrm{d}B_{\mathsf{K}\mathsf{K}}-s_{+}^{(2)}(t)\mathrm{d}t+\left(1+\frac{p-\mathsf{K}}{n}\right)\frac{\theta_{-}^{\prime}\theta_{+}\sum_{\xi}\pi(\xi)\widetilde{f}(\xi)^{\mathsf{q}}}{F+n^{-{\varepsilon}^{\prime}\mathsf{q}+L}}\mathrm{d}t,$

where $s_{+}^{(2)}(t)$ represents another positive term. We can exploit the
rigidity of $\lambda_{\mathsf{K}}$ given by (B.37) and the fact that
$|\lambda_{-,t}-\lambda_{-}|=\operatorname{o}(1)$ (as stated in (i) of Lemma
B.9) to conclude that $\theta_{-}^{\prime}=0$ with high probability, uniformly
in $t\in[T_{1}/2,T_{2}]$. Consequently, we obtain a bound similar to (E.21):

$\sum_{\xi}\frac{\pi(\xi)\widetilde{f}(\xi)^{\mathsf{q}}}{F+n^{-{\varepsilon}^{\prime}\mathsf{q}+L}}\left(\theta_{-}^{\prime}\theta_{+}\mathrm{d}\lambda_{\mathsf{K}}+\theta_{-}^{\prime\prime}\theta_{+}\frac{2\lambda_{\mathsf{K}}}{n}\mathrm{d}t\right)=\mathcal{E}^{(2)}(t)\mathrm{d}t+\operatorname{O}(n^{-1/2+{\varepsilon}})\mathrm{d}B_{\mathsf{K}\mathsf{K}}-s_{+}^{(2)}(t)\mathrm{d}t,$
(E.22)

where $\mathcal{E}^{(2)}(t)$ is a variable satisfying
$\mathcal{E}^{(2)}(t)=\operatorname{O}(n^{-1+2{\varepsilon}})$ with high
probability, uniformly in $t\in[T_{1}/2,T_{2}]$.

Next, we control (E.19). For any $0<l<\mathsf{q}$, we have the trivial bound

$|\widetilde{f}(\xi)|^{\mathsf{q}-l}\lesssim
n^{l{\varepsilon}}\widetilde{f}(\xi)^{\mathsf{q}}+n^{-(\mathsf{q}-l){\varepsilon}}.$
(E.23)

Utilizing (E.23) and the estimates (E.16) and (E.4), we get that

$\sum_{\xi}\frac{\pi(\xi)\mathsf{q}\widetilde{f}(\xi)^{\mathsf{q}-1}}{F+n^{-{\varepsilon}^{\prime}\mathsf{q}+L}}\left(\partial_{\lambda_{\mathsf{K}}}g(\xi)\theta_{-}^{\prime}\theta_{+}\frac{4\lambda_{\mathsf{K}}}{n}+\partial_{\lambda_{1}}g(\xi)\theta_{+}^{\prime}\theta_{-}\frac{4\lambda_{1}}{n}\right)\mathrm{d}t=\operatorname{O}(n^{-1+3\varepsilon})\mathrm{d}t.$
(E.24)

Similarly, the second term of (E.18) can be bounded by

$\mathsf{q}(\mathsf{q}-1)\theta_{-}\theta_{+}\sum_{\xi}\frac{\pi(\xi)\widetilde{f}(\xi)^{\mathsf{q}-2}}{F+n^{-{\varepsilon}^{\prime}\mathsf{q}+L}}\sum_{k}\left(\partial_{\lambda_{k}}g(\xi)\right)^{2}\frac{2\lambda_{k}}{n}\mathrm{d}t=\operatorname{O}(n^{-1+4\varepsilon})\mathrm{d}t.$
(E.25)

Combing (E.21), (E.22), (E.24), and (E.25), we obtain that

$\displaystyle\frac{\mathrm{d}F}{F+n^{-{\varepsilon}^{\prime}\mathsf{q}+L}}=$
$\displaystyle~{}\frac{\mathsf{q}\theta_{-}\theta_{+}}{F+n^{-{\varepsilon}^{\prime}\mathsf{q}+L}}\sum_{\xi}\pi(\xi)\widetilde{f}(\xi)^{\mathsf{q}-1}(\mathrm{d}f(\xi)-\mathrm{d}g(\xi))$
(E.26)
$\displaystyle~{}+\mathcal{E}_{1}(t)\mathrm{d}t+\operatorname{O}(n^{-1/2+{\varepsilon}})\mathrm{d}B_{11}+\operatorname{O}(n^{-1/2+{\varepsilon}})\mathrm{d}B_{\mathsf{K}\mathsf{K}}-(s_{+}^{(1)}(t)+s_{+}^{(2)}(t))\mathrm{d}t,$

where $\mathcal{E}_{1}(t)$ is a variable satisfying
$\mathcal{E}_{1}(t)=\operatorname{O}(n^{-1+4{\varepsilon}})$ with high
probability, uniformly in $t\in[T_{1}/2,T_{2}]$.

Now, to conclude (E.12), it remains to estimate the first term on the RHS of
(E.26). By (4.26), we express the term as

$\displaystyle\theta_{-}\theta_{+}\sum_{\xi}\pi(\xi)\widetilde{f}(\xi)^{\mathsf{q}-1}\mathrm{d}f(\xi)=\theta_{-}\theta_{+}\sum_{\xi}\pi(\xi)\widetilde{f}(\xi)^{\mathsf{q}-1}\sum_{k\not\sim
l}\Upsilon_{kl}(t)\,\xi(k)(1+2\xi(l))\bigl{(}{f(\xi^{k\to
l})-f(\xi)}\bigr{)}\mathrm{d}t.$ (E.27)

On the other hand, using (4.20) and (E.14), we obtain:

$\displaystyle\theta_{-}\theta_{+}\sum_{\xi}\pi(\xi)\widetilde{f}(\xi)^{\mathsf{q}-1}\mathrm{d}g(\xi)=\theta_{-}\theta_{+}\sum_{\xi}\pi(\xi)\widetilde{f}(\xi)^{\mathsf{q}-1}$
(E.28)
$\displaystyle\quad\times\sum_{k}\bigg{[}\partial_{\lambda_{k}}g(\xi)\sum_{l:l\not\sim
k}\frac{\lambda_{k}+\lambda_{l}}{n(\lambda_{k}-\lambda_{l})}\mathrm{d}t+\partial_{\lambda_{k}}g(\xi)\mathrm{d}t+\partial^{2}_{\lambda_{k}}g(\xi)\frac{2\lambda_{k}}{n}\mathrm{d}t+(\partial_{\lambda_{k}}g(\xi))2\sqrt{\lambda}_{k}\frac{\mathrm{d}B_{kk}}{\sqrt{n}}\bigg{]}.$

The diffusion part of the above equation has the same law as

$\theta_{-}\theta_{+}\bigg{[}\sum_{k}\bigg{|}\sum_{\xi}\Big{(}\pi(\xi)\widetilde{f}(\xi)^{\mathsf{q}-1}\left(\partial_{\lambda_{k}}g(\xi)\right)\frac{2\sqrt{\lambda}_{k}}{\sqrt{n}}\Big{)}\bigg{|}^{2}\bigg{]}^{1/2}\mathrm{d}B,$

where $B$ is a standard Brownian motion. Note that the coefficient satisfies
that

$\displaystyle\theta_{-}\theta_{+}\bigg{[}\sum_{k}\bigg{|}\sum_{\xi}\pi(\xi)\widetilde{f}(\xi)^{\mathsf{q}-1}\left(\partial_{\lambda_{k}}g(\xi)\right)\frac{2\sqrt{\lambda}_{k}}{\sqrt{n}}\bigg{|}^{2}\bigg{]}^{1/2}=\operatorname{O}(n^{-1/2+{\varepsilon}})\cdot\theta_{-}\theta_{+}\sum_{\xi}\pi(\xi)|\widetilde{f}(\xi)|^{\mathsf{q}-1}.$
(E.29)

Here, we have used (E.16) and the facts that $\pi(\xi)\asymp 1$ and
$\partial_{\lambda_{k}}g(\xi)$ is non-zero for at most $L$ many $k$’s.
Plugging (E.27)–(E.29) into (E.26) and using (E.23) and (E.16), we can deduce
that

$\displaystyle\frac{\mathrm{d}F}{F+n^{-{\varepsilon}^{\prime}\mathsf{q}+L}}=$
$\displaystyle~{}\sum_{\xi}\frac{\mathsf{q}\theta_{-}\theta_{+}\pi(\xi)\widetilde{f}(\xi)^{\mathsf{q}-1}}{F+n^{-{\varepsilon}^{\prime}\mathsf{q}+L}}\sum_{k\not\sim
l}\Bigg{(}\Upsilon_{kl}\,\xi(k)(1+2\xi(l))\bigl{(}{f(\xi^{k\to
l})-f(\xi)}\bigr{)}-\frac{(\lambda_{k}+\lambda_{l})\partial_{\lambda_{k}}g(\xi)}{n(\lambda_{k}-\lambda_{l})}\Bigg{)}\mathrm{d}t$
$\displaystyle~{}+\mathcal{E}_{2}(t)\mathrm{d}t+\operatorname{O}(n^{-1/2+2{\varepsilon}})\mathrm{d}\widetilde{B}-(s_{+}^{(1)}(t)+s_{+}^{(2)}(t))\mathrm{d}t,$
(E.30)

where $\widetilde{B}$ is another standard Brownian motion and
$\mathcal{E}_{2}(t)$ is a variable satisfying
$\mathcal{E}_{2}(t)=\operatorname{O}(n^{2{\varepsilon}})$ with high
probability, uniformly in $t\in[T_{1}/2,T_{2}]$.

Taking $\eta=n^{-2/3+\mathsf{c}/2}$, we define

$\mathrm{T}_{kl}\equiv\mathrm{T}_{kl}(\eta)\mathrel{\vbox{\hbox{.}\hbox{.}}}=\mathbf{1}(k\not\sim
l)\frac{(\lambda_{k}-\lambda_{l})^{2}}{(\lambda_{k}-\lambda_{l})^{2}+\eta^{2}},\quad\mathrm{T}_{kl}^{c}:=\mathbf{1}(k\not\sim
l)-\mathrm{T}_{kl}.$

Then, we decompose the first term on the RHS of (E.1) as

$\mathsf{q}(F+n^{-{\varepsilon}^{\prime}\mathsf{q}+L})^{-1}\sum_{\xi}\pi(\xi)\widetilde{f}(\xi)^{\mathsf{q}-1}\left(A_{1}+A_{2}+A_{3}\right)\mathrm{d}t,$
(E.31)

where $A_{1}$, $A_{2}$, and $A_{3}$ are defined as

$\displaystyle A_{1}:=\theta_{-}\theta_{+}\sum_{k\not\sim
l}\Upsilon_{kl}\mathrm{T}_{kl}\xi(k)(1+2\xi(l))\bigl{(}{f(\xi^{k\to
l})-f(\xi)}\bigr{)},$ (E.32) $\displaystyle
A_{2}:=-\theta_{-}\theta_{+}\sum_{k\not\sim
l}\mathrm{T}_{kl}\frac{(\lambda_{k}+\lambda_{l})\partial_{\lambda_{k}}g(\xi)}{n(\lambda_{k}-\lambda_{l})},$
(E.33) $\displaystyle A_{3}:=\theta_{-}\theta_{+}\sum_{k\not\sim
l}\mathrm{T}_{kl}^{c}\left[\Upsilon_{kl}\xi(k)(1+2\xi(l))\bigl{(}{f(\xi^{k\to
l})-f(\xi)}\bigr{)}-\frac{(\lambda_{k}+\lambda_{l})\partial_{\lambda_{k}}g(\xi)}{n(\lambda_{k}-\lambda_{l})}\right].$
(E.34)

The remaining part of the proof focuses on controlling the above three terms
one by one. Specifically, we will demonstrate that $A_{1}$ contributes to the
dominant drift term in (E.12).

The term $A_{2}$. We start with the term $A_{2}$:

$\displaystyle A_{2}$ $\displaystyle=\theta_{-}\theta_{+}\sum_{k\not\sim
l}\frac{\partial_{\lambda_{k}}g(\xi)}{n}\frac{\lambda_{l}^{2}-\lambda_{k}^{2}}{(\lambda_{l}-\lambda_{k})^{2}+\eta^{2}}=\theta_{-}\theta_{+}\left[\sum_{k}2c_{n}\lambda_{k}\partial_{\lambda_{k}}g(\xi)\cdot\operatorname{Re}\underline{\mathsf{m}}_{n}(\lambda_{k}+\mathrm{i}\eta)+\sum_{k\not\sim
l}\mathrm{T}_{kl}\frac{\partial_{\lambda_{k}}g(\xi)}{n}\right],$

where in the second step we used the spectral decomposition of
$\underline{\mathsf{m}}_{n}$:

$\operatorname{Re}\underline{\mathsf{m}}_{n}(\lambda_{k}+\mathrm{i}\eta)=\operatorname{Re}\frac{1}{p}\sum_{l}\frac{1}{\lambda_{l}-(\lambda_{k}+\mathrm{i}\eta)}=\frac{1}{p}\sum_{l:l\not\sim
k}\frac{\lambda_{l}-\lambda_{k}}{(\lambda_{l}-\lambda_{k})^{2}+\eta^{2}}.$

Since $\partial_{\lambda_{k}}g(\xi)$ is non-zero for at most $L$ many $k$’s,
using (E.16), the rigidity of $\lambda_{k}$ in (B.37), the trivial bound
$\mathrm{T}_{kl}<1$ and the local law (B.35), we obtain that
$A_{2}=\operatorname{O}(n^{\varepsilon})$ with high probability. Together with
(E.23), it implies that

$\frac{\mathsf{q}}{F+n^{-{\varepsilon}^{\prime}\mathsf{q}+L}}\sum_{\xi}\pi(\xi)\widetilde{f}(\xi)^{\mathsf{q}-1}A_{2}=\operatorname{O}(n^{2{\varepsilon}})\quad\text{w.h.p.}$
(E.35)

The term $A_{3}$. For the term $A_{3}$, we will establish that

$\frac{\mathsf{q}}{F+N^{-{\varepsilon}^{\prime}\mathsf{q}+L}}\sum_{\xi}\pi(\xi)\widetilde{f}(\xi)^{\mathsf{q}-1}A_{3}=\mathcal{E}_{3}(t)-s_{+}^{(3)}(t),$
(E.36)

where $s_{+}^{(3)}(t)\geq 0$ and $\mathcal{E}_{3}(t)$ represents a variable
that satisfies
$\mathcal{E}_{3}(t)=\operatorname{O}(n^{-2/3+4{\varepsilon}+\mathsf{c}})$ with
high probability, uniformly in $t\in[T_{1}/2,T_{2}]$. By utilizing Proposition
E.1, we can express the LHS of (E.36) as

$\displaystyle\frac{\mathsf{q}\theta_{-}\theta_{+}}{F+n^{-{\varepsilon}^{\prime}\mathsf{q}+L}}\sum_{\xi}\pi(\xi)\bigg{[}$
$\displaystyle\frac{1}{2}\sum_{k\not\sim
l}\mathrm{T}_{kl}^{c}\Upsilon_{kl}\xi(k)(1+2\xi(l))\bigl{(}{f(\xi^{k\to
l})-f(\xi)}\bigr{)}\bigl{(}{\widetilde{f}(\xi)^{\mathsf{q}-1}-\widetilde{f}(\xi^{k\to
l})^{\mathsf{q}-1}}\bigr{)}$
$\displaystyle-\widetilde{f}(\xi)^{\mathsf{q}-1}\sum_{k\not\sim
l}\mathrm{T}_{kl}^{c}\frac{(\lambda_{k}+\lambda_{l})\partial_{\lambda_{k}}g(\xi)}{n(\lambda_{k}-\lambda_{l})}\bigg{]}.$
(E.37)

We proceed to control the two terms in (E.37). To simplify the presentation,
for $k\not\sim l$ and a configuration $\xi$, we define the equivalence class
of configurations $[\xi]\equiv[\xi]_{k,l}$ as follows: $\xi^{\prime}\in[\xi]$
if and only if $|\xi^{\prime}|=|\xi|$ and
$\operatorname{supp}(\xi^{\prime}-\xi)\subset\\{\bm{\mathrm{e}}_{k},\bm{\mathrm{e}}_{l}\\}$.
Given $[\xi]$ and $0\leq\ell\leq\xi(k)+\xi(l)$, we define $\xi^{\ell}\in[\xi]$
to be the configuration with $\xi^{\ell}(k)=\ell$. Moreover, for any
$\ell\in[\\![{\xi(k)+\xi(l)}]\\!]$, we define the following functions:

$\displaystyle h_{\xi^{\ell},k,l}$
$\displaystyle:=\pi(\xi^{\ell})\xi^{\ell}(k)(1+2\xi^{\ell}(l))=\pi(\xi^{\ell-1})\xi^{\ell-1}(l)(1+2\xi^{\ell-1}(k)),$
$\displaystyle\quad\widehat{h}_{\xi^{\ell},k,l}$
$\displaystyle:=\pi(\xi^{\ell})\xi^{\ell}(k)=\pi(\xi^{\hat{\ell}})\xi^{\hat{\ell}}(l),\quad\text{with}\quad\hat{\ell}=\xi(k)+\xi(l)-\ell.$

Then, using the definition (4.17) and the anti-symmetry of the coefficient
${(\lambda_{k}+\lambda_{l})}/{(\lambda_{k}-\lambda_{l})}$, we can further
rewrite (E.37) as

$\displaystyle\frac{\mathsf{q}\theta_{-}\theta_{+}}{F+n^{-{\varepsilon}^{\prime}\mathsf{q}+L}}$
$\displaystyle\sum_{k\not\sim
l}\frac{1}{2}\mathrm{T}^{c}_{kl}\sum_{[\xi]}\sum_{\ell=1}^{\xi(k)+\xi(l)}\bigg{[}\Upsilon_{kl}h_{\xi^{\ell},k,l}\big{(}f(\xi^{\ell-1})-f(\xi^{\ell})\big{)}\big{(}\widetilde{f}(\xi^{\ell})^{\mathsf{q}-1}-\widetilde{f}(\xi^{\ell-1})^{\mathsf{q}-1}\big{)}$
$\displaystyle-\widehat{h}_{\xi^{\ell},k,l}\frac{\lambda_{k}+\lambda_{l}}{n(\lambda_{k}-\lambda_{l})}\left(g(\xi^{\ell})\frac{\phi^{\prime}(\lambda_{k},{\varepsilon})}{\phi(\lambda_{k},{\varepsilon})}\widetilde{f}(\xi^{\ell})^{\mathsf{q}-1}-g(\xi^{\hat{\ell}})\frac{\phi^{\prime}(\lambda_{l},{\varepsilon})}{\phi(\lambda_{l},{\varepsilon})}\widetilde{f}(\xi^{\hat{\ell}})^{\mathsf{q}-1}\right)\bigg{]}.$
(E.38)

Using the definitions of $\theta_{\pm}$ and (4.18), we obtain that

$\displaystyle\theta_{-}\theta_{+}g(\xi^{\ell})\frac{\phi^{\prime}(\lambda_{k},{\varepsilon})}{\phi(\lambda_{k},{\varepsilon})}$
$\displaystyle=\theta_{-}\theta_{+}\frac{g(\xi^{\ell})}{\phi(\lambda_{k},{\varepsilon})}\frac{\phi(\lambda_{k},{\varepsilon})-\phi(\lambda_{l},{\varepsilon})}{\lambda_{k}-\lambda_{l}}+\operatorname{O}\left(n^{2{\varepsilon}}(\lambda_{k}-\lambda_{l})\right)$
$\displaystyle=\theta_{-}\theta_{+}g(\xi^{\hat{\ell}})\frac{\phi^{\prime}(\lambda_{l},{\varepsilon})}{\phi(\lambda_{l},{\varepsilon})}+\operatorname{O}\left(n^{2{\varepsilon}}(\lambda_{k}-\lambda_{l})\right).$

By substituting this estimate into (E.38) and using (E.23) again, we deduce
that

$\displaystyle\eqref{hha4}=\frac{\mathsf{q}\theta_{-}\theta_{+}}{F+n^{-{\varepsilon}^{\prime}\mathsf{q}+L}}$
$\displaystyle\sum_{k\not\sim
l}\frac{1}{2}\mathrm{T}^{c}_{kl}\Upsilon_{kl}\sum_{[\xi]}\sum_{\ell=1}^{\xi(k)+\xi(l)}\bigg{[}h_{\xi^{\ell},k,l}\big{(}f(\xi^{\ell-1})-f(\xi^{\ell})\big{)}\big{(}\widetilde{f}(\xi^{\ell})^{\mathsf{q}-1}-\widetilde{f}(\xi^{\ell-1})^{\mathsf{q}-1}\big{)}$
$\displaystyle-2\widehat{h}_{\xi^{\ell},k,l}\frac{g(\xi^{\ell})}{\phi(\lambda_{k},{\varepsilon})}\big{(}\phi(\lambda_{k},{\varepsilon})-\phi(\lambda_{l},{\varepsilon})\big{)}\big{(}\widetilde{f}(\xi^{\ell})^{\mathsf{q}-1}-\widetilde{f}(\xi^{\hat{\ell}})^{\mathsf{q}-1}\big{)}\bigg{]}+\operatorname{O}(n^{-1+3{\varepsilon}}).$

Note
$\widetilde{f}(\xi^{\ell})^{\mathsf{q}-1}-\widetilde{f}(\xi^{\hat{\ell}})^{\mathsf{q}-1}$
can be expressed as a telescoping sum of terms of the form
$\widetilde{f}(\xi^{\ell^{\prime}})^{\mathsf{q}-1}-\widetilde{f}(\xi^{\ell^{\prime}-1})^{\mathsf{q}-1}$
for some $\ell^{\prime}$ between $\ell$ and $\hat{\ell}$. Therefore, there
exist a sequence of deterministic coefficients $c_{\xi,\ell,k,l}$ of order
$\operatorname{O}(1)$ such that we can rewrite the above equation as

$\displaystyle\eqref{hha4}=\frac{\mathsf{q}\theta_{-}\theta_{+}}{F+n^{-{\varepsilon}^{\prime}\mathsf{q}+L}}$
$\displaystyle\sum_{k\not\sim
l}\frac{1}{2}\mathrm{T}^{c}_{kl}\Upsilon_{kl}\sum_{[\xi]}\sum_{\ell=1}^{\xi(k)+\xi(l)}h_{\xi^{\ell},k,l}\Big{[}\big{(}f(\xi^{\ell-1})-f(\xi^{\ell})\big{)}+c_{\xi,\ell,k,l}\frac{g(\xi^{\ell})}{\phi(\lambda_{k},{\varepsilon})}\big{(}\phi(\lambda_{l},{\varepsilon})-\phi(\lambda_{k},{\varepsilon})\big{)}\Big{]}$
$\displaystyle\times\left(\widetilde{f}(\xi^{\ell})^{\mathsf{q}-1}-\widetilde{f}(\xi^{\ell-1})^{\mathsf{q}-1}\right)+\operatorname{O}(n^{-1+3{\varepsilon}})$
$\displaystyle=\frac{\mathsf{q}\theta_{-}\theta_{+}}{F+n^{-{\varepsilon}^{\prime}\mathsf{q}+L}}$
$\displaystyle\sum_{k\not\sim
l}\frac{1}{2}\mathrm{T}^{c}_{kl}\Upsilon_{kl}\sum_{[\xi]}\sum_{\ell=1}^{\xi(k)+\xi(l)}h_{\xi^{\ell},k,l}\Big{[}\big{(}f(\xi^{\ell-1})-f(\xi^{\ell})\big{)}+c_{\xi,\ell,k,l}\big{(}g(\xi^{\ell-1})-g(\xi^{\ell})\big{)}\Big{]}$
$\displaystyle\times\left(\widetilde{f}(\xi^{\ell})^{\mathsf{q}-1}-\widetilde{f}(\xi^{\ell-1})^{\mathsf{q}-1}\right)+\operatorname{O}(n^{-1+3{\varepsilon}}).$
(E.39)

By definition, we have
$f(\xi^{\ell-1})-f(\xi^{\ell})=g(\xi^{\ell-1})-g(\xi^{\ell})+\widetilde{f}(\xi^{\ell-1})-\widetilde{f}(\xi^{\ell}).$
Observe that if $\lvert f(\xi^{\ell-1})-f(\xi^{\ell})\rvert\geq\lvert
g(\xi^{\ell-1})-g(\xi^{\ell})\rvert\log n$, then the corresponding term in
(E.39) is negative, i.e.,

$\Big{[}\big{(}f(\xi^{\ell-1})-f(\xi^{\ell})\big{)}+c_{\xi,\ell,k,l}\big{(}g(\xi^{\ell-1})-g(\xi^{\ell})\big{)}\Big{]}\left(\widetilde{f}(\xi^{\ell})^{\mathsf{q}-1}-\widetilde{f}(\xi^{\ell-1})^{\mathsf{q}-1}\right)<0.$

Otherwise, based on (E.16), we can get that

$\displaystyle\theta_{-}\theta_{+}\Upsilon_{kl}\Big{[}\big{(}f(\xi^{\ell-1})-f(\xi^{\ell})\big{)}+c_{\xi,\ell,k,l}\big{(}g(\xi^{\ell-1})-g(\xi^{\ell})\big{)}\Big{]}\left(\widetilde{f}(\xi^{\ell})^{\mathsf{q}-1}-\widetilde{f}(\xi^{\ell-1})^{\mathsf{q}-1}\right)$
$\displaystyle=\operatorname{O}\left(n^{-1+2{\varepsilon}}(\log
n)^{2}\right)\cdot\theta_{-}\theta_{+}\left(\widetilde{f}(\xi^{\ell})^{\mathsf{q}-2}+\widetilde{f}(\xi^{\ell-1})^{\mathsf{q}-2}\right).$
(E.40)

(We remark that this is the step where it is convenient to define $g_{t}$ in
terms of $\lambda_{k}(t)$ as in (4.17).) By combining these facts, we find
that

$\eqref{eq_finalreducedform}=\mathcal{E}_{3}(t)-s_{+}^{(3)}(t),$

where $s_{+}^{(3)}(t)\geq 0$ and the variable $\mathcal{E}_{3}(t)$ satisfies
the following bound:

$\displaystyle|\mathcal{E}_{3}(t)|$ $\displaystyle\lesssim
n^{-1+3{\varepsilon}}+\frac{\mathsf{q}\theta_{-}\theta_{+}}{F+n^{-{\varepsilon}^{\prime}\mathsf{q}+L}}\sum_{\xi}\pi(\xi)\widetilde{f}(\xi)^{\mathsf{q}-2}\frac{n^{2{\varepsilon}}(\log
n)^{2}}{n}\sum_{k\not\sim
l}\frac{\xi(k)\eta^{2}}{(\lambda_{k}-\lambda_{l})^{2}+\eta^{2}}$
$\displaystyle=n^{-1+3{\varepsilon}}+\frac{\mathsf{q}\theta_{-}\theta_{+}}{F+n^{-{\varepsilon}^{\prime}\mathsf{q}+L}}\sum_{\xi}\pi(\xi)\widetilde{f}(\xi)^{\mathsf{q}-2}n^{2{\varepsilon}}(\log
n)^{2}\sum_{k}\xi(k)\eta\operatorname{Im}\underline{\mathsf{m}}_{n}(\lambda_{k}+\mathrm{i}\eta)$
$\displaystyle\lesssim n^{-1+3{\varepsilon}}+n^{4{\varepsilon}}(\log
n)^{2}\eta\leq n^{-2/3+4{\varepsilon}+\mathsf{c}}\quad\text{w.h.p.}$ (E.41)

Here, in the second step, we used the spectral decomposition of
$\underline{\mathsf{m}}_{n}$, and in the third step, we applied the local law
(B.35) to $\underline{\mathsf{m}}_{n}(\lambda_{k}+\mathrm{i}\eta)$ along with
the rigidity of $\lambda_{k}$ given by (B.37). This concludes (E.36).

The term $A_{1}$. Finally, we prove the following inequality for sufficiently
large $\mathsf{q}$:

$\frac{\mathsf{q}}{F+n^{-{\varepsilon}^{\prime}\mathsf{q}+L}}\sum_{\xi}\pi(\xi)\widetilde{f}(\xi)^{\mathsf{q}-1}A_{1}\leq-C(t)\frac{F}{F+n^{-{\varepsilon}^{\prime}\mathsf{q}+L}}+\operatorname{O}_{\prec}(1),$
(E.42)

where $C(t)$ is a positive variable satisfying $C(t)\gtrsim\eta^{-1/2}$ with
high probability. With the definition of $A_{1}$ in (E.32), we can express the
LHS of (E.42) as follows:

$\displaystyle~{}\frac{\mathsf{q}}{F+n^{-{\varepsilon}^{\prime}\mathsf{q}+L}}\sum_{\xi}\pi(\xi)\widetilde{f}(\xi)^{\mathsf{q}-1}A_{1}$
$\displaystyle=$
$\displaystyle~{}\frac{\mathsf{q}\theta_{-}\theta_{+}}{2n(F+n^{-{\varepsilon}^{\prime}\mathsf{q}+L})}\sum_{\xi}\pi(\xi)\widetilde{f}(\xi)^{\mathsf{q}-1}\sum_{k\not\sim
l}\frac{(\lambda_{k}+\lambda_{l})\xi(k)(1+2\xi(l))}{(\lambda_{k}-\lambda_{l})^{2}+\eta^{2}}\big{(}f(\xi^{k\to
l})-f(\xi)\big{)}.$ (E.43)

Recalling the definitions of $f$ and $\xi^{k\to l}$ in (4.22) and (4.23), we
can express $f(\xi^{k\to l})$ as

$\displaystyle f(\xi^{k\to
l})=\mathbb{E}\bigg{[}\frac{p|\langle{\bf{v}},{\bf{u}}_{l}\rangle|^{2}}{1+2\xi(l)}\prod_{j}\frac{(p|{\bf{v}}^{\top}{\bf{u}}_{j}|^{2})^{\xi(j)-\mathbf{1}_{j=k}}}{a_{\xi(j)-\mathbf{1}_{j=k}}}\bigg{|}\mathcal{F}_{t}\bigg{]},$
(E.44)

where $a_{n}:=(2n-1)!!$ for $n\in\mathbb{N}$. Utilizing the spectral
decomposition of $\mathcal{G}_{1,t}(z)=(W(t)W(t)^{\top}-z)^{-1}$ with
$z=z_{k}:=\lambda_{k}+\mathrm{i}\eta$, we have

$\displaystyle\sum_{l:l\not\sim
k}\frac{(\lambda_{k}+\lambda_{l})|\langle{\bf{v}},{\bf{u}}_{l}\rangle|^{2}}{(\lambda_{l}-\lambda_{k})^{2}+\eta^{2}}=$
$\displaystyle~{}\sum_{l}\frac{2\lambda_{k}|\langle{\bf{v}},{\bf{u}}_{l}\rangle|^{2}}{(\lambda_{l}-\lambda_{k})^{2}+\eta^{2}}-\frac{2\lambda_{k}|\langle{\bf{v}},{\bf{u}}_{k}\rangle|^{2}}{\eta^{2}}+\sum_{l:l\not\sim
k}\frac{(\lambda_{l}-\lambda_{k})|\langle{\bf{v}},{\bf{u}}_{l}\rangle|^{2}}{(\lambda_{l}-\lambda_{k})^{2}+\eta^{2}}$
$\displaystyle=$
$\displaystyle~{}\frac{2\lambda_{k}}{\eta}\operatorname{Im}\left({\bf{v}}^{\top}\mathcal{G}_{1,t}(z_{k}){\bf{v}}\right)+\operatorname{Re}\left({\bf{v}}^{\top}\mathcal{G}_{1,t}(z_{k}){\bf{v}}\right)+\operatorname{O}_{\prec}(n^{-1}\eta^{-2})$
$\displaystyle=$
$\displaystyle~{}\frac{2\lambda_{k}}{\eta}\operatorname{Im}\left({\bf{v}}^{\top}\widetilde{\Pi}_{t}(z_{k}){\bf{v}}\right)+\operatorname{O}_{\prec}\left(\eta^{-1}\Psi_{t}(z_{k})\right).$

Here, in the second step, we applied the delocalization estimate (B.38) to
$|\langle{\bf{v}},{\bf{u}}_{k}\rangle|^{2}$, and in the third step, we applied
the local law (B.35) to ${\bf{v}}^{\top}\mathcal{G}_{1,t}(z_{k}){\bf{v}}$ with
$\widetilde{\Pi}_{t}(z):=-[z(1+m_{t}(z))\Lambda_{t}]^{-1}$. Similarly, by
observing that $\xi(l)$ is nonzero for at most $L$ many $l$’s, and using the
spectral decomposition of $\underline{\mathsf{m}}_{n}$ and the local law
(B.42), we obtain that

$\frac{1}{p}\sum_{l:l\not\sim
k}\frac{(\lambda_{k}+\lambda_{l})(1+2\xi(l))}{(\lambda_{l}-\lambda_{k})^{2}+\eta^{2}}=\frac{2\lambda_{k}}{\eta}\operatorname{Im}m_{1,t}(z_{k})+\operatorname{O}_{\prec}(n^{-1}\eta^{-2}).$

Plugging the above three equations into (E.43) and applying the delocalization
estimate (B.38) to the factors $p|\langle{\bf{v}},{\bf{u}}_{j}\rangle|^{2}$ in
$f(\xi)$ and (E.44), we obtain that

$\displaystyle(\ref{eq_finalparttobeestimated})=$
$\displaystyle~{}\frac{c_{n}\mathsf{q}\theta_{-}\theta_{+}}{F+n^{-{\varepsilon}^{\prime}\mathsf{q}+L}}\sum_{\xi}\pi(\xi)\widetilde{f}(\xi)^{\mathsf{q}-1}\sum_{k}\frac{\lambda_{k}\xi(k)}{\eta}$
$\displaystyle\times\left[\operatorname{Im}\left({\bf{v}}^{\top}\widetilde{\Pi}_{t}(z_{k}){\bf{v}}\right)f(\xi-{\bf
e}_{k})-\operatorname{Im}m_{1,t}(z_{k})f(\xi)+\operatorname{O}_{\prec}\left(\Psi_{t}(z_{k})\right)\right].$
(E.45)

By the induction hypothesis, when $K=L-1$, the estimate (E.6) holds with high
probability for large even integer $\mathsf{p}\geq\mathsf{p}_{0}$, implying
that

$|\widetilde{f}(\xi-\mathbf{e}_{k})|\lesssim
n^{-{\varepsilon}+(L-1)/\mathsf{p}}\quad\text{w.h.p.}\ \Rightarrow\
|\widetilde{f}(\xi-\mathbf{e}_{k})|\prec n^{-{\varepsilon}}.$ (E.46)

On the other hand, by (B.17) and (B.40), we have

$\operatorname{Im}m_{1,t}(z_{k})\asymp\operatorname{Im}m_{t}(z_{k})\gtrsim\min\left\\{\sqrt{\kappa_{k}+\eta},\frac{\eta}{\sqrt{\kappa_{k}+\eta}}\right\\}\gtrsim\sqrt{\eta}\quad\text{w.h.p.},$
(E.47)

where $\kappa_{k}$ is defined as
$\kappa_{k}:=\operatorname{dist}(\lambda_{k},\operatorname{supp}(\varrho_{t}))$,
and we have used the fact that $\kappa_{k}\prec n^{-2/3}$ according to the
rigidity estimate (B.37). With (E.47), for $\eta=n^{-2/3+\mathsf{c}/2}$, we
obtain that

$\frac{\Psi_{t}(z_{k})}{\operatorname{Im}m_{1,t}(z_{k})}\lesssim\left(\frac{1}{n\eta\operatorname{Im}m_{1,t}(z_{k})}\right)^{1/2}+\frac{1}{n\eta\operatorname{Im}m_{1,t}(z_{k})}\lesssim\left(\frac{1}{n\eta^{3/2}}\right)^{1/2}+\frac{1}{n\eta^{3/2}}\lesssim
n^{-3\mathsf{c}/8}.$ (E.48)

With the estimates (E.46) and (E.48), we can simplify (E.45) as follows:

$\displaystyle(\ref{eq_finalparttobeestimated})=\frac{c_{n}\mathsf{q}\theta_{-}\theta_{+}}{F+n^{-{\varepsilon}^{\prime}\mathsf{q}+L}}$
$\displaystyle\sum_{\xi}\pi(\xi)\widetilde{f}(\xi)^{\mathsf{q}-1}\sum_{k}\lambda_{k}\xi(k)\frac{\operatorname{Im}m_{1,t}(z_{k})}{\eta}$
$\displaystyle\times\Big{(}\varphi_{t}(z_{k})g(\xi-{\bf
e}_{k})-f(\xi)+\operatorname{O}_{\prec}\big{(}n^{-{\varepsilon}}+n^{-3\mathsf{c}/8}\big{)}\Big{)},$
(E.49)

where
$\varphi_{t}(z_{k})\equiv\varphi_{t}({\bf{v}},{\bf{v}},z_{k})={\operatorname{Im}({\bf{v}}^{\top}\widetilde{\Pi}_{t}(z_{k}){\bf{v}})}/{\operatorname{Im}m_{1,t}(z_{k})}$
as defined in (4.14). By using (B.18), (C.12), (C.32), and the definition of
$\phi(\lambda_{k},{\varepsilon})$, we can check that

$\displaystyle\varphi_{t}({\bf{v}},{\bf{v}},z_{k})-\phi(\lambda_{k},{\varepsilon})$
$\displaystyle=\left[\varphi_{t}({\bf{v}},{\bf{v}},z_{k})-\varphi_{0}({\bf{v}},{\bf{v}},\lambda_{k})\right]+\left[\varphi_{0}({\bf{v}},{\bf{v}},\lambda_{k})-\phi(\lambda_{k},{\varepsilon})\right]$
$\displaystyle\lesssim\sqrt{t+\eta}+n^{-{\varepsilon}}\leq
2n^{-{\varepsilon}},$ (E.50)

where in the last step, we used the fact that $\sqrt{t+\eta}\lesssim
n^{-1/6+\mathsf{c}/2}+n^{-1/3+\mathsf{c}/4}\leq n^{-{\varepsilon}}$ for
$\mathsf{c}<1/6$ and ${\varepsilon}<1/12$. Plugging (E.50) into (E.1), and in
light of the definition of $g(\xi)$ in (4.17), we find that

$\displaystyle(\ref{eq_finalparttobeestimated})=\frac{c_{n}\mathsf{q}\theta_{-}\theta_{+}}{F+n^{-{\varepsilon}^{\prime}\mathsf{q}+L}}\sum_{\xi}\pi(\xi)\sum_{k}\lambda_{k}\xi(k)\frac{\operatorname{Im}m_{1,t}(z_{k})}{\eta}\left(g(\xi)-f(\xi)+\operatorname{O}_{\prec}\big{(}n^{-{\varepsilon}}+n^{-3\mathsf{c}/8}\big{)}\right)\widetilde{f}(\xi)^{\mathsf{q}-1}$
$\displaystyle=\frac{-c_{n}\mathsf{q}\theta_{-}\theta_{+}}{F+n^{-{\varepsilon}^{\prime}\mathsf{q}+L}}\sum_{\xi}\pi(\xi)\sum_{k}\lambda_{k}\xi(k)\frac{\operatorname{Im}m_{1,t}(z_{k})}{\eta}\left[\widetilde{f}(\xi)^{\mathsf{q}}+\operatorname{O}_{\prec}\left(n^{-\tau}\big{(}\widetilde{f}(\xi)^{\mathsf{q}}+n^{-\mathsf{q}\delta}\big{)}\right)\right]$
$\displaystyle\leq\frac{-c\eta^{-1/2}F}{F+n^{-{\varepsilon}^{\prime}\mathsf{q}+L}}+\operatorname{O}_{\prec}\left(\frac{n^{-\mathsf{q}\delta-\tau+L}\eta^{-1}}{F+n^{-{\varepsilon}^{\prime}\mathsf{q}+L}}\right)=\frac{-c\eta^{-1/2}F}{F+n^{-{\varepsilon}^{\prime}\mathsf{q}+L}}+\operatorname{O}_{\prec}\left(1\right)$

for a constant $c>0$. In the second step, we used (E.23) with ${\varepsilon}$
replaced by a constant
${\varepsilon}^{\prime}<\delta<\min\\{{\varepsilon},3\mathsf{c}/8\\}$, and we
denote $\tau:=\min\\{{\varepsilon},3\mathsf{c}/8\\}-\delta$; in the third
step, we used $\sqrt{\eta}\lesssim\operatorname{Im}m_{1,t}(z_{k})\lesssim 1$
by (E.47) and (B.17); in the last step, we choose $\mathsf{q}$ to be
sufficiently large (depending on $\delta-{\varepsilon}^{\prime}$) such that
$n^{-\mathsf{q}\delta-\tau}\eta^{-1}\leq
n^{-{\varepsilon}^{\prime}\mathsf{q}}$. This proves (E.42) and completes the
proof of Lemma E.2.

##  Appendix F Proof of the results in Section 4.3

###  F.1 Proof of Lemma 4.11

In this subsection, we prove Lemma 4.11 using the strategies outlined in [13,
44]. A key component of the proof relies on the level repulsion estimate. In
the existing literature, the level repulsion estimate for sample covariance
matrices has only been established around the edge [13] or inside the bulk via
the universality of bulk eigenvalue statistics [61]. However, for our specific
applications, we require the level repulsion estimate to hold over the entire
spectrum. Recall the definition of $\mathcal{W}_{t}$ from (4.8), and denote
the eigenvalues of $\mathcal{Q}_{t}=\mathcal{W}_{t}\mathcal{W}_{t}^{\top}$ as
$\\{\bar{\lambda}_{k}(t)\\}$. Given an interval $I\subset\mathbb{R}_{+}$, we
define the counting function

$\overline{\mathsf{N}}_{t}(I)=|\\{i\in[\\![{1,\mathsf{K}}]\\!]:\bar{\lambda}_{i}(t)\in
I\\}|.$

Similarly, we define the counting function ${\mathsf{N}}_{t}$ for the
eigenvalues $\\{\lambda_{k}(t)\\}$ of the matrix DBM defined in (4.4):

${\mathsf{N}}_{t}(I)=|\\{i\in[\\![{1,\mathsf{K}}]\\!]:\lambda_{i}(t)\in
I\\}|.$

Recall that $\gamma_{k}(t)$ represents the quantiles of $\varrho_{t}$ as
defined below (B.33), and $\Delta_{k}$ is defined in (4.29).

######  Proposition F.1 (Level repulsion).

Under the setting of Lemma 4.2, consider a fixed $t\asymp
n^{-1/3+\mathsf{c}}$. There exists a small constant $\delta_{0}\in(0,2/3)$
such that the following statements hold. Define the disjoint union of
intervals $I_{\delta_{0}}^{+}(t)$ as

$I_{\delta_{0}}^{+}(t):=\cup_{k=1}^{q}[a_{2k}(t)-n^{-\delta_{0}},a_{2k-1}(t)+n^{-\delta_{0}}],$

and define
$k_{E}:=\mathrm{arg\,min}\\{|\gamma_{k}(t)-E|:k\in[\\![{\mathsf{K}}]\\!]\\}$
for $E\in I_{\delta_{0}}^{+}$. Then, for every constant
$\delta\in(0,\delta_{0}],$ we can find a constant
$\varepsilon\equiv\varepsilon(\delta)>0$ such that

$\mathbb{P}\left({\mathsf{N}}_{t}\left([E-\Delta_{k_{E}}n^{-\delta},E+\Delta_{k_{E}}n^{-\delta}]\right)\geqslant
2\right)\leqslant n^{-\delta-\varepsilon},\quad\forall E\in
I_{\delta_{0}}^{+}.$ (F.1)

The same estimate also applies to $\overline{\mathsf{N}}_{t}$.

###### Proof.

An analogous level repulsion estimate inside the bulk has been established for
the matrix DBM of Wigner ensembles in [47, Theorem 5.1] and [11, Lemma B.1],
with the key inputs being the local laws and the eigenvalue rigidity estimate.
With the help of our Lemmas B.9 and B.11, we can readily extend the arguments
from [47, 11] to our rectangular matrix DBM $Q(t)$ defined in (4.4).
Specifically, by following the reasoning presented in [47, Theorem 5.1], we
can establish that (F.1) holds for
$E\in\cup_{k=1}^{q}[a_{2k}+\kappa,a_{2k-1}-\kappa]$, where $\kappa$ is an
arbitrarily small positive constant. Furthermore, as explained in the proof of
[11, Lemma B.1], the argument in [47, Theorem 5.1] can be extended a little
bit to show that there exists a sufficiently small constant ${\varepsilon}>0$
so that (F.1) holds for
$E\in\cup_{k=1}^{q}[a_{2k}+n^{-{\varepsilon}},a_{2k-1}-n^{-{\varepsilon}}]$.
Since the extension from Wigner ensembles to sample covariance ensembles is
standard, we omit the complete details here.

To establish (F.1) for $E$ around the spectral edges, specifically for
$E\in\cup_{k=1}^{q}([a_{2k}-n^{-\delta_{0}},a_{2k}+n^{-{\varepsilon}}]\cup[a_{2k-1}-n^{-{\varepsilon}},a_{2k-1}+n^{-\delta_{0}}])$,
we can refer to [13, Proposition 6.3]. Although this result is stated for the
special case with $\Lambda_{0}=I$ in [13], the same argument can be extended
to our more general setting. More precisely, as explained in Lemma 6.4 of
[13], we can show that (F.1) holds for the Gaussian ensemble when $X$ is a
Gaussian random matrix. (The proof of [13, Lemma 6.4] relies on an analysis of
the joint eigenvalue probability density, employing a method developed in the
proof of [14, Theorem 3.2]. While [14, Theorem 3.2] is presented for the
Gaussian $\beta$ ensemble, the same proof carries over almost verbatim to our
setting, where the eigenvalues follow a Laguerre $\beta$ ensemble.) Next, as
explained in the proof of [13, Lemma 6.5], we can establish the level
repulsion estimate (F.1) for $Q(t)$ through a Green’s function comparison
argument between the case with Gaussian $X$ and that with a non-Gaussian $X$
in $W(0)$. Once again, since this extension is standard and straightforward,
we omit the complete details here. This concludes the proof of (F.1) for $E$
near the spectral edges.

By (F.1), we know that the level repulsion estimate applies to the eigenvalues
of $\mathcal{Q}^{G}_{t}:=\mathcal{W}^{G}_{t}(W^{G}_{t})^{\top}$ with
$W^{G}_{t}:=D_{t}^{1/2}V^{\top}X^{G}$. To further prove this estimate for
$\overline{\mathsf{N}}_{t}$, we can once again employ the Green’s function
comparison argument between $\mathcal{Q}^{G}_{t}$ and $\mathcal{Q}_{t}$, as
explained in Lemma 6.5 of [13]. We omit the details here. ∎

With Proposition F.1, we now proceed to complete the proof of Lemma 4.11 by
employing the method developed in [13, 44].

###### Proof of Lemma 4.11.

We only present the proof for $Q(t)$, while the proof for $\mathcal{Q}_{t}$ is
the same. To simplify the presentation, we prove Lemma 4.11 in the special
case where $\theta=\theta(p|\langle{\bf{v}},{\bf{u}}_{i+r}(t)\rangle|^{2})$.
For brevity, we will use the abbreviations
${\bf{u}}_{i}\equiv{\bf{u}}_{i}(t)$, $Q\equiv Q(t)$, and
$\widetilde{\mathcal{G}}_{1}(z)\equiv\widetilde{\mathcal{G}}_{1,t}(z):=(Q(t)-z)^{-1}$
for a fixed $t\asymp n^{-1/3+\mathsf{c}}$. Recall the definition of
$\Delta_{i}$ in (4.29) and the definition of $\eta_{i}$ in (4.30). We now
present the first ingredient of the proof.

######  Lemma F.2.

Under the assumptions of Lemma 4.11, there exists a small constant
${\varepsilon}>0$ such that the following estimate holds for sufficiently
small constants $\delta_{1}\equiv\delta_{1}({\varepsilon})$,
$\delta_{2}\equiv\delta_{2}({\varepsilon},\delta_{1})$, and $\nu>0$:

$\mathbb{E}\theta(p|\langle{\bf{v}},{\bf{u}}_{i+r}(t)\rangle|^{2})=\mathbb{E}\theta\left(\frac{p}{\pi}\int_{I_{i}}\operatorname{Im}\big{[}{\bf{v}}^{\top}\widetilde{\mathcal{G}}_{1}(E+\mathrm{i}\eta_{i}){\bf{v}}\big{]}\chi_{i}(E)\mathrm{d}E\right)+\operatorname{O}(n^{-\nu}),\quad
i\in[\\![{1,\mathsf{K}-r}]\\!],$

where $I_{i}$ is defined in (4.31) and
$\chi_{i}(E):=\mathbf{1}(\lambda_{i+1}<E_{i}^{-}\leqslant\lambda_{i})$. Here,
$E_{i}^{-}=E-n^{\delta_{1}}\eta_{i}$ is defined in (4.30), and we adopt the
convention $\lambda_{\mathsf{K}+1}=0$.

###### Proof.

Our proof follows a similar approach to that of [44, Lemma 3.1], [13, Lemma
7.1], and [26, Lemma 3.2]. The only difference is that we utilize the level
repulsion estimate from Proposition F.1 over the entire spectrum, whereas the
aforementioned works rely on the level repulsion estimate either near the
edges or inside the bulks. We omit the detailed exposition here. ∎

For the second ingredient of the proof, we show that the indicator function
$\chi_{i}(E)$ can be approximated by the smooth functions
$\mathfrak{q}_{i}(x)$ and $\mathfrak{f}_{i}(x)$ defined in (4.28) and (4.32)
up to a negligible error.

######  Lemma F.3.

Under the assumptions of Lemma 4.11, there exists a small constant
${\varepsilon}>0$ such that the following estimate holds for sufficiently
small constants $\delta_{1}\equiv\delta_{1}({\varepsilon})$,
$\delta_{2}\equiv\delta_{2}({\varepsilon},\delta_{1})$, and $\nu>0$:

$\mathbb{E}\theta\left(\frac{p}{\pi}\int_{I_{i}}\operatorname{Im}\big{[}{\bf{v}}^{\top}\widetilde{\mathcal{G}}_{1}(E+\mathrm{i}\eta_{i}){\bf{v}}\big{]}\chi_{i}(E)\mathrm{d}E\right)=\mathbb{E}\theta\left(\frac{p}{\pi}\int_{I_{i}}\operatorname{Im}\big{[}{\bf{v}}^{\top}\widetilde{\mathcal{G}}_{1}(E+\mathrm{i}\eta_{i}){\bf{v}}\big{]}\mathfrak{q}_{i}(\operatorname{Tr}\mathfrak{f}_{i}(Q))\mathrm{d}E\right)+\operatorname{O}(n^{-\nu}),$

where $i\in[\\![{1,\mathsf{K}-r}]\\!]$.

###### Proof.

Using the rigidity of eigenvalues in (B.37), we can obtain that with high
probability,

$\displaystyle
p\int_{I_{i}}\operatorname{Im}\big{[}{\bf{v}}^{\top}\widetilde{\mathcal{G}}_{1}(E+\mathrm{i}\eta_{i}){\bf{v}}\big{]}\chi_{i}(E)\mathrm{d}E$
$\displaystyle=p\int_{I_{i}}\operatorname{Im}\big{[}{\bf{v}}^{\top}\widetilde{\mathcal{G}}_{1}(E+\mathrm{i}\eta_{i}){\bf{v}}\big{]}\mathbf{1}(\mathsf{N}([E^{-}_{i},E^{+}])=i)\mathrm{d}E$
$\displaystyle=p\int_{I_{i}}\operatorname{Im}\big{[}{\bf{v}}^{\top}\widetilde{\mathcal{G}}_{1}(E+\mathrm{i}\eta_{i}){\bf{v}}\big{]}\mathfrak{q}_{i}[\operatorname{Tr}(\chi_{E}(Q))]\mathrm{d}E,$

where $E^{+}$ is defined in (4.30) and we denote
$\chi_{E}(x):=\mathbf{1}(x\in[E_{i}^{-},E^{+}])$. Next, by utilizing
Proposition F.1 and following the proof of (4.11) in [26], we can show that

$\mathbb{E}\theta\left(\frac{p}{\pi}\int_{I_{i}}\operatorname{Im}\big{[}{\bf{v}}^{\top}\widetilde{\mathcal{G}}_{1}(E+\mathrm{i}\eta_{i}){\bf{v}}\big{]}\mathfrak{q}_{i}[\operatorname{Tr}(\chi_{E}(Q))]\mathrm{d}E\right)=\mathbb{E}\theta\left(\frac{p}{\pi}\int_{I_{i}}\operatorname{Im}\big{[}{\bf{v}}^{\top}\widetilde{\mathcal{G}}_{1}(E+\mathrm{i}\eta_{i}){\bf{v}}\big{]}\mathfrak{q}_{i}(\operatorname{Tr}\mathfrak{f}_{i}(Q))\mathrm{d}E\right)+\operatorname{O}(n^{-\nu}).$

Combining the above two equations completes the proof. ∎

We observe that Lemmas F.2 and F.3 together establish Lemma 4.11 when $L=1$.
The proof for the general case of $L>1$ is analogous, as explained in Section
4 of [44]. We omit the details here. ∎

###  F.2 Proof of Lemma 4.12

In this subsection, we prove Lemma 4.12 using a continuous comparison argument
developed in [45, Section 7]. Recall the sequence of interpolating matrices
$W^{s}_{t}\in\mathbb{R}^{p\times n}$ defined in (4.10). For the sake of
simplicity, we denote $Y^{s}\equiv Y^{s}_{t}:=D_{t}^{-1/2}W^{s}_{t}$, i.e.,

$Y^{s}_{i\mu}=\sum_{j=1}^{2p}U^{1}_{ij}\chi^{s}_{j\mu}\mathsf{X}_{j\mu}+\sum_{j=1}^{2p}U^{0}_{ij}(1-\chi^{s}_{j\mu})\mathsf{X}_{j\mu},\quad
i\in\mathcal{I}_{1},\ \mu\in\mathcal{I}_{2},$

where $\\{\chi^{s}_{j\mu}\\}_{j\in[\\![{2p}]\\!],\mu\in[\\![{n}]\\!]}$ are
i.i.d. Bernoulli$(s)$ random variables. Alternatively, we can express $Y^{s}$
as

$Y^{s}=\sum_{j=1}^{2p}\sum_{\mu=1}^{n}Z^{s}_{(j\mu)},\quad\text{with}\quad
Z^{s}_{(j\mu)}:=\mathsf{X}_{j\mu}\Delta^{s}_{(j\mu)},\quad\Delta^{s}_{(j\mu)}:=\chi^{s}_{j\mu}\left(U^{1}\mathbf{e}_{j\mu}\right)+(1-\chi^{s}_{j\mu})\left(U^{0}\mathbf{e}_{j\mu}\right).$
(F.2)

Here, $\mathbf{e}_{j\mu}:=\mathbf{e}_{j}\mathbf{e}_{\mu}^{\top}$ represents a
$2p\times n$ matrix with only one non-zero entry at the $(j,\mu)$-th position.
Furthermore, given a matrix $\Gamma\in\mathbb{R}^{p\times n}$, we define

$Y^{s,\Gamma}_{(j\mu)}:=\sum_{(i,\nu)\neq(j,\mu)}Z_{(i\nu)}^{s}+\Gamma.$ (F.3)

In particular, using this notation, we have
$Y^{s}-Y^{s,\mathbf{0}}_{(j\mu)}=Z_{(j\mu)}^{s}$. Similar to (B.13) and
(B.14), we will use the following linearized resolvent:

$G^{s}_{t}(z):=\begin{pmatrix}-D_{t}^{-1}&Y^{s}\\\
(Y^{s})^{\top}&-z\end{pmatrix}^{-1}=\begin{pmatrix}zD_{t}^{1/2}{\mathcal{G}}_{1,t}^{s}(z)D_{t}^{1/2}&D_{t}Y^{s}\mathcal{G}^{s}_{2,t}\\\
\mathcal{G}^{s}_{2,t}(Y^{s})^{\top}D_{t}&\mathcal{G}^{s}_{2,t}\end{pmatrix},$
(F.4)

where
$\mathcal{G}_{1,t}^{s}:=(D_{t}^{1/2}Y^{s}(Y^{s})^{\top}D_{t}^{1/2}-z)^{-1}$
and $\mathcal{G}_{2,t}^{s}:=((Y^{s})^{\top}D_{t}Y^{s}-z)^{-1}$. For
$s\in\\{0,1\\}$, the local laws for the resolvent $G^{s}$ have already been
established in Lemma B.2 (when $s=0$) and Lemma B.11 (when $s=1$). These local
laws can also be extended to all $s\in(0,1)$, as summarized in the following
lemma.

######  Lemma F.4.

Suppose the assumptions of Lemma 4.12 hold. Define $\widetilde{\Pi}_{t}$ and
$\Psi_{t}$ as in Section B.3 with $\Lambda=D_{0}$ and $\Lambda_{0}=D_{00}$.
Then, there exists a constant $c>0$ such that the following local law holds
uniformly in $z\in\mathbf{D}^{\prime}(c,\tau,t,n)$ for any constant $\tau>0$
and deterministic unit vectors ${\bf{u}},{\bf{v}}\in\mathbb{R}^{p+n}$:

$\left|\mathbf{u}^{\top}G^{s}_{t}(z)\mathbf{v}-\mathbf{u}^{\top}\widetilde{\Pi}_{t}(z)\mathbf{v}\right|\prec\Psi_{t}(z).$
(F.5)

We now prove this lemma using the continuous comparison argument presented in
[45, 68]. In the interest of brevity, we will provide a sketch of the key
points, highlighting the main differences in the argument. For a more
comprehensive understanding, readers can refer to Section 6 of [68].

###### Proof of Lemma F.4.

For any $\eta\geq N^{-1+\tau}$ and a sufficiently small constant
$\delta\in(0,1/100)$, we define

$\eta_{l}:=\eta N^{\delta l}\text{ for }\ l=0,...,K-1,\ \ \ \eta_{K}:=1,$
(F.6)

where $K\equiv K(\eta):=\max\left\\{l\in\mathbb{N}:\eta
N^{\delta(l-1)}<1\right\\}.$ Since $z\mapsto
G_{t}^{s}(z)-\widetilde{\Pi}_{t}(z)$ is Lipschitz continuous in
$\mathbf{D}^{\prime}\equiv\mathbf{D}^{\prime}(c,\tau,t,n)$ with Lipschitz
constant of order $\operatorname{O}(\eta^{-2})=\operatorname{O}(n^{2})$, it
suffices to show that (F.5) holds for all $z$ in a discrete but suitably dense
subset ${\mathbf{S}}\subset\mathbf{D}^{\prime}$. Specifically, we choose
$\mathbf{S}$ to be an $n^{-10}$-net of $\mathbf{D}^{\prime}$ such that
$|\mathbf{S}|\lesssim n^{20}$ and

$E+\mathrm{i}\eta\in\mathbf{S}\Rightarrow
E+\mathrm{i}\eta_{l}\in\mathbf{S}\text{ for }l=1,\ldots,K(\eta).$

Similar to [68, Section 6], our proof relies on an induction argument based on
two scale-dependent properties, denoted as ($\bm{\mathrm{A}}_{d}$) and
($\bm{\mathrm{C}}_{d}$), formulated on the subsets
${\bm{\mathrm{S}}}_{d}\mathrel{\vbox{\hbox{.}\hbox{.}}}=\bigl{\\{}{z\in{\bm{\mathrm{S}}}\mathrel{\vbox{\hbox{.}\hbox{.}}}\operatorname{Im}z\geqslant
n^{-\delta d}}\bigr{\\}},$ $d\in\mathbb{N}$.

  * ($\bm{\mathrm{A}}_{d}$)

For all $z\in{\bm{\mathrm{S}}}_{d}$, we have that for any deterministic unit
vectors ${\bf{u}},{\bf{v}}\in\mathbb{R}^{p+n}$,

$\big{|}\mathbf{u}^{\top}G^{s}_{t}(z)\mathbf{v}-\mathbf{u}^{\top}\widetilde{\Pi}_{t}(z)\mathbf{v}\big{|}\prec
1.$ (F.7)

  * ($\bm{\mathrm{C}}_{d}$)

For all $z\in{\bm{\mathrm{S}}}_{d}$, we have that for any deterministic unit
vectors ${\bf{u}},{\bf{v}}\in\mathbb{R}^{p+n}$,

$\big{|}\mathbf{u}^{\top}G^{s}_{t}(z)\mathbf{v}-\mathbf{u}^{\top}\widetilde{\Pi}_{t}(z)\mathbf{v}\big{|}\prec
n^{24\delta}\Psi(z).$ (F.8)

It is clear that property ($\mathbf{C}_{d}$) implies property
($\mathbf{A}_{d}$) if $\delta$ is chosen sufficiently small (depending on
$\tau$) such that $n^{24\delta}\Psi(z)\ll 1$. On the other hand, if we can
show that property ($\mathbf{A}_{d-1}$) implies property ($\mathbf{C}_{d}$)
for all $d\leqslant\delta^{-1}$, then we can conclude that (F.8) holds for all
$z\in\mathbf{S}$ by induction on $d$. By polarization, it suffices to prove
that when property (${\bm{\mathrm{A}}}_{d-1}$) holds, we have

$\big{|}\mathbf{v}^{\top}G^{s}_{t}(z)\mathbf{v}-\mathbf{v}^{\top}\widetilde{\Pi}_{t}(z)\mathbf{v}\big{|}\prec
n^{24\delta}\Psi(z),$ (F.9)

for all $z\in\mathbf{S}_{d}$ and any deterministic unit vector
$\mathbf{v}\in{\mathbb{C}}^{p+n}$. In fact, we can derive the more general
bound (F.8) by applying (F.9) to the vectors $\mathbf{u}+\mathbf{v}$ and
$\mathbf{u}+\mathrm{i}\mathbf{v}$, respectively.

By utilizing Markov’s inequality, showing (F.9) requires proving that for any
fixed large $\mathsf{q}\in 2\mathbb{N}$ and $d\leqslant\delta^{-1},$ there
exists a constant $C>0$ such that

$|\mathbb{E}F^{{\bf{v}},\,\mathsf{q}}(Y^{s})|\leq
C(n^{24\delta}\Psi)^{\mathsf{q}}$ (F.10)

uniformly for all $z\in\mathbf{S}_{d}$. Here, $F^{{\bf{v}},\mathsf{q}}(Y^{s})$
is defined as

$F^{{\bf{v}},\,\mathsf{q}}(Y^{s})\equiv
F^{{\bf{v}},\,\mathsf{q}}(Y^{s},z):=\big{|}{\bf{v}}^{\top}G_{t}^{s}{\bf{v}}-{\bf{v}}^{\top}\widetilde{\Pi}_{t}{\bf{v}}\big{|}^{\mathsf{q}}.$
(F.11)

To establish (F.10), we follow the arguments presented in Sections 6.1 and 6.2
of [68]. First, we see that (F.5) and consequently (F.10) hold when $s=0$ by
Lemma B.11. Second, similar to [45, Lemma 7.9], we can employ the fundamental
theorem of calculus to obtain that

$\mathbb{E}F^{{\bf{v}},\mathsf{q}}(Y^{1})-\mathbb{E}F^{{\bf{v}},\mathsf{q}}(Y^{0})\;=\;\int_{0}^{1}\mathrm{d}s\,\sum_{j=1}^{2p}\sum_{\mu=1}^{n}\Bigl{[}{\mathbb{E}F^{{\bf{v}},\mathsf{q}}\Bigl{(}{Y^{s,Z^{1}_{(j\mu)}}_{(j\mu)}}\Bigr{)}-\mathbb{E}F^{{\bf{v}},\mathsf{q}}\Bigl{(}{Y^{s,Z^{0}_{(j\mu)}}_{(j\mu)}}\Bigr{)}}\Bigr{]}\,,$
(F.12)

where we use the notation introduced in (F.3) with $\Gamma=Z^{0}_{(j\mu)}$ or
$Z^{1}_{(j\mu)}$. Under the property (${\bm{\mathrm{A}}}_{d-1}$), we claim the
following estimate:

$\sum_{j=1}^{2p}\sum_{\mu=1}^{n}\Bigl{[}{\mathbb{E}F^{{\bf{v}},\mathsf{q}}\Bigl{(}{Y^{s,Z^{1}_{(j\mu)}}_{(j\mu)}}\Bigr{)}-\mathbb{E}F^{{\bf{v}},\mathsf{q}}\Bigl{(}{Y^{s,Z^{0}_{(j\mu)}}_{(j\mu)}}\Bigr{)}}\Bigr{]}\lesssim(n^{24\delta}\Psi)^{\mathsf{q}}+|\mathbb{E}F^{{\bf{v}},\mathsf{q}}(Y^{s})|.$
(F.13)

By substituting this estimate into (F.12) and applying a standard Grönwall’s
argument, we conclude (F.10).

The remainder of the proof focuses on establishing (F.13). To simplify the
presentation, until the end of this proof, we slightly abuse the notation and
define

$G(Y)\equiv G(Y,z):=\begin{pmatrix}-D_{t}^{-1}&Y\\\
Y^{\top}&-z\end{pmatrix}^{-1}$ (F.14)

for given random matrix $Y$. By the definition (F.3), for any given
$\Gamma,\Gamma^{\prime}\in\mathbb{R}^{p\times n}$ and $M\in\mathbb{N}$, we
have the following resolvent expansion:

$\displaystyle
G\left(Y^{s,\Gamma^{\prime}}_{(j\mu)}\right)-G\left(Y^{s,\Gamma}_{(j\mu)}\right)=$
$\displaystyle~{}\sum_{k=1}^{M}G\left(Y^{s,\Gamma}_{(j\mu)}\right)\left[\begin{pmatrix}0&\Gamma-\Gamma^{\prime}\\\
(\Gamma-\Gamma^{\prime})^{\top}&0\end{pmatrix}G\left(Y^{s,\Gamma}_{(j\mu)}\right)\right]^{k}$
(F.15)
$\displaystyle~{}+G\left(Y^{s,\Gamma^{\prime}}_{(j\mu)}\right)\left[\begin{pmatrix}0&\Gamma-\Gamma^{\prime}\\\
(\Gamma-\Gamma^{\prime})^{\top}&0\end{pmatrix}G\left(Y^{s,\Gamma}_{(j\mu)}\right)\right]^{M+1}.$

By employing this expansion with $\Gamma=Z_{(j\mu)}^{s}$ and choosing a
sufficiently large $M$, we can utilize property $({\bf A}_{d-1})$ for
$G_{t}^{s}$ and the rough bound
$\|G(Y^{s,\Gamma^{\prime}}_{(j\mu)})\|\prec\eta^{-1}$ to verify the following
inequalities for any deterministic unit vectors
${\bf{u}},{\bf{v}}\in\mathbb{R}^{p+n}$:

$\left|{\bf{u}}^{\top}G\Bigl{(}{Y^{s,Z^{1}_{(j\mu)}}_{(j\mu)}}\Bigr{)}{\bf{v}}\right|\prec
1,\quad\left|{\bf{u}}^{\top}G\Bigl{(}{Y^{s,Z^{0}_{(j\mu)}}_{(j\mu)}}\Bigr{)}{\bf{v}}\right|\prec
1,\quad\left|{\bf{u}}^{\top}G\Bigl{(}{Y^{s,\mathbf{0}}_{(j\mu)}}\Bigr{)}{\bf{v}}\right|\prec
1.$ (F.16)

To simplify the notation, for fixed $s\in(0,1)$, we further define

$f_{(j\mu)}(\lambda_{0},\lambda_{1})\mathrel{\vbox{\hbox{.}\hbox{.}}}=F^{{\bf{v}},\mathsf{q}}\left(Y^{s,\lambda_{0}U^{0}{\bf
e}_{j\mu}+\lambda_{1}U^{1}{\bf e}_{j\mu}}_{(j\mu)}\right)\,,$ (F.17)

and denote its derivatives as

$f_{(j\mu)}^{(k,l)}=\partial_{\lambda_{0}}^{k}\partial_{\lambda_{1}}^{l}f_{(j\mu)}(\lambda_{0},\lambda_{1}),\quad
k,l\in\mathbb{N}.$

By applying the Taylor expansion and using the estimate (F.16), we derive that

$\displaystyle\mathbb{E}F^{{\bf{v}},\mathsf{q}}\Bigl{(}{Y^{s,Z^{0}_{(j\mu)}}_{(j\mu)}}\Bigr{)}-\mathbb{E}F^{{\bf{v}},\mathsf{q}}\Bigl{(}{Y^{s,\mathbf{0}}_{(j\mu)}}\Bigr{)}=\sum_{\ell=2}^{4\mathsf{q}}\frac{1}{\ell!}\mathbb{E}f_{(j\mu)}^{(\ell,0)}(0,0)\mathbb{E}(\mathsf{X}_{j\mu})^{\ell}+\operatorname{O}_{\prec}(\Psi^{\mathsf{q}}),$
(F.18)
$\displaystyle\mathbb{E}F^{{\bf{v}},\mathsf{q}}\Bigl{(}{Y^{s,Z^{1}_{(j\mu)}}_{(j\mu)}}\Bigr{)}-\mathbb{E}F^{{\bf{v}},\mathsf{q}}\Bigl{(}{Y^{s,\mathbf{0}}_{(j\mu)}}\Bigr{)}\;=\;\sum_{\ell=2}^{4\mathsf{q}}\frac{1}{\ell!}\mathbb{E}f_{(j\mu)}^{(0,\ell)}(0,0)\mathbb{E}(\mathsf{X}_{j\mu})^{\ell}+\operatorname{O}_{\prec}(\Psi^{\mathsf{q}}),$
(F.19)

where we have utilized the mean zero condition $\mathbb{E}\mathsf{X}_{j\mu}=0$
and the independence between $Y^{s,\bm{\mathrm{0}}}_{(j\mu)}$ and
$\mathsf{X}_{j\mu}$. Now, to establish (F.13), we only need to prove that

$\sum_{j,\mu}\mathbb{E}(\mathsf{X}_{j\mu})^{\ell}\left(\mathbb{E}f_{(j\mu)}^{(\ell,0)}(0,0)-\mathbb{E}f_{(j\mu)}^{(0,\ell)}(0,0)\right)\lesssim(n^{24\delta}\Psi)^{\mathsf{q}}+|\mathbb{E}F^{{\bf{v}},\mathsf{q}}(Y^{s})|.$
(F.20)

Similar to (F.18) and (F.19), we can utilize the Taylor expansion and the
estimate (F.16) to derive that

$\displaystyle\mathbb{E}f_{(j\mu)}^{(\mathsf{k},\widetilde{\mathsf{k}})}(0,0)=\mathbb{E}f_{(j\mu)}^{(\mathsf{k},\widetilde{\mathsf{k}})}\left(\chi^{s}_{j\mu}\mathsf{X}_{j\mu},(1-\chi^{s}_{j\mu})\mathsf{X}_{j\mu}\right)$
(F.21)
$\displaystyle-s\sum_{l=2}^{4\mathsf{q}-\mathsf{k}-\widetilde{\mathsf{k}}}\frac{1}{l!}\mathbb{E}f_{(j\mu)}^{(\mathsf{k}+l,\widetilde{\mathsf{k}})}(0,0)\mathbb{E}(\mathsf{X}_{j\mu})^{l}-(1-s)\sum_{\widetilde{l}=2}^{4\mathsf{q}-\mathsf{k}-\widetilde{\mathsf{k}}}\frac{1}{\widetilde{l}!}\mathbb{E}f_{(j\mu)}^{(\mathsf{k},\widetilde{\mathsf{k}}+\widetilde{l})}(0,0)\mathbb{E}(\mathsf{X}_{j\mu})^{\widetilde{l}}+\operatorname{O}_{\prec}(\Psi^{\mathsf{q}-\mathsf{k}-\widetilde{\mathsf{k}}})$

for any fixed $\mathsf{k}$ and $\widetilde{\mathsf{k}}$. In the derivation, we
also used the fact that $\chi^{s}_{j\mu}(1-\chi^{s}_{j\mu})=0$. By repeatedly
applying (F.21), we can further obtain that

$\displaystyle\mathbb{E}f_{(j\mu)}^{(\ell,0)}(0,0)=\sum_{\mathsf{k},\widetilde{\mathsf{k}}:\mathsf{k}+\widetilde{\mathsf{k}}=0}^{4\mathsf{q}-\ell}C_{(j\mu)}^{\mathsf{k},\widetilde{\mathsf{k}}}\sum_{j,\mu}\mathbb{E}f_{(j\mu)}^{(\ell+\mathsf{k},\widetilde{\mathsf{k}})}\left(\chi^{s}_{j\mu}\mathsf{X}_{j\mu},(1-\chi^{s}_{j\mu})\mathsf{X}_{j\mu}\right)+\operatorname{O}_{\prec}(\Psi^{\mathsf{q}-\ell}),$
$\displaystyle\mathbb{E}f_{(j\mu)}^{(0,\ell)}(0,0)=\sum_{\mathsf{k},\widetilde{\mathsf{k}}:\mathsf{k}+\widetilde{\mathsf{k}}=0}^{4\mathsf{q}-\ell}C_{(j\mu)}^{\mathsf{k},\widetilde{\mathsf{k}}}\sum_{j,\mu}\mathbb{E}f_{(j\mu)}^{(\mathsf{k},\ell+\widetilde{\mathsf{k}})}\left(\chi^{s}_{j\mu}\mathsf{X}_{j\mu},(1-\chi^{s}_{j\mu})\mathsf{X}_{j\mu}\right)+\operatorname{O}_{\prec}(\Psi^{\mathsf{q}-\ell}),$

where $C_{(j\mu)}^{\mathsf{k},\widetilde{\mathsf{k}}}$ represents
deterministic coefficients satisfying that $C_{(j\mu)}^{0,0}=1$,

$C_{(j\mu)}^{\mathsf{k},\widetilde{\mathsf{k}}}=\operatorname{O}(n^{-{(\mathsf{k}+\widetilde{\mathsf{k}})}/2}),\quad\text{and}\quad\mathsf{k}=1\
\ \text{or}\ \ \widetilde{\mathsf{k}}=1\implies
C_{(j\mu)}^{\mathsf{k},\widetilde{\mathsf{k}}}=0.$ (F.22)

Now, to conclude (F.20), it suffices to control
$\mathbb{E}f_{(j\mu)}^{(\ell+\mathsf{k},\widetilde{\mathsf{k}})}\left(\chi^{s}_{j\mu}\mathsf{X}_{j\mu},(1-\chi^{s}_{j\mu})\mathsf{X}_{j\mu}\right)-\mathbb{E}f_{(j\mu)}^{(\mathsf{k},\ell+\widetilde{\mathsf{k}})}\left(\chi^{s}_{j\mu}\mathsf{X}_{j\mu},(1-\chi^{s}_{j\mu})\mathsf{X}_{j\mu}\right)$.
For the terms in (F.21), we express them as

$\mathbb{E}f_{(j\mu)}^{(\mathsf{k},\widetilde{\mathsf{k}})}\left(\chi^{s}_{j\mu}\mathsf{X}_{j\mu},(1-\chi^{s}_{j\mu})\mathsf{X}_{j\mu}\right)=\mathbb{E}\Big{[}\Big{(}P^{U^{0}}_{j\mu}\Big{)}^{\mathsf{k}}\Big{(}P^{U^{1}}_{j\mu}\Big{)}^{\widetilde{\mathsf{k}}}F^{{\bf{v}},\mathsf{q}}\Big{]}\left(Y^{s}\right),$
(F.23)

where we adopt the following notation given any $p\times 2p$ matrix $U$:

$P^{U}_{j\mu}:=(U^{\top}\nabla_{Y})_{j\mu}=\sum_{i=1}^{p}{U}_{ij}\frac{\partial}{\partial
Y^{s}_{i\mu}}.$ (F.24)

Note that $P^{U}_{j\mu}$ acting on any resolvent entry yields

$P^{U}_{j\mu}G_{\mathfrak{a}\mathfrak{b}}=-G_{\mathfrak{a}{\bf{u}}_{j}}G_{\mu\mathfrak{b}}-G_{\mathfrak{a}\mu}G_{{\bf{u}}_{j}\mathfrak{b}},\quad\mathfrak{a},\mathfrak{b}\in\mathcal{I},$
(F.25)

where ${\bf{u}}_{j}=U\mathbf{e}_{j}$ represents the $j$-th column vector of
$U$. With the above notations, we observe that to establish (F.20), it is
sufficient to prove that for any $\ell,$ $\mathsf{k}$,
$\widetilde{\mathsf{k}}$ satisfying $\mathsf{k}\neq 1$,
$\widetilde{\mathsf{k}}\neq 1$, $\ell\geq 2$, and
$\ell+\mathsf{k}+\widetilde{\mathsf{k}}\leq 4\mathsf{q},$

$\displaystyle
n^{-(\mathsf{k}+\widetilde{\mathsf{k}}+\ell)/2}\sum_{j,\mu}\mathcal{C}_{(j\mu)}^{\ell}\mathbb{E}\Big{\\{}\Big{[}\Big{(}P^{U^{0}}_{j\mu}\Big{)}^{\ell+\mathsf{k}}\Big{(}P^{U^{1}}_{j\mu}\Big{)}^{\widetilde{\mathsf{k}}}F^{{\bf{v}},\,\mathsf{q}}\Big{]}\left(Y^{s}\right)-\Big{[}\Big{(}P^{U^{0}}_{j\mu}\Big{)}^{\mathsf{k}}\Big{(}P^{U^{1}}_{j\mu}\Big{)}^{\ell+\widetilde{\mathsf{k}}}F^{{\bf{v}},\,\mathsf{q}}\Big{]}\left(Y^{s}\right)\Big{\\}}$
(F.26)
$\displaystyle\lesssim(n^{24\delta}\Psi)^{\mathsf{q}}+|\mathbb{E}F^{{\bf{v}},\mathsf{q}}(Y^{s})|.$

Here, we introduce the notation
$\mathcal{C}_{(j\mu)}^{\ell}:=n^{\ell/2}\mathbb{E}\left(\mathsf{X}_{j\mu}\right)^{\ell}$,
which is of order $\operatorname{O}(1)$ by (2.14). The remaining part of the
proof is dedicated to verifying (F.26). To accomplish this, we divide the
proof into two cases.

Case 1: $\ell+\mathsf{k}+\widetilde{\mathsf{k}}=2.$ In this case, we must have
$\ell=2,\mathsf{k}=\widetilde{\mathsf{k}}=0$. Using the fact that
$U^{0}(U^{0})^{\top}=U^{1}(U^{1})^{\top}=I$, we can easily check that

$\sum_{j,\mu}\left(P^{U^{1}}_{j\mu}\right)^{2}F^{{\bf{v}},\,\mathsf{q}}(Y^{s})=\sum_{j,\mu}\left(P^{U^{0}}_{j\mu}\right)^{2}F^{{\bf{v}},\,\mathsf{q}}(Y^{s})=\sum_{i\mu}\frac{\partial^{2}}{\partial(Y^{s}_{i\mu})^{2}}F^{{\bf{v}},\,\mathsf{q}}(Y^{s}).$

This equation clearly implies (F.26) when
$\ell+\mathsf{k}+\widetilde{\mathsf{k}}=2.$

Case 2: $3\leq\ell+\mathsf{k}+\widetilde{\mathsf{k}}\leq 4\mathsf{q}.$ We
denote $\widetilde{U}=U^{1}-U^{0}.$ It is evident that in order to establish
(F.26), it suffices to prove the following estimate for every
$(k_{1},k_{2})\in\mathbb{N}^{2}$ satisfying $2\leq
k_{1}+k_{2}=\mathsf{k}+\widetilde{\mathsf{k}}+\ell-1\leq 4\mathsf{q}-1$:

$\displaystyle
n^{-(k_{1}+k_{2}+1)/2}\sum_{j,\mu}\mathcal{C}_{(j\mu)}^{\ell}\mathbb{E}\Big{|}\Big{[}\Big{(}P^{U^{0}}_{j\mu}\Big{)}^{k_{1}}\Big{(}P^{U^{1}}_{j\mu}\Big{)}^{k_{2}}P^{\widetilde{U}}_{j\mu}F^{{\bf{v}},\,\mathsf{q}}\Big{]}\Bigl{(}{Y^{s}}\Bigr{)}\Big{|}\lesssim(n^{24\delta}\Psi)^{\mathsf{q}}+|\mathbb{E}F^{{\bf{v}},\mathsf{q}}(Y^{s})|.$
(F.27)

For the proof of (F.27), we need to delve into the detailed structure of the
derivatives. To this end, in view of (F.24) and (F.25), we introduce the
following algebraic objects, which are employed in [45, 68].

######  Definition F.5 (Words).

Given $j\in\mathcal{I}_{1}$ and $\mu\in\mathcal{I}_{2}$, let $\mathcal{W}$ be
the set consisting of words of even length formed from the four letters
$\\{\mathbf{j}_{0},\mathbf{j}_{1},\widetilde{\mathbf{j}},\bm{\mu}\\}$. The
length of a word $w\in{\mathcal{W}}$ is denoted by $2{\bm{l}}(w)\in
2\mathbb{N}$. We use bold symbols to represent the letters in words. For
instance,

$w=\mathbf{t}_{1}\mathbf{s}_{2}\mathbf{t}_{2}\mathbf{s}_{3}\cdots\mathbf{t}_{l}\mathbf{s}_{l+1}$
(F.28)

denotes a word of length $2l$. Define
${\mathcal{W}}_{l}:=\\{w\in\mathcal{W}:{\bm{l}}(w)=l\\}$ as the set of words
of length $2l$ such that the following property holds for a word as in (F.28):

$\mathbf{t}_{i}\mathbf{s}_{i+1}\in\\{\mathbf{j}_{0}\bm{\mu},\bm{\mu}\mathbf{j}_{0},\mathbf{j}_{1}\bm{\mu},\bm{\mu}\mathbf{j}_{1},\widetilde{\mathbf{j}}\bm{\mu},\bm{\mu}\widetilde{\mathbf{j}}\\},\quad
i\in[\\![{l}]\\!].$

Next, we assign a value $[\cdot]$ to each letter as follows:
$[\mathbf{j}_{0}]={\bf{u}}_{j}^{0}$, $[\mathbf{j}_{1}]={\bf{u}}_{j}^{1}$,
$[\widetilde{\mathbf{j}}]=\widetilde{\bf{u}}_{j}$,
$[\bm{\mu}]:=\mathbf{e}_{\mu}$, where ${\bf{u}}_{j}^{0}:=U^{0}\mathbf{e}_{j}$,
${\bf{u}}_{j}^{1}:=U^{1}\mathbf{e}_{j}$,
$\widetilde{\bf{u}}_{j}:=\widetilde{U}\mathbf{e}_{j}/\|\widetilde{U}\mathbf{e}_{j}\|_{2}$.
It is important to distinguish between the abstract letter and its
corresponding value, which is treated as a (generalized) summation index. For
the word $w$ in (F.28) and vectors ${\bf{u}},{\bf{v}}\in\mathbb{C}^{p+n}$, we
assign a random variable $A_{\mathbf{u},\mathbf{v},j,\mu}(w)$ as follows: if
${\bm{l}}(w)=0$, we define

$A_{{\bf{u}},{\bf{v}},j,\mu}(w):=\big{(}G(Y^{s})-\widetilde{\Pi}_{t}\big{)}_{{\bf{u}}{\bf{v}}};$

if ${\bm{l}}(w)\geq 1$, we define

$A_{{\bf{u}},{\bf{v}},j,\mu}(w):=G_{{\bf{u}}[\mathbf{t}_{1}]}(Y^{s})G_{[\mathbf{s}_{2}][\mathbf{t}_{2}]}(Y^{s})\cdots
G_{[\mathbf{s}_{l}][\mathbf{t}_{l}]}(Y^{s})G_{[\mathbf{s}_{l+1}]{\bf{v}}}(Y^{s}).$
(F.29)

Finally, given any $w\in\mathcal{W}_{l}$ of the form (F.28), we denote
$\mathfrak{n}_{0}(w):=\\#\\{i\in[\\![{l}]\\!]:\mathbf{t}_{i}=\mathbf{j}_{0}\
\text{or}\ \mathbf{s}_{i+1}=\mathbf{j}_{0}\\}$,
$\mathfrak{n}_{1}(w):=\\#\\{i\in[\\![{l}]\\!]:\mathbf{t}_{i}=\mathbf{j}_{1}\
\text{or}\ \mathbf{s}_{i+1}=\mathbf{j}_{1}\\}$, and
$\widetilde{\mathfrak{n}}(w):=\\#\\{i\in[\\![{l}]\\!]:\mathbf{t}_{i}=\widetilde{\mathbf{j}}\
\text{or}\ \mathbf{s}_{i+1}=\widetilde{\mathbf{j}}\\}$.

From (F.15) and (F.25), we observe that the above notations are defined such
that for every $l\in\mathbb{N}$, the following equality holds:

$\Big{(}P^{U^{0}}_{j\mu}\Big{)}^{l}\big{(}G(Y^{s})-\widetilde{\Pi}_{t}\big{)}_{{\bf{u}}{\bf{v}}}=(-1)^{l}l!\sum_{w\in\mathcal{W}_{l}:\mathfrak{n}_{0}(w)=l}A_{{\bf{u}},{\bf{v}},j,\mu}(w).$

A similar equation holds if we replace $0$ with $1$. By denoting
$\alpha_{j}:=\|\widetilde{U}\mathbf{e}_{j}\|_{2}$, we have that

$P^{\widetilde{U}}_{j\mu}\big{(}G(Y^{s})-\widetilde{\Pi}_{t}\big{)}_{{\bf{u}}{\bf{v}}}=-\alpha_{j}\sum_{w\in\\{\widetilde{\mathbf{j}}\bm{\mu},\bm{\mu}\widetilde{\mathbf{j}}\\}}A_{{\bf{u}},{\bf{v}},j,\mu}(w).$

Using the above two identities, we can derive that

$\displaystyle\Big{(}P^{U^{0}}_{j\mu}\Big{)}^{k_{1}}\Big{(}P^{U^{1}}_{j\mu}\Big{)}^{k_{2}}P^{\widetilde{U}}_{j\mu}F^{{\bf{v}},\,\mathsf{q}}(Y^{s})$
$\displaystyle=\sum_{\begin{subarray}{c}l_{1},\ldots,l_{\mathsf{q}}\geq 0:\\\
l_{1}+\cdots+l_{\mathsf{q}}=k_{1}+k_{2}+1\end{subarray}}\sum_{w_{1}\in\mathcal{W}_{l_{1}},\ldots,w_{\mathsf{q}}\in\mathcal{W}_{l_{\mathsf{q}}}}^{*}c(w_{1},\ldots,w_{\mathsf{q}})\alpha_{j}$
$\displaystyle\times\prod_{\mathsf{t}=1}^{\mathsf{q}/2}\left[A_{{\bf{v}},{\bf{v}},j,\mu}(w_{\mathsf{t}})\overline{A_{{\bf{v}},{\bf{v}},j,\mu}(w_{\mathsf{t}+\mathsf{q}/2})}\right].$
(F.30)

Here, $c(w_{1},\ldots,w_{\mathsf{q}})$ represents certain deterministic
coefficients of order $\operatorname{O}(1)$, and
$\sum_{w_{1}\in\mathcal{W}_{l_{1}},\ldots,w_{\mathsf{q}}\in\mathcal{W}_{l_{\mathsf{q}}}}^{*}$
denotes the summation over a sequence of words $w_{1},\ldots,w_{\mathsf{q}}$
satisfying the following conditions:

$\sum_{\mathsf{t}=1}^{\mathsf{q}}\mathfrak{n}_{0}(w_{\mathsf{t}})=k_{1},\quad\sum_{\mathsf{t}=1}^{\mathsf{q}}\mathfrak{n}_{1}(w_{\mathsf{t}})=k_{2},\quad\sum_{\mathsf{t}=1}^{\mathsf{q}}\widetilde{\mathfrak{n}}(w_{\mathsf{t}})=1.$
(F.31)

Without loss of generality, suppose there are $\mathsf{p}$ non-zero length
words in (F.30) for some $1\leq\mathsf{p}\leq\mathsf{q}\wedge(k_{1}+k_{2}+1)$,
and these words are $w_{1},\ldots,w_{\mathsf{p}}$. Then, we can identify and
separate the words of length 0 from the product in (F.30), denoted as $w_{0}$.
Now, to show (F.27), it suffices to prove that for any
$1\leq\mathsf{p}\leq\mathsf{q}\wedge(k_{1}+k_{2}+1)$ and a sequence of words
satisfying (F.31), the following estimate holds:

$\displaystyle
n^{-(k_{1}+k_{2}+1)/2}\sum_{j,\mu}\alpha_{j}\mathbb{E}\left|A_{{\bf{v}},{\bf{v}},j,\mu}(w_{0})^{\mathsf{q}-\mathsf{p}}\prod_{r=1}^{\mathsf{p}}A_{{\bf{v}},{\bf{v}},j,\mu}(w_{r})\right|\lesssim(n^{24\delta}\Psi)^{\mathsf{q}}+|\mathbb{E}F^{{\bf{v}},\mathsf{q}}(Y^{s})|.$
(F.32)

The proof of (F.32) follows a similar approach to that of (6.32) in [68]. As
we assume $\mathbb{E}\mathsf{X}^{3}_{i\mu}=0$ according to (ii) of 2.3, we
only need to consider the case when $\ell\geqslant 4$. (It is in this
particular case that the vanishing third-moment condition is relevant and
utilized.) For $j\in\mathcal{I}_{1}$ and $\mu\in\mathcal{I}_{2}$, we define

$\mathcal{R}_{j}:=\big{|}G_{{\bf{v}}{\bf{u}}^{0}_{j}}\big{|}+\big{|}G^{s}_{{\bf{v}}{\bf{u}}^{1}_{j}}\big{|}+\big{|}G_{{\bf{v}}_{\alpha}\widetilde{\bf{u}}_{j}}\big{|},\quad\mathcal{R}_{\mu}:=\big{|}G_{{\bf{v}}\mu}\big{|}.$

Similar to Lemma 6.16 in [68], we can establish the following rough bound
using property ($\mathbf{A}_{d-1}$):

$|A_{{\bf{v}},{\bf{v}},i,\mu}(w)|\prec n^{2\delta({\bm{l}}(w)+1)},$

Additionally, we can get the following bounds: for ${\bm{l}}(w)\geq 2$,

$|A_{{\bf{v}},{\bf{v}},i,\mu}(w)|\prec(\mathcal{R}_{i}^{2}+\mathcal{R}_{\mu}^{2})n^{2\delta({\bm{l}}(w)-1)},$

and for ${\bm{l}}(w)=1$,

$|A_{{\bf{v}},{\bf{v}},i,\mu}(w)|\prec\mathcal{R}_{i}\mathcal{R}_{\mu}.$

With these estimates, we can follow the arguments presented between (6.39) and
(6.42) in [68] to complete the proof of (F.32). (It is worth noting that for
this proof, the condition (4.35), which implies $\alpha_{j}\leq
Cn^{-\mathsf{c}^{\prime}}$, is not necessary.) Since the argument is almost
the same, we omit the details here. This concludes the proof of Lemma F.4. ∎

With the above results, we proceed to complete the proof of Lemma 4.12 by
using the following lemma.

######  Lemma F.6.

Under the assumptions of Lemma 4.12, the following estimate holds for every
$s\in[0,1]$:

$\displaystyle\mathbb{E}\theta\left(\mathcal{U}_{i_{1}}(\mathsf{Q}_{t}^{s}),\ldots,\mathcal{U}_{i_{L}}(\mathsf{Q}_{t}^{s})\right)=\mathbb{E}\theta\left(\overline{\mathcal{U}}_{i_{1}}(\mathsf{Q}_{t}^{s}),\ldots,\overline{\mathcal{U}}_{i_{L}}(\mathsf{Q}_{t}^{s})\right)+\operatorname{O}(n^{-\nu}),$

for a constant $\nu>0$, where
$\overline{\mathcal{U}}_{k}(\mathsf{Q}_{t}^{s})$,