Before we proceed with the proof, let us comment on some apparent difficulties. For this discussion, we assume normalized $F \in \mathcal{P}(\mathbb{R}^d)$ without loss of generality, as we are only concerned with the shape of profiles here. Note that the conservative operator $\mathcal{S}$ in Definition 4.4 associated with $\mathcal{T}$ is actually the composition of the conservative operator $\mathcal{B}$, and the conservative multiplicative operator $\mathcal{M}$ in Definition 2.3, namely, $\mathcal{S} = \mathcal{M} \circ \mathcal{B}$. The former is invariant by translation, whereas the latter is not. However, we have the following inequality for the operator $\mathcal{M}$:

$$ \int_{\mathbb{R}^d} |x|^2 \mathcal{M}[F](x) dx \leq \int_{\mathbb{R}^d} |x|^2 F(dx). \qquad (4.9) $$

This inequality is intuitively clear: by applying the selection operator $\mathcal{M}$, the quadratic moment (that is, the Wasserstein distance to the Dirac distribution at the origin) is non-increasing. In fact, the proof is a straightforward consequence of the following identity:

$$ \int_{\mathbb{R}^d} |x|^2 \mathcal{M}[F](x) dx - \int_{\mathbb{R}^d} |x|^2 F(dx) = \frac{1}{2} \int_{\mathbb{R}^{2d}} (|x|^2 - |y|^2) (e^{-m(x)} - e^{-m(y)}) \frac{F(dx)F(dy)}{\int_{\mathbb{R}^d} e^{-m(z)} F(dz)}. $$

Note that the right-hand-side of this equality is indeed non-positive, provided that $m$ is radially non-decreasing as it is the case for our quadratic choice (1.5). A similar estimate can also be tested on the operator $\mathcal{B}$. Specifically, by a simple change of variable the following identity holds true:

$$ \int_{\mathbb{R}^d} |x|^2 \mathcal{B}[F](x) dx = \mathrm{Var} G + \frac{1}{2} \int_{\mathbb{R}^d} |x|^2 F(dx) + \frac{1}{2} \left| \int_{\mathbb{R}^d} x F(dx) \right|^2. \qquad (4.10) $$

Then, it might seem natural to combine these two relationships in order to control the quadratic moments of the composition $\mathcal{S} = \mathcal{M} \circ \mathcal{B}$. Suppose, for instance, that the center of mass in the last term of (4.10) can be bounded a priori uniformly on $F$. Then, the combination of (4.9) and (4.10) yields the contraction property in Lemma 4.5 due to the reduction by half in the size of the quadratic moment of $F$. Of course, without such a uniform control on the center of mass, the last two terms in (4.10) would contribute together (by Jensen's inequality) with a merely non-expansive factor $\int_{\mathbb{R}^d} |x|^2 F(dx)$ that is not enough for our purpose (we will iterate Lemma 4.5 on the $F_n$ later). This a priori estimate of the center of mass should then be derived mainly from the selection component $\mathcal{M}$, as the reproduction component $\mathcal{B}$ is invariant by translation. However, it cannot be derived solely from $\mathcal{M}$, because the latter can have a dramatic effect on nearly symmetric distributions. This can be seen on the same configuration as in Example 2.4 above: for a sum of nearly symmetrical Dirac masses, the center of mass is close to zero before selection, and close to one of the endpoints after selection (when $h$ is large). Nevertheless, such a configuration is destroyed by $\mathcal{B}$, so that the combination of $\mathcal{M}$ and $\mathcal{B}$ is crucial to guarantee the uniform boundedness of the quadratic moments as shown in the argument below.

**Proof of Lemma 4.5.** Set any $F \in \mathcal{M}_+(\mathbb{R}^d)$. First, note that

$$ \int_{\mathbb{R}^d} |x|^2 \mathcal{S}[F](x) dx = \frac{\int_{\mathbb{R}^{3d}} |x|^2 e^{-m(x)} G\left(x - \frac{x_1 + x_2}{2}\right) \frac{F(dx_1)}{\|F\|_{\mathcal{M}(\mathbb{R}^d)}} \frac{F(dx_2)}{\|F\|_{\mathcal{M}(\mathbb{R}^d)}} dx}{\int_{\mathbb{R}^{3d}} e^{-m(x)} G\left(x - \frac{x_1 + x_2}{2}\right) \frac{F(dx_1)}{\|F\|_{\mathcal{M}(\mathbb{R}^d)}} \frac{F(dx_2)}{\|F\|_{\mathcal{M}(\mathbb{R}^d)}} dx}. $$

Since $G$ is a Gaussian function given by (1.4) and $m$ is a quadratic function given by (1.5), then we can compute explicitly the integrals with respect to $x$ in the numerator and denominator. Specifically, for fixed $x_1, x_2 \in \mathbb{R}^d$ we define $\bar{x} := \frac{x_1+x_2}{2}$ and we obtain

$$ \int_{\mathbb{R}^d} |x|^2 e^{-m(x)} G(x - \bar{x}) dx = \frac{1}{(1 + \alpha)^{d/2}} \exp \left( -\frac{1}{2} \frac{\alpha}{1 + \alpha} |\bar{x}|^2 \right) \left( \frac{1}{1 + \alpha} + \frac{|\bar{x}|^2}{(1 + \alpha)^2} \right), \qquad (4.11) $$

$$ \int_{\mathbb{R}^d} e^{-m(x)} G(x - \bar{x}) dx = \frac{1}{(1 + \alpha)^{d/2}} \exp \left( -\frac{1}{2} \frac{\alpha}{1 + \alpha} |\bar{x}|^2 \right). \qquad (4.12) $$

Therefore, we can write alternatively

$$ \int_{\mathbb{R}^d} |x|^2 \mathcal{S}[F](x) dx = \frac{1}{1+\alpha} + \frac{1}{(1+\alpha)^2} \frac{\int_{\mathbb{R}^{2d}} |\bar{x}|^2 \exp\left(-\frac{1}{2}\frac{\alpha}{1+\alpha}|\bar{x}|^2\right) \frac{F(dx_1)}{\|F\|_{\mathcal{M}(\mathbb{R}^d)}} \frac{F(dx_2)}{\|F\|_{\mathcal{M}(\mathbb{R}^d)}}}{\int_{\mathbb{R}^{2d}} \exp\left(-\frac{1}{2}\frac{\alpha}{1+\alpha}|\bar{x}|^2\right) \frac{F(dx_1)}{\|F\|_{\mathcal{M}(\mathbb{R}^d)}} \frac{F(dx_2)}{\|F\|_{\mathcal{M}(\mathbb{R}^d)}}}. \qquad (4.13) $$

Our ultimate goal is to bound the second term in (4.13). To this end, we provide a lower bound for the denominator and an upper bound for the numerator in two separate steps.
• STEP 1: Consider an arbitrary $R_1 \in \mathbb{R}_+^*$ to be determined later. Then, we obtain

$$ 
\begin{aligned}
& \int_{\mathbb{R}^{2d}} \exp\left(-\frac{1}{2}\frac{\alpha}{1+\alpha}|\bar{x}|^2\right) \frac{F(dx_1)}{\|F\|_{\mathcal{M}(\mathbb{R}^d)}} \frac{F(dx_2)}{\|F\|_{\mathcal{M}(\mathbb{R}^d)}} \\
& \geq \int_{|\bar{x}| \leq R_1} \exp\left(-\frac{1}{2}\frac{\alpha}{1+\alpha}|\bar{x}|^2\right) \frac{F(x_1)}{\|F\|_{L^1(\mathbb{R}^d)}} \frac{F(x_2)}{\|F\|_{L^1(\mathbb{R}^d)}} dx_1 dx_2 \\
& \geq \exp\left(-\frac{1}{2}\frac{\alpha}{1+\alpha}R_1^2\right) \left(1 - \int_{|\bar{x}| > R_1} \frac{F(dx_1)}{\|F\|_{\mathcal{M}(\mathbb{R}^d)}} \frac{F(dx_2)}{\|F\|_{\mathcal{M}(\mathbb{R}^d)}}\right) \\
& \geq \exp\left(-\frac{1}{2}\frac{\alpha}{1+\alpha}R_1^2\right) \left(1 - \frac{1}{2R_1^2} \int_{\mathbb{R}^d} |x|^2 \frac{F(dx)}{\|F\|_{\mathcal{M}(\mathbb{R}^d)}}\right),
\end{aligned}
\quad (4.14) $$

where in last line we have used Young's and Chebyshev's inequalities.

• STEP 2: Consider an arbitrary $R_2 \in \mathbb{R}_+^*$ to be determined later. Then, we obtain

$$ \int_{\mathbb{R}^{2d}} |\bar{x}|^2 \exp\left(-\frac{1}{2}\frac{\alpha}{1+\alpha}|\bar{x}|^2\right) \frac{F(dx_1)}{\|F\|_{\mathcal{M}(\mathbb{R}^d)}} \frac{F(dx_2)}{\|F\|_{\mathcal{M}(\mathbb{R}^d)}} = I_1 + I_2, $$

where each term reads

$$ I_1 := \int_{|\bar{x}| \le R_2} |\bar{x}|^2 \exp\left(-\frac{1}{2}\frac{\alpha}{1+\alpha}|\bar{x}|^2\right) \frac{F(dx_1)}{\|F\|_{\mathcal{M}(\mathbb{R}^d)}} \frac{F(dx_2)}{\|F\|_{\mathcal{M}(\mathbb{R}^d)}}, $$

$$ I_2 := \int_{|\bar{x}| > R_2} |\bar{x}|^2 \exp\left(-\frac{1}{2}\frac{\alpha}{1+\alpha}|\bar{x}|^2\right) \frac{F(dx_1)}{\|F\|_{\mathcal{M}(\mathbb{R}^d)}} \frac{F(dx_2)}{\|F\|_{\mathcal{M}(\mathbb{R}^d)}}. $$

On the one hand, the first term can be readily estimated by

$$ I_1 \le R_2^2 \int_{\mathbb{R}^{2d}} \exp\left(-\frac{1}{2}\frac{\alpha}{1+\alpha}|\bar{x}|^2\right) \frac{F(dx_1)}{\|F\|_{\mathcal{M}(\mathbb{R}^d)}} \frac{F(dx_2)}{\|F\|_{\mathcal{M}(\mathbb{R}^d)}}. \quad (4.15) $$

On the other hand, the second term can be estimated by

$$ 
\begin{aligned}
I_2 &\le \frac{2}{e\gamma}\frac{1+\alpha}{\alpha} \int_{|\bar{x}|>R_2} \exp\left(-\frac{1-\gamma}{2}\frac{\alpha}{1+\alpha}|\bar{x}|^2\right) \frac{F(dx_1)}{\|F\|_{\mathcal{M}(\mathbb{R}^d)}} \frac{F(dx_2)}{\|F\|_{\mathcal{M}(\mathbb{R}^d)}} \\
&\le \frac{2}{e\gamma}\frac{1+\alpha}{\alpha} \exp\left(-\frac{1-\gamma}{2}\frac{\alpha}{1+\alpha}R_2^2\right),
\end{aligned}
\quad (4.16) $$

for any $\gamma \in (0,1)$, where in the first line we have used that $e^{\theta s} \ge \theta e s$ with $s = |\bar{x}|^2$ and $\theta = \frac{\gamma}{2}\frac{\alpha}{1+\alpha}$, and in the second line we have exploited that $F/\|F\|_{\mathcal{M}(\mathbb{R}^d)}$ are probability measures, thus normalized.

• STEP 3: Putting (4.14), (4.15) and (4.16) into (4.13) yields

$$ \int_{\mathbb{R}^d} |x|^2 \mathcal{S}[F](x) dx \le \frac{1}{1+\alpha} + \frac{R_2^2}{(1+\alpha)^2} + \frac{2}{e\gamma} \frac{1}{\alpha(1+\alpha)} \frac{\exp\left(\frac{1}{2}\frac{\alpha}{1+\alpha}R_1^2\right) \exp\left(-\frac{1-\gamma}{2}\frac{\alpha}{1+\alpha}R_2^2\right)}{1 - \frac{1}{2R_1^2} \int_{\mathbb{R}^d} |x|^2 \frac{F(dx)}{\|F\|_{\mathcal{M}(\mathbb{R}^d)}}}, $$

for any $R_1, R_2 \in \mathbb{R}_+^*$ and $\gamma \in (0,1)$. We set the values

$$ R_1^2 = \frac{1}{2(1-\delta)} \int_{\mathbb{R}^d} |x|^2 \frac{F(dx)}{\|F\|_{\mathcal{M}(\mathbb{R}^d)}}, \quad R_2^2 = \frac{R_1^2}{1-\gamma}, $$

with $\gamma, \delta \in (0,1)$. Then, we obtain

$$ \int_{\mathbb{R}^d} |x|^2 \mathcal{S}[F](x) dx \le \frac{1}{1+\alpha} + \frac{2}{e\gamma\delta}\frac{1}{\alpha(1+\alpha)} + \frac{1}{2(1-\delta)(1-\gamma)(1+\alpha)^2} \int_{\mathbb{R}^d} |x|^2 \frac{F(dx)}{\|F\|_{\mathcal{M}(\mathbb{R}^d)}}, $$

thus ending the proof by the arbitrariness of $\gamma, \delta \in (0,1)$. $\square$

By iteration, we then obtain the following contractivity of the variance of generic solutions of (1.1).

**Corollary 4.6** (Propagation of quadratic moments). Assume that $\alpha \in \mathbb{R}_+^*$, set any parameter $\eta \in (\frac{1}{2(1+\alpha)^2}, 1)$ and consider the solution $\{F_n\}_{n \in \mathbb{N}}$ of (1.1) issued at a generic initial datum $F_0 \in \mathcal{M}_+(\mathbb{R}^d)$ with $\int_{\mathbb{R}^d} |x|^2 F_0(dx) < \infty$. Then, there exists $M = M(\alpha, \eta)$ such that

$$ \int_{\mathbb{R}^d} |x|^2 \frac{F_n(x)}{\|F_n\|_{L^1(\mathbb{R}^d)}} dx \le \frac{M}{1-\eta} + \eta^n \int_{\mathbb{R}^d} |x|^2 \frac{F_0(dx)}{\|F_0\|_{\mathcal{M}(\mathbb{R}^d)}}, $$

for any $n \in \mathbb{N}$.
• STEP 1: Consider an arbitrary $R_1 \in \mathbb{R}_+^*$ to be determined later. Then, we obtain

$$ 
\begin{aligned}
& \int_{\mathbb{R}^{2d}} \exp\left(-\frac{1}{2}\frac{\alpha}{1+\alpha}|\bar{x}|^2\right) \frac{F(dx_1)}{\|F\|_{\mathcal{M}(\mathbb{R}^d)}} \frac{F(dx_2)}{\|F\|_{\mathcal{M}(\mathbb{R}^d)}} \\
& \geq \int_{|\bar{x}| \leq R_1} \exp\left(-\frac{1}{2}\frac{\alpha}{1+\alpha}|\bar{x}|^2\right) \frac{F(x_1)}{\|F\|_{L^1(\mathbb{R}^d)}} \frac{F(x_2)}{\|F\|_{L^1(\mathbb{R}^d)}} dx_1 dx_2 \\
& \geq \exp\left(-\frac{1}{2}\frac{\alpha}{1+\alpha}R_1^2\right) \left(1 - \int_{|\bar{x}| > R_1} \frac{F(dx_1)}{\|F\|_{\mathcal{M}(\mathbb{R}^d)}} \frac{F(dx_2)}{\|F\|_{\mathcal{M}(\mathbb{R}^d)}}\right) \\
& \geq \exp\left(-\frac{1}{2}\frac{\alpha}{1+\alpha}R_1^2\right) \left(1 - \frac{1}{2R_1^2} \int_{\mathbb{R}^d} |x|^2 \frac{F(dx)}{\|F\|_{\mathcal{M}(\mathbb{R}^d)}}\right),
\end{aligned}
\quad (4.14) $$

where in last line we have used Young's and Chebyshev's inequalities.

• STEP 2: Consider an arbitrary $R_2 \in \mathbb{R}_+^*$ to be determined later. Then, we obtain

$$ \int_{\mathbb{R}^{2d}} |\bar{x}|^2 \exp\left(-\frac{1}{2}\frac{\alpha}{1+\alpha}|\bar{x}|^2\right) \frac{F(dx_1)}{\|F\|_{\mathcal{M}(\mathbb{R}^d)}} \frac{F(dx_2)}{\|F\|_{\mathcal{M}(\mathbb{R}^d)}} = I_1 + I_2, $$

where each term reads

$$ I_1 := \int_{|\bar{x}| \le R_2} |\bar{x}|^2 \exp\left(-\frac{1}{2}\frac{\alpha}{1+\alpha}|\bar{x}|^2\right) \frac{F(dx_1)}{\|F\|_{\mathcal{M}(\mathbb{R}^d)}} \frac{F(dx_2)}{\|F\|_{\mathcal{M}(\mathbb{R}^d)}}, $$

$$ I_2 := \int_{|\bar{x}| > R_2} |\bar{x}|^2 \exp\left(-\frac{1}{2}\frac{\alpha}{1+\alpha}|\bar{x}|^2\right) \frac{F(dx_1)}{\|F\|_{\mathcal{M}(\mathbb{R}^d)}} \frac{F(dx_2)}{\|F\|_{\mathcal{M}(\mathbb{R}^d)}}. $$

On the one hand, the first term can be readily estimated by

$$ I_1 \le R_2^2 \int_{\mathbb{R}^{2d}} \exp\left(-\frac{1}{2}\frac{\alpha}{1+\alpha}|\bar{x}|^2\right) \frac{F(dx_1)}{\|F\|_{\mathcal{M}(\mathbb{R}^d)}} \frac{F(dx_2)}{\|F\|_{\mathcal{M}(\mathbb{R}^d)}}. \quad (4.15) $$

On the other hand, the second term can be estimated by

$$ 
\begin{aligned}
I_2 &\le \frac{2}{e\gamma}\frac{1+\alpha}{\alpha} \int_{|\bar{x}|>R_2} \exp\left(-\frac{1-\gamma}{2}\frac{\alpha}{1+\alpha}|\bar{x}|^2\right) \frac{F(dx_1)}{\|F\|_{\mathcal{M}(\mathbb{R}^d)}} \frac{F(dx_2)}{\|F\|_{\mathcal{M}(\mathbb{R}^d)}} \\
&\le \frac{2}{e\gamma}\frac{1+\alpha}{\alpha} \exp\left(-\frac{1-\gamma}{2}\frac{\alpha}{1+\alpha}R_2^2\right),
\end{aligned}
\quad (4.16) $$

for any $\gamma \in (0,1)$, where in the first line we have used that $e^{\theta s} \ge \theta e s$ with $s = |\bar{x}|^2$ and $\theta = \frac{\gamma}{2}\frac{\alpha}{1+\alpha}$, and in the second line we have exploited that $F/\|F\|_{\mathcal{M}(\mathbb{R}^d)}$ are probability measures, thus normalized.

• STEP 3: Putting (4.14), (4.15) and (4.16) into (4.13) yields

$$ \int_{\mathbb{R}^d} |x|^2 \mathcal{S}[F](x) dx \le \frac{1}{1+\alpha} + \frac{R_2^2}{(1+\alpha)^2} + \frac{2}{e\gamma} \frac{1}{\alpha(1+\alpha)} \frac{\exp\left(\frac{1}{2}\frac{\alpha}{1+\alpha}R_1^2\right) \exp\left(-\frac{1-\gamma}{2}\frac{\alpha}{1+\alpha}R_2^2\right)}{1 - \frac{1}{2R_1^2} \int_{\mathbb{R}^d} |x|^2 \frac{F(dx)}{\|F\|_{\mathcal{M}(\mathbb{R}^d)}}}, $$

for any $R_1, R_2 \in \mathbb{R}_+^*$ and $\gamma \in (0,1)$. We set the values

$$ R_1^2 = \frac{1}{2(1-\delta)} \int_{\mathbb{R}^d} |x|^2 \frac{F(dx)}{\|F\|_{\mathcal{M}(\mathbb{R}^d)}}, \quad R_2^2 = \frac{R_1^2}{1-\gamma}, $$

with $\gamma, \delta \in (0,1)$. Then, we obtain

$$ \int_{\mathbb{R}^d} |x|^2 \mathcal{S}[F](x) dx \le \frac{1}{1+\alpha} + \frac{2}{e\gamma\delta}\frac{1}{\alpha(1+\alpha)} + \frac{1}{2(1-\delta)(1-\gamma)(1+\alpha)^2} \int_{\mathbb{R}^d} |x|^2 \frac{F(dx)}{\|F\|_{\mathcal{M}(\mathbb{R}^d)}}, $$

thus ending the proof by the arbitrariness of $\gamma, \delta \in (0,1)$. $\square$

By iteration, we then obtain the following contractivity of the variance of generic solutions of (1.1).

**Corollary 4.6** (Propagation of quadratic moments). Assume that $\alpha \in \mathbb{R}_+^*$, set any parameter $\eta \in (\frac{1}{2(1+\alpha)^2}, 1)$ and consider the solution $\{F_n\}_{n \in \mathbb{N}}$ of (1.1) issued at a generic initial datum $F_0 \in \mathcal{M}_+(\mathbb{R}^d)$ with $\int_{\mathbb{R}^d} |x|^2 F_0(dx) < \infty$. Then, there exists $M = M(\alpha, \eta)$ such that

$$ \int_{\mathbb{R}^d} |x|^2 \frac{F_n(x)}{\|F_n\|_{L^1(\mathbb{R}^d)}} dx \le \frac{M}{1-\eta} + \eta^n \int_{\mathbb{R}^d} |x|^2 \frac{F_0(dx)}{\|F_0\|_{\mathcal{M}(\mathbb{R}^d)}}, $$

for any $n \in \mathbb{N}$.
Let us note that, on the one hand the factor $\eta^n$ makes the dependence on the variance of the initial datum $F_0$ negligible as the amount of generations increases. On the other hand, the constant $M$ does not depend on the chosen initial profile $F_0$. This is compatible with our ergodicity property in Theorem 1.1.

**Lemma 4.7** (Control of exponential moments). Assume that $\alpha \in \mathbb{R}_+^*$, set $\theta \in \mathbb{R}_+^*$ so that $\theta < \frac{\alpha}{2}$ and any parameter $\chi > \frac{\alpha}{1+\alpha} \frac{\theta}{2(1+\alpha-2\theta)(\alpha-2\theta)}$. Then, there exists a constant $C = C(\alpha, \theta, \chi)$ such that

$$ \int_{\mathbb{R}^d} e^{\theta|x|^2} S[F](x) dx \le C \left\{ 1 + \exp \left( \chi \int_{\mathbb{R}^d} |x|^2 \frac{F(dx)}{\|F\|_{\mathcal{M}(\mathbb{R}^d)}} \right) \right\}, $$

for any measure $F \in \mathcal{M}_+(\mathbb{R}^d)$.

*Proof.* Since $G$ is a Gaussian function given by (1.4) and $m$ is a quadratic function given by (1.5), then we can compute explicitly the integrals with respect to $x$ in the numerator and denominator like formulas (4.11)-(4.12) in the proof of Lemma 4.5 and we obtain

$$ \int_{\mathbb{R}^d} e^{\theta|x|^2} S[F](x) dx = \left( \frac{1+\alpha}{1+\alpha-2\theta} \right)^{d/2} \frac{\int_{\mathbb{R}^{2d}} \exp\left(-\frac{1}{2}\frac{\alpha-2\theta}{1+\alpha-2\theta}|\bar{x}|^2\right) \frac{F(dx_1)}{\|F\|_{\mathcal{M}(\mathbb{R}^d)}} \frac{F(dx_2)}{\|F\|_{\mathcal{M}(\mathbb{R}^d)}}}{\int_{\mathbb{R}^{2d}} \exp\left(-\frac{1}{2}\frac{\alpha}{1+\alpha}|\bar{x}|^2\right) \frac{F(dx_1)}{\|F\|_{\mathcal{M}(\mathbb{R}^d)}} \frac{F(dx_2)}{\|F\|_{\mathcal{M}(\mathbb{R}^d)}}}, \quad (4.17) $$

where $\bar{x} = \frac{x_1+x_2}{2}$. Again, we provide bounds for the numerator and denominator in (4.17). Regarding the denominator, we recover estimate (4.14) for any $R_1 \in \mathbb{R}_+^*$. Thus, we focus on the bound of the numerator

$$ \int_{\mathbb{R}^{2d}} \exp\left(-\frac{1}{2}\frac{\alpha-2\theta}{1+\alpha-2\theta}|\bar{x}|^2\right) \frac{F(dx_1)}{\|F\|_{\mathcal{M}(\mathbb{R}^d)}} \frac{F(dx_2)}{\|F\|_{\mathcal{M}(\mathbb{R}^d)}} = I_1 + I_2, $$

where each term reads

$$ I_1 := \int_{|\bar{x}| \le R_2} \exp\left(-\frac{1}{2}\frac{\alpha-2\theta}{1+\alpha-2\theta}|\bar{x}|^2\right) \frac{F(dx_1)}{\|F\|_{\mathcal{M}(\mathbb{R}^d)}} \frac{F(dx_2)}{\|F\|_{\mathcal{M}(\mathbb{R}^d)}}, $$

$$ I_2 := \int_{|\bar{x}| > R_2} \exp\left(-\frac{1}{2}\frac{\alpha-2\theta}{1+\alpha-2\theta}|\bar{x}|^2\right) \frac{F(dx_1)}{\|F\|_{\mathcal{M}(\mathbb{R}^d)}} \frac{F(dx_2)}{\|F\|_{\mathcal{M}(\mathbb{R}^d)}}. $$

On the one hand, note that

$$ I_1 \le \exp\left(\frac{\theta R_2^2}{(1+\alpha-2\theta)^2}\right) \int_{\mathbb{R}^{2d}} \exp\left(-\frac{1}{2}\frac{\alpha}{1+\alpha}|\bar{x}|^2\right) \frac{F(dx_1)}{\|F\|_{\mathcal{M}(\mathbb{R}^d)}} \frac{F(dx_2)}{\|F\|_{\mathcal{M}(\mathbb{R}^d)}}, \quad (4.18) $$

where we have used the mean value theorem applied to the function $r \in \mathbb{R}_+^* \mapsto \frac{r}{1+r}$ to derive the relation

$$ \frac{\alpha}{1+\alpha} - \frac{\alpha-2\theta}{1+\alpha-2\theta} \le \frac{2\theta}{(1+\alpha-2\theta)^2}. $$

On the other hand,

$$ I_2 \le \exp\left(-\frac{1}{2}\frac{\alpha-2\theta}{1+\alpha-2\theta}R_2^2\right). \quad (4.19) $$

Putting (4.14), (4.18) and (4.19) into (4.17) we obtain

$$ \int_{\mathbb{R}^d} e^{\theta|x|^2} S[F](x) dx \le \left( \frac{1+\alpha}{1+\alpha-2\theta} \right)^{d/2} \left\{ \exp\left( \frac{\theta R_2^2}{(1+\alpha-2\theta)^2} \right) + \frac{\exp\left( -\frac{1}{2} \frac{\alpha-2\theta}{1+\alpha-2\theta} R_2^2 \right) \exp\left( \frac{1}{2} \frac{\alpha}{1+\alpha} R_1^2 \right)}{1 - \frac{1}{2R_1^2} \int_{\mathbb{R}^d} |x|^2 \frac{F(dx)}{\|F\|_{\mathcal{M}(\mathbb{R}^d)}}} \right\}, $$

for any $R_1, R_2 \in \mathbb{R}_+^*$. As in the proof of Lemma 4.5 we choose them as follows

$$ R_1^2 = \frac{1}{2(1-\delta)} \int_{\mathbb{R}^d} |x|^2 \frac{F(dx)}{\|F\|_{\mathcal{M}(\mathbb{R}^d)}}, \quad R_2^2 = \frac{\alpha}{1+\alpha} \frac{1+\alpha-2\theta}{\alpha-2\theta} R_1^2, $$

where $\delta \in (0, 1)$. Then, we obtain

$$ \int_{\mathbb{R}^d} e^{\theta|x|^2} S[F](x) dx \le \left( \frac{1+\alpha}{1+\alpha-2\theta} \right)^{d/2} \left\{ \frac{1}{\delta} + \exp \left( \frac{\alpha}{1+\alpha} \frac{\theta}{(1+\alpha-2\theta)(\alpha-2\theta)} \frac{1}{2(1-\delta)} \int_{\mathbb{R}^d} |x|^2 \frac{F(dx)}{\|F\|_{\mathcal{M}(\mathbb{R}^d)}} \right) \right\}, $$

and we end the proof by arbitrariness of $\delta \in (0, 1)$. $\square$
Let us note that, on the one hand the factor $\eta^n$ makes the dependence on the variance of the initial datum $F_0$ negligible as the amount of generations increases. On the other hand, the constant $M$ does not depend on the chosen initial profile $F_0$. This is compatible with our ergodicity property in Theorem 1.1.

**Lemma 4.7** (Control of exponential moments). Assume that $\alpha \in \mathbb{R}_+^*$, set $\theta \in \mathbb{R}_+^*$ so that $\theta < \frac{\alpha}{2}$ and any parameter $\chi > \frac{\alpha}{1+\alpha} \frac{\theta}{2(1+\alpha-2\theta)(\alpha-2\theta)}$. Then, there exists a constant $C = C(\alpha, \theta, \chi)$ such that

$$ \int_{\mathbb{R}^d} e^{\theta|x|^2} S[F](x) dx \le C \left\{ 1 + \exp \left( \chi \int_{\mathbb{R}^d} |x|^2 \frac{F(dx)}{\|F\|_{\mathcal{M}(\mathbb{R}^d)}} \right) \right\}, $$

for any measure $F \in \mathcal{M}_+(\mathbb{R}^d)$.

*Proof.* Since $G$ is a Gaussian function given by (1.4) and $m$ is a quadratic function given by (1.5), then we can compute explicitly the integrals with respect to $x$ in the numerator and denominator like formulas (4.11)-(4.12) in the proof of Lemma 4.5 and we obtain

$$ \int_{\mathbb{R}^d} e^{\theta|x|^2} S[F](x) dx = \left( \frac{1+\alpha}{1+\alpha-2\theta} \right)^{d/2} \frac{\int_{\mathbb{R}^{2d}} \exp\left(-\frac{1}{2}\frac{\alpha-2\theta}{1+\alpha-2\theta}|\bar{x}|^2\right) \frac{F(dx_1)}{\|F\|_{\mathcal{M}(\mathbb{R}^d)}} \frac{F(dx_2)}{\|F\|_{\mathcal{M}(\mathbb{R}^d)}}}{\int_{\mathbb{R}^{2d}} \exp\left(-\frac{1}{2}\frac{\alpha}{1+\alpha}|\bar{x}|^2\right) \frac{F(dx_1)}{\|F\|_{\mathcal{M}(\mathbb{R}^d)}} \frac{F(dx_2)}{\|F\|_{\mathcal{M}(\mathbb{R}^d)}}}, \quad (4.17) $$

where $\bar{x} = \frac{x_1+x_2}{2}$. Again, we provide bounds for the numerator and denominator in (4.17). Regarding the denominator, we recover estimate (4.14) for any $R_1 \in \mathbb{R}_+^*$. Thus, we focus on the bound of the numerator

$$ \int_{\mathbb{R}^{2d}} \exp\left(-\frac{1}{2}\frac{\alpha-2\theta}{1+\alpha-2\theta}|\bar{x}|^2\right) \frac{F(dx_1)}{\|F\|_{\mathcal{M}(\mathbb{R}^d)}} \frac{F(dx_2)}{\|F\|_{\mathcal{M}(\mathbb{R}^d)}} = I_1 + I_2, $$

where each term reads

$$ I_1 := \int_{|\bar{x}| \le R_2} \exp\left(-\frac{1}{2}\frac{\alpha-2\theta}{1+\alpha-2\theta}|\bar{x}|^2\right) \frac{F(dx_1)}{\|F\|_{\mathcal{M}(\mathbb{R}^d)}} \frac{F(dx_2)}{\|F\|_{\mathcal{M}(\mathbb{R}^d)}}, $$

$$ I_2 := \int_{|\bar{x}| > R_2} \exp\left(-\frac{1}{2}\frac{\alpha-2\theta}{1+\alpha-2\theta}|\bar{x}|^2\right) \frac{F(dx_1)}{\|F\|_{\mathcal{M}(\mathbb{R}^d)}} \frac{F(dx_2)}{\|F\|_{\mathcal{M}(\mathbb{R}^d)}}. $$

On the one hand, note that

$$ I_1 \le \exp\left(\frac{\theta R_2^2}{(1+\alpha-2\theta)^2}\right) \int_{\mathbb{R}^{2d}} \exp\left(-\frac{1}{2}\frac{\alpha}{1+\alpha}|\bar{x}|^2\right) \frac{F(dx_1)}{\|F\|_{\mathcal{M}(\mathbb{R}^d)}} \frac{F(dx_2)}{\|F\|_{\mathcal{M}(\mathbb{R}^d)}}, \quad (4.18) $$

where we have used the mean value theorem applied to the function $r \in \mathbb{R}_+^* \mapsto \frac{r}{1+r}$ to derive the relation

$$ \frac{\alpha}{1+\alpha} - \frac{\alpha-2\theta}{1+\alpha-2\theta} \le \frac{2\theta}{(1+\alpha-2\theta)^2}. $$

On the other hand,

$$ I_2 \le \exp\left(-\frac{1}{2}\frac{\alpha-2\theta}{1+\alpha-2\theta}R_2^2\right). \quad (4.19) $$

Putting (4.14), (4.18) and (4.19) into (4.17) we obtain

$$ \int_{\mathbb{R}^d} e^{\theta|x|^2} S[F](x) dx \le \left( \frac{1+\alpha}{1+\alpha-2\theta} \right)^{d/2} \left\{ \exp\left( \frac{\theta R_2^2}{(1+\alpha-2\theta)^2} \right) + \frac{\exp\left( -\frac{1}{2}\frac{\alpha-2\theta}{1+\alpha-2\theta}R_2^2 \right) \exp\left( \frac{1}{2}\frac{\alpha}{1+\alpha}R_1^2 \right)}{1 - \frac{1}{2R_1^2} \int_{\mathbb{R}^d} |x|^2 \frac{F(dx)}{\|F\|_{\mathcal{M}(\mathbb{R}^d)}}} \right\}, $$

for any $R_1, R_2 \in \mathbb{R}_+^*$. As in the proof of Lemma 4.5 we choose them as follows

$$ R_1^2 = \frac{1}{2(1-\delta)} \int_{\mathbb{R}^d} |x|^2 \frac{F(dx)}{\|F\|_{\mathcal{M}(\mathbb{R}^d)}}, \quad R_2^2 = \frac{\alpha}{1+\alpha} \frac{1+\alpha-2\theta}{\alpha-2\theta} R_1^2, $$

where $\delta \in (0, 1)$. Then, we obtain

$$ \int_{\mathbb{R}^d} e^{\theta|x|^2} S[F](x) dx \le \left( \frac{1+\alpha}{1+\alpha-2\theta} \right)^{d/2} \left\{ \frac{1}{\delta} + \exp\left( \frac{\alpha}{1+\alpha} \frac{\theta}{(1+\alpha-2\theta)(\alpha-2\theta)} \frac{1}{2(1-\delta)} \int_{\mathbb{R}^d} |x|^2 \frac{F(dx)}{\|F\|_{\mathcal{M}(\mathbb{R}^d)}} \right) \right\}, $$

and we end the proof by arbitrariness of $\delta \in (0, 1)$. $\square$
Therefore, small exponential moments can be controlled by quadratic moments. We can then couple Lemma 4.7 and Corollary 4.6 to obtain the propagation of exponential moments in the time-discrete problem (1.1).

**Corollary 4.8 (Propagation of exponential moments).** Assume that $\alpha \in \mathbb{R}_+^*$, set $\theta \in \mathbb{R}_+^*$ with $\theta < \frac{\alpha}{2}$ and any parameter $\eta \in (\frac{1}{2(1+\alpha)^2}, 1)$. Consider the solution $\{F_n\}_{n \in \mathbb{N}}$ of (1.1) issued at a generic initial datum $F_0 \in \mathcal{M}_+(\mathbb{R}^d)$ with $\int_{\mathbb{R}^d} |x|^2 F_0(dx) < \infty$. Then, there exist $C = C(\alpha, \theta, \eta)$ and $C' = C'(\alpha, \theta, \eta)$ such that

$$ \int_{\mathbb{R}^d} e^{\theta|x|^2} \frac{F_n(x)}{\|F_n\|_{L^1(\mathbb{R}^d)}} dx \le C \left\{ 1 + \exp \left( C' \left( 1 + \eta^{n-1} \int_{\mathbb{R}^d} |x|^2 \frac{F_0(dx)}{\|F_0\|_{\mathcal{M}(\mathbb{R}^d)}} \right) \right) \right\}, $$

for any $n \in \mathbb{N}$.

## 5. REFORMULATING THE RECURSION

In this section we find an appropriate reformulation of the recursion for the solutions $\{F_n\}_{n \in \mathbb{N}}$ of the time-discrete problem (1.1), as already discussed in Section 2.4. Note that by iterating $\mathcal{T}$ n-times we find that $F_n$ depends on the initial datum $F_0 \in \mathcal{M}_+(\mathbb{R}^d)$ via a high-dimensional integral parametrized by variables indexed on the binary tree $\mathbb{T}_*^n$ (representing the traits of the ancestors in the pedigree chart). Specifically, for any initial datum $F_0 \in \mathcal{M}_+(\mathbb{R}^d)$ we obtain:

### *Iteration 1:*

By the explicit definition of $\mathcal{T}$ in (1.2) we readily obtain

$$ F_1(x) = \frac{1}{\|F_0\|_{\mathcal{M}(\mathbb{R}^d)}} \int_{\mathbb{R}^{2d}} e^{-m(x)} G \left( x - \frac{x_1 + x_2}{2} \right) F_0(x_1) F_0(x_2) dx_1, \quad (5.1) $$

where we denote $\mathbf{x}_1 = (x_1, x_2) \in \mathbb{R}^{2d}$.

### *Iteration 2:*

By (5.1) with $F_0$ replaced by $F_1$ and using $F_1 = \mathcal{T}[F_0]$ in the right hand side we get

$$ F_2(x) = \frac{1}{\|F_1\|_{L^1(\mathbb{R}^d)} \|F_0\|_{\mathcal{M}(\mathbb{R}^d)}^2} \int_{\mathbb{R}^{6d}} e^{-(m(x)+m(x_1)+m(x_2))} G \left( x - \frac{x_1+x_2}{2} \right) \\ \times G \left( x_1 - \frac{x_{11}+x_{12}}{2} \right) G \left( x_2 - \frac{x_{21}+x_{22}}{2} \right) F_0(x_{11}) F_0(x_{12}) F_0(x_{21}) F_0(x_{22}) dx_2, \quad (5.2) $$

where we denote $\mathbf{x}_2 = (x_1, x_2, x_{11}, x_{12}, x_{21}, x_{22}) \in \mathbb{R}^{6d}$.

### *Iteration n* $\in \mathbb{N}$:

By a clear inductive process, we can iterate the operator $\mathcal{T}$ as many times $n \in \mathbb{N}$ as needed to recover an explicit dependence of $F_n$ on the initial datum $F_0$. Of course, this generates a high-dimensional integral involving $2(2^n - 1)$ variables that we can label along the vertices of the perfect binary tree $\mathbb{T}^n$ according to the *universal address system* notation in Section 3.1. Namely, we obtain

$$ F_n(x) = \frac{1}{\prod_{m=0}^{n-1} \|F_m\|_{L^1(\mathbb{R}^d)}^{2^{n-1-m}}} \times \int_{\mathbb{R}^{2(2^{n-1})d}} \exp \left( -\sum_{i \in \mathbb{T}^n} m(x_i) \right) \prod_{i \in \mathbb{T}^n} G \left( x_i - \frac{x_{i1} + x_{i2}}{2} \right) \prod_{i \in L^n} F_0(x_i) dx_n, \quad (5.3) $$

where we have used the tree-indexed notation $\mathbf{x}_n = (x_i)_{i \in \mathbb{T}_*^n} \in \mathbb{R}^{2(2^{n-1})d}$ in Remark 3.1 and $x_0 := x$.

The ultimate goal of this section is to find a suitable change of variables as discussed in Section 2.4. We shall see that the effect of such a change of variable on the above high-dimensional integral in (5.3) will reveal that actual dependence of $\{F_n\}_{n \in \mathbb{N}}$ on the shape of $F_0$. We will see that such a dependence is indeed weak and the initial datum is rapidly “forgotten” across the different levels of the tree $\mathbb{T}^n$. This ergodicity property will be crucial in our analysis and will be exploited later in Section 6 to derive our main result in Theorem 1.1.
Therefore, small exponential moments can be controlled by quadratic moments. We can then couple Lemma 4.7 and Corollary 4.6 to obtain the propagation of exponential moments in the time-discrete problem (1.1).

**Corollary 4.8 (Propagation of exponential moments).** Assume that $\alpha \in \mathbb{R}_+^*$, set $\theta \in \mathbb{R}_+^*$ with $\theta < \frac{\alpha}{2}$ and any parameter $\eta \in (\frac{1}{2(1+\alpha)^2}, 1)$. Consider the solution $\{F_n\}_{n \in \mathbb{N}}$ of (1.1) issued at a generic initial datum $F_0 \in \mathcal{M}_+(\mathbb{R}^d)$ with $\int_{\mathbb{R}^d} |x|^2 F_0(dx) < \infty$. Then, there exist $C = C(\alpha, \theta, \eta)$ and $C' = C'(\alpha, \theta, \eta)$ such that

$$ \int_{\mathbb{R}^d} e^{\theta|x|^2} \frac{F_n(x)}{\|F_n\|_{L^1(\mathbb{R}^d)}} dx \le C \left\{ 1 + \exp \left( C' \left( 1 + \eta^{n-1} \int_{\mathbb{R}^d} |x|^2 \frac{F_0(dx)}{\|F_0\|_{\mathcal{M}(\mathbb{R}^d)}} \right) \right) \right\}, $$

for any $n \in \mathbb{N}$.

## 5. REFORMULATING THE RECURSION

In this section we find an appropriate reformulation of the recursion for the solutions $\{F_n\}_{n \in \mathbb{N}}$ of the time-discrete problem (1.1), as already discussed in Section 2.4. Note that by iterating $\mathcal{T}$ n-times we find that $F_n$ depends on the initial datum $F_0 \in \mathcal{M}_+(\mathbb{R}^d)$ via a high-dimensional integral parametrized by variables indexed on the binary tree $\mathbb{T}_*^n$ (representing the traits of the ancestors in the pedigree chart). Specifically, for any initial datum $F_0 \in \mathcal{M}_+(\mathbb{R}^d)$ we obtain:

### *Iteration 1:*

By the explicit definition of $\mathcal{T}$ in (1.2) we readily obtain

$$ F_1(x) = \frac{1}{\|F_0\|_{\mathcal{M}(\mathbb{R}^d)}} \int_{\mathbb{R}^{2d}} e^{-m(x)} G \left( x - \frac{x_1 + x_2}{2} \right) F_0(x_1) F_0(x_2) dx_1, \quad (5.1) $$

where we denote $\mathbf{x}_1 = (x_1, x_2) \in \mathbb{R}^{2d}$.

### *Iteration 2:*

By (5.1) with $F_0$ replaced by $F_1$ and using $F_1 = \mathcal{T}[F_0]$ in the right hand side we get

$$ F_2(x) = \frac{1}{\|F_1\|_{L^1(\mathbb{R}^d)} \|F_0\|_{\mathcal{M}(\mathbb{R}^d)}^2} \int_{\mathbb{R}^{6d}} e^{-(m(x)+m(x_1)+m(x_2))} G \left( x - \frac{x_1+x_2}{2} \right) \\ \times G \left( x_1 - \frac{x_{11}+x_{12}}{2} \right) G \left( x_2 - \frac{x_{21}+x_{22}}{2} \right) F_0(x_{11}) F_0(x_{12}) F_0(x_{21}) F_0(x_{22}) dx_2, \quad (5.2) $$

where we denote $\mathbf{x}_2 = (x_1, x_2, x_{11}, x_{12}, x_{21}, x_{22}) \in \mathbb{R}^{6d}$.

### *Iteration n* ∈ $\mathbb{N}$:

By a clear inductive process, we can iterate the operator $\mathcal{T}$ as many times $n \in \mathbb{N}$ as needed to recover an explicit dependence of $F_n$ on the initial datum $F_0$. Of course, this generates a high-dimensional integral involving $2(2^n - 1)$ variables that we can label along the vertices of the perfect binary tree $\mathbb{T}^n$ according to the *universal address system* notation in Section 3.1. Namely, we obtain

$$ F_n(x) = \frac{1}{\prod_{m=0}^{n-1} \|F_m\|_{L^1(\mathbb{R}^d)}^{2^{n-1-m}}} \\ \times \int_{\mathbb{R}^{2(2^{n-1})d}} \exp \left( -\sum_{i \in \mathbb{T}^n} m(x_i) \right) \prod_{i \in \mathbb{T}^n} G \left( x_i - \frac{x_{i1} + x_{i2}}{2} \right) \prod_{i \in L^n} F_0(x_i) dx_n, \quad (5.3) $$

where we have used the tree-indexed notation $\mathbf{x}_n = (x_i)_{i \in \mathbb{T}_*^n} \in \mathbb{R}^{2(2^{n-1})d}$ in Remark 3.1 and $x_0 := x$.

The ultimate goal of this section is to find a suitable change of variables as discussed in Section 2.4. We shall see that the effect of such a change of variable on the above high-dimensional integral in (5.3) will reveal that actual dependence of $\{F_n\}_{n \in \mathbb{N}}$ on the shape of $F_0$. We will see that such a dependence is indeed weak and the initial datum is rapidly “forgotten” across the different levels of the tree $\mathbb{T}^n$. This ergodicity property will be crucial in our analysis and will be exploited later in Section 6 to derive our main result in Theorem 1.1.
5.1. **A change of variables across the binary tree.** The goal of this section is to appropriately reformulate the high-dimensional integral (5.3) by exploiting the special form of the Gaussian mixing kernel $G$ in (1.4) and the quadratic selection function $m$ in (1.5). Namely, we shall derive a suitable change of variables according to a $n$-step backwards process starting at the leaf variables $x_j$ with $j \in \mathbb{L}^n$ and ending at the root variable $x_0 = x$. More specifically, at each level $m$ in the tree we shall change the reference frame of the variables $x_i$ indexed with $i \in \mathbb{L}_m^n$ so that the quadratic form involving those variables in the exponential of (5.3) gets appropriately centered at its minimum. By doing so we find that the minimum is located at a contracted value of the variable $x_{c(i)} \in \mathbb{L}_{m-1}^n$ at the level $m-1$ below, which is indexed by its children $c(i)$ in the binary tree (recall notation in Section 3.1). By appropriately propagating the backwards process across the tree, and by tracking the accumulated contraction of variables, we shall see that the eventual dependence on the root variable $x$ is weakly gradually forgotten as $n \to \infty$. This suggests an apparent form of ergodicity, which we will be crucial later in Section 6.

First, for the sake of clarity, we illustrate our method in the simpler case $n = 2$. Later, we address the general case with $n > 2$ driven by the same strategy. As proposed in Section 2.4, we define the following choice for the rescaled profiles that will be used along the change of variables.

**Definition 5.1** (Rescaled distributions). Consider any $F \in \mathcal{M}_+(\mathbb{R}^d)$. Then, we define

$$ \bar{F}(dx) := e^{m(x)} \frac{F(dx)}{\bar{F}_{\alpha=0}(x)}, \quad x \in \mathbb{R}^d. $$

Here, we remind that $\bar{F}_{\alpha=0} = G_{0,2}$ is the Gaussian eigenfunction corresponding to (1.12) with $\alpha = 0$. We also recall that there is some freedom in the normalization, and in particular the term $e^m$ is not mandatory but it is convenient for an easier sorting of the various terms, as we anticipated in Section 2.4. By substituting (1.4) and (1.5) into (5.1) and writing the integral in the right hand side in terms of the rescaled function $\bar{F}_0$ in Definition 5.1, we obtain that

$$ \begin{aligned} (e^m F_2)(x) = & \frac{1}{(2\pi)^{3d/2}} \frac{1}{(4\pi)^{2d}} \frac{1}{\|F_1\|_{L^1(\mathbb{R}^d)} \|F_0\|_{\mathcal{M}(\mathbb{R}^d)}^2} \\ & \times \int_{\mathbb{R}^{6d}} \exp \left( -\frac{\alpha}{2} (|x_1|^2 + |x_2|^2) - \frac{1}{2} \left| x - \frac{x_1+x_2}{2} \right|^2 \right) \\ & \times \exp \left( -\frac{\alpha_0}{2} (|x_{11}|^2 + |x_{12}|^2 + |x_{21}|^2 + |x_{22}|^2) - \frac{1}{2} \left| x_1 - \frac{x_{11}+x_{12}}{2} \right|^2 - \frac{1}{2} \left| x_2 - \frac{x_{21}+x_{22}}{2} \right|^2 \right) \\ & \times \bar{F}_0(x_{11}) \bar{F}_0(x_{12}) \bar{F}_0(x_{21}) \bar{F}_0(x_{22}) \, dx_2. \end{aligned} \quad (5.4) $$

Here, the parameter $\alpha_0 \in \mathbb{R}_+^*$ has been defined by

$$ \alpha_0 := \frac{1}{2} + \alpha, $$

so that the quadratic terms in the first term of the third line of (5.4) correct the rescaling $\bar{F}_0$ of $F_0$. The new formulation will be obtained by an appropriate change of variables, so that the quadratic forms inside the above exponential get appropriately centered. To do so, notice that this formula involves as many variables $x_i$ as indices $i$ in the perfect binary tree $\mathbb{T}^2$. Indeed, we have sorted the different terms in (5.4) in such a way that the second line only involves variables $x_1$ and $x_2$ in the first level $\mathbb{L}_1^2$ of the tree, whilst the third line contains the terms involving the variables $x_{11}, x_{12}, x_{21}$ and $x_{22}$ in the second level of the tree $\mathbb{L}_2^2$ (i.e., the leaves $\mathbb{L}^2$). We will divide the method into two steps. First, we will address the change of the variables indexed by the leaves $\mathbb{L}^2$. Second, we will perform the change of variables indexed by the first level of the tree $\mathbb{L}_1^2$.

• **STEP 1:** Change of variables for $x_{11}, x_{12}, x_{21}$ and $x_{22}$.

Consider some coefficient $k_1 > 0$ to be determined later and define the change of variables $x_{11} \to y_{11}$, $x_{12} \to y_{12}$, $x_{21} \to y_{21}$ and $x_{22} \to y_{22}$ given by

$$ x_{11} = k_1 x_1 + y_{11}, \qquad x_{12} = k_1 x_1 + y_{12}, $$

$$ x_{21} = k_1 x_2 + y_{21}, \qquad x_{22} = k_1 x_2 + y_{22}. $$
5.1. **A change of variables across the binary tree.** The goal of this section is to appropriately reformulate the high-dimensional integral (5.3) by exploiting the special form of the Gaussian mixing kernel $G$ in (1.4) and the quadratic selection function $m$ in (1.5). Namely, we shall derive a suitable change of variables according to a $n$-step backwards process starting at the leaf variables $x_j$ with $j \in \mathbb{L}^n$ and ending at the root variable $x_0 = x$. More specifically, at each level $m$ in the tree we shall change the reference frame of the variables $x_i$ indexed with $i \in \mathbb{L}_m^n$ so that the quadratic form involving those variables in the exponential of (5.3) gets appropriately centered at its minimum. By doing so we find that the minimum is located at a contracted value of the variable $x_{c(i)} \in \mathbb{L}_{m-1}^n$ at the level $m-1$ below, which is indexed by its children $c(i)$ in the binary tree (recall notation in Section 3.1). By appropriately propagating the backwards process across the tree, and by tracking the accumulated contraction of variables, we shall see that the eventual dependence on the root variable $x$ is weakly gradually forgotten as $n \to \infty$. This suggests an apparent form of ergodicity, which we will be crucial later in Section 6.

First, for the sake of clarity, we illustrate our method in the simpler case $n = 2$. Later, we address the general case with $n > 2$ driven by the same strategy. As proposed in Section 2.4, we define the following choice for the rescaled profiles that will be used along the change of variables.

**Definition 5.1** (Rescaled distributions). Consider any $F \in \mathcal{M}_+(\mathbb{R}^d)$. Then, we define

$$ \bar{F}(dx) := e^{m(x)} \frac{F(dx)}{\bar{F}_{\alpha=0}(x)}, \quad x \in \mathbb{R}^d. $$

Here, we remind that $F_{\alpha=0} = G_{0,2}$ is the Gaussian eigenfunction corresponding to (1.12) with $\alpha = 0$. We also recall that there is some freedom in the normalization, and in particular the term $e^m$ is not mandatory but it is convenient for an easier sorting of the various terms, as we anticipated in Section 2.4. By substituting (1.4) and (1.5) into (5.1) and writing the integral in the right hand side in terms of the rescaled function $\bar{F}_0$ in Definition 5.1, we obtain that

$$ \begin{aligned} (e^m F_2)(x) = & \frac{1}{(2\pi)^{3d/2}} \frac{1}{(4\pi)^{2d}} \frac{1}{\|F_1\|_{L^1(\mathbb{R}^d)} \|F_0\|_{\mathcal{M}(\mathbb{R}^d)}^2} \\ & \times \int_{\mathbb{R}^{6d}} \exp \left( -\frac{\alpha}{2} (|x_1|^2 + |x_2|^2) - \frac{1}{2} \left| x - \frac{x_1+x_2}{2} \right|^2 \right) \\ & \times \exp \left( -\frac{\alpha_0}{2} (|x_{11}|^2 + |x_{12}|^2 + |x_{21}|^2 + |x_{22}|^2) - \frac{1}{2} \left| x_1 - \frac{x_{11}+x_{12}}{2} \right|^2 - \frac{1}{2} \left| x_2 - \frac{x_{21}+x_{22}}{2} \right|^2 \right) \\ & \times \bar{F}_0(x_{11}) \bar{F}_0(x_{12}) \bar{F}_0(x_{21}) \bar{F}_0(x_{22}) \, dx_2. \end{aligned} \quad (5.4) $$

Here, the parameter $\alpha_0 \in \mathbb{R}_+^*$ has been defined by

$$ \alpha_0 := \frac{1}{2} + \alpha, $$

so that the quadratic terms in the first term of the third line of (5.4) correct the rescaling $\bar{F}_0$ of $F_0$. The new formulation will be obtained by an appropriate change of variables, so that the quadratic forms inside the above exponential get appropriately centered. To do so, notice that this formula involves as many variables $x_i$ as indices $i$ in the perfect binary tree $\mathbb{T}^2$. Indeed, we have sorted the different terms in (5.4) in such a way that the second line only involves variables $x_1$ and $x_2$ in the first level $\mathbb{L}_1^2$ of the tree, whilst the third line contains the terms involving the variables $x_{11}, x_{12}, x_{21}$ and $x_{22}$ in the second level of the tree $\mathbb{L}_2^2$ (i.e., the leaves $\mathbb{L}^2$). We will divide the method into two steps. First, we will address the change of the variables indexed by the leaves $\mathbb{L}^2$. Second, we will perform the change of variables indexed by the first level of the tree $\mathbb{L}_1^2$.

• **STEP 1:** Change of variables for $x_{11}, x_{12}, x_{21}$ and $x_{22}$.

Consider some coefficient $k_1 > 0$ to be determined later and define the change of variables $x_{11} \to y_{11}$, $x_{12} \to y_{12}$, $x_{21} \to y_{21}$ and $x_{22} \to y_{22}$ given by

$$ x_{11} = k_1 x_1 + y_{11}, \qquad x_{12} = k_1 x_1 + y_{12}, $$

$$ x_{21} = k_1 x_2 + y_{21}, \qquad x_{22} = k_1 x_2 + y_{22}. $$
On the one hand, using the change of variables in the terms inside the exponential of (5.4) which involve $x_{11}$ and $x_{12}$, we obtain

$$ 
\begin{aligned}
& \frac{\alpha_0}{2}|x_{11}|^2 + \frac{\alpha_0}{2}|x_{12}|^2 + \frac{1}{2}\left|x_1 - \frac{x_{11} + x_{12}}{2}\right|^2 \\
&= \frac{\alpha_0}{2}|k_1x_1 + y_{11}|^2 + \frac{\alpha_0}{2}|k_1x_1 + y_{12}|^2 + \frac{1}{2}\left|(1-k_1)x_1 - \frac{y_{11} + y_{12}}{2}\right|^2 \\
&= \frac{1}{2}\left(2\alpha_0k_1^2 + (1-k_1)^2\right)|x_1|^2 + \frac{\alpha_0}{2}|y_{11}|^2 + \frac{\alpha_0}{2}|y_{12}|^2 + \frac{1}{8}|y_{11} + y_{12}|^2 \\
&\quad + \left(\alpha_0k_1 - \frac{1-k_1}{2}\right)x_1 \cdot y_{11} + \left(\alpha_0k_1 - \frac{1-k_1}{2}\right)x_1 \cdot y_{12}.
\end{aligned}
$$

On the other hand, using the change of variables in the terms which involve $x_{21}$ and $x_{22}$, we get

$$ 
\begin{aligned}
& \frac{\alpha_0}{2}|x_{21}|^2 + \frac{\alpha_0}{2}|x_{22}|^2 + \frac{1}{2}\left|x_2 - \frac{x_{21} + x_{22}}{2}\right|^2 \\
&= \frac{\alpha_0}{2}|k_1x_2 + y_{21}|^2 + \frac{\alpha_0}{2}|k_1x_2 + y_{22}|^2 + \frac{1}{2}\left|(1-k_1)x_2 - \frac{y_{21} + y_{22}}{2}\right|^2 \\
&= \frac{1}{2}\left(2\alpha_0k_1^2 + (1-k_1)^2\right)|x_2|^2 + \frac{\alpha_0}{2}|y_{21}|^2 + \frac{\alpha_0}{2}|y_{22}|^2 + \frac{1}{8}|y_{21} + y_{22}|^2 \\
&\quad + \left(\alpha_0k_1 - \frac{1-k_1}{2}\right)x_2 \cdot y_{21} + \left(\alpha_0k_1 - \frac{1-k_1}{2}\right)x_2 \cdot y_{22}.
\end{aligned}
$$

Notice that one can eliminate the crossed terms in both expressions by choosing

$$ k_1 := \frac{1}{1 + 2\alpha_0} = \frac{1}{2(1 + \alpha)}. $$

In that case, adding both terms we obtain that

$$ 
\begin{aligned}
& \frac{\alpha_0}{2}|x_{11}|^2 + \frac{\alpha_0}{2}|x_{12}|^2 + \frac{1}{2}\left|x_1 - \frac{x_{11} + x_{12}}{2}\right|^2 + \frac{\alpha_0}{2}|x_{21}|^2 + \frac{\alpha_0}{2}|x_{22}|^2 + \frac{1}{2}\left|x_2 - \frac{x_{21} + x_{22}}{2}\right|^2 \\
&= \frac{1}{2}(1-k_1)|x_1|^2 + \frac{1}{4k_1}|y_{11}|^2 + \frac{1}{4k_1}|y_{12}|^2 - \frac{1}{8}|y_{11} - y_{12}|^2 \\
&\quad + \frac{1}{2}(1-k_1)|x_2|^2 + \frac{1}{4k_1}|y_{21}|^2 + \frac{1}{4k_1}|y_{22}|^2 - \frac{1}{8}|y_{21} - y_{22}|^2.
\end{aligned}
$$

Putting everything together into (5.4) yields

$$ 
\begin{aligned}
(e^m F_2)(x) = & \frac{1}{(2\pi)^{3d/2}} \frac{1}{(4\pi)^{2d}} \frac{1}{\|F_1\|_{L^1(\mathbb{R}^d)} \|F_0\|_{\mathcal{M}(\mathbb{R}^d)}^2} \\
& \times \int_{\mathbb{R}^{6d}} \exp \left( -\frac{\alpha_1}{2} (|x_1|^2 + |x_2|^2) - \frac{1}{2} \left| x - \frac{x_1 + x_2}{2} \right|^2 \right) \\
& \times \exp \left( -\frac{1}{4k_1} (|y_{11}|^2 + |y_{12}|^2 + |y_{21}|^2 + |y_{22}|^2) + \frac{1}{8} (|y_{11} - y_{12}|^2 + |y_{21} - y_{22}|^2) \right) \\
& \times \bar{F}_0(k_1 x_1 + y_{11}) \bar{F}_0(k_1 x_1 + y_{12}) \bar{F}_0(k_1 x_2 + y_{21}) \bar{F}_0(k_1 x_2 + y_{22}) dx_1 dx_2 dy_{11} dy_{12} dy_{21} dy_{22},
\end{aligned}
\quad (5.5) $$

where the parameter $\alpha_1 \in \mathbb{R}_+^*$ has been defined by

$$ \alpha_1 := 1 - k_1 + \alpha, $$

in order to recombine the initial terms in the second line of (5.4) involving variables indexed by the first level $L_1^2$ with the new remainders that have appeared from the previous step.

• STEP 2: Change of variables for $x_1$ and $x_2$.

Consider some coefficient $k_2 > 0$ to be determined later and define the change of variables $x_1 \to y_1$ and $x_2 \to y_2$ given by

$$ x_1 = k_2x + y_1, \quad x_2 = k_2x + y_2. $$
On the one hand, using the change of variables in the terms inside the exponential of (5.4) which involve $x_{11}$ and $x_{12}$, we obtain

$$ 
\begin{aligned}
& \frac{\alpha_0}{2}|x_{11}|^2 + \frac{\alpha_0}{2}|x_{12}|^2 + \frac{1}{2}\left|x_1 - \frac{x_{11} + x_{12}}{2}\right|^2 \\
&= \frac{\alpha_0}{2}|k_1x_1 + y_{11}|^2 + \frac{\alpha_0}{2}|k_1x_1 + y_{12}|^2 + \frac{1}{2}\left|(1-k_1)x_1 - \frac{y_{11} + y_{12}}{2}\right|^2 \\
&= \frac{1}{2}(2\alpha_0 k_1^2 + (1-k_1)^2)|x_1|^2 + \frac{\alpha_0}{2}|y_{11}|^2 + \frac{\alpha_0}{2}|y_{12}|^2 + \frac{1}{8}|y_{11} + y_{12}|^2 \\
&\quad + \left(\alpha_0 k_1 - \frac{1-k_1}{2}\right)x_1 \cdot y_{11} + \left(\alpha_0 k_1 - \frac{1-k_1}{2}\right)x_1 \cdot y_{12}.
\end{aligned}
$$

On the other hand, using the change of variables in the terms which involve $x_{21}$ and $x_{22}$, we get

$$ 
\begin{aligned}
& \frac{\alpha_0}{2}|x_{21}|^2 + \frac{\alpha_0}{2}|x_{22}|^2 + \frac{1}{2}\left|x_2 - \frac{x_{21} + x_{22}}{2}\right|^2 \\
&= \frac{\alpha_0}{2}|k_1x_2 + y_{21}|^2 + \frac{\alpha_0}{2}|k_1x_2 + y_{22}|^2 + \frac{1}{2}\left|(1-k_1)x_2 - \frac{y_{21} + y_{22}}{2}\right|^2 \\
&= \frac{1}{2}(2\alpha_0 k_1^2 + (1-k_1)^2)|x_2|^2 + \frac{\alpha_0}{2}|y_{21}|^2 + \frac{\alpha_0}{2}|y_{22}|^2 + \frac{1}{8}|y_{21} + y_{22}|^2 \\
&\quad + \left(\alpha_0 k_1 - \frac{1-k_1}{2}\right)x_2 \cdot y_{21} + \left(\alpha_0 k_1 - \frac{1-k_1}{2}\right)x_2 \cdot y_{22}.
\end{aligned}
$$

Notice that one can eliminate the crossed terms in both expressions by choosing

$$ k_1 := \frac{1}{1 + 2\alpha_0} = \frac{1}{2(1 + \alpha)}. $$

In that case, adding both terms we obtain that

$$ 
\begin{aligned}
& \frac{\alpha_0}{2}|x_{11}|^2 + \frac{\alpha_0}{2}|x_{12}|^2 + \frac{1}{2}\left|x_1 - \frac{x_{11} + x_{12}}{2}\right|^2 + \frac{\alpha_0}{2}|x_{21}|^2 + \frac{\alpha_0}{2}|x_{22}|^2 + \frac{1}{2}\left|x_2 - \frac{x_{21} + x_{22}}{2}\right|^2 \\
&= \frac{1}{2}(1 - k_1)|x_1|^2 + \frac{1}{4k_1}|y_{11}|^2 + \frac{1}{4k_1}|y_{12}|^2 - \frac{1}{8}|y_{11} - y_{12}|^2 \\
&\quad + \frac{1}{2}(1 - k_1)|x_2|^2 + \frac{1}{4k_1}|y_{21}|^2 + \frac{1}{4k_1}|y_{22}|^2 - \frac{1}{8}|y_{21} - y_{22}|^2.
\end{aligned}
$$

Putting everything together into (5.4) yields

$$ 
\begin{aligned}
(e^m F_2)(x) = & \frac{1}{(2\pi)^{3d/2}} \frac{1}{(4\pi)^{2d}} \frac{1}{\|F_1\|_{L^1(\mathbb{R}^d)} \|F_0\|_{\mathcal{M}(\mathbb{R}^d)}^2} \\
& \times \int_{\mathbb{R}^{6d}} \exp \left( -\frac{\alpha_1}{2} (|x_1|^2 + |x_2|^2) - \frac{1}{2} \left| x - \frac{x_1 + x_2}{2} \right|^2 \right) \\
& \times \exp \left( -\frac{1}{4k_1} (|y_{11}|^2 + |y_{12}|^2 + |y_{21}|^2 + |y_{22}|^2) + \frac{1}{8} (|y_{11} - y_{12}|^2 + |y_{21} - y_{22}|^2) \right) \\
& \times \bar{F}_0(k_1 x_1 + y_{11}) \bar{F}_0(k_1 x_1 + y_{12}) \bar{F}_0(k_1 x_2 + y_{21}) \bar{F}_0(k_1 x_2 + y_{22}) dx_1 dx_2 dy_{11} dy_{12} dy_{21} dy_{22},
\end{aligned}
\quad (5.5) $$

where the parameter $\alpha_1 \in \mathbb{R}_+^*$ has been defined by

$$ \alpha_1 := 1 - k_1 + \alpha, $$

in order to recombine the initial terms in the second line of (5.4) involving variables indexed by the first level $L_1^2$ with the new remainders that have appeared from the previous step.

• STEP 2: Change of variables for $x_1$ and $x_2$.

Consider some coefficient $k_2 > 0$ to be determined later and define the change of variables $x_1 \to y_1$ and $x_2 \to y_2$ given by

$$ x_1 = k_2x + y_1, \quad x_2 = k_2x + y_2. $$
This time, using the change of variables in the terms of the exponential in the second line of (5.5) yields

$$ 
\begin{aligned}
& \frac{\alpha_1}{2}|x_1|^2 + \frac{\alpha_1}{2}|x_2|^2 + \frac{1}{2}\left|x - \frac{x_1+x_2}{2}\right|^2 \\
&= \frac{\alpha_1}{2}|k_2x + y_1|^2 + \frac{\alpha_1}{2}|k_2x + y_2|^2 + \frac{1}{2}\left|(1-k_2)x - \frac{y_1+y_2}{2}\right|^2 \\
&= \frac{1}{2}(2\alpha_1k_2^2 + (1-k_2)^2)|x|^2 + \frac{\alpha_1}{2}|y_1|^2 + \frac{\alpha_1}{2}|y_2|^2 + \frac{1}{8}|y_1+y_2|^2 \\
&\quad + \left(\alpha_1k_2 - \frac{1-k_2}{2}\right)x \cdot y_1 + \left(\alpha_1k_2 - \frac{1-k_2}{2}\right)x \cdot y_2.
\end{aligned}
$$

Again, we can cancel the crossed term by choosing

$$ k_2 := \frac{1}{1 + 2\alpha_1} = \frac{1}{3 + 2\alpha - 2k_1}. $$

Namely, we obtain that

$$ 
\begin{aligned}
& \frac{\alpha_1}{2}|x_1|^2 + \frac{\alpha_1}{2}|x_2|^2 + \frac{1}{2}\left|x - \frac{x_1+x_2}{2}\right|^2 \\
&= \frac{1}{2}(1-k_2)|x|^2 + \frac{1}{4k_2}|y_1|^2 + \frac{1}{4k_2}|y_2|^2 - \frac{1}{8}|y_1 - y_2|^2.
\end{aligned}
$$

Then, putting everything together into (5.5) yields

$$ 
\begin{aligned}
F_2(x) = & \frac{1}{(2\pi)^{3d/2}} \frac{1}{(4\pi)^{2d}} \frac{1}{\|F_1\|_{L^1(\mathbb{R}^d)} \|F_0\|_{\mathcal{M}(\mathbb{R}^d)}^2} e^{-\frac{1+\alpha-k_2}{2}|x|^2} \\
& \times \int_{\mathbb{R}^{6d}} \exp \left( -\frac{1}{4k_2} (|y_1|^2 + |y_2|^2) + \frac{1}{8} |y_1 - y_2|^2 \right) \\
& \times \exp \left( -\frac{1}{4k_1} (|y_{11}|^2 + |y_{12}|^2 + |y_{21}|^2 + |y_{22}|^2) + \frac{1}{8} (|y_{11} - y_{12}|^2 + |y_{21} - y_{22}|^2) \right) \\
& \times \bar{F}_0(k_1 k_2 x + k_1 y_1 + y_{11}) \bar{F}_0(k_1 k_2 x + k_1 y_1 + y_{12}) \\
& \times \bar{F}_0(k_1 k_2 x + k_1 y_2 + y_{21}) \bar{F}_0(k_1 k_2 x + k_1 y_2 + y_{22}) dy_2,
\end{aligned}
\quad (5.6) $$

where we denote again $\mathbf{y}_2 = (y_1, y_2, y_{11}, y_{12}, y_{21}, y_{22}) \in \mathbb{R}^{6d}$.

Before stating the main result for general $n \in \mathbb{N}$, we collect some natural notation according to the preceding computations, that will be useful here on. First, we define the following sequences of coefficients.

**Definition 5.2 (Coefficients).** Consider any $\alpha \in \mathbb{R}_+$.

(1) The coefficients $\{k_n\}_{n \in \mathbb{N}}$ are defined by the recursive formula

$$ 
\begin{aligned}
k_1 &:= \frac{1}{2(1+\alpha)}, \\
k_n &:= \frac{1}{3+2\alpha-2k_{n-1}}, \quad \text{for } n \ge 2.
\end{aligned}
\quad (5.7) $$

(2) The coefficients $\{\kappa_n\}_{n \in \mathbb{N}}$ are defined by the recursive formula

$$ 
\begin{aligned}
\kappa_0 &:= 1, \\
\kappa_n &:= k_1 \cdots k_n, \quad \text{for } n \ge 1.
\end{aligned}
\quad (5.8) $$

We note that for $n = 2$ the above coefficients $k_1, k_2$ reduce to those appearing in the previous reformulation (5.6) of the recursion. Since it will be used later, we study the asymptotic behavior of such sequences of coefficients as $n \to \infty$.

**Lemma 5.3 (Asymptotic behavior of the coefficients).** Consider the sequences $\{k_n\}_{n \in \mathbb{N}}$ and $\{\kappa_n\}_{n \in \mathbb{N}}$ in Definition 5.2. Then, the following properties hold true:

(i) (Coefficients $k_n$) The coefficients $\{k_n\}_{n \in \mathbb{N}}$ are positive numbers and $k_n \searrow k_\alpha$ as $n \to \infty$, where $k_\alpha \in \mathbb{R}_+^*$ is the smallest root of the equation

$$ \frac{1}{3 + 2\alpha - 2k_{\alpha}} = k_{\alpha}. \quad (5.9) $$
This time, using the change of variables in the terms of the exponential in the second line of (5.5) yields

$$ 
\begin{aligned}
& \frac{\alpha_1}{2} |x_1|^2 + \frac{\alpha_1}{2} |x_2|^2 + \frac{1}{2} \left| x - \frac{x_1 + x_2}{2} \right|^2 \\
&= \frac{\alpha_1}{2} |k_2 x + y_1|^2 + \frac{\alpha_1}{2} |k_2 x + y_2|^2 + \frac{1}{2} \left| (1-k_2)x - \frac{y_1+y_2}{2} \right|^2 \\
&= \frac{1}{2} (2\alpha_1 k_2^2 + (1-k_2)^2) |x|^2 + \frac{\alpha_1}{2} |y_1|^2 + \frac{\alpha_1}{2} |y_2|^2 + \frac{1}{8} |y_1+y_2|^2 \\
&\quad + \left( \alpha_1 k_2 - \frac{1-k_2}{2} \right) x \cdot y_1 + \left( \alpha_1 k_2 - \frac{1-k_2}{2} \right) x \cdot y_2.
\end{aligned}
$$

Again, we can cancel the crossed term by choosing

$$ k_2 := \frac{1}{1 + 2\alpha_1} = \frac{1}{3 + 2\alpha - 2k_1}. $$

Namely, we obtain that

$$ 
\begin{aligned}
& \frac{\alpha_1}{2} |x_1|^2 + \frac{\alpha_1}{2} |x_2|^2 + \frac{1}{2} \left| x - \frac{x_1 + x_2}{2} \right|^2 \\
&= \frac{1}{2} (1 - k_2) |x|^2 + \frac{1}{4k_2} |y_1|^2 + \frac{1}{4k_2} |y_2|^2 - \frac{1}{8} |y_1 - y_2|^2.
\end{aligned}
$$

Then, putting everything together into (5.5) yields

$$ 
\begin{aligned}
F_2(x) = & \frac{1}{(2\pi)^{3d/2}} \frac{1}{(4\pi)^{2d}} \frac{1}{\|F_1\|_{L^1(\mathbb{R}^d)} \|F_0\|_{\mathcal{M}(\mathbb{R}^d)}^2} e^{-\frac{1+\alpha-k_2}{2} |x|^2} \\
& \times \int_{\mathbb{R}^{6d}} \exp \left( -\frac{1}{4k_2} (|y_1|^2 + |y_2|^2) + \frac{1}{8} |y_1 - y_2|^2 \right) \\
& \times \exp \left( -\frac{1}{4k_1} (|y_{11}|^2 + |y_{12}|^2 + |y_{21}|^2 + |y_{22}|^2) + \frac{1}{8} (|y_{11} - y_{12}|^2 + |y_{21} - y_{22}|^2) \right) \\
& \times \bar{F}_0(k_1 k_2 x + k_1 y_1 + y_{11}) \bar{F}_0(k_1 k_2 x + k_1 y_1 + y_{12}) \\
& \times \bar{F}_0(k_1 k_2 x + k_1 y_2 + y_{21}) \bar{F}_0(k_1 k_2 x + k_1 y_2 + y_{22}) dy_2,
\end{aligned}
\quad (5.6) $$

where we denote again $\mathbf{y}_2 = (y_1, y_2, y_{11}, y_{12}, y_{21}, y_{22}) \in \mathbb{R}^{6d}$.

Before stating the main result for general $n \in \mathbb{N}$, we collect some natural notation according to the preceding computations, that will be useful here on. First, we define the following sequences of coefficients.

**Definition 5.2 (Coefficients).** Consider any $\alpha \in \mathbb{R}_+$.

(1) The coefficients $\{k_n\}_{n \in \mathbb{N}}$ are defined by the recursive formula

$$ 
\begin{aligned}
k_1 &:= \frac{1}{2(1+\alpha)}, \\
k_n &:= \frac{1}{3+2\alpha-2k_{n-1}}, \quad \text{for } n \ge 2.
\end{aligned}
\quad (5.7) $$

(2) The coefficients $\{\kappa_n\}_{n \in \mathbb{N}}$ are defined by the recursive formula

$$ 
\begin{aligned}
\kappa_0 &:= 1, \\
\kappa_n &:= k_1 \cdots k_n, \quad \text{for } n \ge 1.
\end{aligned}
\quad (5.8) $$

We note that for $n = 2$ the above coefficients $k_1, k_2$ reduce to those appearing in the previous reformulation (5.6) of the recursion. Since it will be used later, we study the asymptotic behavior of such sequences of coefficients as $n \to \infty$.

**Lemma 5.3 (Asymptotic behavior of the coefficients).** Consider the sequences $\{k_n\}_{n \in \mathbb{N}}$ and $\{\kappa_n\}_{n \in \mathbb{N}}$ in Definition 5.2. Then, the following properties hold true:

(i) (Coefficients $k_n$) The coefficients $\{k_n\}_{n \in \mathbb{N}}$ are positive numbers and $k_n \searrow k_\alpha$ as $n \to \infty$, where $k_\alpha \in \mathbb{R}_+^*$ is the smallest root of the equation

$$ \frac{1}{3 + 2\alpha - 2k_{\alpha}} = k_{\alpha}. \quad (5.9) $$
Specifically, $\mathbf{k}_\alpha$ is given explicitly by formula (1.15) and it is related to the coefficient $\mathbf{r}_\alpha$ in formula (3.8) of Lemma 3.8 by $\mathbf{r}_\alpha = 2\mathbf{k}_\alpha^2$. In addition,

$$k_n - \mathbf{k}_\alpha \le C \mathbf{r}_\alpha^{n-1}, \qquad (5.10)$$

for every $n \in \mathbb{N}$, where $\mathbf{r}_\alpha$ is and $C \in \mathbb{R}_+^*$ depends only on $\alpha$.

(ii) (Coefficients $\kappa_n$) The coefficients $\{\kappa_n\}_{n \in \mathbb{N}}$ are positive numbers and decay geometrically with maximal rate $\mathbf{k}_\alpha$. Specifically, for every $\varepsilon \in \mathbb{R}_+^*$ there exists $C_\varepsilon \in \mathbb{R}_+^*$ such that

$$\kappa_n \le C_\varepsilon (\mathbf{k}_\alpha + \varepsilon)^n, \qquad (5.11)$$

for every $n \in \mathbb{N}$.

*Proof.* On the one hand, the properties of $\{\kappa_n\}_{n \in \mathbb{N}}$ readily follow from those of $\{k_n\}_{n \in \mathbb{N}}$ and the relation (5.8) in Definition 5.2. Specifically, for any given $q \in \mathbb{N}$ and any $n \ge q$ we get the decomposition

$$\kappa_n = \left( \prod_{m=0}^{q} k_m \right) \left( \prod_{m=q+1}^{n} k_m \right) \le k_1^q k_{q+1}^{n-q} = \frac{k_1^q}{k_{q+1}^q} k_{q+1}^n \le \frac{k_1^q}{k_{q+1}^q} (C \mathbf{r}_\alpha^{q+1} + \mathbf{k}_\alpha)^n.$$

Therefore, taking $q$ sufficiently large so that $C \mathbf{r}_\alpha^{q+1} \le \varepsilon$ we conclude (5.11). We then focus on the properties of $\{k_n\}_{n \in \mathbb{N}}$ and we use a similar strategy like in Lemma 3.8.

• STEP 1: Monotone convergence of $\{k_n\}_{n \in \mathbb{N}}$.

Define the following function

$$f(x) := \frac{1}{3 + 2\alpha - 2x}, \quad x \in (-\infty, \frac{3}{2} + \alpha).$$

Then, by (5.7) $\{k_n\}_{n \in \mathbb{N}}$ obeys the following recursive relation

$$k_n = f(k_{n-1}), \quad n > 1.$$

Let us consider the fixed points $x_- < x_+$ of $f$, i.e,

$$x_{\pm} = \frac{(3 + 2\alpha) \pm \sqrt{(1 + 2\alpha)^2 + 8\alpha}}{4}.$$

By inspection, it is clear that $0 < \mathbf{k}_\alpha = x_- < x_+ < \frac{3}{2} + \alpha$ and

$$x_- < f(x) < x, \quad \text{if } x \in (x_-, x_+),$$

where $\mathbf{k}_\alpha$ is given in (1.15). Since $k_1 \in (x_-, x_+)$, then we conclude that $\{k_n\}_{n \in \mathbb{N}}$ is a well defined, positive and monotonically decreasing sequence contained in the compact interval $[x_-, x_+]$. Therefore, it must converge towards some limit $\ell$, which is a solution of $f(x) = x$ i.e., $\ell \in \{x_-, x_+\}$. Since the full sequence $\{k_n\}_{n \in \mathbb{N}}$ is below $x_+$ and decreasing, we then conclude that $\ell = x_- = \mathbf{k}_\alpha$. In particular, we obtain

$$\mathbf{k}_\alpha < k_n \le k_1, \qquad (5.12)$$

for any $n \in \mathbb{N}$.

• STEP 2: Convergence rates.

As for Lemma 3.8, we shall use the special algebraic structure of the function $f$, which is the restriction to $(-\infty, \frac{3}{2} + \alpha)$ of the Möbius function $M: \mathbb{R} \setminus \{\frac{3}{2} + \alpha\} \to \mathbb{R}$ given by

$$M(x) := \frac{1}{3 + 2\alpha - 2x}, \quad x \in \mathbb{R} \setminus \left\{\frac{3}{2} + \alpha\right\}.$$

By definition of $\{k_n\}_{n \in \mathbb{N}}$ in (5.7) we obtain

$$k_n - x_+ = \frac{1}{3 + 2\alpha - 2k_{n-1}} - \frac{1}{3 + 2\alpha - 2x_+} = \frac{2}{(3 + 2\alpha - 2k_{n-1})(3 + 2\alpha - 2x_+)}(k_{n-1} - x_+).$$

Similarly, we obtain

$$k_n - x_- = \frac{1}{3 + 2\alpha - 2k_{n-1}} - \frac{1}{3 + 2\alpha - 2x_-} = \frac{2}{(3 + 2\alpha - 2k_{n-1})(3 + 2\alpha - 2x_-)}(k_{n-1} - x_-).$$

By dividing both expressions and iterating such an identity we get

$$\frac{k_n - x_-}{k_n - x_+} = \left( \frac{3 + 2\alpha - 2x_+}{3 + 2\alpha - 2x_-} \right)^{n-1} \frac{k_1 - x_-}{k_1 - x_+}, \qquad (5.13)$$
Specifically, $\mathbf{k}_\alpha$ is given explicitly by formula (1.15) and it is related to the coefficient $\mathbf{r}_\alpha$ in formula (3.8) of Lemma 3.8 by $\mathbf{r}_\alpha = 2\mathbf{k}_\alpha^2$. In addition,

$$k_n - \mathbf{k}_\alpha \le C \mathbf{r}_\alpha^{n-1}, \qquad (5.10)$$

for every $n \in \mathbb{N}$, where $\mathbf{r}_\alpha$ is and $C \in \mathbb{R}_+^*$ depends only on $\alpha$.

(ii) (Coefficients $\kappa_n$) The coefficients $\{\kappa_n\}_{n \in \mathbb{N}}$ are positive numbers and decay geometrically with maximal rate $\mathbf{k}_\alpha$. Specifically, for every $\varepsilon \in \mathbb{R}_+^*$ there exists $C_\varepsilon \in \mathbb{R}_+^*$ such that

$$\kappa_n \le C_\varepsilon (\mathbf{k}_\alpha + \varepsilon)^n, \qquad (5.11)$$

for every $n \in \mathbb{N}$.

*Proof.* On the one hand, the properties of $\{\kappa_n\}_{n \in \mathbb{N}}$ readily follow from those of $\{k_n\}_{n \in \mathbb{N}}$ and the relation (5.8) in Definition 5.2. Specifically, for any given $q \in \mathbb{N}$ and any $n \ge q$ we get the decomposition

$$\kappa_n = \left( \prod_{m=0}^{q} k_m \right) \left( \prod_{m=q+1}^{n} k_m \right) \le k_1^q k_{q+1}^{n-q} = \frac{k_1^q}{k_{q+1}^q} k_{q+1}^n \le \frac{k_1^q}{k_{q+1}^q} (C \mathbf{r}_\alpha^{q+1} + \mathbf{k}_\alpha)^n.$$

Therefore, taking $q$ sufficiently large so that $C \mathbf{r}_\alpha^{q+1} \le \varepsilon$ we conclude (5.11). We then focus on the properties of $\{k_n\}_{n \in \mathbb{N}}$ and we use a similar strategy like in Lemma 3.8.

• STEP 1: Monotone convergence of $\{k_n\}_{n \in \mathbb{N}}$.

Define the following function

$$f(x) := \frac{1}{3 + 2\alpha - 2x}, \quad x \in (-\infty, \frac{3}{2} + \alpha).$$

Then, by (5.7) $\{k_n\}_{n \in \mathbb{N}}$ obeys the following recursive relation

$$k_n = f(k_{n-1}), \quad n > 1.$$

Let us consider the fixed points $x_- < x_+$ of $f$, i.e,

$$x_{\pm} = \frac{(3 + 2\alpha) \pm \sqrt{(1 + 2\alpha)^2 + 8\alpha}}{4}.$$

By inspection, it is clear that $0 < \mathbf{k}_\alpha = x_- < x_+ < \frac{3}{2} + \alpha$ and

$$x_- < f(x) < x, \quad \text{if } x \in (x_-, x_+),$$

where $\mathbf{k}_\alpha$ is given in (1.15). Since $k_1 \in (x_-, x_+)$, then we conclude that $\{k_n\}_{n \in \mathbb{N}}$ is a well defined, positive and monotonically decreasing sequence contained in the compact interval $[x_-, x_+]$. Therefore, it must converge towards some limit $\ell$, which is a solution of $f(x) = x$ i.e., $\ell \in \{x_-, x_+\}$. Since the full sequence $\{k_n\}_{n \in \mathbb{N}}$ is below $x_+$ and decreasing, we then conclude that $\ell = x_- = \mathbf{k}_\alpha$. In particular, we obtain

$$\mathbf{k}_\alpha < k_n \le k_1, \qquad (5.12)$$

for any $n \in \mathbb{N}$.

• STEP 2: Convergence rates.

As for Lemma 3.8, we shall use the special algebraic structure of the function $f$, which is the restriction to $(-\infty, \frac{3}{2} + \alpha)$ of the Möbius function $M: \mathbb{R} \setminus \{\frac{3}{2} + \alpha\} \to \mathbb{R}$ given by

$$M(x) := \frac{1}{3 + 2\alpha - 2x}, \quad x \in \mathbb{R} \setminus \left\{\frac{3}{2} + \alpha\right\}.$$

By definition of $\{k_n\}_{n \in \mathbb{N}}$ in (5.7) we obtain

$$k_n - x_+ = \frac{1}{3 + 2\alpha - 2k_{n-1}} - \frac{1}{3 + 2\alpha - 2x_+} = \frac{2}{(3 + 2\alpha - 2k_{n-1})(3 + 2\alpha - 2x_+)}(k_{n-1} - x_+).$$

Similarly, we obtain

$$k_n - x_- = \frac{1}{3 + 2\alpha - 2k_{n-1}} - \frac{1}{3 + 2\alpha - 2x_-} = \frac{2}{(3 + 2\alpha - 2k_{n-1})(3 + 2\alpha - 2x_-)}(k_{n-1} - x_-).$$

By dividing both expressions and iterating such an identity we get

$$\frac{k_n - x_-}{k_n - x_+} = \left( \frac{3 + 2\alpha - 2x_+}{3 + 2\alpha - 2x_-} \right)^{n-1} \frac{k_1 - x_-}{k_1 - x_+}, \qquad (5.13)$$
<figure><img src="image_9.png" /><figcaption>FIGURE 7. Unique path joining the leaf $j = 121$ to the root in the perfect binary tree $T^3$. The corresponding lineage map takes the form $\Phi_n^{121}(x; \mathbf{y}_3) = \kappa_3 x + \kappa_2 y_1 + \kappa_1 y_{12} + y_{121}$.</figcaption></figure>

for every $n \in \mathbb{N}$. In fact, the basis can be restated as follows

$$ \frac{3 + 2\alpha - 2x_+}{3 + 2\alpha - 2x_-} = \frac{3 + 2\alpha - \sqrt{(1 + 2\alpha)^2 + 8\alpha}}{3 + 2\alpha + \sqrt{(1 + 2\alpha)^2 + 8\alpha}} = \frac{8}{\left(3 + 2\alpha + \sqrt{(1 + 2\alpha)^2 + 8\alpha}\right)^2} \equiv \mathbf{r}_{\alpha} < 1, $$

where $\mathbf{r}_{\alpha}$ is determined by (3.8) in Lemma 3.8. Therefore, using (5.12) and (5.13) we conclude that

$$ k_n - \mathbf{k}_{\alpha} = \frac{k_1 - x_{-}}{x_{+} - k_1}(x_{+} - k_n) \mathbf{r}_{\alpha}^{n-1} \leq \frac{k_1 - x_{-}}{x_{+} - k_1}(x_{+} - x_{-}) \mathbf{r}_{\alpha}^{n-1}, $$

for any $n \in \mathbb{N}$, thus ending the proof.

<figure><img src="image_10.png" /></figure>

**Definition 5.4 (Quadratic forms).** Consider the sequence $\{k_n\}_{n \in \mathbb{N}}$ in Definition 5.2. Then, we define the quadratic form $Q_n = Q_n(\mathbf{y}_n)$ by

$$ Q_n(\mathbf{y}_n) := \sum_{m=0}^{n-1} \sum_{i \in \mathbb{L}_m^n} \left( \frac{1}{4k_{n-m}} (|y_{i1}|^2 + |y_{i2}|^2) - \frac{1}{8} |y_{i1} - y_{i2}|^2 \right), \quad \mathbf{y}_n \in \mathbb{R}^{2(2^n-1)d}, \quad (5.14) $$

where we are using the tree-indexed notation $\mathbf{y}_n = (y_i)_{i \in \mathcal{T}_*^n} \in \mathbb{R}^{2(2^n-1)d}$ in Remark 3.1.

Again, note that when $n = 2$ the above quadratic form $Q_2(\mathbf{y}_2)$ reduces to the one inside the exponential of (5.6). We now show the following uniform control of the quadratic forms $Q_n$ as $n \to \infty$.

**Lemma 5.5 (Uniform positive definite quadratic forms).** Consider the quadratic form $Q_n = Q_n(\mathbf{y}_n)$ in Definition 5.4 and define the couple of coefficients $\beta_{\min}, \beta_{\max} \in \mathbb{R}_+^*$ by

$$ \beta_{\min} := \frac{1 + 2\alpha}{4}, \quad \beta_{\max} := \frac{1}{4 \mathbf{k}_{\alpha}}, $$

where $\mathbf{k}_{\alpha}$ is given by formula (1.15). Then, we obtain that

$$ \beta_{\min} \| \mathbf{y}_n \|^2 \leq Q_n(\mathbf{y}_n) \leq \beta_{\max} \| \mathbf{y}_n \|^2, $$

where we are using the tree-indexed notation $\mathbf{y}_n = (y_i)_{i \in \mathcal{T}_*^n} \in \mathbb{R}^{2(2^n-1)d}$ in Remark 3.1.

*Proof.* By virtue of Young's inequality, we obtain the following lower and upper bound for $Q_n$

$$ \sum_{m=0}^{n-1} \sum_{i \in \mathbb{L}_m^n} \frac{1}{4} \left( \frac{1}{k_{n-m}} - 1 \right) (|y_{i1}|^2 + |y_{i2}|^2) \leq Q_n(\mathbf{y}_n) \leq \sum_{m=0}^{n-1} \sum_{i \in \mathbb{L}_m^n} \frac{1}{4} \frac{1}{k_{n-m}} (|y_{i1}|^2 + |y_{i2}|^2), $$

for every $\mathbf{y}_n = (y_i)_{i \in \mathcal{T}_*^n} \in \mathbb{R}^{2(2^n-1)d}$. Then, the result follows from Lemma 5.3, which guarantees the uniform control $\mathbf{k}_{\alpha} < k_{m-n} \leq k_1$ for every $m = 0, \dots, n-1$. $\square$

**Definition 5.6 (Lineage maps).** We define the lineage maps $\Phi_n^j = \Phi_n^j(x; \mathbf{y}_n)$ associated to the leaf $j \in \mathbb{L}^n$ (see Figure 7) as follows

$$ \Phi_n^j(x; \mathbf{y}_n) := \kappa_n x + \sum_{m=0}^{n-1} \kappa_m y_{cm}(j), \quad \mathbf{y}_n \in \mathbb{R}^{2(2^n-1)d}, \quad (5.15) $$

where $x \in \mathbb{R}^d$ represents the root value, $\mathbf{y}_n = (y_i)_{i \in \mathcal{T}_*^n} \in \mathbb{R}^{2(2^n-1)d}$ is represented according to the tree-indexed notation in Remark 3.1, $\{\kappa_n\}_{n \in \mathbb{N}}$ is given in Definition 5.2, and $c^m = c \circ \dots \circ c$ is the m times iterated map $c: \mathcal{T}_*^n \to \hat{T}^n$ which, to any vertex $i \in \mathcal{T}_*^n$, it associates its child $c(i) \in \hat{T}^n$ (cf. Section 3.1).
<figure><img src="image_9.png" /><figcaption>FIGURE 7. Unique path joining the leaf $j = 121$ to the root in the perfect binary tree $T^3$. The corresponding lineage map takes the form $\Phi_n^{121}(x; \mathbf{y}_3) = \kappa_3 x + \kappa_2 y_1 + \kappa_1 y_{12} + y_{121}$.</figcaption></figure>

for every $n \in \mathbb{N}$. In fact, the basis can be restated as follows

$$ \frac{3 + 2\alpha - 2x_+}{3 + 2\alpha - 2x_-} = \frac{3 + 2\alpha - \sqrt{(1 + 2\alpha)^2 + 8\alpha}}{3 + 2\alpha + \sqrt{(1 + 2\alpha)^2 + 8\alpha}} = \frac{8}{\left(3 + 2\alpha + \sqrt{(1 + 2\alpha)^2 + 8\alpha}\right)^2} \equiv \mathbf{r}_{\alpha} < 1, $$

where $\mathbf{r}_{\alpha}$ is determined by (3.8) in Lemma 3.8. Therefore, using (5.12) and (5.13) we conclude that

$$ k_n - \mathbf{k}_{\alpha} = \frac{k_1 - x_{-}}{x_{+} - k_1}(x_{+} - k_n) \mathbf{r}_{\alpha}^{n-1} \leq \frac{k_1 - x_{-}}{x_{+} - k_1}(x_{+} - x_{-}) \mathbf{r}_{\alpha}^{n-1}, $$

for any $n \in \mathbb{N}$, thus ending the proof.

<figure><img src="image_10.png" /></figure>

**Definition 5.4 (Quadratic forms).** Consider the sequence $\{k_n\}_{n \in \mathbb{N}}$ in Definition 5.2. Then, we define the quadratic form $Q_n = Q_n(\mathbf{y}_n)$ by

$$ Q_n(\mathbf{y}_n) := \sum_{m=0}^{n-1} \sum_{i \in \mathbb{L}_m^n} \left( \frac{1}{4k_{n-m}} (|y_{i1}|^2 + |y_{i2}|^2) - \frac{1}{8} |y_{i1} - y_{i2}|^2 \right), \quad \mathbf{y}_n \in \mathbb{R}^{2(2^n-1)d}, \quad (5.14) $$

where we are using the tree-indexed notation $\mathbf{y}_n = (y_i)_{i \in \mathbb{T}_*^n} \in \mathbb{R}^{2(2^n-1)d}$ in Remark 3.1.

Again, note that when $n = 2$ the above quadratic form $Q_2(\mathbf{y}_2)$ reduces to the one inside the exponential of (5.6). We now show the following uniform control of the quadratic forms $Q_n$ as $n \to \infty$.

**Lemma 5.5 (Uniform positive definite quadratic forms).** Consider the quadratic form $Q_n = Q_n(\mathbf{y}_n)$ in Definition 5.4 and define the couple of coefficients $\beta_{\min}, \beta_{\max} \in \mathbb{R}_+^*$ by

$$ \beta_{\min} := \frac{1 + 2\alpha}{4}, \quad \beta_{\max} := \frac{1}{4 \mathbf{k}_{\alpha}}, $$

where $\mathbf{k}_{\alpha}$ is given by formula (1.15). Then, we obtain that

$$ \beta_{\min} \| \mathbf{y}_n \|^2 \leq Q_n(\mathbf{y}_n) \leq \beta_{\max} \| \mathbf{y}_n \|^2, $$

where we are using the tree-indexed notation $\mathbf{y}_n = (y_i)_{i \in \mathbb{T}_*^n} \in \mathbb{R}^{2(2^n-1)d}$ in Remark 3.1.

*Proof.* By virtue of Young's inequality, we obtain the following lower and upper bound for $Q_n$

$$ \sum_{m=0}^{n-1} \sum_{i \in \mathbb{L}_m^n} \frac{1}{4} \left( \frac{1}{k_{n-m}} - 1 \right) (|y_{i1}|^2 + |y_{i2}|^2) \leq Q_n(\mathbf{y}_n) \leq \sum_{m=0}^{n-1} \sum_{i \in \mathbb{L}_m^n} \frac{1}{4} \frac{1}{k_{n-m}} (|y_{i1}|^2 + |y_{i2}|^2), $$

for every $\mathbf{y}_n = (y_i)_{i \in \mathbb{T}_*^n} \in \mathbb{R}^{2(2^n-1)d}$. Then, the result follows from Lemma 5.3, which guarantees the uniform control $\mathbf{k}_{\alpha} < k_{m-n} \leq k_1$ for every $m = 0, \dots, n-1$. $\square$

**Definition 5.6 (Lineage maps).** We define the lineage maps $\Phi_n^j = \Phi_n^j(x; \mathbf{y}_n)$ associated to the leaf $j \in \mathbb{L}^n$ (see Figure 7) as follows

$$ \Phi_n^j(x; \mathbf{y}_n) := \kappa_n x + \sum_{m=0}^{n-1} \kappa_m y_{cm}(j), \quad \mathbf{y}_n \in \mathbb{R}^{2(2^n-1)d}, \quad (5.15) $$

where $x \in \mathbb{R}^d$ represents the root value, $\mathbf{y}_n = (y_i)_{i \in \mathbb{T}_*^n} \in \mathbb{R}^{2(2^n-1)d}$ is represented according to the tree-indexed notation in Remark 3.1, $\{\kappa_n\}_{n \in \mathbb{N}}$ is given in Definition 5.2, and $c^m = c \circ \dots \circ c$ is the m times iterated map $c: \mathbb{T}_*^n \to \hat{\mathbb{T}}^n$ which, to any vertex $i \in \mathbb{T}_*^n$, it associates its child $c(i) \in \hat{\mathbb{T}}^n$ (cf. Section 3.1).
We are now ready to state the main result of this section extending formula (5.6) to any $n \in \mathbb{N}$. The starting point is again formula (5.3). In the particular case of Gaussian mixing kernel (1.4) and quadratic selection function (1.5), it takes the explicit form

$$ (e^m F_n)(x) = \frac{1}{(2\pi)^{(2^n-1)d/2}} \frac{1}{(4\pi)^{2^n-1d}} \frac{1}{\prod_{m=0}^{n-1} \|F_m\|_{L^1(\mathbb{R}^d)}^{2^n-1-m}} \\ \times \int_{\mathbb{R}^{2(2^n-1)d}} \exp \left( -\sum_{m=0}^{n-2} \sum_{i \in \mathbb{L}_m^n} \left[ \frac{\alpha}{2} (|x_{i1}|^2 + |x_{i2}|^2) + \frac{1}{2} \left| x_i - \frac{x_{i1} + x_{i2}}{2} \right|^2 \right] \right) \\ \times \exp \left( -\sum_{i \in \mathbb{L}_{n-1}^n} \left[ \frac{1}{2} (\alpha + \frac{1}{2}) (|x_{i1}|^2 + |x_{i2}|^2) + \frac{1}{2} \left| x_i - \frac{x_{i1} + x_{i2}}{2} \right|^2 \right] \right) \\ \times \prod_{j \in \mathbb{L}^n} \bar{F}_0(x_j) \, dx_n \tag{5.16} $$

for any $x \in \mathbb{R}^d$ and $n \in \mathbb{N}$. Above, we have used the tree-indexed notation $\mathbf{x}_n = (x_i)_{i \in \mathbb{T}_*^n} \in \mathbb{R}^{2(2^n-1)d}$ in Remark 3.1, we have set $x_0 := x$ and we have considered the rescaled distribution $\bar{F}_0$ according to Definition 5.1. Again, notice that the exponential terms in the first term of the third line of (5.16) have been introduced to appropriately correct the rescaled distribution $\bar{F}_0$ of $F_0$. As a consequence, the factor $(4\pi)^{-2^{n-1}d}$ arises from the normalization by $\mathbf{F}_{\alpha=0}(x_j)$ whilst the factor $(2\pi)^{-(2^n-1)d/2}$ comes from the repeated products of the Gaussian mixing kernel $G$. The main result then reads as follows

**Proposition 5.7** (Reformulation of the iterations I). Assume that $\alpha \in \mathbb{R}_+$ and set any initial datum $F_0 \in \mathcal{M}_+(\mathbb{R}^d)$. Hence, the solution $\{F_n\}_{n \in \mathbb{N}}$ to the time-discrete problem (1.1) admits the following form

$$ F_n(x) = \frac{e^{-\frac{1+\alpha-k_n}{2}|x|^2}}{(4\pi)^{2^n-1}d(2\pi)^{(2^n-1)d/2}} \frac{1}{\prod_{m=0}^{n-1} \|F_m\|_{L^1(\mathbb{R}^d)}^{2^n-1-m}} \int_{\mathbb{R}^{2(2^n-1)d}} e^{-Q_n(\mathbf{y}_n)} \prod_{j \in \mathbb{L}^n} \bar{F}_0(\Phi_n^j(x; \mathbf{y}_n)) \, d\mathbf{y}_n, \tag{5.17} $$

for every $x \in \mathbb{R}^d$ and $n \in \mathbb{N}$, where we are using the tree-indexed notation $\mathbf{y}_n = (y_i)_{i \in \mathbb{T}_*^n} \in \mathbb{R}^{2(2^n-1)d}$ in Remark 3.1 along with further notation from Definitions 5.1, 5.2, 5.4 and 5.6.

As we anticipated for the particular case $n = 2$ at the beginning of this section, the proof of Proposition 5.7 will be based on the following change of variables

$$ x_{i1} = k_{n-m}x_i + y_{i1}, \quad x_{i2} = k_{n-m}x_i + y_{i2}, $$

for every $i \in \mathbb{L}_m^n$, which we will apply in a backwards way starting at indices $i \in \mathbb{L}_{n-1}^n$ and ending at $i = \emptyset$. Again, the objective of such a change of variables is to appropriately move the reference frame so that the above quadratic forms $\frac{1}{2} |x_i - \frac{x_{i1}+x_{i2}}{2}|^2$ inside the exponential of (5.16) get appropriately centered and no crossed terms remain when the selection parts $\frac{\alpha}{2}(|x_{i1}|^2 + |x_{i2}|^2)$ are taken into account. For simplicity, along the proof we shall restrict to tracking how the various terms inside the exponentials in (5.16) get modified after the change of variables.

### Proof of Proposition 5.7.

* **STEP 1:** Change of variables $x_{i1} \to y_{i1}$ and $x_{i2} \to y_{i2}$ for $i \in \mathbb{L}_{n-1}^n$. We consider the change of variables given by

$$ x_{i1} = k_1 x_i + y_{i1}, \quad x_{i2} = k_1 x_i + y_{i2}, $$
We are now ready to state the main result of this section extending formula (5.6) to any $n \in \mathbb{N}$. The starting point is again formula (5.3). In the particular case of Gaussian mixing kernel (1.4) and quadratic selection function (1.5), it takes the explicit form

$$ (e^m F_n)(x) = \frac{1}{(2\pi)^{(2^n-1)d/2}} \frac{1}{(4\pi)^{2^n-1d}} \frac{1}{\prod_{m=0}^{n-1} \|F_m\|_{L^1(\mathbb{R}^d)}^{2^n-1-m}} \\ \times \int_{\mathbb{R}^{2(2^n-1)d}} \exp \left( -\sum_{m=0}^{n-2} \sum_{i \in \mathbb{L}_m^n} \left[ \frac{\alpha}{2} (|x_{i1}|^2 + |x_{i2}|^2) + \frac{1}{2} \left| x_i - \frac{x_{i1} + x_{i2}}{2} \right|^2 \right] \right) \\ \times \exp \left( -\sum_{i \in \mathbb{L}_{n-1}^n} \left[ \frac{1}{2} (\alpha + \frac{1}{2}) (|x_{i1}|^2 + |x_{i2}|^2) + \frac{1}{2} \left| x_i - \frac{x_{i1} + x_{i2}}{2} \right|^2 \right] \right) \\ \times \prod_{j \in \mathbb{L}^n} \bar{F}_0(x_j) \, dx_n, \qquad (5.16) $$

for any $x \in \mathbb{R}^d$ and $n \in \mathbb{N}$. Above, we have used the tree-indexed notation $\mathbf{x}_n = (x_i)_{i \in \mathbb{T}_*^n} \in \mathbb{R}^{2(2^n-1)d}$ in Remark 3.1, we have set $x_0 := x$ and we have considered the rescaled distribution $\bar{F}_0$ according to Definition 5.1. Again, notice that the exponential terms in the first term of the third line of (5.16) have been introduced to appropriately correct the rescaled distribution $\bar{F}_0$ of $F_0$. As a consequence, the factor $(4\pi)^{-2^{n-1}d}$ arises from the normalization by $\mathbf{F}_{\alpha=0}(x_j)$ whilst the factor $(2\pi)^{-(2^n-1)d/2}$ comes from the repeated products of the Gaussian mixing kernel $G$. The main result then reads as follows

**Proposition 5.7** (Reformulation of the iterations I). Assume that $\alpha \in \mathbb{R}_+$ and set any initial datum $F_0 \in \mathcal{M}_+(\mathbb{R}^d)$. Hence, the solution $\{F_n\}_{n \in \mathbb{N}}$ to the time-discrete problem (1.1) admits the following form

$$ F_n(x) = \frac{e^{-\frac{1+\alpha-k_n}{2}|x|^2}}{(4\pi)^{2^n-1-d}(2\pi)^{(2^n-1)d/2}} \frac{1}{\prod_{m=0}^{n-1} \|F_m\|_{L^1(\mathbb{R}^d)}^{2^n-1-m}} \int_{\mathbb{R}^{2(2^n-1)d}} e^{-Q_n(\mathbf{y}_n)} \prod_{j \in \mathbb{L}^n} \bar{F}_0(\Phi_n^j(x; \mathbf{y}_n)) \, d\mathbf{y}_n, \quad (5.17) $$

for every $x \in \mathbb{R}^d$ and $n \in \mathbb{N}$, where we are using the tree-indexed notation $\mathbf{y}_n = (y_i)_{i \in \mathbb{T}_*^n} \in \mathbb{R}^{2(2^n-1)d}$ in Remark 3.1 along with further notation from Definitions 5.1, 5.2, 5.4 and 5.6.

As we anticipated for the particular case $n = 2$ at the beginning of this section, the proof of Proposition 5.7 will be based on the following change of variables

$$ x_{i1} = k_{n-m}x_i + y_{i1}, \quad x_{i2} = k_{n-m}x_i + y_{i2}, $$

for every $i \in \mathbb{L}_m^n$, which we will apply in a backwards way starting at indices $i \in \mathbb{L}_{n-1}^n$ and ending at $i = \emptyset$. Again, the objective of such a change of variables is to appropriately move the reference frame so that the above quadratic forms $\frac{1}{2} |x_i - \frac{x_{i1}+x_{i2}}{2}|^2$ inside the exponential of (5.16) get appropriately centered and no crossed terms remain when the selection parts $\frac{\alpha}{2}(|x_{i1}|^2 + |x_{i2}|^2)$ are taken into account. For simplicity, along the proof we shall restrict to tracking how the various terms inside the exponentials in (5.16) get modified after the change of variables.

### Proof of Proposition 5.7.

* **STEP 1:** Change of variables $x_{i1} \to y_{i1}$ and $x_{i2} \to y_{i2}$ for $i \in \mathbb{L}_{n-1}^n$. We consider the change of variables given by

$$ x_{i1} = k_1 x_i + y_{i1}, \quad x_{i2} = k_1 x_i + y_{i2}, $$
where $k_1$ is given by Definition 5.2. The collection of all the terms for $i \in \mathbb{L}_{n-1}^n$ in (5.16) then read

$$
\begin{align*}
& \frac{\alpha_0}{2} |x_{i1}|^2 + \frac{\alpha_0}{2} |x_{i2}|^2 + \frac{1}{2} \left| x_i - \frac{x_{i1} + x_{i2}}{2} \right|^2 \\
&= \frac{\alpha_0}{2} |k_1 x_i + y_{i1}|^2 + \frac{\alpha_0}{2} |k_1 x_i + y_{i2}|^2 + \frac{1}{2} \left| (1-k_1)x_i - \frac{y_{i1} + y_{i2}}{2} \right|^2 \\
&= \frac{1}{2} (2\alpha_0 k_1^2 + (1-k_1)^2) |x_i|^2 + \frac{\alpha_0}{2} |y_{i1}|^2 + \frac{\alpha_0}{2} |y_{i2}|^2 + \frac{1}{8} |y_{i1} + y_{i2}|^2 \\
&\quad + \left( \alpha_0 k_1 - \frac{1-k_1}{2} \right) x_i \cdot y_{i1} + \left( \alpha_0 k_1 - \frac{1-k_1}{2} \right) x_i \cdot y_{i2} \\
&= \frac{1}{2} (2\alpha_0 k_1^2 + (1-k_1)^2) |x_i|^2 + \frac{\alpha_0}{2} |y_{i1}|^2 + \frac{\alpha_0}{2} |y_{i2}|^2 + \frac{1}{8} |y_{i1} + y_{i2}|^2 \\
&= \frac{1}{2} (1-k_1) |x_i|^2 + \frac{1}{4k_1} |y_{i1}|^2 + \frac{1}{4k_1} |y_{i2}|^2 - \frac{1}{8} |y_{i1} - y_{i2}|^2,
\end{align*}
$$

where we have defined the coefficient $\alpha_0 \in \mathbb{R}_+^*$ by

$$
\alpha_0 := \alpha + \frac{1}{2},
$$

to recombine terms in the second and third lines of (5.16), and we have used the following relation
between $k_1$ and $\alpha_0$ in order to cancel the crossed terms, i.e.,

$$
k_1 = \frac{1}{1 + 2\alpha_0}.
$$

• STEP 2: Change of variables $x_{i1} \to y_{i1}$ and $x_{i2} \to y_{i2}$ for $i \in \mathbb{L}_{n-2}^n$.
Again, we consider the change of variables given by

$x_{i1} = k_2 x_i + y_{i1}, \quad x_{i2} = k_2 x_i + y_{i2},$

where $k_2$ is given by Definition 5.2. Putting together the terms for $i \in \mathbb{L}_{n-2}^n$ in (5.16) and the above
$x_i$-dependent remainder in the above expression in STEP 1 yields the following term under the above
change of variables

$$
\begin{align*}
& \frac{\alpha_1}{2} |x_{i1}|^2 + \frac{\alpha_1}{2} |x_{i2}|^2 + \frac{1}{2} \left| x_i - \frac{x_{i1} + x_{i2}}{2} \right|^2 \\
&= \frac{\alpha_1}{2} |k_2 x_i + y_{i1}|^2 + \frac{\alpha_1}{2} |k_2 x_i + y_{i2}|^2 + \frac{1}{2} \left| (1-k_2)x_i - \frac{y_{i1} + y_{i2}}{2} \right|^2 \\
&= \frac{1}{2} (2\alpha_1 k_2^2 + (1-k_2)^2) |x_i|^2 + \frac{\alpha_1}{2} |y_{i1}|^2 + \frac{\alpha_1}{2} |y_{i2}|^2 + \frac{1}{8} |y_{i1} + y_{i2}|^2 \\
&\quad + \left( \alpha_1 k_2 - \frac{1-k_2}{2} \right) x_i \cdot y_{i1} + \left( \alpha_1 k_2 - \frac{1-k_2}{2} \right) x_i \cdot y_{i2} \\
&= \frac{1}{2} (2\alpha_1 k_2^2 + (1-k_2)^2) |x_i|^2 + \frac{\alpha_1}{2} |y_{i1}|^2 + \frac{\alpha_1}{2} |y_{i2}|^2 + \frac{1}{8} |y_{i1} + y_{i2}|^2 \\
&= \frac{1}{2} (1-k_2) |x_i|^2 + \frac{1}{4k_2} |y_{i1}|^2 + \frac{1}{4k_2} |y_{i2}|^2 - \frac{1}{8} |y_{i1} - y_{i2}|^2,
\end{align*}
$$

where we have defined again coefficient the $\alpha_1 \in \mathbb{R}_+^*$ by

$$
\alpha_1 := 1 + \alpha - k_1,
$$

in order to absorb the above-mentioned $x_i$-dependent remainder. In addition, note that we have used
again the following relation between $k_2$ and $\alpha_1$ to cancel the crossed term, i.e.,

$$
k_2 = \frac{1}{1 + 2\alpha_1}.
$$

• STEP 3: Change of variables $x_{i1} \to y_{i1}$ and $x_{i2} \to y_{i2}$ for $i \in \mathbb{L}_{n-3}^n$.
We consider the change of variables given by

$x_{i1} = k_3 x_i + y_{i1},$

$x_{i2} = k_3 x_i + y_{i2},$
where $k_1$ is given by Definition 5.2. The collection of all the terms for $i \in \mathbb{L}_{n-1}^n$ in (5.16) then read

$$
\begin{align*}
& \frac{\alpha_0}{2} |x_{i1}|^2 + \frac{\alpha_0}{2} |x_{i2}|^2 + \frac{1}{2} \left| x_i - \frac{x_{i1} + x_{i2}}{2} \right|^2 \\
&= \frac{\alpha_0}{2} |k_1 x_i + y_{i1}|^2 + \frac{\alpha_0}{2} |k_1 x_i + y_{i2}|^2 + \frac{1}{2} \left| (1-k_1)x_i - \frac{y_{i1} + y_{i2}}{2} \right|^2 \\
&= \frac{1}{2} (2\alpha_0 k_1^2 + (1-k_1)^2) |x_i|^2 + \frac{\alpha_0}{2} |y_{i1}|^2 + \frac{\alpha_0}{2} |y_{i2}|^2 + \frac{1}{8} |y_{i1} + y_{i2}|^2 \\
&\quad + \left( \alpha_0 k_1 - \frac{1-k_1}{2} \right) x_i \cdot y_{i1} + \left( \alpha_0 k_1 - \frac{1-k_1}{2} \right) x_i \cdot y_{i2} \\
&= \frac{1}{2} (2\alpha_0 k_1^2 + (1-k_1)^2) |x_i|^2 + \frac{\alpha_0}{2} |y_{i1}|^2 + \frac{\alpha_0}{2} |y_{i2}|^2 + \frac{1}{8} |y_{i1} + y_{i2}|^2 \\
&= \frac{1}{2} (1-k_1) |x_i|^2 + \frac{1}{4k_1} |y_{i1}|^2 + \frac{1}{4k_1} |y_{i2}|^2 - \frac{1}{8} |y_{i1} - y_{i2}|^2,
\end{align*}
$$

where we have defined the coefficient $\alpha_0 \in \mathbb{R}_+^*$ by

$$
\alpha_0 := \alpha + \frac{1}{2},
$$

to recombine terms in the second and third lines of (5.16), and we have used the following relation
between $k_1$ and $\alpha_0$ in order to cancel the crossed terms, i.e.,

$$
k_1 = \frac{1}{1 + 2\alpha_0}.
$$

• STEP 2: Change of variables $x_{i1} \rightarrow y_{i1}$ and $x_{i2} \rightarrow y_{i2}$ for $i \in \mathbb{L}_{n-2}^n$.
Again, we consider the change of variables given by

$x_{i1} = k_2 x_i + y_{i1}, \quad x_{i2} = k_2 x_i + y_{i2},$

where $k_2$ is given by Definition 5.2. Putting together the terms for $i \in \mathbb{L}_{n-2}^n$ in (5.16) and the above
$x_i$-dependent remainder in the above expression in STEP 1 yields the following term under the above
change of variables

$$
\begin{align*}
& \frac{\alpha_1}{2} |x_{i1}|^2 + \frac{\alpha_1}{2} |x_{i2}|^2 + \frac{1}{2} \left| x_i - \frac{x_{i1} + x_{i2}}{2} \right|^2 \\
&= \frac{\alpha_1}{2} |k_2 x_i + y_{i1}|^2 + \frac{\alpha_1}{2} |k_2 x_i + y_{i2}|^2 + \frac{1}{2} \left| (1-k_2)x_i - \frac{y_{i1} + y_{i2}}{2} \right|^2 \\
&= \frac{1}{2} (2\alpha_1 k_2^2 + (1-k_2)^2) |x_i|^2 + \frac{\alpha_1}{2} |y_{i1}|^2 + \frac{\alpha_1}{2} |y_{i2}|^2 + \frac{1}{8} |y_{i1} + y_{i2}|^2 \\
&\quad + \left( \alpha_1 k_2 - \frac{1-k_2}{2} \right) x_i \cdot y_{i1} + \left( \alpha_1 k_2 - \frac{1-k_2}{2} \right) x_i \cdot y_{i2} \\
&= \frac{1}{2} (2\alpha_1 k_2^2 + (1-k_2)^2) |x_i|^2 + \frac{\alpha_1}{2} |y_{i1}|^2 + \frac{\alpha_1}{2} |y_{i2}|^2 + \frac{1}{8} |y_{i1} + y_{i2}|^2 \\
&= \frac{1}{2} (1-k_2) |x_i|^2 + \frac{1}{4k_2} |y_{i1}|^2 + \frac{1}{4k_2} |y_{i2}|^2 - \frac{1}{8} |y_{i1} - y_{i2}|^2,
\end{align*}
$$

where we have defined again coefficient the $\alpha_1 \in \mathbb{R}_+^*$ by

$$
\alpha_1 := 1 + \alpha - k_1,
$$

in order to absorb the above-mentioned $x_i$-dependent remainder. In addition, note that we have used
again the following relation between $k_2$ and $\alpha_1$ to cancel the crossed term, i.e.,

$$
k_2 = \frac{1}{1 + 2\alpha_1}.
$$

• STEP 3: Change of variables $x_{i1} \rightarrow y_{i1}$ and $x_{i2} \rightarrow y_{i2}$ for $i \in \mathbb{L}_{n-3}^n$.
We consider the change of variables given by

$x_{i1} = k_3 x_i + y_{i1},$

$x_{i2} = k_3 x_i + y_{i2},$
where $k_3$ is given by Definition 5.2. Putting together the terms for $i \in \mathbb{L}_{n-3}^n$ in (5.16) and the above $x_i$-dependent remainder in the above expression in STEP 2 yields the following term under the above change of variables

$$
\begin{align*}
& \frac{\alpha_2}{2} |x_{i1}|^2 + \frac{\alpha_2}{2} |x_{i2}|^2 + \frac{1}{2} \left| x_i - \frac{x_{i1} + x_{i2}}{2} \right|^2 \\
&= \frac{\alpha_2}{2} |k_3 x_i + y_{i1}|^2 + \frac{\alpha_2}{2} |k_3 x_i + y_{i2}|^2 + \frac{1}{2} \left| (1-k_3)x_i - \frac{y_{i1} + y_{i2}}{2} \right|^2 \\
&= \frac{1}{2} (2\alpha_2 k_3^2 + (1-k_3)^2) |x_i|^2 + \frac{\alpha_2}{2} |y_{i1}|^2 + \frac{\alpha_2}{2} |y_{i2}|^2 + \frac{1}{8} |y_{i1} + y_{i2}|^2 \\
&\quad + \left( \alpha_2 k_3 - \frac{1-k_3}{2} \right) x_i \cdot y_{i1} + \left( \alpha_2 k_3 - \frac{1-k_3}{2} \right) x_i \cdot y_{i2} \\
&= \frac{1}{2} (2\alpha_2 k_3^2 + (1-k_3)^2) |x_i|^2 + \frac{\alpha_2}{2} |y_{i1}|^2 + \frac{\alpha_2}{2} |y_{i2}|^2 + \frac{1}{8} |y_{i1} + y_{i2}|^2 \\
&= \frac{1}{2} (1-k_3) |x_i|^2 + \frac{1}{4k_3} |y_{i1}|^2 + \frac{1}{4k_3} |y_{i2}|^2 - \frac{1}{8} |y_{i1} - y_{i2}|^2,
\end{align*}
$$

where we have defined again coefficient the $\alpha_2 \in \mathbb{R}_+^*$

$$
\alpha_2 := 1 + \alpha - k_2,
$$

in order to absorb the above-mentioned $x_i$-dependent remainder. In addition, note that we have used
again the following relation between $k_3$ and $\alpha_2$ to cancel the crossed term, i.e.,

$$
k_3 = \frac{1}{1 + 2\alpha_2}.
$$

Following a similar recursive process, we readily identify the quadratic form $Q_n$ in Definition 5.4 in the
exponential of the integrand in the final expression (5.17).

• STEP 4: Identifying the lineage maps $\Phi_n^j$ and the exponential factor.

On the one hand, the hierarchy of changes of variables immediately leads to

$$
x_j = \Phi_n^j(x; \mathbf{y}_n),
$$

for any leaf $j \in \mathbb{L}^n$, thanks to Definition 5.2 for the coefficients $\{\kappa_n\}_{n \in \mathbb{N}}$ and Definition 5.6 for the lineage maps $\Phi_n^j$. On the other hand, since the Jacobian determinant of each change of variables is 1, no further factors appear during the iterative process except for the exponential $x$-dependent remainder $e^{-\frac{1+\alpha-k_n}{2}|x|^2}$ of the last step of the recurrence, which is not absorbed in the quadratic form $Q_n$. $\square$

5.2. **Probabilistic reinterpretation.** For general purposes, in this part we reformulate the result in Proposition 5.7 in appropriate probabilistic terms. More precisely, this reformulation will not only provide shorter formulas, but it will also allow identifying a problematic key point when studying the asymptotics of the high-dimensional integral (5.17) in Section 6, namely, the presence of non-negligible correlations between the various factors indexed with indexes over leaves $j \in \mathbb{L}^n$ of the tree, which do not dissipate even for long time $n \to \infty$. We shall use the following notation.

**Definition 5.8 (High-dimensional normal distributions).** Consider any integer $n \in \mathbb{N}$ and

$$
\Sigma_n^m := \begin{pmatrix} \frac{(2 - k_{n-m})k_{n-m}}{1 - k_{n-m}} I_d & \frac{k_{n-m}^2}{1 - k_{n-m}} I_d \\ -\frac{k_{n-m}^2}{1 - k_{n-m}} I_d & \frac{(2 - k_{n-m})k_{n-m}}{1 - k_{n-m}} I_d \end{pmatrix}, \quad (5.18)
$$

for every $m = 0, \dots, n-1$. Then, we will define the random vector $\mathbf{Y}_n = (Y_i)_{i \in \mathbb{T}_*^n}$ distributed according to the following multivariate normal distribution

$$
G_n(\mathbf{y}_n) := \prod_{m=0}^{n-1} \prod_{i \in \mathbb{L}_m^n} \frac{1}{(2\pi)^d \sqrt{\det(\Sigma_n^m)}} \exp \left( -\frac{1}{2} (y_{i1}^T, y_{i2}^T)(\Sigma_n^m)^{-1} \begin{pmatrix} y_{i1} \\ y_{i2} \end{pmatrix} \right), \quad \mathbf{y}_n \in \mathbb{R}^{2(2^n-1)d}, \quad (5.19)
$$

where again we are using the tree-indexed notation $\mathbf{y}_n = (y_i)_{i \in \mathbb{T}_*^n} \in \mathbb{R}^{2(2^n-1)d}$ in Remark 3.1.

Notice that for every index $i \in \mathbb{L}_m^n$ with $m = 0, \dots, n-1$ the pair $(Y_{i1}, Y_{i2})$ is normally distributed according to $\mathcal{N}(0, \Sigma_n^m)$ and $Y_{i1}$ and $Y_{i2}$ are certainly correlated. Nevertheless, the vector $(Y_{i1}, Y_{i2})$ is independent on any other random vector $Y_j$ for indices $j \neq i1$ and $j \neq i2$. We are now ready to introduce the following probabilistic reformulation of Proposition 5.7.
where $k_3$ is given by Definition 5.2. Putting together the terms for $i \in \mathbb{L}_{n-3}^n$ in (5.16) and the above $x_i$-dependent remainder in the above expression in STEP 2 yields the following term under the above change of variables

$$ 
\begin{aligned}
& \frac{\alpha_2}{2}|x_{i1}|^2 + \frac{\alpha_2}{2}|x_{i2}|^2 + \frac{1}{2}\left|x_i - \frac{x_{i1} + x_{i2}}{2}\right|^2 \\
& = \frac{\alpha_2}{2}|k_3x_i + y_{i1}|^2 + \frac{\alpha_2}{2}|k_3x_i + y_{i2}|^2 + \frac{1}{2}\left|(1-k_3)x_i - \frac{y_{i1} + y_{i2}}{2}\right|^2 \\
& = \frac{1}{2}(2\alpha_2k_3^2 + (1-k_3)^2)|x_i|^2 + \frac{\alpha_2}{2}|y_{i1}|^2 + \frac{\alpha_2}{2}|y_{i2}|^2 + \frac{1}{8}|y_{i1} + y_{i2}|^2 \\
& \quad + \left(\alpha_2k_3 - \frac{1-k_3}{2}\right)x_i \cdot y_{i1} + \left(\alpha_2k_3 - \frac{1-k_3}{2}\right)x_i \cdot y_{i2} \\
& = \frac{1}{2}(2\alpha_2k_3^2 + (1-k_3)^2)|x_i|^2 + \frac{\alpha_2}{2}|y_{i1}|^2 + \frac{\alpha_2}{2}|y_{i2}|^2 + \frac{1}{8}|y_{i1} + y_{i2}|^2 \\
& = \frac{1}{2}(1-k_3)|x_i|^2 + \frac{1}{4k_3}|y_{i1}|^2 + \frac{1}{4k_3}|y_{i2}|^2 - \frac{1}{8}|y_{i1} - y_{i2}|^2,
\end{aligned}
$$

where we have defined again coefficient the $\alpha_2 \in \mathbb{R}_+^*$

$$ \alpha_2 := 1 + \alpha - k_2, $$

in order to absorb the above-mentioned $x_i$-dependent remainder. In addition, note that we have used again the following relation between $k_3$ and $\alpha_2$ to cancel the crossed term, i.e.,

$$ k_3 = \frac{1}{1 + 2\alpha_2}. $$

Following a similar recursive process, we readily identify the quadratic form $Q_n$ in Definition 5.4 in the exponential of the integrand in the final expression (5.17).

• **STEP 4:** Identifying the lineage maps $\Phi_n^j$ and the exponential factor.

On the one hand, the hierarchy of changes of variables immediately leads to

$$ x_j = \Phi_n^j(x; \mathbf{y}_n), $$

for any leaf $j \in \mathbb{L}^n$, thanks to Definition 5.2 for the coefficients $\{\kappa_n\}_{n \in \mathbb{N}}$ and Definition 5.6 for the lineage maps $\Phi_n^j$. On the other hand, since the Jacobian determinant of each change of variables is 1, no further factors appear during the iterative process except for the exponential $x$-dependent remainder $e^{-\frac{1+\alpha-k_n}{2}|x|^2}$ of the last step of the recurrence, which is not absorbed in the quadratic form $Q_n$. $\square$

5.2. **Probabilistic reinterpretation.** For general purposes, in this part we reformulate the result in Proposition 5.7 in appropriate probabilistic terms. More precisely, this reformulation will not only provide shorter formulas, but it will also allow identifying a problematic key point when studying the asymptotics of the high-dimensional integral (5.17) in Section 6, namely, the presence of non-negligible correlations between the various factors indexed with indexes over leaves $j \in \mathbb{L}^n$ of the tree, which do not dissipate even for long time $n \to \infty$. We shall use the following notation.

**Definition 5.8 (High-dimensional normal distributions).** Consider any integer $n \in \mathbb{N}$ and

$$ \Sigma_n^m := \begin{pmatrix} \frac{(2 - k_{n-m})k_{n-m}}{1 - k_{n-m}} I_d & \frac{k_{n-m}^2}{1 - k_{n-m}} I_d \\ -\frac{k_{n-m}^2}{1 - k_{n-m}} I_d & \frac{(2 - k_{n-m})k_{n-m}}{1 - k_{n-m}} I_d \end{pmatrix}, \quad (5.18) $$

for every $m = 0, \dots, n-1$. Then, we will define the random vector $\mathbf{Y}_n = (Y_i)_{i \in \mathbb{T}_*^n}$ distributed according to the following multivariate normal distribution

$$ G_n(\mathbf{y}_n) := \prod_{m=0}^{n-1} \prod_{i \in \mathbb{L}_m^n} \frac{1}{(2\pi)^d \sqrt{\det(\Sigma_n^m)}} \exp \left( -\frac{1}{2} (y_{i1}^T, y_{i2}^T)(\Sigma_n^m)^{-1} \begin{pmatrix} y_{i1} \\ y_{i2} \end{pmatrix} \right), \quad \mathbf{y}_n \in \mathbb{R}^{2(2^n-1)d}, \quad (5.19) $$

where again we are using the tree-indexed notation $\mathbf{y}_n = (y_i)_{i \in \mathbb{T}_*^n} \in \mathbb{R}^{2(2^n-1)d}$ in Remark 3.1.

Notice that for every index $i \in \mathbb{L}_m^n$ with $m = 0, \dots, n-1$ the pair $(Y_{i1}, Y_{i2})$ is normally distributed according to $N(0, \Sigma_n^m)$ and $Y_{i1}$ and $Y_{i2}$ are certainly correlated. Nevertheless, the vector $(Y_{i1}, Y_{i2})$ is independent on any other random vector $Y_j$ for indices $j \neq i1$ and $j \neq i2$. We are now ready to introduce the following probabilistic reformulation of Proposition 5.7.
**Proposition 5.9** (Reformulation of the iterations II). Assume that $\alpha \in \mathbb{R}_+^*$ and set any initial datum $F_0 \in \mathcal{M}_+(\mathbb{R}^d)$. Hence, the solution $\{F_n\}_{n \in \mathbb{N}}$ to the time-discrete problem (1.1) admits the following form

$$F_n(x) = \frac{e^{-\frac{1+\alpha-k_n}{2}|x|^2}}{(2\pi)^{d/2} 2^{2n-1d}} \left( \prod_{m=0}^{n-1} \frac{\left( \frac{4k_{n-m}^2}{1-k_{n-m}} \right)^{2m-1d}}{\|F_m\|_{L^1(\mathbb{R}^d)}^{2n-1-m}} \right) \mathbb{E} \left[ \prod_{j \in \mathbb{L}^n} \bar{F}_0(\Phi_n^j(x; \mathbf{Y}_n)) \right], \quad (5.20)$$

for every $x \in \mathbb{R}^d$ and $n \in \mathbb{N}$, where we are using the high-dimensional random vector $\mathbf{Y}_n = (Y_i)_{i \in \mathbb{T}_*^n}$ in Definition 5.8 along with further notation from Definitions 5.1, 5.2, 5.4 and 5.6.

*Proof.* By Proposition 5.7 we obtain

$$F_n(x) = \frac{e^{-\frac{1+\alpha-k_n}{2}|x|^2}}{(4\pi)^{2n-1d}(2\pi)^{(2n-1)d/2}} \frac{1}{\prod_{m=0}^{n-1} \|F_m\|_{L^1(\mathbb{R}^d)}^{2n-1-m}} \int_{\mathbb{R}^{2(2n-1)d}} e^{-Q_n(\mathbf{y}_n)} \prod_{j \in \mathbb{L}^n} \bar{F}_0(\Phi_n^j(x; \mathbf{y}_n)) \, d\mathbf{y}_n \\ = \frac{(2\pi)^{(2n-1)d/2}}{(4\pi)^{2n-1d}} e^{-\frac{1+\alpha-k_n}{2}|x|^2} \left( \prod_{m=0}^{n-1} \frac{\left( \frac{4k_{n-m}^2}{1-k_{n-m}} \right)^{2m-1d}}{\|F_m\|_{L^1(\mathbb{R}^d)}^{2n-1-m}} \right) \int_{\mathbb{R}^{2(2n-1)d}} G_n(\mathbf{y}_n) \prod_{j \in \mathbb{L}^n} \bar{F}_0(\Phi_n^j(x; \mathbf{y}_n)) \, d\mathbf{y}_n,$$

for every $x \in \mathbb{R}^d$ and each $n \in \mathbb{N}$, where in the last line we have used that the inverse and determinant of the covariance matrices $\Sigma_n^m$ in (5.18) take the form

$$(\Sigma_n^m)^{-1} := \begin{pmatrix} \frac{2 - k_{n-m}}{4k_{n-m}} I_d & \frac{1}{4} I_d \\ \frac{1}{4} I_d & \frac{2 - k_{n-m}}{4k_{n-m}} I_d \end{pmatrix}, \quad \det(\Sigma_n^m) = \left( \frac{4k_{n-m}^2}{1 - k_{n-m}} \right)^d.$$

Then, the result follows by applying the law of the unconscious statistician (LOTUS) to relate the preceding integral to the expectation of the random variable $\prod_{j \in \mathbb{L}^n} \bar{F}_0(\Phi_n^j(x; \mathbf{Y}_n))$. $\square$

As illustrated in the overall map in Figure 2, the reformulation of the iterations in Propositions 5.7 and 5.9 will be fundamental in order to characterize the long-time behavior of the solution $\{F_n\}_{n \in \mathbb{N}}$ to the time-discrete problem (1.1). More specifically, we need to unravel the asymptotic behavior of the high-dimensional integral encoded in the expectation term (5.20). As discussed in Sections 2.3 and 2.4, one is able to control it in the log-Lipschitz norm under stringent constraints on the initial data, which are not fully satisfactory. In view of the structure of the expectation term, alternatively one may expect to be able to interchange expectations with products and control each factor separately. However, as we show below, this naive idea is bound to fail since the involved random variables are correlated, and they do not uncorrelate even in the limit $n \to \infty$. This complicated structure imposes severe problems when handling formula (5.20), and one needs more powerfull ideas, which we introduce later in Lemma 6.3.

**Remark 5.10** (Non-negligible correlations). Since $\mathbf{Y}_n$ in Definition 5.8 is normally distributed and the lineage map $\Phi_n^j(x; \mathbf{Y}_n)$ in Definition 5.6 are affine transformations of $\mathbf{y}_n$, then all $(\Phi_n^j(x; \mathbf{Y}_n))_{j \in \mathbb{L}^n}$ are also normally distributed. Unfortunately, the components of $(\Phi_n^j(x; \mathbf{Y}_n))_{j \in \mathbb{L}^n}$ are correlated since, in particular, all the components $(i1, i2)$ are correlated in view of the structure (5.18) of their covariance matrices. Specifically, a straightforward computation shows that the covariance matrix

$$K_{ij} = \mathbb{E} \left[ (\Phi_n^i(x; \mathbf{Y}_n) - \mathbb{E} \Phi_n^i(x; \mathbf{Y}_n)) \otimes (\Phi_n^j(x; \mathbf{Y}_n) - \mathbb{E} \Phi_n^j(x; \mathbf{Y}_n)) \right],$$

between any couple $(i, j) \in \mathbb{L}^n \times \mathbb{L}^n$ of components (see Table 1 for the notation $\otimes$ of the Kronecker product) takes the following explicit form

$$K_{ij} = \begin{cases} \left( \sum_{q=0}^{n-1} \kappa_q^2 \frac{(2 - k_{q+1})k_{q+1}}{1 - k_{q+1}} \right) I_d, & \text{if } |i \wedge j| = n, \\ \left( -\kappa_{n-i \wedge j}^2 \frac{k_{n-i \wedge j}^2}{1 - k_{n-i \wedge j}} + \sum_{q=n-i \wedge j}^{n-1} \kappa_q^2 \frac{(2 - k_{q+1})k_{q+1}}{1 - k_{q+1}} \right) I_d, & \text{if } 0 < |i \wedge j| < n, \\ \left( -\kappa_{n-1}^2 \frac{k_n^2}{1 - k_n} \right) I_d, & \text{if } |i \wedge j| = 0. \end{cases}$$

Here, $|i \wedge j|$ represents the level of the highest common descendant $i \wedge j$ of $i$ and $j$ (cf. Section 3.1). Note that correlations only disappears for leaves $(i, j)$ such that $n - |i \wedge j| \to \infty$, that is, when the length of the
**Proposition 5.9** (Reformulation of the iterations II). Assume that $\alpha \in \mathbb{R}_+^*$ and set any initial datum $F_0 \in \mathcal{M}_+(\mathbb{R}^d)$. Hence, the solution $\{F_n\}_{n \in \mathbb{N}}$ to the time-discrete problem (1.1) admits the following form

$$F_n(x) = \frac{e^{-\frac{1+\alpha-k_n}{2}|x|^2}}{(2\pi)^{d/2} 2^{2n-1d}} \left( \prod_{m=0}^{n-1} \frac{\left( \frac{4k_{n-m}^2}{1-k_{n-m}} \right)^{2^{m-1}d}}{\|F_m\|_{L^1(\mathbb{R}^d)}^{2^{n-1-m}}} \right) \mathbb{E} \left[ \prod_{j \in \mathbb{L}^n} \bar{F}_0(\Phi_n^j(x; \mathbf{Y}_n)) \right], \quad (5.20)$$

for every $x \in \mathbb{R}^d$ and $n \in \mathbb{N}$, where we are using the high-dimensional random vector $\mathbf{Y}_n = (Y_i)_{i \in \mathbb{T}_*^n}$ in Definition 5.8 along with further notation from Definitions 5.1, 5.2, 5.4 and 5.6.

*Proof.* By Proposition 5.7 we obtain

$$F_n(x) = \frac{e^{-\frac{1+\alpha-k_n}{2}|x|^2}}{(4\pi)^{2n-1d}(2\pi)^{(2n-1)d/2}} \frac{1}{\prod_{m=0}^{n-1} \|F_m\|_{L^1(\mathbb{R}^d)}^{2^{n-1-m}}} \int_{\mathbb{R}^{2(2^n-1)d}} e^{-Q_n(\mathbf{y}_n)} \prod_{j \in \mathbb{L}^n} \bar{F}_0(\Phi_n^j(x; \mathbf{y}_n)) \, d\mathbf{y}_n \\ = \frac{(2\pi)^{(2n-1)d/2}}{(4\pi)^{2n-1d}} e^{-\frac{1+\alpha-k_n}{2}|x|^2} \left( \prod_{m=0}^{n-1} \frac{\left( \frac{4k_{n-m}^2}{1-k_{n-m}} \right)^{2^{m-1}d}}{\|F_m\|_{L^1(\mathbb{R}^d)}^{2^{n-1-m}}} \right) \int_{\mathbb{R}^{2(2^n-1)d}} G_n(\mathbf{y}_n) \prod_{j \in \mathbb{L}^n} \bar{F}_0(\Phi_n^j(x; \mathbf{y}_n)) \, d\mathbf{y}_n,$$

for every $x \in \mathbb{R}^d$ and each $n \in \mathbb{N}$, where in the last line we have used that the inverse and determinant of the covariance matrices $\Sigma_n^m$ in (5.18) take the form

$$(\Sigma_n^m)^{-1} := \begin{pmatrix} \frac{2 - k_{n-m}}{4k_{n-m}} I_d & \frac{1}{4} I_d \\ \frac{1}{4} I_d & \frac{2 - k_{n-m}}{4k_{n-m}} I_d \end{pmatrix}, \quad \det(\Sigma_n^m) = \left( \frac{4k_{n-m}^2}{1 - k_{n-m}} \right)^d.$$

Then, the result follows by applying the law of the unconscious statistician (LOTUS) to relate the preceding integral to the expectation of the random variable $\prod_{j \in \mathbb{L}^n} \bar{F}_0(\Phi_n^j(x; \mathbf{Y}_n))$. $\square$

As illustrated in the overall map in Figure 2, the reformulation of the iterations in Propositions 5.7 and 5.9 will be fundamental in order to characterize the long-time behavior of the solution $\{F_n\}_{n \in \mathbb{N}}$ to the time-discrete problem (1.1). More specifically, we need to unravel the asymptotic behavior of the high-dimensional integral encoded in the expectation term (5.20). As discussed in Sections 2.3 and 2.4, one is able to control it in the log-Lipschitz norm under stringent constraints on the initial data, which are not fully satisfactory. In view of the structure of the expectation term, alternatively one may expect to be able to interchange expectations with products and control each factor separately. However, as we show below, this naive idea is bound to fail since the involved random variables are correlated, and they do not uncorrelate even in the limit $n \to \infty$. This complicated structure imposes severe problems when handling formula (5.20), and one needs more powerfull ideas, which we introduce later in Lemma 6.3.

**Remark 5.10** (Non-negligible correlations). Since $\mathbf{Y}_n$ in Definition 5.8 is normally distributed and the lineage map $\Phi_n^j(x; \mathbf{Y}_n)$ in Definition 5.6 are affine transformations of $\mathbf{y}_n$, then all $(\Phi_n^j(x; \mathbf{Y}_n))_{j \in \mathbb{L}^n}$ are also normally distributed. Unfortunately, the components of $(\Phi_n^j(x; \mathbf{Y}_n))_{j \in \mathbb{L}^n}$ are correlated since, in particular, all the components $(i1, i2)$ are correlated in view of the structure (5.18) of their covariance matrices. Specifically, a straightforward computation shows that the covariance matrix

$$K_{ij} = \mathbb{E} \left[ (\Phi_n^i(x; \mathbf{Y}_n) - \mathbb{E} \Phi_n^i(x; \mathbf{Y}_n)) \otimes (\Phi_n^j(x; \mathbf{Y}_n) - \mathbb{E} \Phi_n^j(x; \mathbf{Y}_n)) \right],$$

between any couple $(i, j) \in \mathbb{L}^n \times \mathbb{L}^n$ of components (see Table 1 for the notation $\otimes$ of the Kronecker product) takes the following explicit form

$$K_{ij} = \begin{cases} \left( \sum_{q=0}^{n-1} \kappa_q^2 \frac{(2-k_{q+1})k_{q+1}}{1-k_{q+1}} \right) I_d, & \text{if } |i \wedge j| = n, \\ \left( -\kappa_{n-i \wedge j-1}^2 \frac{k_{n-i \wedge j}^2}{1-k_{n-i \wedge j}^2} + \sum_{q=n-i \wedge j}^{n-1} \kappa_q^2 \frac{(2-k_{q+1})k_{q+1}}{1-k_{q+1}} \right) I_d, & \text{if } 0 < |i \wedge j| < n, \\ \left( -\kappa_{n-1}^2 \frac{k_n^2}{1-k_n} \right) I_d, & \text{if } |i \wedge j| = 0. \end{cases}$$

Here, $|i \wedge j|$ represents the level of the highest common descendant $i \wedge j$ of $i$ and $j$ (cf. Section 3.1). Note that correlations only disappears for leaves $(i, j)$ such that $n - |i \wedge j| \to \infty$, that is, when the length of the
path from i (or j) to the highest common descendant i ∧ j diverges as the height n of the tree increases. Unfortunately, correlations are non-negligible in all other cases even for n → ∞ .

6. UNIQUENESS OF EQUILIBRIA AND QUANTITATIVE ERGODICITY PROPERTY

In this part we shall prove our main results, namely, Theorem 1.1 and Corollary 1.2. First, we will derive a local convergence result as n → ∞ of the profiles { F$_{n}$ }$_{n}$$_{∈}$ N solving (1.1) with generic initial datum F$_{0}$ ∈ M $_{+}$( R $^{d}$) towards the Gaussian profile F α in (1.12). This weaker local analogue of Theorem 1.1 will be enough to derive the uniqueness of the solution ( λ $_{α}$, F $_{α}$) of the eigenproblem (1.11) in Corollary 1.2. Second, we will extend our local result into our main global result in Theorem 1.1 by quantifying the relaxation of the growth rate ‖ F$_{n}$ ‖$_{L}$ $^{1}$( R $^{d}$) / ‖ F$_{n}$$_{−}$$_{1}$ ‖$_{L}$ $^{1}$( R $^{d}$) and the normalized profiles F$_{n}$ / ‖ F$_{n}$ ‖$_{L}$ $^{1}$( R $^{d}$) towards the unique solution ( λ $_{α}$, F $_{α}$) of the eigenproblem (1.11). Our strategy will exploit the high-dimensional iterative formula (5.20) in Proposition 5.9.

6.1. Local convergence result. In this part we shall prove our local convergence result. To this end, we will show that the dependency of the expectation term in (5.20) on the variable x decays to zero as n → ∞ uniformly over compact sets. To simplify notation, we shall define the following operator.

Definition 6.1 (Expectation operators) . Under the notation in Definitions 5.6 and 5.8, we define

$$
\mathcal{E}_n[F](x) := \mathbb{E} \left[ \prod_{j \in \mathcal{L}^n} F(\Phi_n^j(x; \mathbf{Y}_n)) \right], \quad x \in \mathbb{R}^d, \qquad (6.1)
$$

for any n ∈ N and any non-negative measurable function F : R d → R .

Note that E $_{n}$[ F ]( x ) possibly takes extended values in ( −∞ , + ∞ ] if the expectation does not exist. However, if F is a bounded function (and it will be often the case), we obtain that E $_{n}$[ F ]( x ) always exists and is finite. The following technical lemma for the tail of the incomplete Gamma function will be required later to determine the asymptotic behavior of the above operator E $_{n}$.

Lemma 6.2. Consider the (upper) incomplete Gamma function with parameter a ∈ R ∗ $_{+}$, i.e.,

$$
\Gamma(a, x) := \int_{x}^{+\infty} s^{a-1} e^{-s} ds, \quad x \in \mathbb{R}_{+}^{*}. \qquad (6.2)
$$

Then, the following estimate holds true

$$
\Gamma(a+1, x) \leq 2e^{-x/2}, \tag{6.3}
$$

for every x, a ∈ R ∗ + verifying the constraints

$$
x > e, \quad 2a < \frac{x}{\log(x)}. \tag{6.4}
$$

Proof. Let us restate the incomplete Gamma function as follows

$$
\Gamma(a+1, x) = \int_{x}^{+\infty} \exp(-\psi_{a}(s)s) ds,
$$

where ψ$_{a}$ ( s ) is the function

$$
\psi_a(s) := 1 - a \frac{\log(s)}{s}, \quad s > 0.
$$

Notice that ψ$_{a}$ is monotonically increasing for s ∈ [ e, ∞ ) . Hence, we obtain that

$$
\psi_a(s) \geq \psi_a(x) \geq \frac{1}{2},
$$

for every s ≥ x , where we have used (6.4) on each of the above inequalities. Consequently, we get

$$
\Gamma(a+1, x) \leq \int_{x}^{+\infty} e^{-\frac{s}{2}} ds,
$$

and this ends the proof.

Therefore, we obtain the following weak dependency of E $_{n}$[ F ] on the variable x as n → ∞ .
path from i (or j) to the highest common descendant i ∧ j diverges as the height n of the tree increases. Unfortunately, correlations are non-negligible in all other cases even for $n \to \infty$.

## 6. UNIQUENESS OF EQUILIBRIA AND QUANTITATIVE ERGODICITY PROPERTY

In this part we shall prove our main results, namely, Theorem 1.1 and Corollary 1.2. First, we will derive a local convergence result as $n \to \infty$ of the profiles $\{F_n\}_{n \in \mathbb{N}}$ solving (1.1) with generic initial datum $F_0 \in \mathcal{M}_+(\mathbb{R}^d)$ towards the Gaussian profile $F_\alpha$ in (1.12). This weaker local analogue of Theorem 1.1 will be enough to derive the uniqueness of the solution $(\lambda_\alpha, F_\alpha)$ of the eigenproblem (1.11) in Corollary 1.2. Second, we will extend our local result into our main global result in Theorem 1.1 by quantifying the relaxation of the growth rate $\|F_n\|_{L^1(\mathbb{R}^d)} / \|F_{n-1}\|_{L^1(\mathbb{R}^d)}$ and the normalized profiles $F_n / \|F_n\|_{L^1(\mathbb{R}^d)}$ towards the unique solution $(\lambda_\alpha, F_\alpha)$ of the eigenproblem (1.11). Our strategy will exploit the high-dimensional iterative formula (5.20) in Proposition 5.9.

**6.1. Local convergence result.** In this part we shall prove our local convergence result. To this end, we will show that the dependency of the expectation term in (5.20) on the variable $x$ decays to zero as $n \to \infty$ uniformly over compact sets. To simplify notation, we shall define the following operator.

**Definition 6.1 (Expectation operators).** Under the notation in Definitions 5.6 and 5.8, we define

$$ \mathcal{E}_n[F](x) := \mathbb{E} \left[ \prod_{j \in \mathbb{L}^n} F(\Phi_n^j(x; \mathbf{Y}_n)) \right], \quad x \in \mathbb{R}^d, \qquad (6.1) $$

for any $n \in \mathbb{N}$ and any non-negative measurable function $F: \mathbb{R}^d \to \mathbb{R}$.

Note that $\mathcal{E}_n[F](x)$ possibly takes extended values in $(-\infty, +\infty]$ if the expectation does not exist. However, if $F$ is a bounded function (and it will be often the case), we obtain that $\mathcal{E}_n[F](x)$ always exists and is finite. The following technical lemma for the tail of the incomplete Gamma function will be required later to determine the asymptotic behavior of the above operator $\mathcal{E}_n$.

**Lemma 6.2.** Consider the (upper) incomplete Gamma function with parameter $a \in \mathbb{R}_+^*$, i.e.,

$$ \Gamma(a, x) := \int_x^{+\infty} s^{a-1} e^{-s} ds, \quad x \in \mathbb{R}_+^*. \qquad (6.2) $$

Then, the following estimate holds true

$$ \Gamma(a + 1, x) \le 2e^{-x/2}, \qquad (6.3) $$

for every $x, a \in \mathbb{R}_+^*$ verifying the constraints

$$ x > e, \quad 2a < \frac{x}{\log(x)}. \qquad (6.4) $$

**Proof.** Let us restate the incomplete Gamma function as follows

$$ \Gamma(a + 1, x) = \int_x^{+\infty} \exp(-\psi_a(s)) ds, $$

where $\psi_a(s)$ is the function

$$ \psi_a(s) := 1 - a \frac{\log(s)}{s}, \quad s > 0. $$

Notice that $\psi_a$ is monotonically increasing for $s \in [e, \infty)$. Hence, we obtain that

$$ \psi_a(s) \ge \psi_a(x) \ge \frac{1}{2}, $$

for every $s \ge x$, where we have used (6.4) on each of the above inequalities. Consequently, we get

$$ \Gamma(a + 1, x) \le \int_x^{+\infty} e^{-\frac{s}{2}} ds, $$

and this ends the proof.

Therefore, we obtain the following weak dependency of $\mathcal{E}_n[F]$ on the variable $x$ as $n \to \infty$.
**Lemma 6.3** (Asymptotics of $\mathcal{E}_n$). Assume that $\alpha \in \mathbb{R}_+^*$ and consider any $F \in C^1(\mathbb{R}^d)$ such that

$$ 
\begin{gathered}
0 < F(x) \le C, \\
|\nabla \log F(x)| \le D(1 + |x|),
\end{gathered}
\qquad (6.5) $$

for every $x \in \mathbb{R}^d$ and appropriate $C, D \in \mathbb{R}_+^*$. Then, for every $\varepsilon \in \mathbb{R}_+^*$, there exists a large enough constant $C_\varepsilon \in \mathbb{R}_+^*$ such that the following property holds true

$$ |\nabla \log \mathcal{E}_n[F](x)| \le C_\varepsilon \left( (2\mathbf{k}_\alpha + \varepsilon)^n + (\mathbf{r}_\alpha + \varepsilon)^n |x| + e^{-\frac{(2+\varepsilon)^n}{\varepsilon}} e^{C_\varepsilon(2\mathbf{r}_\alpha + \varepsilon)^n |x|^2} \right), \qquad (6.6) $$

for every $x \in \mathbb{R}^d$ and any $n \in \mathbb{N}$, where $\mathcal{E}_n$ is the operator in Definition 6.1.

*Proof.* Before proving the main result, note that the hypothesis (6.5) guarantee that $F$ is bounded below by a Gaussian profile. Indeed, by the fundamental theorem of calculus we obtain the following relation

$$ \log F(x) = \log F(0) + \int_0^1 \nabla(\log F)(\theta x) \cdot x \, d\theta, $$

for any $x \in \mathbb{R}^d$. Using the linear growth assumption in (6.5) for the log-derivative $\nabla \log F$ we obtain that

$$ \log F(x) \ge \log F(0) - D|x| \int_0^1 (1 + \theta|x|) \, d\theta \ge \log F(0) - \frac{D}{2} - D|x|^2, $$

for any $x \in \mathbb{R}^d$, where in the last inequality we have integrated with respect to $\theta$ and we have used Young's inequality to interpolate $|x|$ by $|x|^2$. Therefore, we obtain the following Gaussian lower bound

$$ F(x) \ge c e^{-\beta|x|^2}, \qquad (6.7) $$

for each $x \in \mathbb{R}^d$ and some $c, \beta \in \mathbb{R}_+^*$ depending only on $F(0)$ and $D$.

For simplicity of notation we define

$$ E_n := \mathcal{E}_n[F], $$

for any $n \in \mathbb{N}$. In the sequel, we shall study the behavior of $\nabla \log E_n$ as $n \to \infty$.

• STEP 1: Control of $\sum_{j \in \mathbb{L}^n} |\Phi_n^j(x; \mathbf{y}_n)|$.

By definition (5.15) of the lineage map $\Phi_n^j$, we obtain that

$$ \sum_{j \in \mathbb{L}^n} |\Phi_n^j(x; \mathbf{y}_n)| = \sum_{j \in \mathbb{L}^n} \left| \kappa_n x + \sum_{m=0}^{n-1} \kappa_m y_{c^m}(j) \right| \le 2^n \kappa_n |x| + \sum_{m=0}^{n-1} 2^m \kappa_m \sum_{i \in \mathbb{L}_{n-m}^n} |y_i|, $$

where we are using the tree-indexed notation $\mathbf{y}_n = (y_i)_{i \in \mathbb{T}_*^n} \in \mathbb{R}^{2(2^n-1)d}$ in Remark 3.1 and $c: \mathbb{T}_*^n \to \hat{\mathbb{T}}^n$ is defined in Section 3.1. Note that in the last inequality we have used that each index of the level $\mathbb{L}_{n-m}^n$ in the tree contributes to the amount of $2^m$ leaves. Then, we obtain

$$ 
\begin{align}
\sum_{j \in \mathbb{L}^n} |\Phi_n^j(x; \mathbf{y}_n)| &\le 2^n \kappa_n |x| + \left( \sum_{m=0}^{n-1} \sum_{i \in \mathbb{L}_{n-m}^n} (2^m \kappa_m)^2 \right)^{1/2} \left( \sum_{m=0}^{n-1} \sum_{i \in \mathbb{L}_{n-m}^n} |y_i|^2 \right)^{1/2} \tag{6.8} \\
&= 2^n \kappa_n |x| + \left( \sum_{m=0}^{n-1} 2^{n-m} (2^m \kappa_m)^2 \right)^{1/2} \| \mathbf{y}_n \| \nonumber \\
&\le 2^n \kappa_n |x| + \sqrt{2} \, 2^{n/2} \| \mathbf{y}_n \|, \nonumber
\end{align}
$$

for any $\mathbf{y}_n \in \mathbb{R}^{2(2^n-1)d}$, where in the first line we have used the Cauchy-Schwarz inequality on the joint sum over indices $m$ and $i$, in the second line we have used the definition $\|\mathbf{y}_n\| = (\sum_{i \in \mathbb{T}_*^n} |y_i|^2)^{1/2}$ of the $\ell_2$ norm (cf. Table 1), and in the last line we have used that

$$ \sum_{m=0}^{n-1} 2^{n-m} (2^m \kappa_m)^2 \le 2^n \sum_{m=0}^{n-1} 2^m k_1^{2m} \le 2^n \sum_{m=0}^{\infty} 2^m k_1^{2m} = \frac{2^n}{1 - 2k_1^2} \le 2^{n+1}. $$

Above, we have used that $k_n \le k_1$ for all $n \in \mathbb{N}$ by virtue of Lemma 5.3(i) and also that

$$ \frac{1}{1 - 2k_1^2} = \frac{1 + 4\alpha + 2\alpha^2}{2(1 + \alpha)^2} = \frac{2(1 + \alpha^2)}{(1 + \alpha)^2 + 2\alpha + \alpha^2} < 2. $$
**Lemma 6.3** (Asymptotics of $\mathcal{E}_n$). Assume that $\alpha \in \mathbb{R}_+^*$ and consider any $F \in C^1(\mathbb{R}^d)$ such that

$$ 
\begin{gathered}
0 < F(x) \le C, \\
|\nabla \log F(x)| \le D(1 + |x|),
\end{gathered}
\qquad (6.5) 
$$

for every $x \in \mathbb{R}^d$ and appropriate $C, D \in \mathbb{R}_+^*$. Then, for every $\varepsilon \in \mathbb{R}_+^*$, there exists a large enough constant $C_\varepsilon \in \mathbb{R}_+^*$ such that the following property holds true

$$ |\nabla \log \mathcal{E}_n[F](x)| \le C_\varepsilon \left( (2\mathbf{k}_\alpha + \varepsilon)^n + (\mathbf{r}_\alpha + \varepsilon)^n |x| + e^{-\frac{(2+\varepsilon)^n}{\varepsilon}} e^{C_\varepsilon(2\mathbf{r}_\alpha + \varepsilon)^n |x|^2} \right), \quad (6.6) $$

for every $x \in \mathbb{R}^d$ and any $n \in \mathbb{N}$, where $\mathcal{E}_n$ is the operator in Definition 6.1.

*Proof.* Before proving the main result, note that the hypothesis (6.5) guarantee that $F$ is bounded below by a Gaussian profile. Indeed, by the fundamental theorem of calculus we obtain the following relation

$$ \log F(x) = \log F(0) + \int_0^1 \nabla(\log F)(\theta x) \cdot x \, d\theta, $$

for any $x \in \mathbb{R}^d$. Using the linear growth assumption in (6.5) for the log-derivative $\nabla \log F$ we obtain that

$$ \log F(x) \ge \log F(0) - D|x| \int_0^1 (1 + \theta|x|) \, d\theta \ge \log F(0) - \frac{D}{2} - D|x|^2, $$

for any $x \in \mathbb{R}^d$, where in the last inequality we have integrated with respect to $\theta$ and we have used Young's inequality to interpolate $|x|$ by $|x|^2$. Therefore, we obtain the following Gaussian lower bound

$$ F(x) \ge c e^{-\beta|x|^2}, \qquad (6.7) $$

for each $x \in \mathbb{R}^d$ and some $c, \beta \in \mathbb{R}_+^*$ depending only on $F(0)$ and $D$.

For simplicity of notation we define

$$ E_n := \mathcal{E}_n[F], $$

for any $n \in \mathbb{N}$. In the sequel, we shall study the behavior of $\nabla \log E_n$ as $n \to \infty$.

• STEP 1: Control of $\sum_{j \in \mathbb{L}^n} |\Phi_n^j(x; \mathbf{y}_n)|$.

By definition (5.15) of the lineage map $\Phi_n^j$, we obtain that

$$ \sum_{j \in \mathbb{L}^n} |\Phi_n^j(x; \mathbf{y}_n)| = \sum_{j \in \mathbb{L}^n} \left| \kappa_n x + \sum_{m=0}^{n-1} \kappa_m y_{c^m}(j) \right| \le 2^n \kappa_n |x| + \sum_{m=0}^{n-1} 2^m \kappa_m \sum_{i \in \mathbb{L}_{n-m}^n} |y_i|, $$

where we are using the tree-indexed notation $\mathbf{y}_n = (y_i)_{i \in \mathbb{T}_*^n} \in \mathbb{R}^{2(2^n-1)d}$ in Remark 3.1 and $c: \mathbb{T}_*^n \to \hat{\mathbb{T}}^n$ is defined in Section 3.1. Note that in the last inequality we have used that each index of the level $\mathbb{L}_{n-m}^n$ in the tree contributes to the amount of $2^m$ leaves. Then, we obtain

$$ 
\begin{align} 
\sum_{j \in \mathbb{L}^n} |\Phi_n^j(x; \mathbf{y}_n)| & \le 2^n \kappa_n |x| + \left( \sum_{m=0}^{n-1} \sum_{i \in \mathbb{L}_{n-m}^n} (2^m \kappa_m)^2 \right)^{1/2} \left( \sum_{m=0}^{n-1} \sum_{i \in \mathbb{L}_{n-m}^n} |y_i|^2 \right)^{1/2} \tag{6.8} \\
& = 2^n \kappa_n |x| + \left( \sum_{m=0}^{n-1} 2^{n-m} (2^m \kappa_m)^2 \right)^{1/2} \| \mathbf{y}_n \| \\
& \le 2^n \kappa_n |x| + \sqrt{2} \, 2^{n/2} \| \mathbf{y}_n \|, \nonumber
\end{align} 
$$

for any $\mathbf{y}_n \in \mathbb{R}^{2(2^n-1)d}$, where in the first line we have used the Cauchy-Schwarz inequality on the joint sum over indices $m$ and $i$, in the second line we have used the definition $\|\mathbf{y}_n\| = (\sum_{i \in \mathbb{T}_*^n} |y_i|^2)^{1/2}$ of the $\ell_2$ norm (cf. Table 1), and in the last line we have used that

$$ \sum_{m=0}^{n-1} 2^{n-m} (2^m \kappa_m)^2 \le 2^n \sum_{m=0}^{n-1} 2^m k_1^{2m} \le 2^n \sum_{m=0}^{\infty} 2^m k_1^{2m} = \frac{2^n}{1 - 2k_1^2} \le 2^{n+1}. $$

Above, we have used that $k_n \le k_1$ for all $n \in \mathbb{N}$ by virtue of Lemma 5.3(i) and also that

$$ \frac{1}{1 - 2k_1^2} = \frac{1 + 4\alpha + 2\alpha^2}{2(1 + \alpha)^2} = \frac{2(1 + \alpha^2)}{(1 + \alpha)^2 + 2\alpha + \alpha^2} < 2. $$
• STEP 2: Derivative of $\log E_n$.

Taking derivatives on (6.1), using the second property in (6.5) along with (6.8) entails

$$ 
\begin{aligned}
|\nabla \log E_n(x)| &\le \kappa_n \frac{\mathbb{E}\left[\left(\sum_{j \in \mathcal{L}^n} |\nabla \log F(\Phi_n^j(x; \mathbf{Y}_n))| \prod_{j \in \mathcal{L}^n} F(\Phi_n^j(x; \mathbf{Y}_n))\right]}{\mathbb{E}\left[\prod_{j \in \mathcal{L}^n} F(\Phi_n^j(x; \mathbf{Y}_n))\right]} \\
&\lesssim \kappa_n \frac{\mathbb{E}\left[\left(\sum_{j \in \mathcal{L}^n} (1 + |\Phi_n^j(x; \mathbf{Y}_n)|) \prod_{j \in \mathcal{L}^n} F(\Phi_n^j(x; \mathbf{Y}_n))\right]}{\mathbb{E}\left[\prod_{j \in \mathcal{L}^n} F(\Phi_n^j(x; \mathbf{Y}_n))\right]} \\
&\lesssim 2^n \kappa_n + 2^n \kappa_n^2 |x| + 2^{n/2} \kappa_n R_n(x),
\end{aligned} 
\quad (6.9) $$

for every $x \in \mathbb{R}^d$, where the remainder $R_n$ reads

$$ R_n(x) := \frac{\mathbb{E}\left[\|\mathbf{Y}_n\| \prod_{j \in \mathcal{L}^n} F(\Phi_n^j(x; \mathbf{Y}_n))\right]}{\mathbb{E}\left[\prod_{j \in \mathcal{L}^n} F(\Phi_n^j(x; \mathbf{Y}_n))\right]}, \quad x \in \mathbb{R}^d. $$

This estimation provides a first insight about the ergodicity property. Indeed the correlations, as measured by $\kappa_n$, decreases faster than $2^{-n}$ (cf. Lemma 5.3), meaning that the contribution $2^n \kappa_n$ decays fast to zero. It is also the case of the other contributions, locally in $x$, but the last one, namely $2^{n/2} \kappa_n R_n(x)$, requires more work. Consider any $\sqrt{2} < q < \frac{1}{\sqrt{2} \kappa_\alpha}$, which exists because $2 \mathbf{k}_\alpha < 1$ for $\alpha \in \mathbb{R}_+^*$, and split $R_n(x) = R_n^1(x) + R_n^2(x)$ in terms of the functions

$$ R_n^1(x) := \frac{1}{Z_n(x)} \int_{\|\mathbf{y}_n\| \le q^n} \|\mathbf{y}_n\| \prod_{j \in \mathcal{L}^n} F(\Phi_n^j(x; \mathbf{y}_n)) e^{-Q_n(\mathbf{y}_n)} d\mathbf{y}_n, \quad (6.10) $$

$$ R_n^2(x) := \frac{1}{Z_n(x)} \int_{\|\mathbf{y}_n\| > q^n} \|\mathbf{y}_n\| \prod_{j \in \mathcal{L}^n} F(\Phi_n^j(x; \mathbf{y}_n)) e^{-Q_n(\mathbf{y}_n)} d\mathbf{y}_n, \quad (6.11) $$

and the normalization factor

$$ Z_n(x) := \int_{\mathbb{R}^{2(2^{n-1})d}} \prod_{j \in \mathcal{L}^n} F(\Phi_n^j(x; \mathbf{y}_n)) e^{-Q_n(\mathbf{y}_n)} d\mathbf{y}_n. \quad (6.12) $$

• STEP 3: Lower bound for $Z_n$.

Using the above lower bound by a Gaussian in (6.7) along with the upper bound for $Q_n$ in Lemma 5.5 implies the following estimate of $Z_n$ in (6.12):

$$ 
\begin{aligned}
Z_n(x) &\ge c^{2^n} \int_{\mathbb{R}^{2(2^{n-1})d}} \exp \left( -\beta \sum_{j \in \mathcal{L}^n} |\Phi_n^j(x; \mathbf{y}_n)|^2 \right) \exp(-\beta_{\max} \|\mathbf{y}_n\|^2) d\mathbf{y}_n \\
&\ge c^{2^n} e^{-2^{2n+1} \kappa_n^2 \beta \|x\|^2} \int_{\mathbb{R}^{2(2^{n-1})d}} \exp \left( -(2^{n+2} \beta + \beta_{\max}) \|\mathbf{y}_n\|^2 \right) d\mathbf{y}_n \\
&= c^{2^n} e^{-2^{2n+1} \kappa_n^2 \beta |x|^2} \left| \mathbb{S}^{2(2^{n-1})d-1} \right| \int_0^{+\infty} r^{2(2^{n-1})d-1} \exp \left( -(2^{n+2} \beta + \beta_{\max}) r^2 \right) dr \\
&= \frac{c^{2^n} e^{-2^{2n+1} \kappa_n^2 \beta |x|^2}}{2(2^{n+2} \beta + \beta_{\max})^{(2^n-1)d}} \left| \mathbb{S}^{2(2^{n-1})d-1} \right| \int_0^{+\infty} s^{(2^n-1)d-1} e^{-s} ds \\
&= \frac{c^{2^n} e^{-2^{2n+1} \kappa_n^2 \beta |x|^2}}{2(2^{n+2} \beta + \beta_{\max})^{(2^n-1)d}} \left| \mathbb{S}^{2(2^{n-1})d-1} \right| \Gamma((2^n-1)d),
\end{aligned}
 $$

where in the last step we have used the Gamma function, that is defined by

$$ \Gamma(a) := \int_{0}^{+\infty} s^{a-1} e^{-s} ds, \quad a > 0. $$

Recall that area of the hypersphere is $|\mathbb{S}^{d-1}| = \frac{2\pi^{d/2}}{\Gamma(\frac{d}{2})}$. Then, we obtain

$$ Z_n(x) \ge \frac{c^{2^n} \pi^{(2^n-1)d}}{(2^{n+2}\beta + \beta_{\max})^{(2^n-1)d}} e^{-2^{2n+1}\kappa_n^2 \beta |x|^2}, \quad (6.13) $$

for every $x \in \mathbb{R}^d$.
• STEP 2: Derivative of $\log E_n$.

Taking derivatives on (6.1), using the second property in (6.5) along with (6.8) entails

$$ 
\begin{aligned}
|\nabla \log E_n(x)| &\le \kappa_n \frac{\mathbb{E}\left[\left(\sum_{j \in \mathcal{L}^n} |\nabla \log F(\Phi_n^j(x; \mathbf{Y}_n))| \prod_{j \in \mathcal{L}^n} F(\Phi_n^j(x; \mathbf{Y}_n))\right]}{\mathbb{E}\left[\prod_{j \in \mathcal{L}^n} F(\Phi_n^j(x; \mathbf{Y}_n))\right]} \\
&\lesssim \kappa_n \frac{\mathbb{E}\left[\left(\sum_{j \in \mathcal{L}^n} (1 + |\Phi_n^j(x; \mathbf{Y}_n)|) \prod_{j \in \mathcal{L}^n} F(\Phi_n^j(x; \mathbf{Y}_n))\right]}{\mathbb{E}\left[\prod_{j \in \mathcal{L}^n} F(\Phi_n^j(x; \mathbf{Y}_n))\right]} \\
&\lesssim 2^n \kappa_n + 2^n \kappa_n^2 |x| + 2^{n/2} \kappa_n R_n(x),
\end{aligned}
\quad (6.9) $$

for every $x \in \mathbb{R}^d$, where the remainder $R_n$ reads

$$ R_n(x) := \frac{\mathbb{E}\left[\|\mathbf{Y}_n\| \prod_{j \in \mathcal{L}^n} F(\Phi_n^j(x; \mathbf{Y}_n))\right]}{\mathbb{E}\left[\prod_{j \in \mathcal{L}^n} F(\Phi_n^j(x; \mathbf{Y}_n))\right]}, \quad x \in \mathbb{R}^d. $$

This estimation provides a first insight about the ergodicity property. Indeed the correlations, as measured by $\kappa_n$, decreases faster than $2^{-n}$ (cf. Lemma 5.3), meaning that the contribution $2^n \kappa_n$ decays fast to zero. It is also the case of the other contributions, locally in $x$, but the last one, namely $2^{n/2} \kappa_n R_n(x)$, requires more work. Consider any $\sqrt{2} < q < \frac{1}{\sqrt{2} \kappa_\alpha}$, which exists because $2 k_\alpha < 1$ for $\alpha \in \mathbb{R}_+^*$, and split $R_n(x) = R_n^1(x) + R_n^2(x)$ in terms of the functions

$$ R_n^1(x) := \frac{1}{Z_n(x)} \int_{\|\mathbf{y}_n\| \le q^n} \|\mathbf{y}_n\| \prod_{j \in \mathcal{L}^n} F(\Phi_n^j(x; \mathbf{y}_n)) e^{-Q_n(\mathbf{y}_n)} d\mathbf{y}_n, \quad (6.10) $$

$$ R_n^2(x) := \frac{1}{Z_n(x)} \int_{\|\mathbf{y}_n\| > q^n} \|\mathbf{y}_n\| \prod_{j \in \mathcal{L}^n} F(\Phi_n^j(x; \mathbf{y}_n)) e^{-Q_n(\mathbf{y}_n)} d\mathbf{y}_n, \quad (6.11) $$

and the normalization factor

$$ Z_n(x) := \int_{\mathbb{R}^{2(2^{n-1})d}} \prod_{j \in \mathcal{L}^n} F(\Phi_n^j(x; \mathbf{y}_n)) e^{-Q_n(\mathbf{y}_n)} d\mathbf{y}_n. \quad (6.12) $$

• STEP 3: Lower bound for $Z_n$.

Using the above lower bound by a Gaussian in (6.7) along with the upper bound for $Q_n$ in Lemma 5.5 implies the following estimate of $Z_n$ in (6.12):

$$ 
\begin{aligned}
Z_n(x) &\ge c^{2^n} \int_{\mathbb{R}^{2(2^{n-1})d}} \exp \left( -\beta \sum_{j \in \mathcal{L}^n} |\Phi_n^j(x; \mathbf{y}_n)|^2 \right) \exp(-\beta_{\max} \|\mathbf{y}_n\|^2) d\mathbf{y}_n \\
&\ge c^{2^n} e^{-2^{2n+1} \kappa_n^2 \beta \|x\|^2} \int_{\mathbb{R}^{2(2^{n-1})d}} \exp \left( -(2^{n+2} \beta + \beta_{\max}) \|\mathbf{y}_n\|^2 \right) d\mathbf{y}_n \\
&= c^{2^n} e^{-2^{2n+1} \kappa_n^2 \beta |x|^2} \left| \mathbb{S}^{2(2^{n-1})d-1} \right| \int_0^{+\infty} r^{2(2^{n-1})d-1} \exp \left( -(2^{n+2} \beta + \beta_{\max}) r^2 \right) dr \\
&= \frac{c^{2^n} e^{-2^{2n+1} \kappa_n^2 \beta |x|^2}}{2(2^{n+2} \beta + \beta_{\max})^{(2^n-1)d}} \left| \mathbb{S}^{2(2^{n-1})d-1} \right| \int_0^{+\infty} s^{(2^n-1)d-1} e^{-s} ds \\
&= \frac{c^{2^n} e^{-2^{2n+1} \kappa_n^2 \beta |x|^2}}{2(2^{n+2} \beta + \beta_{\max})^{(2^n-1)d}} \left| \mathbb{S}^{2(2^{n-1})d-1} \right| \Gamma((2^n-1)d),
\end{aligned}
$$

where in the last step we have used the Gamma function, that is defined by

$$ \Gamma(a) := \int_{0}^{+\infty} s^{a-1} e^{-s} ds, \quad a > 0. $$

Recall that area of the hypersphere is $|\mathbb{S}^{d-1}| = \frac{2\pi^{d/2}}{\Gamma(\frac{d}{2})}$. Then, we obtain

$$ Z_n(x) \ge \frac{c^{2^n} \pi^{(2^n-1)d}}{(2^{n+2}\beta + \beta_{\max})^{(2^n-1)d}} e^{-2^{2n+1}\kappa_n^2 \beta |x|^2}, \quad (6.13) $$

for every $x \in \mathbb{R}^d$.
• STEP 4: Upper bound for $R_n^1$ and $R_n^2$.
On the one hand, notice that

$$R_n^1(x) \le q^n, \quad (6.14)$$

for all $x \in \mathbb{R}^d$ by definition of $R_n^1$ in (6.10). On the other hand, using the lower bound of $Q_n$ in Lemma 5.5 we can estimate $R_n^2$ in (6.11) by

$$\begin{aligned}
R_n^2(x) &\le \frac{C^{2n}}{q^n Z_n(x)} \int_{\|\mathbf{y}_n\| \ge q^n} \|\mathbf{y}_n\|^2 e^{-\beta_{\min} \|\mathbf{y}_n\|} d\mathbf{y}_n \\
&= \frac{C^{2n}}{q^n Z_n(x)} \left|\mathbb{S}^{2(2^n-1)d-1}\right| \int_{q^n}^{+\infty} r^{2(2^n-1)d+1} e^{-\beta_{\min} r^2} dr \\
&= \frac{C^{2n}}{q^n Z_n(x)} \left|\mathbb{S}^{2(2^n-1)d-1}\right| \frac{1}{2\beta_{\min}^{(2^n-1)d+1}} \int_{\beta_{\min} q^{2n}}^{+\infty} s^{(2^n-1)d} e^{-s} ds \\
&= \frac{C^{2n} \pi^{(2^n-1)d}}{q^n \beta_{\min}^{(2^n-1)d+1} Z_n(x)} \frac{\Gamma((2^n-1)d+1, \beta_{\min} q^{2n})}{\Gamma((2^n-1)d)},
\end{aligned} \quad (6.15)$$

for all $x \in \mathbb{R}^d$, where we have used again the incomplete Gamma function (6.2) in Lemma 6.2. Putting the lower estimate (6.13) for $Z_n$ into (6.15) yields

$$R_n^2(x) \lesssim e^{2^{2n+1} \kappa_n^2 \beta |x|^2} \frac{C^{2n} (2^{n+2} \beta + \beta_{\max})^{(2^n-1)d}}{q^n c^{2n} \beta_{\min}^{(2^n-1)d}} \frac{\Gamma((2^n-1)d+1, \beta_{\min} q^{2n})}{\Gamma((2^n-1)d)}, \quad (6.16)$$

for all $x \in \mathbb{R}^d$. Now, take $a = (2^n - 1)d$ and $x = \beta_{\min} q^{2n}$ and notice that the constraints (6.4) in Lemma 6.2 hold true for large enough $n \in \mathbb{N}$ as long as we choose $q > \sqrt{2}$. Hence, we obtain the estimate

$$\Gamma((2^n - 1)d + 1, \beta_{\min} q^{2n}) \le 2e^{-\frac{\beta_{\min}}{2} q^{2n}},$$

for all $n \ge n_* = n_*(\beta_{\min}, d, q)$. Using Stirling's formula for the remaining Gamma function in (6.16), we conclude the estimate

$$R_n^2(x) \lesssim e^{2^{2n+1} \kappa_n^2 \beta |x|^2} \frac{C^{2n} (2^{n+2} \beta + \beta_{\max})^{(2^n-1)d}}{q^n c^{2n} \beta_{\min}^{(2^n-1)d}} \frac{e^{-\frac{\beta_{\min}}{2} q^{2n}}}{(2^n - 1)^{1/2} \left(\frac{(2^n-1)d}{e}\right)^{(2^n-1)d}}, \quad (6.17)$$

for every $x \in \mathbb{R}^d$ and each $n \ge n_*$.

• STEP 5: Final conclusion.

Let us put the preceding estimates (6.14) and (6.17) into (6.9) to achieve

$$|\nabla \log E_n(x)| \lesssim 2^n \kappa_n + 2^{n/2} q^n \kappa_n + 2^n \kappa_n^2 |x| + \frac{M^{2n}}{e^{\frac{\beta_{\min}}{2} q^{2n}}} e^{2^{2n+1} \kappa_n^2 \beta |x|^2}, \quad (6.18)$$

for every $x \in \mathbb{R}^d$ and each $n \ge n_*$, where $M = M(C, c, \beta, \beta_{\min}, \beta_{\max}, d, q)$ is a universal constant that we have introduced to absorb all the double exponential factors of the form $a^{2n}$ with $a > 0$ in the fourth term along with any further lower order factor. On the one hand, choosing $q > \sqrt{2}$ sufficiently close to $\sqrt{2}$ and using (5.11) in Lemma 5.3 we obtain the following asymptotics for the factors in the first three terms of (6.18):

$$\begin{aligned}
2^n \kappa_n &\le C_\varepsilon (2 \mathbf{k}_\alpha + \varepsilon)^n, \\
2^{n/2} q^n \kappa_n &\le C_\varepsilon (2 \mathbf{k}_\alpha + \varepsilon)^n, \\
2^n \kappa_n^2 &\le C_\varepsilon (2 (\mathbf{k}_\alpha + \varepsilon)^2)^n \le C_\varepsilon (\mathbf{r}_\alpha + \varepsilon)^n,
\end{aligned}$$

for arbitrarily small $\varepsilon \in \mathbb{R}_+^*$ and sufficiently large $C_\varepsilon \in \mathbb{R}_+^*$. Note that in the last identity we have used the relation $\mathbf{r}_\alpha = 2 \mathbf{k}_\alpha^2$. On the other hand, since $q$ has been taken larger than $\sqrt{2}$ and arbitrarily close to $\sqrt{2}$, then the double exponential in the denominator of the last term in (6.18) kills the double exponential on the numerator independently on the explicit value of $M$. Therefore, we obtain

$$\frac{M^{2n}}{e^{\frac{\beta_{\min}}{2} q^{2n}}} \le C_\varepsilon e^{-\frac{(2+\varepsilon)^n}{\varepsilon}},$$

by appropriately increasing $C_\varepsilon$ if necessary. Similarly, we infer the following control for the remaining $x$-dependent exponential factor in the fourth term of (6.18)

$$e^{2^{2n+1} \kappa_n^2 \beta |x|^2} \le e^{C_\varepsilon (4 (\mathbf{k}_\alpha + \varepsilon)^2)^n |x|^2} \le e^{C_\varepsilon (2 \mathbf{r}_\alpha + \varepsilon)^n |x|^2},$$
• STEP 4: Upper bound for $R_n^1$ and $R_n^2$.
On the one hand, notice that

$$R_n^1(x) \le q^n, \quad (6.14)$$

for all $x \in \mathbb{R}^d$ by definition of $R_n^1$ in (6.10). On the other hand, using the lower bound of $Q_n$ in Lemma 5.5 we can estimate $R_n^2$ in (6.11) by

$$\begin{align*}
R_n^2(x) &\le \frac{C^{2n}}{q^n Z_n(x)} \int_{\|\mathbf{y}_n\| \ge q^n} \|\mathbf{y}_n\|^2 e^{-\beta_{\min} \|\mathbf{y}_n\|} d\mathbf{y}_n \\
&= \frac{C^{2n}}{q^n Z_n(x)} \left|\mathbb{S}^{2(2^n-1)d-1}\right| \int_{q^n}^{+\infty} r^{2(2^n-1)d+1} e^{-\beta_{\min} r^2} dr \\
&= \frac{C^{2n}}{q^n Z_n(x)} \left|\mathbb{S}^{2(2^n-1)d-1}\right| \frac{1}{2\beta_{\min}^{(2^n-1)d+1}} \int_{\beta_{\min} q^{2n}}^{+\infty} s^{(2^n-1)d} e^{-s} ds \\
&= \frac{C^{2n} \pi^{(2^n-1)d}}{q^n \beta_{\min}^{(2^n-1)d+1} Z_n(x)} \frac{\Gamma((2^n-1)d+1, \beta_{\min} q^{2n})}{\Gamma((2^n-1)d)}, \quad (6.15)
\end{align*}$$

for all $x \in \mathbb{R}^d$, where we have used again the incomplete Gamma function (6.2) in Lemma 6.2. Putting the lower estimate (6.13) for $Z_n$ into (6.15) yields

$$R_n^2(x) \lesssim e^{2^{2n+1} \kappa_n^2 \beta |x|^2} \frac{C^{2n} (2^{n+2} \beta + \beta_{\max})^{(2^n-1)d}}{q^n c^{2n} \beta_{\min}^{(2^n-1)d}} \frac{\Gamma((2^n-1)d+1, \beta_{\min} q^{2n})}{\Gamma((2^n-1)d)}, \quad (6.16)$$

for all $x \in \mathbb{R}^d$. Now, take $a = (2^n - 1)d$ and $x = \beta_{\min} q^{2n}$ and notice that the constraints (6.4) in Lemma 6.2 hold true for large enough $n \in \mathbb{N}$ as long as we choose $q > \sqrt{2}$. Hence, we obtain the estimate

$$\Gamma((2^n - 1)d + 1, \beta_{\min} q^{2n}) \le 2e^{-\frac{\beta_{\min}}{2} q^{2n}},$$

for all $n \ge n_* = n_*(\beta_{\min}, d, q)$. Using Stirling's formula for the remaining Gamma function in (6.16), we conclude the estimate

$$R_n^2(x) \lesssim e^{2^{2n+1} \kappa_n^2 \beta |x|^2} \frac{C^{2n} (2^{n+2} \beta + \beta_{\max})^{(2^n-1)d}}{q^n c^{2n} \beta_{\min}^{(2^n-1)d}} \frac{e^{-\frac{\beta_{\min}}{2} q^{2n}}}{(2^n - 1)^{1/2} \left(\frac{(2^n-1)d}{e}\right)^{(2^n-1)d}}, \quad (6.17)$$

for every $x \in \mathbb{R}^d$ and each $n \ge n_*$.

• STEP 5: Final conclusion.

Let us put the preceding estimates (6.14) and (6.17) into (6.9) to achieve

$$|\nabla \log E_n(x)| \lesssim 2^n \kappa_n + 2^{n/2} q^n \kappa_n + 2^n \kappa_n^2 |x| + \frac{M^{2n}}{e^{\frac{\beta_{\min}}{2} q^{2n}}} e^{2^{2n+1} \kappa_n^2 \beta |x|^2}, \quad (6.18)$$

for every $x \in \mathbb{R}^d$ and each $n \ge n_*$, where $M = M(C, c, \beta, \beta_{\min}, \beta_{\max}, d, q)$ is a universal constant that we have introduced to absorb all the double exponential factors of the form $a^{2n}$ with $a > 0$ in the fourth term along with any further lower order factor. On the one hand, choosing $q > \sqrt{2}$ sufficiently close to $\sqrt{2}$ and using (5.11) in Lemma 5.3 we obtain the following asymptotics for the factors in the first three terms of (6.18):

$$\begin{align*}
2^n \kappa_n &\le C_\varepsilon (2 \mathbf{k}_\alpha + \varepsilon)^n, \\
2^{n/2} q^n \kappa_n &\le C_\varepsilon (2 \mathbf{k}_\alpha + \varepsilon)^n, \\
2^n \kappa_n^2 &\le C_\varepsilon (2 (\mathbf{k}_\alpha + \varepsilon)^2)^n \le C_\varepsilon (\mathbf{r}_\alpha + \varepsilon)^n,
\end{align*}$$

for arbitrarily small $\varepsilon \in \mathbb{R}_+$ and sufficiently large $C_\varepsilon \in \mathbb{R}_+^*$. Note that in the last identity we have used the relation $\mathbf{r}_\alpha = 2 \mathbf{k}_\alpha^2$. On the other hand, since $q$ has been taken larger than $\sqrt{2}$ and arbitrarily close to $\sqrt{2}$, then the double exponential in the denominator of the last term in (6.18) kills the double exponential on the numerator independently on the explicit value of $M$. Therefore, we obtain

$$\frac{M^{2n}}{e^{\frac{\beta_{\min}}{2} q^{2n}}} \le C_\varepsilon e^{-\frac{(2+\varepsilon)^n}{\varepsilon}},$$

by appropriately increasing $C_\varepsilon$ if necessary. Similarly, we infer the following control for the remaining $x$-dependent exponential factor in the fourth term of (6.18)

$$e^{2^{2n+1} \kappa_n^2 \beta |x|^2} \le e^{C_\varepsilon (4 (\mathbf{k}_\alpha + \varepsilon)^2)^n |x|^2} \le e^{C_\varepsilon (2 \mathbf{r}_\alpha + \varepsilon)^n |x|^2},$$
which yields the announced estimate (6.6).

**Remark 6.4.** Under the assumptions in Lemma 6.3, notice that the above result implies in particular that for every $\epsilon \in \mathbb{R}_+^*$ and $R \in \mathbb{R}_+^*$, there exists a large enough constant $C_{\epsilon,R} \in \mathbb{R}_+^*$ such that

$$ \sup_{|x| \le R} |\nabla \log \mathcal{E}_n[F](x)| \le C_{\epsilon,R} (2 \mathbf{k}_\alpha)^n, $$

for any $n \in \mathbb{N}$. In particular, $\nabla \log \mathcal{E}_n[F]$ converge to zero as $n \to \infty$ uniformly on compact sets.

The above observation is the cornerstone to prove the following local version of our main Theorem 1.1.

**Corollary 6.5 (Local convergence result).** Assume that $\alpha \in \mathbb{R}_+^*$, set any initial datum $F_0 \in \mathcal{M}_+(\mathbb{R}^d)$ and consider the solution $\{F_n\}_{n \in \mathbb{N}}$ to the time-discrete problem (1.1). Then, for every $\epsilon \in \mathbb{R}_+^*$ and $R \in \mathbb{R}_+^*$, there exists a large enough constant $C_{\epsilon,R} \in \mathbb{R}_+^*$ such that

$$ \sup_{|x| \le R} \left| \frac{F_n(x)}{\|F_n\|_{L^1(\mathbb{R}^d)}} - \mathbf{F}_\alpha(x) \right| \le C_{\epsilon,R} (2 \mathbf{k}_\alpha)^n, \quad (6.19) $$

for any $n \in \mathbb{N}$, where $\mathbf{F}_\alpha$ is the Gaussian profile (1.12). In particular, $F_n/\|F_n\|_{L^1(\mathbb{R}^d)}$ converge to $\mathbf{F}_\alpha$ as $n \to \infty$ uniformly on compact sets.

In the proof of Corollary 6.5 we will require a careful control of the asymptotic behavior of the ratios $F_n(0)/\|F_n\|_{L^1(\mathbb{R}^d)}$ as $n \to \infty$, which we provide below. First, note that if $\{F_n\}_{n \in \mathbb{N}}$ is the solution (3.4) of the time-discrete problem (1.1) issued at a Gaussian initial datum $F_0$ like in Proposition 3.7, we have

$$ \frac{F_n(0)}{\|F_n\|_{L^1(\mathbb{R}^d)}} = G_{\mu_n, \sigma_n^2}(0) \to \mathbf{F}_\alpha(0) = \frac{1}{(2\pi \sigma_\alpha^2)^{d/2}}, $$

as $n \to \infty$ explicitly, where $\mu_n$ and $\sigma_n^2$ are determined by the recursive relations (3.5) and we have used that $\mu_n \to 0$ and $\sigma_n^2 \to \sigma_\alpha^2$ according to Lemma 3.8. In the following lemma, we provide a similar quantitative result with convergence rates for generic initial data $F_0 \in \mathcal{M}_+(\mathbb{R}^d)$.

**Lemma 6.6.** Assume that $\alpha \in \mathbb{R}_+^*$, set any initial datum $F_0 \in \mathcal{M}_+(\mathbb{R}^d)$ and consider the solution $\{F_n\}_{n \in \mathbb{N}}$ to the time-discrete problem (1.1). For any $\epsilon \in \mathbb{R}_+^*$ there exists $C_\epsilon \in \mathbb{R}_+^*$ sufficiently large with

$$ \left| \frac{F_n(0)}{\|F_n\|_{L^1(\mathbb{R}^d)}} - \mathbf{F}_\alpha(0) \right| \le C_\epsilon (2 \mathbf{k}_\alpha + \epsilon)^n, \quad (6.20) $$

for any $n \in \mathbb{N}$.

*Proof.* Let us fix any $\epsilon \in \mathbb{R}_+^*$ arbitrarily small so that $2\mathbf{k}_\alpha + \epsilon < 1$. Again, this is possible because $2\mathbf{k}_\alpha < 1$ for $\alpha \in \mathbb{R}_+^*$. For simplicity of notation, we shall define the following sequence of coefficients:

$$ c_n := \frac{F_n(0)}{\|F_n\|_{L^1(\mathbb{R}^d)}}, \quad n \in \mathbb{N}. $$

Our goal then reduces to studying the asymptotic behavior of $\{c_n\}_{n \in \mathbb{N}}$ and obtaining quantitative convergence rates. Note that by formula (5.20) in Proposition 5.9 for the recursion, we obtain

$$ \frac{1}{c_n} = \frac{\int_{\mathbb{R}^d} e^{-\frac{1+\alpha-k_n}{2}|x|^2} \mathcal{E}_n[\bar{F}_0](x) dx}{\mathcal{E}_n[\bar{F}_0](0)} = \int_{|x| \le R_n} e^{-\frac{1+\alpha-k_n}{2}|x|^2} \frac{\mathcal{E}_n[\bar{F}_0](x)}{\mathcal{E}_n[\bar{F}_0](0)} dx + \frac{1}{c_n} \int_{|x| > R_n} \frac{F_n(x)}{\|F_n\|_{L^1(\mathbb{R}^d)}} dx, $$

where we set the sequence of radii $\{R_n\}_{n \in \mathbb{N}}$ as follows

$$ R_n := \frac{1}{(2 \mathbf{k}_\alpha)^{\frac{n}{2m+1}}}, \quad n \in \mathbb{N}, \quad (6.21) $$

for a fixed value $m \in \mathbb{N}$, which we take large enough so that

$$ (2 \mathbf{k}_\alpha)^{\frac{2m}{2m+1}} \le 2 \mathbf{k}_\alpha + \frac{\epsilon}{2}. \quad (6.22) $$

Whilst our choice of $R_n$ is not justified at first glance, we claim that it has been taken as to minimize the decay rate on the error terms $\tilde{E}_n$ and $E_n$ below. By solving the implicit equation on $c_n$ we infer

$$ c_n = \frac{1 - \tilde{E}_n}{(2\pi \sigma_\alpha^2)^{d/2} + E_n}, \quad (6.23) $$
which yields the announced estimate (6.6).

**Remark 6.4.** Under the assumptions in Lemma 6.3, notice that the above result implies in particular that for every $\varepsilon \in \mathbb{R}_+^*$ and $R \in \mathbb{R}_+^*$, there exists a large enough constant $C_{\varepsilon,R} \in \mathbb{R}_+^*$ such that

$$ \sup_{|x| \le R} |\nabla \log \mathcal{E}_n[F](x)| \le C_{\varepsilon,R} (2 \mathbf{k}_\alpha)^n, $$

for any $n \in \mathbb{N}$. In particular, $\nabla \log \mathcal{E}_n[F]$ converge to zero as $n \to \infty$ uniformly on compact sets.

The above observation is the cornerstone to prove the following local version of our main Theorem 1.1.

**Corollary 6.5 (Local convergence result).** Assume that $\alpha \in \mathbb{R}_+^*$, set any initial datum $F_0 \in \mathcal{M}_+(\mathbb{R}^d)$ and consider the solution $\{F_n\}_{n \in \mathbb{N}}$ to the time-discrete problem (1.1). Then, for every $\varepsilon \in \mathbb{R}_+^*$ and $R \in \mathbb{R}_+^*$, there exists a large enough constant $C_{\varepsilon,R} \in \mathbb{R}_+^*$ such that

$$ \sup_{|x| \le R} \left| \frac{F_n(x)}{\|F_n\|_{L^1(\mathbb{R}^d)}} - \mathbf{F}_\alpha(x) \right| \le C_{\varepsilon,R} (2 \mathbf{k}_\alpha)^n, \quad (6.19) $$

for any $n \in \mathbb{N}$, where $\mathbf{F}_\alpha$ is the Gaussian profile (1.12). In particular, $F_n/\|F_n\|_{L^1(\mathbb{R}^d)}$ converge to $\mathbf{F}_\alpha$ as $n \to \infty$ uniformly on compact sets.

In the proof of Corollary 6.5 we will require a careful control of the asymptotic behavior of the ratios $F_n(0)/\|F_n\|_{L^1(\mathbb{R}^d)}$ as $n \to \infty$, which we provide below. First, note that if $\{F_n\}_{n \in \mathbb{N}}$ is the solution (3.4) of the time-discrete problem (1.1) issued at a Gaussian initial datum $F_0$ like in Proposition 3.7, we have

$$ \frac{F_n(0)}{\|F_n\|_{L^1(\mathbb{R}^d)}} = G_{\mu_n, \sigma_n^2}(0) \to \mathbf{F}_\alpha(0) = \frac{1}{(2\pi \sigma_\alpha^2)^{d/2}}, $$

as $n \to \infty$ explicitly, where $\mu_n$ and $\sigma_n^2$ are determined by the recursive relations (3.5) and we have used that $\mu_n \to 0$ and $\sigma_n^2 \to \sigma_\alpha^2$ according to Lemma 3.8. In the following lemma, we provide a similar quantitative result with convergence rates for generic initial data $F_0 \in \mathcal{M}_+(\mathbb{R}^d)$.

**Lemma 6.6.** Assume that $\alpha \in \mathbb{R}_+^*$, set any initial datum $F_0 \in \mathcal{M}_+(\mathbb{R}^d)$ and consider the solution $\{F_n\}_{n \in \mathbb{N}}$ to the time-discrete problem (1.1). For any $\varepsilon \in \mathbb{R}_+^*$ there exists $C_\varepsilon \in \mathbb{R}_+^*$ sufficiently large with

$$ \left| \frac{F_n(0)}{\|F_n\|_{L^1(\mathbb{R}^d)}} - \mathbf{F}_\alpha(0) \right| \le C_\varepsilon (2 \mathbf{k}_\alpha + \varepsilon)^n, \quad (6.20) $$

for any $n \in \mathbb{N}$.

*Proof.* Let us fix any $\varepsilon \in \mathbb{R}_+^*$ arbitrarily small so that $2\mathbf{k}_\alpha + \varepsilon < 1$. Again, this is possible because $2\mathbf{k}_\alpha < 1$ for $\alpha \in \mathbb{R}_+^*$. For simplicity of notation, we shall define the following sequence of coefficients:

$$ c_n := \frac{F_n(0)}{\|F_n\|_{L^1(\mathbb{R}^d)}}, \quad n \in \mathbb{N}. $$

Our goal then reduces to studying the asymptotic behavior of $\{c_n\}_{n \in \mathbb{N}}$ and obtaining quantitative convergence rates. Note that by formula (5.20) in Proposition 5.9 for the recursion, we obtain

$$ \frac{1}{c_n} = \frac{\int_{\mathbb{R}^d} e^{-\frac{1+\alpha-k_n}{2}|x|^2} \mathcal{E}_n[\bar{F}_0](x) dx}{\mathcal{E}_n[\bar{F}_0](0)} = \int_{|x| \le R_n} e^{-\frac{1+\alpha-k_n}{2}|x|^2} \frac{\mathcal{E}_n[\bar{F}_0](x)}{\mathcal{E}_n[\bar{F}_0](0)} dx + \frac{1}{c_n} \int_{|x| > R_n} \frac{F_n(x)}{\|F_n\|_{L^1(\mathbb{R}^d)}} dx, $$

where we set the sequence of radii $\{R_n\}_{n \in \mathbb{N}}$ as follows

$$ R_n := \frac{1}{(2 \mathbf{k}_\alpha)^{\frac{n}{2m+1}}}, \quad n \in \mathbb{N}, \quad (6.21) $$

for a fixed value $m \in \mathbb{N}$, which we take large enough so that

$$ (2 \mathbf{k}_\alpha)^{\frac{2m}{2m+1}} \le 2 \mathbf{k}_\alpha + \frac{\varepsilon}{2}. \quad (6.22) $$

Whilst our choice of $R_n$ is not justified at first glance, we claim that it has been taken as to minimize the decay rate on the error terms $\tilde{E}_n$ and $E_n$ below. By solving the implicit equation on $c_n$ we infer

$$ c_n = \frac{1 - \tilde{E}_n}{(2\pi \sigma_\alpha^2)^{d/2} + E_n}, \quad (6.23) $$
for every $n \in \mathbb{N}$, where the error terms $\tilde{E}_n$ and $E_n$ take the form

$$ 
\begin{aligned}
\tilde{E}_n &:= \int_{|x|>R_n} \frac{F_n(x)}{\|F_n\|_{L^1(\mathbb{R}^d)}} dx, \\
E_n &:= \int_{|x|\le R_n} e^{-\frac{1+\alpha-k_n}{2}|x|^2} \frac{\mathcal{E}_n[\bar{F}_0](x)}{\mathcal{E}_n[\bar{F}_0](0)} dx - (2\pi \boldsymbol{\sigma}_\alpha^2)^{d/2},
\end{aligned}
\quad (6.24) $$

Additionally, we can split the second error term as $E_n = E_{n,1} - E_{n,2} + E_{n,3}$ with

$$ 
\begin{aligned}
E_{n,1} &:= \left( \frac{2\pi}{1+\alpha-k_n} \right)^{d/2} - \left( \frac{2\pi}{1+\alpha-\mathbf{k}_\alpha} \right)^{d/2}, \\
E_{n,2} &:= \int_{|x|>R_n} e^{-\frac{1+\alpha-k_n}{2}|x|^2} dx, \\
E_{n,3} &:= \int_{|x|\le R_n} e^{-\frac{1+\alpha-k_n}{2}|x|^2} \left( \frac{\mathcal{E}_n[\bar{F}_0](x)}{\mathcal{E}_n[\bar{F}_0](0)} - 1 \right) dx,
\end{aligned}
\quad (6.25) $$

where we have used the relationship $\boldsymbol{\sigma}_\alpha^2(1 + \alpha - \mathbf{k}_\alpha) = 1$, which results from (1.13) and (1.15), and therefore the decomposition of $E_n$ becomes clear since $E_{n,1}$ can be reformulated as

$$ E_{n,1} = \int_{\mathbb{R}^d} e^{-\frac{1+\alpha-k_n}{2}|x|^2} dx - (2\pi \boldsymbol{\sigma}_\alpha^2)^{d/2}. $$

On the one hand, given any fixed value $\theta \in \mathbb{R}^+$ with $\theta < \frac{\alpha}{2}$ (for instance $\theta = \frac{\alpha}{4}$) the propagation of exponential moments in Corollary 4.8 implies that

$$ E_{\theta} := \sup_{n \in \mathbb{N}} \int_{\mathbb{R}^d} e^{\theta|x|^2} \frac{F_n(x)}{\|F_n\|_{L^1(\mathbb{R}^d)}} dx < \infty. $$

By expanding the exponential in power series, one obtains in particular uniformly bounded moments of any order, and in particular, of order $2m$, namely

$$ M_m := \sup_{n \in \mathbb{N}} \int_{\mathbb{R}^d} |x|^{2m} \frac{F_n(x)}{\|F_n\|_{L^1(\mathbb{R}^d)}} dx \le \frac{E_\theta m!}{\theta^m} < \infty. $$

Hence, we infer the following control on the first error term in (6.24)

$$ \tilde{E}_n \le \frac{1}{R_n^{2m}} \int_{\mathbb{R}^d} |x|^{2m} \frac{F_n(x)}{\|F_n\|_{L^1(\mathbb{R}^d)}} dx \lesssim \frac{1}{R_n^{2m}}, \quad (6.26) $$

for any $n \in \mathbb{N}$. On the other hand, by Lemma 5.3 we have that $E_{n,1}$ in (6.25) can be bounded by

$$ |E_{n,1}| \lesssim k_n - \mathbf{k}_\alpha \lesssim \mathbf{r}_\alpha^n. \quad (6.27) $$

By direct calculation we also obtain

$$ E_{n,2} = |\mathbb{S}^{d-1}| \int_{R_n}^{+\infty} r^{d-1} e^{-\frac{1+\alpha-k_n}{2}r^2} dr \lesssim \Gamma\left(\frac{d}{2}, \frac{1+\alpha-k_n}{2}R_n^2\right) \lesssim \exp\left(-\frac{1+\alpha-k_n}{4(2\mathbf{k}_\alpha)^{\frac{2n}{2m+1}}}\right), \quad (6.28) $$

where $\Gamma(a, x)$ is the incomplete Gamma function (6.2) and we have used estimate (6.3) in Lemma 6.2. Note that the constraint (6.4) is trivially satisfied since $R_n \to \infty$ by our choice (6.21). Finally, we control the error term $E_{n,3}$ in (6.25) under the addition assumption that $\bar{F}_0$ satisfies the hypothesis (6.5) in Lemma 6.3. Whilst this is not always true, we show at the end of the proof that we can always assume so without loss of generality by replacing the argument on $\bar{F}_0$ by an advance enough time step $\bar{F}_m$ so that selection has properly shaped the Gaussian tails. Under this condition, Lemma 6.3 implies that given any $\varepsilon' \in \mathbb{R}_+^*$ there is $C_{\varepsilon'} \in \mathbb{R}_+^*$ so that, for $|x| \le R_n$, we have

$$ 
\begin{aligned}
& \left| \int_0^1 (\nabla \log \mathcal{E}_n)[\bar{F}_0](\theta x) \cdot x \, d\theta \right| \\
& \lesssim (2\mathbf{k}_\alpha + \varepsilon')^n R_n + (\mathbf{r}_\alpha + \varepsilon')^n R_n^2 + e^{-\frac{(2+\varepsilon')^n}{\varepsilon'} \mathcal{E}_{\varepsilon'}(2\mathbf{r}_\alpha + \varepsilon')^n R_n^2} R_n \\
& \lesssim \left( \frac{2\mathbf{k}_\alpha}{(2\mathbf{k}_\alpha)^{\frac{1}{2m+1}}} + \frac{\varepsilon}{2} \right)^n + \left( \frac{\mathbf{r}_\alpha}{(2\mathbf{k}_\alpha)^{\frac{2}{2m+1}}} + \frac{\varepsilon}{2} \right)^n + e^{-\frac{(2+\varepsilon')^n}{\varepsilon'} \exp \left[ C_{\varepsilon'} \left( \frac{2\mathbf{r}_\alpha}{(2\mathbf{k}_\alpha)^{\frac{2}{2m+1}}} + \frac{\varepsilon}{2} \right)^n \right]} (2\mathbf{k}_\alpha)^{\frac{n}{2m+1}} \\
& \lesssim \left( (2\mathbf{k}_\alpha)^{\frac{2m}{2m+1}} + \frac{\varepsilon}{2} \right)^n \lesssim (2\mathbf{k}_\alpha + \varepsilon')^n.
\end{aligned}
$$

In the second inequality above we have taken $\varepsilon'$ small enough compared with $\varepsilon$. In the third inequality we have used the relation $\mathbf{r}_\alpha = 2\mathbf{k}_\alpha^2$, which guarantees that all the contributions in the third line are
for every $n \in \mathbb{N}$, where the error terms $\tilde{E}_n$ and $E_n$ take the form

$$ 
\begin{aligned}
\tilde{E}_n &:= \int_{|x|>R_n} \frac{F_n(x)}{\|F_n\|_{L^1(\mathbb{R}^d)}} dx, \\
E_n &:= \int_{|x|\le R_n} e^{-\frac{1+\alpha-k_n}{2}|x|^2} \frac{\mathcal{E}_n[\bar{F}_0](x)}{\mathcal{E}_n[\bar{F}_0](0)} dx - (2\pi \boldsymbol{\sigma}_\alpha^2)^{d/2},
\end{aligned}
\quad (6.24) $$

Additionally, we can split the second error term as $E_n = E_{n,1} - E_{n,2} + E_{n,3}$ with

$$ 
\begin{aligned}
E_{n,1} &:= \left( \frac{2\pi}{1+\alpha-k_n} \right)^{d/2} - \left( \frac{2\pi}{1+\alpha-\mathbf{k}_\alpha} \right)^{d/2}, \\
E_{n,2} &:= \int_{|x|>R_n} e^{-\frac{1+\alpha-k_n}{2}|x|^2} dx, \\
E_{n,3} &:= \int_{|x|\le R_n} e^{-\frac{1+\alpha-k_n}{2}|x|^2} \left( \frac{\mathcal{E}_n[\bar{F}_0](x)}{\mathcal{E}_n[\bar{F}_0](0)} - 1 \right) dx,
\end{aligned}
\quad (6.25) $$

where we have used the relationship $\boldsymbol{\sigma}_\alpha^2(1 + \alpha - \mathbf{k}_\alpha) = 1$, which results from (1.13) and (1.15), and therefore the decomposition of $E_n$ becomes clear since $E_{n,1}$ can be reformulated as

$$ E_{n,1} = \int_{\mathbb{R}^d} e^{-\frac{1+\alpha-k_n}{2}|x|^2} dx - (2\pi \boldsymbol{\sigma}_\alpha^2)^{d/2}. $$

On the one hand, given any fixed value $\theta \in \mathbb{R}^+$ with $\theta < \frac{\alpha}{2}$ (for instance $\theta = \frac{\alpha}{4}$) the propagation of exponential moments in Corollary 4.8 implies that

$$ E_{\theta} := \sup_{n \in \mathbb{N}} \int_{\mathbb{R}^d} e^{\theta|x|^2} \frac{F_n(x)}{\|F_n\|_{L^1(\mathbb{R}^d)}} dx < \infty. $$

By expanding the exponential in power series, one obtains in particular uniformly bounded moments of any order, and in particular, of order $2m$, namely

$$ M_m := \sup_{n \in \mathbb{N}} \int_{\mathbb{R}^d} |x|^{2m} \frac{F_n(x)}{\|F_n\|_{L^1(\mathbb{R}^d)}} dx \le \frac{E_\theta m!}{\theta^m} < \infty. $$

Hence, we infer the following control on the first error term in (6.24)

$$ \tilde{E}_n \le \frac{1}{R_n^{2m}} \int_{\mathbb{R}^d} |x|^{2m} \frac{F_n(x)}{\|F_n\|_{L^1(\mathbb{R}^d)}} dx \lesssim \frac{1}{R_n^{2m}}, \quad (6.26) $$

for any $n \in \mathbb{N}$. On the other hand, by Lemma 5.3 we have that $E_{n,1}$ in (6.25) can be bounded by

$$ |E_{n,1}| \lesssim k_n - \mathbf{k}_\alpha \lesssim \mathbf{r}_\alpha^n. \quad (6.27) $$

By direct calculation we also obtain

$$ E_{n,2} = |\mathbb{S}^{d-1}| \int_{R_n}^{+\infty} r^{d-1} e^{-\frac{1+\alpha-k_n}{2}r^2} dr \lesssim \Gamma\left(\frac{d}{2}, \frac{1+\alpha-k_n}{2}R_n^2\right) \lesssim \exp\left(-\frac{1+\alpha-k_n}{4(2\mathbf{k}_\alpha)^{\frac{2n}{2m+1}}}\right), \quad (6.28) $$

where $\Gamma(a, x)$ is the incomplete Gamma function (6.2) and we have used estimate (6.3) in Lemma 6.2. Note that the constraint (6.4) is trivially satisfied since $R_n \to \infty$ by our choice (6.21). Finally, we control the error term $E_{n,3}$ in (6.25) under the addition assumption that $\bar{F}_0$ satisfies the hypothesis (6.5) in Lemma 6.3. Whilst this is not always true, we show at the end of the proof that we can always assume so without loss of generality by replacing the argument on $\bar{F}_0$ by an advance enough time step $\bar{F}_m$ so that selection has properly shaped the Gaussian tails. Under this condition, Lemma 6.3 implies that given any $\varepsilon' \in \mathbb{R}_+^*$ there is $C_{\varepsilon'} \in \mathbb{R}_+^*$ so that, for $|x| \le R_n$, we have

$$ 
\begin{aligned}
& \left| \int_0^1 (\nabla \log \mathcal{E}_n)[\bar{F}_0](\theta x) \cdot x \, d\theta \right| \\
& \lesssim (2\mathbf{k}_\alpha + \varepsilon')^n R_n + (\mathbf{r}_\alpha + \varepsilon')^n R_n^2 + e^{-\frac{(2+\varepsilon')^n}{\varepsilon'} \mathcal{E}_{\varepsilon'}(2\mathbf{r}_\alpha + \varepsilon')^n R_n^2} R_n \\
& \lesssim \left( \frac{2\mathbf{k}_\alpha}{(2\mathbf{k}_\alpha)^{\frac{1}{2m+1}}} + \frac{\varepsilon}{2} \right)^n + \left( \frac{\mathbf{r}_\alpha}{(2\mathbf{k}_\alpha)^{\frac{2}{2m+1}}} + \frac{\varepsilon}{2} \right)^n + e^{-\frac{(2+\varepsilon')^n}{\varepsilon'} \exp \left[ C_{\varepsilon'} \left( \frac{2\mathbf{r}_\alpha}{(2\mathbf{k}_\alpha)^{\frac{2}{2m+1}}} + \frac{\varepsilon}{2} \right)^n \right]} \frac{1}{(2\mathbf{k}_\alpha)^{\frac{n}{2m+1}}} \\
& \lesssim \left( (2\mathbf{k}_\alpha)^{\frac{2m}{2m+1}} + \frac{\varepsilon}{2} \right)^n \lesssim (2\mathbf{k}_\alpha + \varepsilon')^n.
\end{aligned}
$$

In the second inequality above we have taken $\varepsilon'$ small enough compared with $\varepsilon$. In the third inequality we have used the relation $\mathbf{r}_\alpha = 2\mathbf{k}_\alpha^2$, which guarantees that all the contributions in the third line are
dominated by the first term $((2\mathbf{k}_\alpha)^{\frac{2m}{2m+1}+\frac{\varepsilon}{2}})^n$. Finally, in the last inequality we have used our choice of $m \in \mathbb{N}$ large enough so that (6.22) holds. Therefore, the mean value theorem implies

$$ \left| \frac{\mathcal{E}_n[\bar{F}_0](x)}{\mathcal{E}_n[\bar{F}_0](0)} - 1 \right| = \left| \exp \left( \int_0^1 (\nabla \log \mathcal{E}_n[\bar{F}_0])(\theta x) \cdot x \, d\theta \right) - 1 \right| \lesssim (2\mathbf{k}_\alpha + \varepsilon)^n. $$

Hence, $E_{n,3}$ in (6.25) can be controlled by

$$ |E_{n,3}| \le \int_{|x| \le R_n} e^{-\frac{1+\alpha-k_n}{2}|x|^2} \left| \frac{\mathcal{E}_n[\bar{F}_0](x)}{\mathcal{E}_n[\bar{F}_0](0)} - 1 \right| dx \lesssim (2\mathbf{k}_\alpha + \varepsilon)^n. \quad (6.29) $$

Putting (6.27), (6.28) and (6.29) together yields

$$ |E_n| \lesssim (2\mathbf{k}_\alpha + \varepsilon)^n. \quad (6.30) $$

Thereby, using the mean value theorem on (6.23) along with the bounds (6.26) and (6.30) concludes that

$$ \left| c_n - \frac{1}{(2\pi \sigma_{\alpha}^2)^{d/2}} \right| \lesssim \tilde{E}_n + |E_n| \lesssim (2\mathbf{k}_{\alpha} + \varepsilon)^n. $$

To end the proof, we show that we can always assume that the hypothesis (6.5) in Lemma 6.3 is satisfied without loss of generality. Indeed, note that after two iterations $\bar{F}_2$ already satisfies (6.5)$_2$ thanks to Corollary 4.2, but only two iterations are not enough to guarantee (6.5)$_1$, in general. Nevertheless, applying Corollary 4.2 leads to the following upper bound on the normalized profiles:

$$ \bar{F}_m(x) \le C_m e^{\frac{1}{2}\left(\alpha + \frac{1}{2} - \frac{1}{\bar{\sigma}_m^2}\right)|x|^2}, $$

for appropriate $C_m \in \mathbb{R}_+^*$, with variances $\{\bar{\sigma}_m^2\}_{m \in \mathbb{N}}$ defined by the recurrence (4.5), i.e.,

$$ \frac{1}{\bar{\sigma}_{m+1}^2} = \alpha + \frac{1}{1 + \frac{\bar{\sigma}_m^2}{2}}, \quad m \in \mathbb{N}, $$

and with initial datum $\bar{\sigma}_1^2 = 1/\alpha$. The boundedness condition (6.5)$_2$ is then satisfied by $\bar{F}_m$ if we can show that the prefactor $\alpha + \frac{1}{2} - \frac{1}{\bar{\sigma}_m^2}$ in the above exponential bound becomes non-positive for large enough $m$. This is where the precise choice of normalization in Definition 5.1 plays a role. Specifically, recall that so defined $\bar{\sigma}_m^2 \to \sigma_\alpha^2$ and we have precise relaxation rates by Lemma 3.8. Hence,

$$ \alpha + \frac{1}{2} - \frac{1}{\bar{\sigma}_m^2} = \left( \alpha + \frac{1}{2} - \frac{1}{\sigma_\alpha^2} \right) + \left( \frac{1}{\sigma_\alpha^2} - \frac{1}{\bar{\sigma}_m^2} \right) \le \left( \alpha + \frac{1}{2} - \frac{1}{\sigma_\alpha^2} \right) + C_v r_\alpha^m, $$

for any $m \in \mathbb{N}$. Since we chose a normalization $\bar{F}_n$ of the profiles by a Gaussian $G_{0,\sigma^2}$ with variance $\sigma^2$ strictly larger than the variance $\sigma_\alpha^2$ of the equilibrium $F_\alpha$ (more specifically $\frac{1}{\sigma_\alpha^2} = \frac{1}{\sigma^2} = \alpha + \frac{1}{2}$), we can guarantee that the right hand side above is non-positive if $m \ge n_\alpha$ for a sufficiently large $n_\alpha \ge 2$ (depending only on $\alpha$). This justifies that $\bar{F}_{n_\alpha}$ satisfies both conditions in (6.5). $\square$

We are now in position to prove the above local convergence result in Corollary 6.5

*Proof of Corollary 6.5.* By formula (5.20) in Proposition 5.9 and Definition 6.1 we obtain

$$ \nabla \log F_n(x) = -(1 + \alpha - k_n)x + \nabla \log \mathcal{E}_n[\bar{F}_0](x), $$

for any $x \in \mathbb{R}^d$ and every $n \in \mathbb{N}$. Using the mean value theorem we then achieve

$$ \begin{align} \log \frac{F_n(x)}{F_n(0)} - \log \frac{F_\alpha(x)}{F_\alpha(0)} &= \int_0^1 (\nabla \log F_n)(\theta x) \cdot x \, d\theta + \frac{1}{2\sigma_\alpha^2} |x|^2 \nonumber \\ &= -\frac{1+\alpha-k_n}{2} |x| + \int_0^1 \nabla \log \mathcal{E}_n[\bar{F}_0](\theta x) \cdot x \, d\theta + \frac{1}{2\sigma_\alpha^2} |x|^2 \tag{6.31} \\ &= \frac{k_n-k_\alpha}{2} |x|^2 + \int_0^1 \nabla \log \mathcal{E}_n[\bar{F}_0](\theta x) \cdot x \, d\theta, \nonumber \end{align} $$

where we have used the relation $1/\sigma_\alpha^2 = 1 + \alpha - k_\alpha$. On the one hand, the first term in the right hand side of (6.31) converges to zero thanks to Lemma 5.3. On the other hand, for the second term we shall apply Lemma 6.3. To do so we need to guarantee again that $\bar{F}_0$ verifies the hypothesis (6.5) of such a lemma. Recall that we can always assume the condition (6.5) without loss of generality (cf. last step in the proof of Lemma 6.6). Therefore, we can apply Remark 6.4 with $F = \bar{F}_0$ and obtain that $\nabla \log \mathcal{E}_n[\bar{F}_0]$
dominated by the first term $((2\mathbf{k}_\alpha)^{\frac{2m}{2m+1} + \frac{\epsilon}{2}})^n$. Finally, in the last inequality we have used our choice of $m \in \mathbb{N}$ large enough so that (6.22) holds. Therefore, the mean value theorem implies

$$ \left| \frac{\mathcal{E}_n[\bar{F}_0](x)}{\mathcal{E}_n[\bar{F}_0](0)} - 1 \right| = \left| \exp \left( \int_0^1 (\nabla \log \mathcal{E}_n[\bar{F}_0])(\theta x) \cdot x \, d\theta \right) - 1 \right| \lesssim (2\mathbf{k}_\alpha + \epsilon)^n. $$

Hence, $E_{n,3}$ in (6.25) can be controlled by

$$ |E_{n,3}| \le \int_{|x| \le R_n} e^{-\frac{1+\alpha-k_n}{2}|x|^2} \left| \frac{\mathcal{E}_n[\bar{F}_0](x)}{\mathcal{E}_n[\bar{F}_0](0)} - 1 \right| dx \lesssim (2\mathbf{k}_\alpha + \epsilon)^n. \quad (6.29) $$

Putting (6.27), (6.28) and (6.29) together yields

$$ |E_n| \lesssim (2\mathbf{k}_\alpha + \epsilon)^n. \quad (6.30) $$

Thereby, using the mean value theorem on (6.23) along with the bounds (6.26) and (6.30) concludes that

$$ \left| c_n - \frac{1}{(2\pi \sigma_{\alpha}^2)^{d/2}} \right| \lesssim \tilde{E}_n + |E_n| \lesssim (2 \mathbf{k}_{\alpha} + \epsilon)^n. $$

To end the proof, we show that we can always assume that the hypothesis (6.5) in Lemma 6.3 is satisfied without loss of generality. Indeed, note that after two iterations $\bar{F}_2$ already satisfies (6.5)$_2$ thanks to Corollary 4.2, but only two iterations are not enough to guarantee (6.5)$_1$, in general. Nevertheless, applying Corollary 4.2 leads to the following upper bound on the normalized profiles:

$$ \bar{F}_m(x) \le C_m e^{\frac{1}{2}\left(\alpha + \frac{1}{2} - \frac{1}{\bar{\sigma}_m^2}\right)|x|^2}, $$

for appropriate $C_m \in \mathbb{R}_+^*$, with variances $\{\bar{\sigma}_m^2\}_{m \in \mathbb{N}}$ defined by the recurrence (4.5), i.e.,

$$ \frac{1}{\bar{\sigma}_{m+1}^2} = \alpha + \frac{1}{1 + \frac{\bar{\sigma}_m^2}{2}}, \quad m \in \mathbb{N}, $$

and with initial datum $\bar{\sigma}_1^2 = 1/\alpha$. The boundedness condition (6.5)$_2$ is then satisfied by $\bar{F}_m$ if we can show that the prefactor $\alpha + \frac{1}{2} - \frac{1}{\bar{\sigma}_m^2}$ in the above exponential bound becomes non-positive for large enough $m$. This is where the precise choice of normalization in Definition 5.1 plays a role. Specifically, recall that so defined $\bar{\sigma}_m^2 \to \sigma_\alpha^2$ and we have precise relaxation rates by Lemma 3.8. Hence,

$$ \alpha + \frac{1}{2} - \frac{1}{\bar{\sigma}_m^2} = \left( \alpha + \frac{1}{2} - \frac{1}{\sigma_\alpha^2} \right) + \left( \frac{1}{\sigma_\alpha^2} - \frac{1}{\bar{\sigma}_m^2} \right) \le \left( \alpha + \frac{1}{2} - \frac{1}{\sigma_\alpha^2} \right) + C_v r_\alpha^m, $$

for any $m \in \mathbb{N}$. Since we chose a normalization $\bar{F}_n$ of the profiles by a Gaussian $G_{0,\sigma^2}$ with variance $\sigma^2$ strictly larger than the variance $\sigma_\alpha^2$ of the equilibrium $F_\alpha$ (more specifically $\frac{1}{\sigma_\alpha^2} = \frac{1}{\sigma^2} = \alpha + \frac{1}{2}$), we can guarantee that the right hand side above is non-positive if $m \ge n_\alpha$ for a sufficiently large $n_\alpha \ge 2$ (depending only on $\alpha$). This justifies that $\bar{F}_{n_\alpha}$ satisfies both conditions in (6.5). $\square$

We are now in position to prove the above local convergence result in Corollary 6.5

*Proof of Corollary 6.5.* By formula (5.20) in Proposition 5.9 and Definition 6.1 we obtain

$$ \nabla \log F_n(x) = -(1 + \alpha - k_n)x + \nabla \log \mathcal{E}_n[\bar{F}_0](x), $$

for any $x \in \mathbb{R}^d$ and every $n \in \mathbb{N}$. Using the mean value theorem we then achieve

$$ \begin{align} \log \frac{F_n(x)}{F_n(0)} - \log \frac{F_\alpha(x)}{F_\alpha(0)} &= \int_0^1 (\nabla \log F_n)(\theta x) \cdot x \, d\theta + \frac{1}{2 \sigma_\alpha^2} |x|^2 \nonumber \\ &= -\frac{1+\alpha-k_n}{2} |x| + \int_0^1 \nabla \log \mathcal{E}_n[\bar{F}_0](\theta x) \cdot x \, d\theta + \frac{1}{2 \sigma_\alpha^2} |x|^2 \tag{6.31} \\ &= \frac{k_n-k_\alpha}{2} |x|^2 + \int_0^1 \nabla \log \mathcal{E}_n[\bar{F}_0](\theta x) \cdot x \, d\theta, \nonumber \end{align} $$

where we have used the relation $1/\sigma_\alpha^2 = 1 + \alpha - k_\alpha$. On the one hand, the first term in the right hand side of (6.31) converges to zero thanks to Lemma 5.3. On the other hand, for the second term we shall apply Lemma 6.3. To do so we need to guarantee again that $\bar{F}_0$ verifies the hypothesis (6.5) of such a lemma. Recall that we can always assume the condition (6.5) without loss of generality (cf. last step in the proof of Lemma 6.6). Therefore, we can apply Remark 6.4 with $F = \bar{F}_0$ and obtain that $\nabla \log \mathcal{E}_n[\bar{F}_0]$
converges to zero uniformly over compact sets with explicit convergence rates. Putting it into (6.31) we obtain more explicitly

$$ \sup_{|x| \le R} \left| \log \frac{F_n(x)}{F_n(0)} - \log \frac{\mathbf{F}_\alpha(x)}{\mathbf{F}_\alpha(0)} \right| \le C_{\varepsilon,R} (2 \mathbf{k}_\alpha)^n, $$

for any $R \in \mathbb{R}_+^*$, any $n \in \mathbb{N}$ and a large enough $C_{\varepsilon,R} \in \mathbb{R}_+^*$. This, together with the above control on the asymptotics of $\|F_n\|_{L^1(\mathbb{R}^d)}/F_n(0)$ in Lemma 6.6, allow proving (6.19). $\square$

In particular, note that the above Corollary 6.5 is enough to prove the uniqueness of solutions to the eigenproblem (1.11), as stated in Corollary 1.2.

*Proof of Corollary 1.2.* Let $(\lambda, F)$ be any solution of the eigenproblem (1.11). Then, the ansatz (1.10) defines a solution $F_n = \lambda^n F$ of the time-discrete problem (1.1). By Corollary 6.5 we obtain that

$$ \lim_{n \to \infty} \sup_{|x| \le R} \left| \frac{\lambda^n F(x)}{\lambda^n} - \mathbf{F}_\alpha(x) \right| = 0, $$

for any $R \in \mathbb{R}_+^*$. Hence, $F \equiv \mathbf{F}_\alpha$, and thus $\lambda = \lambda_\alpha$. $\square$

Whilst the above local convergence result is enough to identify asymptotically the profile $\mathbf{F}_\alpha$, a global result with quantitative convergence rates is still missing. In particular, note that the constants $C_{\varepsilon,R}$ above blow up when $R \to \infty$ as we see explicitly in Lemma 6.3. A second drawback of this we are unable to characterize the long-time behavior of the mass $\|F_n\|_{L^1(\mathbb{R}^n)}$ in terms of the eigenvalue $\lambda_\alpha$ via this method. In the following section, we give an answer to both questions by better exploiting the previous fundamental Lemma 6.3 and using the propagation of quadratic and exponential moments in Section 4.2.

## 6.2. Global convergence result.
We are now ready to prove our main result. Let us emphasize that our final convergence result in Theorem 1.1 is presented using the Kullback-Leibler divergence, which is a very different metric from the log-Lipschitz type metrics used in the previous Section 6.1 for the local convergence results. As anticipated in Remark 1.4, this decision does not obey aesthetic reasons only, but we actually need to move from uniform norms (like the log-Lipschitz norm) to averaged norms (like the Kullback-Leibler divergence) in order to address the deficiency encountered in Lemma 6.3. Specifically, recall that for generic initial data $F_0 \in \mathcal{M}_+(\mathbb{R}^d)$, the log-Lipschitz norms of the high-dimensional integral $\varepsilon_n[\bar{F}_0]$ cannot be controlled uniformly due to additional exponentially growing terms.

*Proof of Theorem 1.1.*

*STEP 1:* Convergence of the profiles $F_n/\|F_n\|_{L^1(\mathbb{R}^d)}$.

Since $\mathbf{F}_\alpha$ is Gaussian (thus strongly log-concave), then the logarithmic-Sobolev inequality holds true and therefore we obtain the following control of the relative entropy by the relative Fisher information

$$ \mathcal{D}_{\text{KL}} \left( \frac{F_n}{\|F_n\|} \middle| \mathbf{F}_{\alpha} \right) \le \frac{\sigma_{\alpha}^2}{2} \int_{\mathbb{R}^d} \left| \nabla \log \left( \frac{F_n(x)}{\mathbf{F}_{\alpha}(x)} \right) \right|^2 \frac{F_n(x)}{\|F_n\|_{L^1(\mathbb{R}^d)}} dx, $$

see Corollary 5.7.2 and Section 9.3.1 in [2] for details. By the reformulation of $F_n$ in formula (5.20) of Proposition 5.9 and using the notation in Definition 6.1 we have that

$$ \nabla \log \left( \frac{F_n(x)}{\mathbf{F}_\alpha(x)} \right) = (k_n - \mathbf{k}_\alpha) x + \nabla \log \varepsilon_n[\bar{F}_0](x). $$

Therefore, we obtain the following upper bound

$$ \mathcal{D}_{\text{KL}} \left( \frac{F_n}{\|F_n\|} \middle| \mathbf{F}_{\alpha} \right) \lesssim D_{n,1} + D_{n,2}, \qquad (6.32) $$

where each factor takes the form

$$ D_{n,1} := |k_n - \mathbf{k}_\alpha|^2 \int_{\mathbb{R}^d} |x|^2 \frac{F_n(x)}{\|F_n\|_{L^1(\mathbb{R}^d)}} dx, $$

$$ D_{n,2} := \int_{\mathbb{R}^d} \left| \nabla \log \varepsilon_n[\bar{F}_0](x) \right|^2 \frac{F_n(x)}{\|F_n\|_{L^1(\mathbb{R}^d)}} dx. $$

On the one hand, for $D_{n,1}$ we can use the propagation of quadratic moments in Corollary 4.6 together with the explicit relaxation rates of $\{k_n\}_{n \in \mathbb{N}}$ to $\mathbf{k}_\alpha$ in Lemma 5.3 to show that

$$ D_{n,1} \lesssim \mathbf{r}_\alpha^{2n}. \qquad (6.33) $$
converges to zero uniformly over compact sets with explicit convergence rates. Putting it into (6.31) we obtain more explicitly

$$ \sup_{|x| \le R} \left| \log \frac{F_n(x)}{F_n(0)} - \log \frac{\mathbf{F}_\alpha(x)}{\mathbf{F}_\alpha(0)} \right| \le C_{\varepsilon,R} (2 \mathbf{k}_\alpha)^n, $$

for any $R \in \mathbb{R}_+^*$, any $n \in \mathbb{N}$ and a large enough $C_{\varepsilon,R} \in \mathbb{R}_+^*$. This, together with the above control on the asymptotics of $\|F_n\|_{L^1(\mathbb{R}^d)}/F_n(0)$ in Lemma 6.6, allow proving (6.19). $\square$

In particular, note that the above Corollary 6.5 is enough to prove the uniqueness of solutions to the eigenproblem (1.11), as stated in Corollary 1.2.

*Proof of Corollary 1.2.* Let $(\lambda, F)$ be any solution of the eigenproblem (1.11). Then, the ansatz (1.10) defines a solution $F_n = \lambda^n F$ of the time-discrete problem (1.1). By Corollary 6.5 we obtain that

$$ \lim_{n \to \infty} \sup_{|x| \le R} \left| \frac{\lambda^n F(x)}{\lambda^n} - \mathbf{F}_\alpha(x) \right| = 0, $$

for any $R \in \mathbb{R}_+^*$. Hence, $F \equiv \mathbf{F}_\alpha$, and thus $\lambda = \lambda_\alpha$. $\square$

Whilst the above local convergence result is enough to identify asymptotically the profile $\mathbf{F}_\alpha$, a global result with quantitative convergence rates is still missing. In particular, note that the constants $C_{\varepsilon,R}$ above blow up when $R \to \infty$ as we see explicitly in Lemma 6.3. A second drawback of this we are unable to characterize the long-time behavior of the mass $\|F_n\|_{L^1(\mathbb{R}^n)}$ in terms of the eigenvalue $\lambda_\alpha$ via this method. In the following section, we give an answer to both questions by better exploiting the previous fundamental Lemma 6.3 and using the propagation of quadratic and exponential moments in Section 4.2.

## 6.2. Global convergence result.
We are now ready to prove our main result. Let us emphasize that our final convergence result in Theorem 1.1 is presented using the Kullback-Leibler divergence, which is a very different metric from the log-Lipschitz type metrics used in the previous Section 6.1 for the local convergence results. As anticipated in Remark 1.4, this decision does not obey aesthetic reasons only, but we actually need to move from uniform norms (like the log-Lipschitz norm) to averaged norms (like the Kullback-Leibler divergence) in order to address the deficiency encountered in Lemma 6.3. Specifically, recall that for generic initial data $F_0 \in \mathcal{M}_+(\mathbb{R}^d)$, the log-Lipschitz norms of the high-dimensional integral $\varepsilon_n[\bar{F}_0]$ cannot be controlled uniformly due to additional exponentially growing terms.

*Proof of Theorem 1.1.*

• **STEP 1:** Convergence of the profiles $F_n/\|F_n\|_{L^1(\mathbb{R}^d)}$.

Since $\mathbf{F}_\alpha$ is Gaussian (thus strongly log-concave), then the logarithmic-Sobolev inequality holds true and therefore we obtain the following control of the relative entropy by the relative Fisher information

$$ \mathcal{D}_{\text{KL}} \left( \frac{F_n}{\|F_n\|} \middle| \mathbf{F}_{\alpha} \right) \le \frac{\sigma_{\alpha}^2}{2} \int_{\mathbb{R}^d} \left| \nabla \log \left( \frac{F_n(x)}{\mathbf{F}_{\alpha}(x)} \right) \right|^2 \frac{F_n(x)}{\|F_n\|_{L^1(\mathbb{R}^d)}} dx, $$

see Corollary 5.7.2 and Section 9.3.1 in [2] for details. By the reformulation of $F_n$ in formula (5.20) of Proposition 5.9 and using the notation in Definition 6.1 we have that

$$ \nabla \log \left( \frac{F_n(x)}{\mathbf{F}_\alpha(x)} \right) = (k_n - \mathbf{k}_\alpha) x + \nabla \log \varepsilon_n[\bar{F}_0](x). $$

Therefore, we obtain the following upper bound

$$ \mathcal{D}_{\text{KL}} \left( \frac{F_n}{\|F_n\|} \middle| \mathbf{F}_{\alpha} \right) \lesssim D_{n,1} + D_{n,2}, \qquad (6.32) $$

where each factor takes the form

$$ D_{n,1} := |k_n - \mathbf{k}_\alpha|^2 \int_{\mathbb{R}^d} |x|^2 \frac{F_n(x)}{\|F_n\|_{L^1(\mathbb{R}^d)}} dx, $$

$$ D_{n,2} := \int_{\mathbb{R}^d} \left| \nabla \log \varepsilon_n[\bar{F}_0](x) \right|^2 \frac{F_n(x)}{\|F_n\|_{L^1(\mathbb{R}^d)}} dx. $$

On the one hand, for $D_{n,1}$ we can use the propagation of quadratic moments in Corollary 4.6 together with the explicit relaxation rates of $\{k_n\}_{n \in \mathbb{N}}$ to $\mathbf{k}_\alpha$ in Lemma 5.3 to show that

$$ D_{n,1} \lesssim \mathbf{r}_{\alpha}^{2n}. \qquad (6.33) $$
On the other hand, for $D_{n,2}$ we shall refine the argument in the proof of Corollary 6.5. Specifically, recall again that we can assume that $\bar{F}_0$ satisfies the hypothesis (6.5) without loss of generality (cf. proof of Corollary 6.5). Then, applying Lemma 6.3 in order to control $\nabla \log \mathcal{E}_n[\bar{F}_0]$ with fixed $\varepsilon \in \mathbb{R}_+^*$ yields

$$D_{n,2} \lesssim ((2 \mathbf{k}_\alpha)^2 + \varepsilon)^n \int_{\mathbb{R}^d} \left(1 + |x|^2 + e^{2C_\varepsilon(2 \mathbf{r}_\alpha + \varepsilon)^n |x|^2}\right) \frac{F_n(x)}{\|F_n\|_{L^1(\mathbb{R}^d)}} dx,$$

for some sufficiently large $C_\varepsilon \in \mathbb{R}_+^*$. Above we have used that $(2 \mathbf{k}_\alpha)^n$ is the decay rate that controls all the others in formula (6.6). Take $n \ge n_\varepsilon$ sufficiently large such that

$$\theta_\varepsilon := \sup_{n \ge n_\varepsilon} 2 C_\varepsilon (2 \mathbf{r}_\alpha + \varepsilon)^n < \frac{\alpha}{2}.$$

Then, Corollaries 4.6 and 4.8 imply that the above quadratic and exponential moments are uniformly bounded with respect to $n$ for $n \ge n_\varepsilon$. Therefore, we have

$$D_{n,2} \lesssim ((2 \mathbf{k}_\alpha)^2 + \varepsilon)^n. \qquad (6.34)$$

Finally, putting the estimates (6.33) and (6.34) into the split (6.32), and noticing that $\mathbf{r}_\alpha^2 < (2 \mathbf{k}_\alpha)^2$ thanks to the relation $\mathbf{r}_\alpha = 2 \mathbf{k}_\alpha^2$, ends this part of the proof.

• STEP 2: Convergence of the growth rates $\|F_n\|_{L^1(\mathbb{R}^d)}/\|F_{n-1}\|_{L^1(\mathbb{R}^d)}$.

For simplicity of notation we shall define the following sequence of coefficients:

$$\lambda_n := \frac{\|F_n\|_{L^1(\mathbb{R}^d)}}{\|F_{n-1}\|_{L^1(\mathbb{R}^d)}}, \quad n \in \mathbb{N}.$$

Our goal then reduces to studying the asymptotic behavior of $\{\lambda_n\}_{n \in \mathbb{N}}$ and obtaining quantitative convergence rates. By definition of operator $\mathcal{T}$ in (1.2) we obtain that

$$\lambda_n = \int_{\mathbb{R}^d} e^{-m(x)} \int_{\mathbb{R}^{2d}} G\left(x - \frac{x_1 + x_2}{2}\right) \frac{F_{n-1}(x_1)}{\|F_{n-1}\|_{L^1(\mathbb{R}^d)}} \frac{F_{n-1}(x_2)}{\|F_{n-1}\|_{L^1(\mathbb{R}^d)}} dx_1 dx_2 dx.$$

By direct computation of the integral with respect to $x$ as it was done in (4.11) we find that

$$\lambda_n = \int_{\mathbb{R}^{2d}} H(x_1, x_2) \frac{F_{n-1}(x_1)}{\|F_{n-1}\|_{L^1(\mathbb{R}^d)}} \frac{F_{n-1}(x_2)}{\|F_{n-1}\|_{L^1(\mathbb{R}^d)}} dx_1 dx_2, \qquad (6.35)$$

where the function $H(x_1, x_2)$ takes the form

$$H(x_1, x_2) := \frac{1}{(1+\alpha)^{d/2}} \exp \left( -\frac{\alpha}{(1+\alpha)} \left| \frac{x_1+x_2}{2} \right|^2 \right), \quad (x_1, x_2) \in \mathbb{R}^d. \qquad (6.36)$$

In addition, recalling the explicit form of the eigen-pair $(\lambda_\alpha, \mathbf{F}_\alpha)$ in (1.12) implies

$$\lambda_\alpha = \int_{\mathbb{R}^{2d}} H(x_1, x_2) \mathbf{F}_\alpha(x_1) \mathbf{F}_\alpha(x_2) dx_1 dx_2. \qquad (6.37)$$

Therefore, taking the difference of (6.35) and (6.37) and noticing that $H$ in (6.36) is bounded yields

$$|\lambda_n - \lambda_\alpha| \lesssim \left\| \frac{F_n}{\|F_n\|_{L^1(\mathbb{R}^d)}} \otimes \frac{F_n}{\|F_n\|_{L^1(\mathbb{R}^d)}} - \mathbf{F}_\alpha \otimes \mathbf{F}_\alpha \right\|_{L^1(\mathbb{R}^{2d})}, \qquad (6.38)$$

for any $n \in \mathbb{N}$, where $P \otimes Q$ denotes the Kronecker product of two measures $P, Q \in \mathcal{M}(\mathbb{R}^d)$, namely

$$\int_{\mathbb{R}^{2d}} \varphi(x, y) (P \otimes Q)(dx, dy) = \int_{\mathbb{R}^d} \left( \int_{\mathbb{R}^d} \varphi(x, y) P(dx) \right) Q(dy),$$

for all $\varphi \in C_c(\mathbb{R}^d)$ (cf. Table 1). We do not have a direct convergence result of $F_n/\|F_n\|_{L^1(\mathbb{R}^d)}$ in $L^1$ norms, but we do have convergence of the Kullback-Leibler divergence thanks to STEP 1. By Pinsker's inequality, the latter metric controls the former, namely

$$\begin{aligned} & \left\| \frac{F_n}{\|F_n\|_{L^1(\mathbb{R}^d)}} \otimes \frac{F_n}{\|F_n\|_{L^1(\mathbb{R}^d)}} - \mathbf{F}_\alpha \otimes \mathbf{F}_\alpha \right\|_{L^1(\mathbb{R}^{2d})} \\ & \le \frac{1}{\sqrt{2}} \sqrt{\mathcal{D}_{\text{KL}} \left( \frac{F_n}{\|F_n\|_{L^1(\mathbb{R}^d)}} \otimes \frac{F_n}{\|F_n\|_{L^1(\mathbb{R}^d)}} \right) \| \mathbf{F}_\alpha \otimes \mathbf{F}_\alpha \|} \qquad (6.39) \\ & = \frac{1}{\sqrt{2}} \sqrt{2 \mathcal{D}_{\text{KL}} \left( \frac{F_n}{\|F_n\|_{L^1(\mathbb{R}^d)}} \otimes \| \mathbf{F}_\alpha \| \right)}, \end{aligned}$$
On the other hand, for $D_{n,2}$ we shall refine the argument in the proof of Corollary 6.5. Specifically, recall again that we can assume that $\bar{F}_0$ satisfies the hypothesis (6.5) without loss of generality (cf. proof of Corollary 6.5). Then, applying Lemma 6.3 in order to control $\nabla \log \mathcal{E}_n[\bar{F}_0]$ with fixed $\varepsilon \in \mathbb{R}_+^*$ yields

$$D_{n,2} \lesssim ((2 \mathbf{k}_\alpha)^2 + \varepsilon)^n \int_{\mathbb{R}^d} \left(1 + |x|^2 + e^{2C_\varepsilon(2 \mathbf{r}_\alpha + \varepsilon)^n |x|^2}\right) \frac{F_n(x)}{\|F_n\|_{L^1(\mathbb{R}^d)}} dx,$$

for some sufficiently large $C_\varepsilon \in \mathbb{R}_+^*$. Above we have used that $(2 \mathbf{k}_\alpha)^n$ is the decay rate that controls all the others in formula (6.6). Take $n \ge n_\varepsilon$ sufficiently large such that

$$\theta_\varepsilon := \sup_{n \ge n_\varepsilon} 2 C_\varepsilon (2 \mathbf{r}_\alpha + \varepsilon)^n < \frac{\alpha}{2}.$$

Then, Corollaries 4.6 and 4.8 imply that the above quadratic and exponential moments are uniformly bounded with respect to $n$ for $n \ge n_\varepsilon$. Therefore, we have

$$D_{n,2} \lesssim ((2 \mathbf{k}_\alpha)^2 + \varepsilon)^n. \qquad (6.34)$$

Finally, putting the estimates (6.33) and (6.34) into the split (6.32), and noticing that $\mathbf{r}_\alpha^2 < (2 \mathbf{k}_\alpha)^2$ thanks to the relation $\mathbf{r}_\alpha = 2 \mathbf{k}_\alpha^2$, ends this part of the proof.

• STEP 2: Convergence of the growth rates $\|F_n\|_{L^1(\mathbb{R}^d)}/\|F_{n-1}\|_{L^1(\mathbb{R}^d)}$.

For simplicity of notation we shall define the following sequence of coefficients:

$$\lambda_n := \frac{\|F_n\|_{L^1(\mathbb{R}^d)}}{\|F_{n-1}\|_{L^1(\mathbb{R}^d)}}, \quad n \in \mathbb{N}.$$

Our goal then reduces to studying the asymptotic behavior of $\{\lambda_n\}_{n \in \mathbb{N}}$ and obtaining quantitative convergence rates. By definition of operator $\mathcal{T}$ in (1.2) we obtain that

$$\lambda_n = \int_{\mathbb{R}^d} e^{-m(x)} \int_{\mathbb{R}^{2d}} G\left(x - \frac{x_1 + x_2}{2}\right) \frac{F_{n-1}(x_1)}{\|F_{n-1}\|_{L^1(\mathbb{R}^d)}} \frac{F_{n-1}(x_2)}{\|F_{n-1}\|_{L^1(\mathbb{R}^d)}} dx_1 dx_2 dx.$$

By direct computation of the integral with respect to $x$ as it was done in (4.11) we find that

$$\lambda_n = \int_{\mathbb{R}^{2d}} H(x_1, x_2) \frac{F_{n-1}(x_1)}{\|F_{n-1}\|_{L^1(\mathbb{R}^d)}} \frac{F_{n-1}(x_2)}{\|F_{n-1}\|_{L^1(\mathbb{R}^d)}} dx_1 dx_2, \qquad (6.35)$$

where the function $H(x_1, x_2)$ takes the form

$$H(x_1, x_2) := \frac{1}{(1+\alpha)^{d/2}} \exp \left( -\frac{\alpha}{(1+\alpha)} \left| \frac{x_1+x_2}{2} \right|^2 \right), \quad (x_1, x_2) \in \mathbb{R}^d. \qquad (6.36)$$

In addition, recalling the explicit form of the eigen-pair $(\lambda_\alpha, \mathbf{F}_\alpha)$ in (1.12) implies

$$\lambda_\alpha = \int_{\mathbb{R}^{2d}} H(x_1, x_2) \mathbf{F}_\alpha(x_1) \mathbf{F}_\alpha(x_2) dx_1 dx_2. \qquad (6.37)$$

Therefore, taking the difference of (6.35) and (6.37) and noticing that $H$ in (6.36) is bounded yields

$$|\lambda_n - \lambda_\alpha| \lesssim \left\| \frac{F_n}{\|F_n\|_{L^1(\mathbb{R}^d)}} \otimes \frac{F_n}{\|F_n\|_{L^1(\mathbb{R}^d)}} - \mathbf{F}_\alpha \otimes \mathbf{F}_\alpha \right\|_{L^1(\mathbb{R}^{2d})}, \qquad (6.38)$$

for any $n \in \mathbb{N}$, where $P \otimes Q$ denotes the Kronecker product of two measures $P, Q \in \mathcal{M}(\mathbb{R}^d)$, namely

$$\int_{\mathbb{R}^{2d}} \varphi(x, y) (P \otimes Q)(dx, dy) = \int_{\mathbb{R}^d} \left( \int_{\mathbb{R}^d} \varphi(x, y) P(dx) \right) Q(dy),$$

for all $\varphi \in C_c(\mathbb{R}^d)$ (cf. Table 1). We do not have a direct convergence result of $F_n/\|F_n\|_{L^1(\mathbb{R}^d)}$ in $L^1$ norms, but we do have convergence of the Kullback-Leibler divergence thanks to STEP 1. By Pinsker's inequality, the latter metric controls the former, namely

$$\begin{aligned} & \left\| \frac{F_n}{\|F_n\|_{L^1(\mathbb{R}^d)}} \otimes \frac{F_n}{\|F_n\|_{L^1(\mathbb{R}^d)}} - \mathbf{F}_\alpha \otimes \mathbf{F}_\alpha \right\|_{L^1(\mathbb{R}^{2d})} \\ & \le \frac{1}{\sqrt{2}} \sqrt{\mathcal{D}_{\text{KL}} \left( \frac{F_n}{\|F_n\|_{L^1(\mathbb{R}^d)}} \otimes \frac{F_n}{\|F_n\|_{L^1(\mathbb{R}^d)}} \right) \| \mathbf{F}_\alpha \otimes \mathbf{F}_\alpha \|} \qquad (6.39) \\ & = \frac{1}{\sqrt{2}} \sqrt{2 \mathcal{D}_{\text{KL}} \left( \frac{F_n}{\|F_n\|_{L^1(\mathbb{R}^d)}} \otimes \| \mathbf{F}_\alpha \| \right)}, \end{aligned}$$
where in last step we have used the tensorization property of the Kullback-Leibler divergence. The result then follows from (6.38)-(6.39) and the above explicit convergence rates of the normalized profiles $F_n/||F_n||_{L^1(\mathbb{R}^d)}$ in STEP 1. $\square$

## 7. NUMERICAL EXPERIMENTS

For our numerical simulation, we have restricted to one-dimensional traits (i.e., $d = 1$) and we have considered a step function of the following form as initial datum:

$$F_0 = \frac{1}{Z}(30 \mathbb{1}_{[-7,-3]} + 20 \mathbb{1}_{[7.5,12.5]} + 50 \mathbb{1}_{[30,40]} + 30 \mathbb{1}_{[52,5,57,5]}), \quad (7.1)$$

where $Z \in \mathbb{R}_+^*$ is a normalizing factor so that $F_0 \in \mathcal{P}(\mathbb{R})$. Let $\{F_n\}_{n \in \mathbb{N}}$ be the solution to the time-discrete problem (1.1) starting at $F_0$. Our numerical simulation is performed with PYTHON on the finite computation domain $[-15, 60]$ for the trait variable $x$, which contains the support of the previous initial datum $F_0$ and all the mass of each $F_n$ (except for a negligible Gaussian tail). In our simulation, we set $\Delta x = 0.001$ as the step for the discretization of our computational domain. In particular, with such a step we compute numerical integrals with respect to $x$ according to the left rectangle rule.

In the following, we explore numerically two different scenarios: weak selection and strong selection. Specifically, we obtain numerical approximations for the profiles $F_n$ in both regimes and we illustrate numerically the convergence of the growth rates $||F_n||_{L^1(\mathbb{R}^d)}/||F_{n-1}||_{L^1(\mathbb{R}^d)}$ towards the eigenvalue $\lambda_\alpha$, and the relaxation of the normalized profiles $F_n/||F_n||_{L^1(\mathbb{R}^d)}$ towards the eigenfunction $F_\alpha$. As a consequence, we derive numerical approximations of the convergence rates to be compared with the theoretical results in this paper. More specifically, we note that the theoretical convergence rates in Theorem 1.1 are sharp, except a mismatch for the rates of growth of mass, which was discussed previously in Remark 3.10 (see also Figure 6). Indeed, we actually attain numerically the same convergence rates as in Corollary 3.9, which were sharp for Gaussian initial data.

### 7.1. Weak selection.
In this part, we consider a small value $\alpha = 0.015$. This leads to the following numerical values of the features of the equilibrium:

$$\lambda_\alpha \approx 0.9857, \quad \sigma_\alpha^2 \approx 1.8897.$$

characterizing the eigenpair $(\lambda_\alpha, F_\alpha)$, according to (1.12) and (1.13). Note that since the selection parameter has been taken very small, then the eigenvalue and the variance of the eigenfunction are close to those at linkage equilibrium, namely, $\lambda_{\alpha=0} = 1$ and $\sigma_{\alpha=0}^2 = 2$ (see Remark 3.5 and Figure 4).

In Figure 8 we observe the relaxation of the normalized profiles $F_n/||F_n||_{L^1(\mathbb{R}^d)}$ towards the eigenfunction $F_\alpha$ along the time iterations $n = 0, 1, 2, 3, 4, 7, 150$. We remark on the strong contraction of the variance during the first few iterations leading to a well-identified Gaussian-shaped profile at time $n = 3$. At time $n = 7$ the variance of the profile is already close to that of the equilibrium $F_\alpha$, but its mean is still substantially shifted to the right. Thereafter, we observe the gradual motion of the profiles towards the left at a much lower speed. After 150 iterations, the approximate mean and variance of $F_n/||F_n||_{L^1(\mathbb{R}^d)}$ are given by $\approx 0.0474$ and $\approx 1.8897$ respectively, where the latter agrees with the above exact value $\sigma_\alpha^2$ up to 5 digits.

In Figure 9 we represent the convergence of the growth rates $||F_n||_{L^1(\mathbb{R}^d)}/||F_{n-1}||_{L^1(\mathbb{R}^d)}$ towards the eigenvalue $\lambda_\alpha$. After 150 iterations we obtain that the growth rate takes the approximate value $\approx 0.9857$, which again agrees with the exact value $\lambda_\alpha$ above up to 5 digits.

In Figure 10 we have represented the errors

$$\varepsilon_n^{\text{prof}} := \mathcal{D}_{\text{KL}} \left( \frac{F_n}{\|F_n\|_{L^1(\mathbb{R}^d)}} \middle| \mathbf{F}_\alpha \right), \quad \varepsilon_n^{\text{mass}} := \left| \frac{\|F_n\|_{L^1(\mathbb{R}^d)}}{\|F_{n-1}\|_{L^1(\mathbb{R}^d)}} - \lambda_\alpha \right|,$$

in semi-log plots so that the horizontal axis appears in the natural scale and represents each time iteration, and the vertical axis contains the logarithm of the errors. As observed, both plots reduce to essentially straight lines, suggesting exponential relaxation. Interestingly, the numerical rate of convergence in the Kullback-Leibler divergence measured in Figure 10a is approximately 0.9441. It coincides numerically with the rate of relaxation among the subclass of Gaussian solutions (see Corollary 3.9), namely, $\lambda_\alpha^4 \approx 0.9441$, which in turns is identical to the theoretical rate $(2k_\alpha)^2$ obtained in Theorem 1.1. These results are in perfect agreement with the other convergence results in Figure 10b: the numerical rate of convergence in the rate of growth of mass is approximately 0.9442, which is close to the one among the class of Gaussian solutions $\lambda_\alpha^4 \approx 0.9441$, in contrast with the theoretical upper bound obtained in Theorem 1.1, namely $2k_\alpha \approx 0.9717$. Moreover, the numerical rate of convergence of the variance of the normalized profiles $F_n/||F_n||_{L^1(\mathbb{R}^d)}$ is much faster: it is approximately $\approx 0.4721 < 0.5$, again close to
where in last step we have used the tensorization property of the Kullback-Leibler divergence. The result then follows from (6.38)-(6.39) and the above explicit convergence rates of the normalized profiles $F_n/||F_n||_{L^1(\mathbb{R}^d)}$ in STEP 1. $\square$

## 7. NUMERICAL EXPERIMENTS

For our numerical simulation, we have restricted to one-dimensional traits (i.e., $d=1$) and we have considered a step function of the following form as initial datum:

$$F_0 = \frac{1}{Z}(30 \mathbb{1}_{[-7,-3]} + 20 \mathbb{1}_{[7.5,12.5]} + 50 \mathbb{1}_{[30,40]} + 30 \mathbb{1}_{[52,5,57,5]}), \quad (7.1)$$

where $Z \in \mathbb{R}_+^*$ is a normalizing factor so that $F_0 \in \mathcal{P}(\mathbb{R})$. Let $\{F_n\}_{n \in \mathbb{N}}$ be the solution to the time-discrete problem (1.1) starting at $F_0$. Our numerical simulation is performed with PYTHON on the finite computation domain $[-15, 60]$ for the trait variable $x$, which contains the support of the previous initial datum $F_0$ and all the mass of each $F_n$ (except for a negligible Gaussian tail). In our simulation, we set $\Delta x = 0.001$ as the step for the discretization of our computational domain. In particular, with such a step we compute numerical integrals with respect to $x$ according to the left rectangle rule.

In the following, we explore numerically two different scenarios: weak selection and strong selection. Specifically, we obtain numerical approximations for the profiles $F_n$ in both regimes and we illustrate numerically the convergence of the growth rates $||F_n||_{L^1(\mathbb{R}^d)}/||F_{n-1}||_{L^1(\mathbb{R}^d)}$ towards the eigenvalue $\lambda_\alpha$, and the relaxation of the normalized profiles $F_n/||F_n||_{L^1(\mathbb{R}^d)}$ towards the eigenfunction $F_\alpha$. As a consequence, we derive numerical approximations of the convergence rates to be compared with the theoretical results in this paper. More specifically, we note that the theoretical convergence rates in Theorem 1.1 are sharp, except a mismatch for the rates of growth of mass, which was discussed previously in Remark 3.10 (see also Figure 6). Indeed, we actually attain numerically the same convergence rates as in Corollary 3.9, which were sharp for Gaussian initial data.

### 7.1. Weak selection.
In this part, we consider a small value $\alpha = 0.015$. This leads to the following numerical values of the features of the equilibrium:

$$\lambda_{\alpha} \approx 0.9857, \quad \sigma_{\alpha}^{2} \approx 1.8897.$$

characterizing the eigenpair $(\lambda_\alpha, F_\alpha)$, according to (1.12) and (1.13). Note that since the selection parameter has been taken very small, then the eigenvalue and the variance of the eigenfunction are close to those at linkage equilibrium, namely, $\lambda_{\alpha=0} = 1$ and $\sigma_{\alpha=0}^2 = 2$ (see Remark 3.5 and Figure 4).

In Figure 8 we observe the relaxation of the normalized profiles $F_n/||F_n||_{L^1(\mathbb{R}^d)}$ towards the eigenfunction $F_\alpha$ along the time iterations $n = 0, 1, 2, 3, 4, 7, 150$. We remark on the strong contraction of the variance during the first few iterations leading to a well-identified Gaussian-shaped profile at time $n = 3$. At time $n = 7$ the variance of the profile is already close to that of the equilibrium $F_\alpha$, but its mean is still substantially shifted to the right. Thereafter, we observe the gradual motion of the profiles towards the left at a much lower speed. After 150 iterations, the approximate mean and variance of $F_n/||F_n||_{L^1(\mathbb{R}^d)}$ are given by $\approx 0.0474$ and $\approx 1.8897$ respectively, where the latter agrees with the above exact value $\sigma_\alpha^2$ up to 5 digits.

In Figure 9 we represent the convergence of the growth rates $||F_n||_{L^1(\mathbb{R}^d)}/||F_{n-1}||_{L^1(\mathbb{R}^d)}$ towards the eigenvalue $\lambda_\alpha$. After 150 iterations we obtain that the growth rate takes the approximate value $\approx 0.9857$, which again agrees with the exact value $\lambda_\alpha$ above up to 5 digits.

In Figure 10 we have represented the errors

$$\varepsilon_n^{\text{prof}} := \mathcal{D}_{\text{KL}} \left( \frac{F_n}{||F_n||_{L^1(\mathbb{R}^d)}} \middle| \mathbf{F}_\alpha \right), \quad \varepsilon_n^{\text{mass}} := \left| \frac{||F_n||_{L^1(\mathbb{R}^d)}}{||F_{n-1}||_{L^1(\mathbb{R}^d)}} - \lambda_\alpha \right|,$$

in semi-log plots so that the horizontal axis appears in the natural scale and represents each time iteration, and the vertical axis contains the logarithm of the errors. As observed, both plots reduce to essentially straight lines, suggesting exponential relaxation. Interestingly, the numerical rate of convergence in the Kullback-Leibler divergence measured in Figure 10a is approximately 0.9441. It coincides numerically with the rate of relaxation among the subclass of Gaussian solutions (see Corollary 3.9), namely, $\lambda_\alpha^4 \approx 0.9441$, which in turns is identical to the theoretical rate $(2k_\alpha)^2$ obtained in Theorem 1.1. These results are in perfect agreement with the other convergence results in Figure 10b: the numerical rate of convergence in the rate of growth of mass is approximately 0.9442, which is close to the one among the class of Gaussian solutions $\lambda_\alpha^4 \approx 0.9441$, in contrast with the theoretical upper bound obtained in Theorem 1.1, namely $2k_\alpha \approx 0.9717$. Moreover, the numerical rate of convergence of the variance of the normalized profiles $F_n/||F_n||_{L^1(\mathbb{R}^d)}$ is much faster: it is approximately $\approx 0.4721 < 0.5$, again close to
<figure><img src="image_11.png" /><figcaption>FIGURE 8. Relaxation of the normalized profiles $F_n/||F_n||_{L^1(\mathbb{R}^d)}$ towards the eigenfunction $F_\alpha$ along the time iterations $n = 0, 1, 2, 3, 4, 7$ and $150$ for a step function (7.1) as initial datum and weak selection parameter $\alpha = 0.015$. The vertical dotted line represents the location of the mean of the equilibrium profile $F_\alpha$.</figcaption></figure>

the expected value for Gaussian solutions, namely, $r_\alpha \approx 0.472$. These numerical results illustrate the two-step process arising in the relaxation dynamics: convergence towards a Gaussian profile occurs much faster than relaxation of the mean towards the origin due to (weak) selection. Alternatively speaking, the equilibrium variance builds up much faster than the center of the distribution gets to the origin.

7.2. **Strong selection.** In this part, we consider a larger value of $\alpha$. We opted for $\alpha = 0.4$. Indeed, if $\alpha$ is taken too large then there is no visual difference between the solution and the equilibrium configuration after one single iteration. This time, the numerical values of the features of the equilibrium are:

$$ \lambda_{\alpha} \approx 0.7944, \quad \sigma_{\alpha}^2 \approx 0.9221. $$
<figure><img src="image_11.png" /><figcaption>FIGURE 8. Relaxation of the normalized profiles $F_n/||F_n||_{L^1(\mathbb{R}^d)}$ towards the eigenfunction $F_\alpha$ along the time iterations $n = 0, 1, 2, 3, 4, 7$ and $150$ for a step function (7.1) as initial datum and weak selection parameter $\alpha = 0.015$. The vertical dotted line represents the location of the mean of the equilibrium profile $F_\alpha$.</figcaption></figure>

the expected value for Gaussian solutions, namely, $r_\alpha \approx 0.472$. These numerical results illustrate the two-step process arising in the relaxation dynamics: convergence towards a Gaussian profile occurs much faster than relaxation of the mean towards the origin due to (weak) selection. Alternatively speaking, the equilibrium variance builds up much faster than the center of the distribution gets to the origin.

7.2. **Strong selection.** In this part, we consider a larger value of $\alpha$. We opted for $\alpha = 0.4$. Indeed, if $\alpha$ is taken too large then there is no visual difference between the solution and the equilibrium configuration after one single iteration. This time, the numerical values of the features of the equilibrium are:

$$ \lambda_{\alpha} \approx 0.7944, \quad \sigma_{\alpha}^2 \approx 0.9221. $$
<figure><img src="image_12.png" /><figcaption>FIGURE 9. Relaxation of the growth rates $\|F_n\|_{L^1(\mathbb{R}^d)}/\|F_{n-1}\|_{L^1(\mathbb{R}^d)}$ towards the eigenvalue $\lambda_\alpha$ along time iterations $0 \le n \le 150$ for a step function (7.1) as initial datum and weak selection parameter $\alpha = 0.015$.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>FIGURE 10. Numerical computation of the rates of convergence of the normalized profiles and the growth rates: (a) semi-log plot of the errors $\varepsilon_n^{\text{prof}} := \mathcal{D}_{\text{KL}} (F_n/\|F_n\|_{L^1} \|\mathbf{F}_\alpha)$, and (b) semi-log plot of the errors $\varepsilon_n^{\text{mass}} := \|F_n\|_{L^1}/\|F_{n-1}\|_{L^1} - \lambda_\alpha$.</figcaption></figure>

The latter substantially differs from the value at linkage equilibrium (i.e., $\sigma_{\alpha=0}^2 = 2$). In Figure 11a we note that the solution resembles a Gaussian distribution after a couple of iterations. After 15 iterations we obtain that the normalized profile $F_n/\|F_n\|_{L^1(\mathbb{R}^d)}$ has approximate mean and variance respectively given by $\approx 0.0017$ and $\approx 0.9221$. In Figure 11b we observe the convergence of the growth rates $\|F_n\|_{L^1(\mathbb{R}^d)}/\|F_{n-1}\|_{L^1(\mathbb{R}^d)}$ towards the eigenvalue $\lambda_\alpha$. After 15 iterations the growth rate becomes $\approx 0.7944$, which again agrees with the exact value $\lambda_\alpha$ above up to 5 digits.

Similar computations as in the previous Figure 10 allow finding numerical approximations for the rate of convergence of the normalized profiles and the growth rates. Specifically, we obtain an approximation 0.3966 for the rate of convergence of the Kullback-Leibler divergence. Once again, such a value is close to the rate of relaxation among the subclass of Gaussian solutions, namely, $\lambda_\alpha^4 \approx 0.3983$, which in turns agrees with the theoretical rate $(2k_\alpha)^2$ obtained in Theorem 1.1. Similarly, the numerical rate of convergence of the growth rates is approximately 0.3901, which is close to the one among the class of Gaussian solutions $\lambda_\alpha^4 \approx 0.3983$, in contrast with the upper bound $2k_\alpha \approx 0.6311$ in Theorem 1.1.

8. CONCLUSIONS AND PERSPECTIVES

In this paper, we have proven asynchronous exponential growth in a quantitative genetics model for the evolution of the distribution of traits in a population governed by sexual reproduction and multiplicative effect of selection. Our model assumes time-discrete non-overlapping generations, which rule out an eventual mixing with previous generations of ancestors. In addition, our non-linear sexual reproduction operator is set in agreement with Fisher's infinitesimal model, and we have chosen selection to act on
<figure><img src="image_12.png" /><figcaption>FIGURE 9. Relaxation of the growth rates $\|F_n\|_{L^1(\mathbb{R}^d)}/\|F_{n-1}\|_{L^1(\mathbb{R}^d)}$ towards the eigenvalue $\lambda_\alpha$ along time iterations $0 \le n \le 150$ for a step function (7.1) as initial datum and weak selection parameter $\alpha = 0.015$.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>FIGURE 10. Numerical computation of the rates of convergence of the normalized profiles and the growth rates: (a) semi-log plot of the errors $\varepsilon_n^{\text{prof}} := \mathcal{D}_{\text{KL}} (F_n/\|F_n\|_{L^1} \|\mathbf{F}_\alpha)$, and (b) semi-log plot of the errors $\varepsilon_n^{\text{mass}} := \|F_n\|_{L^1}/\|F_{n-1}\|_{L^1} - \lambda_\alpha$.</figcaption></figure>

The latter substantially differs from the value at linkage equilibrium (i.e., $\sigma^2_{\alpha=0} = 2$). In Figure 11a we note that the solution resembles a Gaussian distribution after a couple of iterations. After 15 iterations we obtain that the normalized profile $F_n/\|F_n\|_{L^1(\mathbb{R}^d)}$ has approximate mean and variance respectively given by $\approx 0.0017$ and $\approx 0.9221$. In Figure 11b we observe the convergence of the growth rates $\|F_n\|_{L^1(\mathbb{R}^d)}/\|F_{n-1}\|_{L^1(\mathbb{R}^d)}$ towards the eigenvalue $\lambda_\alpha$. After 15 iterations the growth rate becomes $\approx 0.7944$, which again agrees with the exact value $\lambda_\alpha$ above up to 5 digits.

Similar computations as in the previous Figure 10 allow finding numerical approximations for the rate of convergence of the normalized profiles and the growth rates. Specifically, we obtain an approximation 0.3966 for the rate of convergence of the Kullback-Leibler divergence. Once again, such a value is close to the rate of relaxation among the subclass of Gaussian solutions, namely, $\lambda_\alpha^4 \approx 0.3983$, which in turns agrees with the theoretical rate $(2k_\alpha)^2$ obtained in Theorem 1.1. Similarly, the numerical rate of convergence of the growth rates is approximately 0.3901, which is close to the one among the class of Gaussian solutions $\lambda_\alpha^4 \approx 0.3983$, in contrast with the upper bound $2k_\alpha \approx 0.6311$ in Theorem 1.1.

## 8. CONCLUSIONS AND PERSPECTIVES

In this paper, we have proven asynchronous exponential growth in a quantitative genetics model for the evolution of the distribution of traits in a population governed by sexual reproduction and multiplicative effect of selection. Our model assumes time-discrete non-overlapping generations, which rule out an eventual mixing with previous generations of ancestors. In addition, our non-linear sexual reproduction operator is set in agreement with Fisher's infinitesimal model, and we have chosen selection to act on