# Ergodicity of the Fisher infinitesimal model with quadratic selection

Vincent Calvez, Thomas Lepoutre, David Poyato

► To cite this version:

Vincent Calvez, Thomas Lepoutre, David Poyato. Ergodicity of the Fisher infinitesimal model with quadratic selection. Nonlinear Analysis: Theory, Methods and Applications, 2024, 238, pp.113392. 10.1016/j.na.2023.113392. hal-03274734v2

HAL Id: hal-03274734

https://hal.science/hal-03274734v2

Submitted on 15 Jan 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

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

Distributed under a Creative Commons Attribution 4.0 International License
# Ergodicity of the Fisher infinitesimal model with quadratic selection

Vincent Calvez, Thomas Lepoutre, David Poyato

► To cite this version:

Vincent Calvez, Thomas Lepoutre, David Poyato. Ergodicity of the Fisher infinitesimal model with quadratic selection. Nonlinear Analysis: Theory, Methods and Applications, 2024, 238, pp.113392. 10.1016/j.na.2023.113392. hal-03274734v2

HAL Id: hal-03274734

https://hal.science/hal-03274734v2

Submitted on 15 Jan 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

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

Distributed under a Creative Commons Attribution 4.0 International License
ERGODICITY OF THE FISHER INFINITESIMAL MODEL WITH QUADRATIC SELECTION

VINCENT CALVEZ, THOMAS LEPOUTRE, AND DAVID POYATO

ABSTRACT. We study the convergence towards a unique equilibrium distribution of the solutions to a time-discrete model with non-overlapping generations arising in quantitative genetics. The model describes the dynamics of a phenotypic distribution with respect to a multi-dimensional trait, which is shaped by selection and Fisher's infinitesimal model of sexual reproduction. We extend some previous works devoted to the time-continuous analogues, that followed a perturbative approach in the regime of weak selection, by exploiting the contractivity of the infinitesimal model operator in the Wasserstein metric. Here, we tackle the case of quadratic selection by a global approach. We establish uniqueness of the equilibrium distribution and exponential convergence of the renormalized profile. Our technique relies on an accurate control of the propagation of information across the large binary trees of ancestors (the pedigree chart), and reveals an ergodicity property, meaning that the shape of the initial datum is quickly forgotten across generations. We combine this information with appropriate estimates for the emergence of Gaussian tails and propagation of quadratic and exponential moments to derive quantitative convergence rates. Our result can be interpreted as a generalization of the Krein-Rutman theorem in a genuinely non-linear, and non-monotone setting.

# 1. INTRODUCTION

Fisher's infinitesimal model (also known as the *polygenic model*) is a widely used statistical model in quantitative genetics initially proposed by R. FISHER [15]. It assumes that the genetic component of a quantitative phenotypical trait is affected by an infinite number of loci with infinitesimal and additive allelic effects and claims that the genetic component of descendants' traits is normally distributed around the mean value of parents' traits, with a constant (genetic) variance across generations. This model allowed reconciling Mendelian inheritance and the continuous trait variations documented by F. GALTON via a Central Limit Theorem. More specifically, by taking limits when the number of underlying loci tends to infinity on a model with Mendelian inheritance, N. BARTON, A. ETHERIDGE and A. VÉBER [5] recently proved rigorously the validity of Fisher's infinitesimal model under various evolutionary processes (e.g., natural selection).

In this paper we study a time-discrete evolution problem for the distribution of a phenotypical trait $x \in \mathbb{R}^d$ in a population undergoing sexual reproduction and the effect of natural selection. Specifically, starting at any initial configuration $F_0 \in \mathcal{M}_+(\mathbb{R}^d)$ of trait distribution, we analyze the long term dynamics of trait distributions $\{F_n\}_{n \in \mathbb{N}}$ across successive generations $n \in \mathbb{N}$, which solves the following recursion

$$F_n = \mathcal{T}[F_{n-1}], \qquad (1.1)$$

for any $n \in \mathbb{N}$. The operator $\mathcal{T}$ encodes the balanced effect of sexual reproduction in the population (under the infinitesimal model) and natural selection. Specifically, $\mathcal{T}$ is defined by

$$\mathcal{T}[F] := e^{-m}\mathcal{B}[F], \qquad (1.2)$$

for any trait distribution $F \in \mathcal{M}_+(\mathbb{R}^d)$. On the one hand, $m = m(x) \ge 0$ is called the selection function and represents the mortality effect of a trait-dependent natural selection on the population, so that $e^{-m}$ stands for the survival probability in the next generation. On the other hand, the operator $\mathcal{B}$ is chosen

<footnote>Date: January 11, 2024.</footnote>

<footnote n="2020">2020 Mathematics Subject Classification. 35B40; 35P30; 35Q92; 47G20; 92D15.</footnote>

<footnote>Key words and phrases. Integro-differential equations, asymptotic behavior, nonlinear spectral theory, quantitative genetics.</footnote>

<footnote>Acknowledgment. VC and DP have received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation program (grant agreement No 865711). DP has also received funding from the European Union's Horizon Europe research and innovation program under the Marie Skłodowska-Curie grant agreement No 101064402, and partially from the State Research Agency (SRA) of the Spanish Ministry of Science and Innovation and European Regional Development Fund (ERDF), project PID2022-137228OB-I00, and by Modeling Nature Research Unit, project QUAL21-011.</footnote>
ERGODICITY OF THE FISHER INFINITESIMAL MODEL WITH QUADRATIC SELECTION

VINCENT CALVEZ, THOMAS LEPOUTRE, AND DAVID POYATO

ABSTRACT. We study the convergence towards a unique equilibrium distribution of the solutions to a time-discrete model with non-overlapping generations arising in quantitative genetics. The model describes the dynamics of a phenotypic distribution with respect to a multi-dimensional trait, which is shaped by selection and Fisher's infinitesimal model of sexual reproduction. We extend some previous works devoted to the time-continuous analogues, that followed a perturbative approach in the regime of weak selection, by exploiting the contractivity of the infinitesimal model operator in the Wasserstein metric. Here, we tackle the case of quadratic selection by a global approach. We establish uniqueness of the equilibrium distribution and exponential convergence of the renormalized profile. Our technique relies on an accurate control of the propagation of information across the large binary trees of ancestors (the pedigree chart), and reveals an ergodicity property, meaning that the shape of the initial datum is quickly forgotten across generations. We combine this information with appropriate estimates for the emergence of Gaussian tails and propagation of quadratic and exponential moments to derive quantitative convergence rates. Our result can be interpreted as a generalization of the Krein-Rutman theorem in a genuinely non-linear, and non-monotone setting.

# 1. INTRODUCTION

Fisher's infinitesimal model (also known as the *polygenic model*) is a widely used statistical model in quantitative genetics initially proposed by R. FISHER [15]. It assumes that the genetic component of a quantitative phenotypical trait is affected by an infinite number of loci with infinitesimal and additive allelic effects and claims that the genetic component of descendants' traits is normally distributed around the mean value of parents' traits, with a constant (genetic) variance across generations. This model allowed reconciling Mendelian inheritance and the continuous trait variations documented by F. GALTON via a Central Limit Theorem. More specifically, by taking limits when the number of underlying loci tends to infinity on a model with Mendelian inheritance, N. BARTON, A. ETHERIDGE and A. VÉBER [5] recently proved rigorously the validity of Fisher's infinitesimal model under various evolutionary processes (e.g., natural selection).

In this paper we study a time-discrete evolution problem for the distribution of a phenotypical trait $x \in \mathbb{R}^d$ in a population undergoing sexual reproduction and the effect of natural selection. Specifically, starting at any initial configuration $F_0 \in \mathcal{M}_+(\mathbb{R}^d)$ of trait distribution, we analyze the long term dynamics of trait distributions $\{F_n\}_{n \in \mathbb{N}}$ across successive generations $n \in \mathbb{N}$, which solves the following recursion

$$F_n = \mathcal{T}[F_{n-1}], \qquad (1.1)$$

for any $n \in \mathbb{N}$. The operator $\mathcal{T}$ encodes the balanced effect of sexual reproduction in the population (under the infinitesimal model) and natural selection. Specifically, $\mathcal{T}$ is defined by

$$\mathcal{T}[F] := e^{-m}\mathcal{B}[F], \qquad (1.2)$$

for any trait distribution $F \in \mathcal{M}_+(\mathbb{R}^d)$. On the one hand, $m = m(x) \ge 0$ is called the selection function and represents the mortality effect of a trait-dependent natural selection on the population, so that $e^{-m}$ stands for the survival probability in the next generation. On the other hand, the operator $\mathcal{B}$ is chosen

<footnote>Date: January 11, 2024.</footnote>

<footnote n="2020">2020 Mathematics Subject Classification. 35B40; 35P30; 35Q92; 47G20; 92D15.</footnote>

<footnote>Key words and phrases. Integro-differential equations, asymptotic behavior, nonlinear spectral theory, quantitative genetics.</footnote>

<footnote>Acknowledgment. VC and DP have received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation program (grant agreement No 865711). DP has also received funding from the European Union's Horizon Europe research and innovation program under the Marie Skłodowska-Curie grant agreement No 101064402, and partially from the State Research Agency (SRA) of the Spanish Ministry of Science and Innovation and European Regional Development Fund (ERDF), project PID2022-137228OB-I00, and by Modeling Nature Research Unit, project QUAL21-011.</footnote>
to be Fisher's infinitesimal operator and it takes the form

$$ \mathcal{B}[F](x) := \int_{\mathbb{R}^{2d}} G \left( x - \frac{x_1 + x_2}{2} \right) \frac{F(x_1)F(x_2)}{\int_{\mathbb{R}^d} F(x') dx'} dx_1 dx_2, \quad x \in \mathbb{R}^d, \qquad (1.3) $$

for any trait distribution $F \in \mathcal{M}_+(\mathbb{R}^d)$. Here $G = G(x)$ is a probability density (the mixing kernel) and the factor $G(x - \frac{x_1+x_2}{2})$ represents the transition probability that two given individuals with trait values $x_1, x_2 \in \mathbb{R}^d$ will mate and yield a descendant with trait value $x \in \mathbb{R}^d$. In other words the resulting trait distributes around the mean value $\frac{x_1+x_2}{2}$ of parents' trait with law $G$. By definition, $\mathcal{B}[F]$ quantifies the number of births after all possible matches of any couple of individuals according to the trait distribution $F$. Altogether, $\mathcal{T}[F]$ quantifies the amount of offspring of a population distributed according to $F$ having resisted the effect of selection.

The above sexual reproduction operator $\mathcal{B}$ has recently pulled the attention of both the applied and more theoretical communities, cf. [5, 10, 18, 19, 20, 21]. In this paper we shall restrict to Gaussian mixing kernel and quadratic selection function, i.e.,

$$ G(x) := \frac{1}{(2\pi)^{d/2}} e^{-\frac{|x|^2}{2}}, \quad x \in \mathbb{R}^d, \qquad (1.4) $$

$$ m(x) := \frac{\alpha}{2} |x|^2, \quad x \in \mathbb{R}^d, \qquad (1.5) $$

where $\alpha \in \mathbb{R}_+$ is a fixed parameter. Thereby, the trait of offsprings is normally distributed around the mean value of the trait of parents by assumption (1.4), thus reducing to the standard infinitesimal model when the assumptions of the Central Limit Theorem are met [5]. For simplicity, we set a Gaussian $G$ with unit variance, but any value of the genetic variance could also be considered (see nondimensionalization in Appendix A).

Before introducing our results, we shall relate the previous time-discrete problem (1.1) to analogous time-continuous quantitative genetics models of evolutionary dynamics that have been studied in the literature. Meanwhile, we will anticipate the major difficulties that can be faced when analyzing the long-time dynamics of (1.1). To this end, we consider, the following type of integro-differential equations for the evolutionary dynamics of a trait distribution $f(t, x)$:

$$ \left\{ \begin{array}{ll} \partial_t f = -m(x)f + \mathcal{R}[f], & t \ge 0, x \in \mathbb{R}^d, \\ f(0,x) = f_0(x), & x \in \mathbb{R}^d, \end{array} \right. \qquad (1.6) $$

where, $m = m(x)$ is again the mortality rate and $\mathcal{R} = \mathcal{R}[f]$ is the reproduction operator. Hence, $\mathcal{R}[f](x)$ determines the amount of births with trait value $x \in \mathbb{R}^d$ per unit of time. As for (1.1), the resulting dynamics of the population becomes a consequence of a balance between selection encoded by the trait-dependent mortality and the diversity generated by the growth term $\mathcal{R}$ across generations.

Many studies consider a linear reproduction operator $\mathcal{R}$, associated with a probability density $K = K(x)$ characterizing the mutational effects at birth, of the form

$$ \mathcal{R}[f](x) := \int_{\mathbb{R}^d} K(x - y)f(t, y) dy, \quad x \in \mathbb{R}^d. \qquad (1.7) $$

The factor $K(x - y)$ determines the probability that an individual with trait value $y \in \mathbb{R}^d$ produces a descendant with trait value $x \in \mathbb{R}^d$ (possibly deviating from $y$). This class of linear reproduction operators is well-suited for an asexual mode of reproduction. This includes parabolic equations in the limit of small variance of $K$. In particular, we refer to a series of works about the long-time asymptotics in the regime of small variance, initiated in [3, 4, 14], including an additional density-dependent competition term that makes the analysis non standard (see e.g. [11] and references therein for the well-posedness of the constrained Hamilton-Jacobi equation derived in the limit).

Recently, inspired by the infinitesimal model, the case of sexual reproduction has been addressed by invoking the preceding nonlinear version $\mathcal{B}$ in (1.3) as reproduction operator $\mathcal{R} = \mathcal{B}$. Several asymptotic regimes have been addressed: large reproduction rate [18, 20, 21], small variance asymptotics [10, 19]. In the latter case, the limiting problem keeps the non-local nature of the problem, being of a finite-difference type, rather than a Hamilton-Jacobi PDE. Before we continue the discussion about the state-of-the-art, let us emphasize that there is no restriction on the parameter $\alpha$ in the present work.

Since both $f \mapsto mf$ and $f \mapsto \mathcal{B}[f]$ are 1-homogeneous operators, we can seek for special steady solutions of (1.3) through the following ansatz:

$$ f(t, x) = e^{\lambda t} F(x), \quad (t, x) \in \mathbb{R}_+ \times \mathbb{R}^d, \qquad (1.8) $$
to be Fisher's infinitesimal operator and it takes the form

$$ \mathcal{B}[F](x) := \int_{\mathbb{R}^{2d}} G\left(x - \frac{x_1 + x_2}{2}\right) \frac{F(x_1)F(x_2)}{\int_{\mathbb{R}^d} F(x') dx'} dx_1 dx_2, \quad x \in \mathbb{R}^d, \qquad (1.3) $$

for any trait distribution $F \in \mathcal{M}_+(\mathbb{R}^d)$. Here $G = G(x)$ is a probability density (the mixing kernel) and the factor $G(x - \frac{x_1+x_2}{2})$ represents the transition probability that two given individuals with trait values $x_1, x_2 \in \mathbb{R}^d$ will mate and yield a descendant with trait value $x \in \mathbb{R}^d$. In other words the resulting trait distributes around the mean value $\frac{x_1+x_2}{2}$ of parents' trait with law $G$. By definition, $\mathcal{B}[F]$ quantifies the number of births after all possible matches of any couple of individuals according to the trait distribution $F$. Altogether, $\mathcal{T}[F]$ quantifies the amount of offspring of a population distributed according to $F$ having resisted the effect of selection.

The above sexual reproduction operator $\mathcal{B}$ has recently pulled the attention of both the applied and more theoretical communities, cf. [5, 10, 18, 19, 20, 21]. In this paper we shall restrict to Gaussian mixing kernel and quadratic selection function, i.e.,

$$ G(x) := \frac{1}{(2\pi)^{d/2}} e^{-\frac{|x|^2}{2}}, \quad x \in \mathbb{R}^d, \qquad (1.4) $$

$$ m(x) := \frac{\alpha}{2} |x|^2, \quad x \in \mathbb{R}^d, \qquad (1.5) $$

where $\alpha \in \mathbb{R}_+$ is a fixed parameter. Thereby, the trait of offsprings is normally distributed around the mean value of the trait of parents by assumption (1.4), thus reducing to the standard infinitesimal model when the assumptions of the Central Limit Theorem are met [5]. For simplicity, we set a Gaussian $G$ with unit variance, but any value of the genetic variance could also be considered (see nondimensionalization in Appendix A).

Before introducing our results, we shall relate the previous time-discrete problem (1.1) to analogous time-continuous quantitative genetics models of evolutionary dynamics that have been studied in the literature. Meanwhile, we will anticipate the major difficulties that can be faced when analyzing the long-time dynamics of (1.1). To this end, we consider, the following type of integro-differential equations for the evolutionary dynamics of a trait distribution $f(t, x)$:

$$ \left\{ \begin{array}{ll} \partial_t f = -m(x)f + \mathcal{R}[f], & t \ge 0, x \in \mathbb{R}^d, \\ f(0, x) = f_0(x), & x \in \mathbb{R}^d, \end{array} \right. \qquad (1.6) $$

where, $m = m(x)$ is again the mortality rate and $\mathcal{R} = \mathcal{R}[f]$ is the reproduction operator. Hence, $\mathcal{R}[f](x)$ determines the amount of births with trait value $x \in \mathbb{R}^d$ per unit of time. As for (1.1), the resulting dynamics of the population becomes a consequence of a balance between selection encoded by the trait-dependent mortality and the diversity generated by the growth term $\mathcal{R}$ across generations.

Many studies consider a linear reproduction operator $\mathcal{R}$, associated with a probability density $K = K(x)$ characterizing the mutational effects at birth, of the form

$$ \mathcal{R}[f](x) := \int_{\mathbb{R}^d} K(x - y)f(t, y) dy, \quad x \in \mathbb{R}^d. \qquad (1.7) $$

The factor $K(x - y)$ determines the probability that an individual with trait value $y \in \mathbb{R}^d$ produces a descendant with trait value $x \in \mathbb{R}^d$ (possibly deviating from $y$). This class of linear reproduction operators is well-suited for an asexual mode of reproduction. This includes parabolic equations in the limit of small variance of $K$. In particular, we refer to a series of works about the long-time asymptotics in the regime of small variance, initiated in [3, 4, 14], including an additional density-dependent competition term that makes the analysis non standard (see e.g. [11] and references therein for the well-posedness of the constrained Hamilton-Jacobi equation derived in the limit).

Recently, inspired by the infinitesimal model, the case of sexual reproduction has been addressed by invoking the preceding nonlinear version $\mathcal{B}$ in (1.3) as reproduction operator $\mathcal{R} = \mathcal{B}$. Several asymptotic regimes have been addressed: large reproduction rate [18, 20, 21], small variance asymptotics [10, 19]. In the latter case, the limiting problem keeps the non-local nature of the problem, being of a finite-difference type, rather than a Hamilton-Jacobi PDE. Before we continue the discussion about the state-of-the-art, let us emphasize that there is no restriction on the parameter $\alpha$ in the present work.

Since both $f \mapsto mf$ and $f \mapsto \mathcal{B}[f]$ are 1-homogeneous operators, we can seek for special steady solutions of (1.3) through the following ansatz:

$$ f(t, x) = e^{\lambda t} F(x), \quad (t, x) \in \mathbb{R}_+ \times \mathbb{R}^d, \qquad (1.8) $$
The parameter $\lambda \in \mathbb{R}$ represents the rate at which the number of individuals grows (if $\lambda \ge 0$) or decreases (if $\lambda \le 0$), and $F = F(x) \ge 0$ is an unknown probability density. By imposing such an ansatz on (1.6)-(1.3), the following generalized spectral problem arises for the couple $(\lambda, F)$:

$$ \left\{ \begin{array}{ll} \lambda F(x) + m(x) F(x) = \mathcal{B}[F](x), & x \in \mathbb{R}^d, \\ \int_{\mathbb{R}^d} F(x') dx' = 1. & \end{array} \right. \qquad (1.9) $$

Note that the operator $\mathcal{B}$ is genuinely non-linear so that methods based on the Krein-Rutman theory or maximum principles cannot be applied straightforwardly, see [6, 13] and the references therein for the linear case. Further, usual extensions of the Krein-Rutman theory to 1-homogeneous operators [16] cannot be applied neither because $\mathcal{B}$ is not monotone. To date, the main strategy behind the existence of solutions of (1.9) relies on a suitable application of Schauder fixed-point theorem to the operator $F \mapsto (\lambda + m)^{-1}\mathcal{B}[F]$ over an appropriate cone of $L^1(\mathbb{R}^d)$ that is conserved by the nonlinear operator, see [8]. However, uniqueness cannot be achieved by this method. Moreover, it has been proven in [10, Corollary 1.5] that several equilibrium states $(\lambda, F)$ can co-exist in the presence of multiple local minima of $m$ (provided the variance of $G$ is small enough). That is, the generalized eigenproblem (1.9) does not admit a unique positive eigenfunction, in contrast with general conclusions of the Krein-Rutman theory.

Recently, G. RAOUL addressed the long-time dynamics of (1.6) in 1D, with $\mathcal{R} = \mathcal{B}$ and a trait dependent fecundity rate [21]. He obtained local uniqueness and exponential relaxation under the assumption of weak and localized (compactly supported) selection effects. To this end, he controlled locally in space the Wasserstein distance between the solution and the stationary state, using the uniform contraction property of $\mathcal{B}$ in the space of probability measures sharing the same center of mass. Unfortunately, the Wasserstein metric is not fully compatible with multiplicative operators, such as trait-dependent fecundity (see the discussion in [21, Section 3.4] and Section 2 below). This can be circumvented under the additional assumption that the trait density is locally uniformly bounded below, following [7]. Obviously, this cannot hold globally for integrable densities, hence G. RAOUL developed estimates of the distribution's tail to complete the contraction estimates. Also, a lot of attention has to be paid to the dynamics of the center of the distribution which is essentially driven by selection. Indeed, in the case of flat selection ($m \equiv 0$), the problem is invariant by translation, so that local uniqueness cannot hold. To conclude this discussion, let us emphasize that global uniqueness and the asymptotic behavior of generic solutions to the evolution problem (1.6)-(1.3) is still open. Below, we provide a first result in this direction, for the time-discrete problem (1.1), though.

We remark that the time-discrete version (1.1) that we propose in this paper can be partially regarded as a discretization in time of the above time-continuous problem (1.6)-(1.3), with non-overlapping generations (see Appendix A for further details). As for the time-continuous problem, we could seek special solutions to the time-discrete problem (1.1) in the following form

$$ F_n(x) = \lambda^n F(x), \quad (n, x) \in \mathbb{N} \times \mathbb{R}^d. \qquad (1.10) $$

Again, $\lambda \in \mathbb{R}_+^*$ is the rate of growth (if $\lambda \ge 1$) or decrease (if $\lambda \le 1$) of individuals, and $F = F(x)$ is an unknown probability density. This yields the following generalized eigenproblem for the couple $(\lambda, F)$:

$$ \left\{ \begin{array}{ll} \lambda F(x) = \mathcal{T}[F](x), & x \in \mathbb{R}^d, \\ \int_{\mathbb{R}^d} F(x') dx' = 1. & \end{array} \right. \qquad (1.11) $$

In this paper we aim to address the following questions:

(Q1) Does the eigenproblem (1.11) have a unique solution $(\lambda_\alpha, F_\alpha)$, with $\lambda_\alpha \in \mathbb{R}$ and $F_\alpha$ being a probability measure, for each $\alpha \in \mathbb{R}_+^*$?

(Q2) Consider any generic initial datum $F_0 \in \mathcal{M}_+(\mathbb{R}^d)$ and its associated solution $\{F_n\}_{n \in \mathbb{N}}$ of the time-discrete problem (1.1). Do the renormalized profiles $F_n/\|F_n\|_{L^1(\mathbb{R}^d)}$ converge to the unique steady profile $F_\alpha$ solving (1.11) when $n \to \infty$?

We shall prove that the answer to both questions is affirmative. It stands to reason that similar existence and local uniqueness results like in [8, 10] could be extended to the new eigenproblem (1.11) by applying Schauder and Banach fixed point theorems. Nevertheless, in this paper we introduce a novel method that unravels an ergodicity property of the operator $\mathcal{T}$, leading to quantitative estimates for the relaxation of profiles $\{F_n\}_{n \in \mathbb{N}}$ towards $F_\alpha$. First, we prove that an explicit (Gaussian) solution to the eigenproblem (1.11) exists. Second, by computing $n$ iterations of the operator $\mathcal{T}$, that is $F_n = \mathcal{T}^n[F_0]$, we notice that information of $F_n$ at the trait value $x$ is propagated from the initial datum $F_0$ across $2^n$ ancestors over a binary tree with height $n$ and rooted at $x$ (the pedigree chart). Interestingly, an appropriate reformulation of $\mathcal{T}^n$ in the case of Gaussian mixing (1.4) and quadratic selection (1.5) shows that the dependence of the solution $\{F_n\}_{n \in \mathbb{N}}$ on the initial datum $F_0$ is rapidly lost across the different
The parameter $\lambda \in \mathbb{R}$ represents the rate at which the number of individuals grows (if $\lambda \ge 0$) or decreases (if $\lambda \le 0$), and $F = F(x) \ge 0$ is an unknown probability density. By imposing such an ansatz on (1.6)-(1.3), the following generalized spectral problem arises for the couple $(\lambda, F)$:

$$ \left\{ \begin{array}{ll} \lambda F(x) + m(x)F(x) = \mathcal{B}[F](x), & x \in \mathbb{R}^d, \\ \int_{\mathbb{R}^d} F(x') dx' = 1. & \end{array} \right. \qquad (1.9) $$

Note that the operator $\mathcal{B}$ is genuinely non-linear so that methods based on the Krein-Rutman theory or maximum principles cannot be applied straightforwardly, see [6, 13] and the references therein for the linear case. Further, usual extensions of the Krein-Rutman theory to 1-homogeneous operators [16] cannot be applied neither because $\mathcal{B}$ is not monotone. To date, the main strategy behind the existence of solutions of (1.9) relies on a suitable application of Schauder fixed-point theorem to the operator $F \mapsto (\lambda + m)^{-1}\mathcal{B}[F]$ over an appropriate cone of $L^1(\mathbb{R}^d)$ that is conserved by the nonlinear operator, see [8]. However, uniqueness cannot be achieved by this method. Moreover, it has been proven in [10, Corollary 1.5] that several equilibrium states $(\lambda, F)$ can co-exist in the presence of multiple local minima of $m$ (provided the variance of $G$ is small enough). That is, the generalized eigenproblem (1.9) does not admit a unique positive eigenfunction, in contrast with general conclusions of the Krein-Rutman theory.

Recently, G. RAOUL addressed the long-time dynamics of (1.6) in 1D, with $\mathcal{R} = \mathcal{B}$ and a trait dependent fecundity rate [21]. He obtained local uniqueness and exponential relaxation under the assumption of weak and localized (compactly supported) selection effects. To this end, he controlled locally in space the Wasserstein distance between the solution and the stationary state, using the uniform contraction property of $\mathcal{B}$ in the space of probability measures sharing the same center of mass. Unfortunately, the Wasserstein metric is not fully compatible with multiplicative operators, such as trait-dependent fecundity (see the discussion in [21, Section 3.4] and Section 2 below). This can be circumvented under the additional assumption that the trait density is locally uniformly bounded below, following [7]. Obviously, this cannot hold globally for integrable densities, hence G. RAOUL developed estimates of the distribution's tail to complete the contraction estimates. Also, a lot of attention has to be paid to the dynamics of the center of the distribution which is essentially driven by selection. Indeed, in the case of flat selection ($m \equiv 0$), the problem is invariant by translation, so that local uniqueness cannot hold. To conclude this discussion, let us emphasize that global uniqueness and the asymptotic behavior of generic solutions to the evolution problem (1.6)-(1.3) is still open. Below, we provide a first result in this direction, for the time-discrete problem (1.1), though.

We remark that the time-discrete version (1.1) that we propose in this paper can be partially regarded as a discretization in time of the above time-continuous problem (1.6)-(1.3), with non-overlapping generations (see Appendix A for further details). As for the time-continuous problem, we could seek special solutions to the time-discrete problem (1.1) in the following form

$$ F_n(x) = \lambda^n F(x), \quad (n, x) \in \mathbb{N} \times \mathbb{R}^d. \qquad (1.10) $$

Again, $\lambda \in \mathbb{R}_+^*$ is the rate of growth (if $\lambda \ge 1$) or decrease (if $\lambda \le 1$) of individuals, and $F = F(x)$ is an unknown probability density. This yields the following generalized eigenproblem for the couple $(\lambda, F)$:

$$ \left\{ \begin{array}{ll} \lambda F(x) = \mathcal{T}[F](x), & x \in \mathbb{R}^d, \\ \int_{\mathbb{R}^d} F(x') dx' = 1. & \end{array} \right. \qquad (1.11) $$

In this paper we aim to address the following questions:

(Q1) Does the eigenproblem (1.11) have a unique solution $(\lambda_\alpha, F_\alpha)$, with $\lambda_\alpha \in \mathbb{R}$ and $F_\alpha$ being a probability measure, for each $\alpha \in \mathbb{R}_+^*$?

(Q2) Consider any generic initial datum $F_0 \in \mathcal{M}_+(\mathbb{R}^d)$ and its associated solution $\{F_n\}_{n \in \mathbb{N}}$ of the time-discrete problem (1.1). Do the renormalized profiles $F_n/\|F_n\|_{L^1(\mathbb{R}^d)}$ converge to the unique steady profile $F_\alpha$ solving (1.11) when $n \to \infty$?

We shall prove that the answer to both questions is affirmative. It stands to reason that similar existence and local uniqueness results like in [8, 10] could be extended to the new eigenproblem (1.11) by applying Schauder and Banach fixed point theorems. Nevertheless, in this paper we introduce a novel method that unravels an ergodicity property of the operator $\mathcal{T}$, leading to quantitative estimates for the relaxation of profiles $\{F_n\}_{n \in \mathbb{N}}$ towards $F_\alpha$. First, we prove that an explicit (Gaussian) solution to the eigenproblem (1.11) exists. Second, by computing $n$ iterations of the operator $\mathcal{T}$, that is $F_n = \mathcal{T}^n[F_0]$, we notice that information of $F_n$ at the trait value $x$ is propagated from the initial datum $F_0$ across $2^n$ ancestors over a binary tree with height $n$ and rooted at $x$ (the pedigree chart). Interestingly, an appropriate reformulation of $\mathcal{T}^n$ in the case of Gaussian mixing (1.4) and quadratic selection (1.5) shows that the dependence of the solution $\{F_n\}_{n \in \mathbb{N}}$ on the initial datum $F_0$ is rapidly lost across the different
levels of the tree. More specifically, a strong convergence of generic solutions {$F_n$} of the time-discrete problem (1.1) towards the steady Gaussian profile $F_\alpha$ solving (1.11) is achieved locally with respect to $x$. Third, we prove an appropriate propagation of quadratic and exponential moments, leading to uniform tightness of the family {$F_n$}$_{n \in \mathbb{N}}$. Finally, we glue all the information together and conclude the final global convergence result in question (Q2) in relative entropy. We refer to Section 2 for a more detailed sketch of our strategy of proof. Specifically, we obtain our main result:

**Theorem 1.1.** Assume that $\alpha \in \mathbb{R}_+^*$ and set any initial non-negative measure $F_0 \in \mathcal{M}_+(\mathbb{R}^d)$. The solution {$F_n$}$_{n \in \mathbb{N}}$ to the time-discrete problem (1.1) verifies that the growth rates $\|F_n\|_{L^1(\mathbb{R}^d)}/\|F_{n-1}\|_{L^1(\mathbb{R}^d)}$ relax towards $\lambda_\alpha$ and the normalized profiles $F_n/\|F_n\|_{L^1(\mathbb{R}^d)}$ relax towards $F_\alpha$ with

$$ \lambda_{\alpha} := \left(1 + \alpha \left(1 + \frac{\sigma_{\alpha}^2}{2}\right)\right)^{-\frac{d}{2}}, \quad F_{\alpha} := G_{0, \sigma_{\alpha}^2}, \qquad (1.12) $$

and the variance $\sigma_{\alpha}^2 \in \mathbb{R}_+^*$ is the unique positive root of the equation

$$ \frac{1}{\sigma_{\alpha}^2} = \alpha + \frac{1}{1 + \frac{\sigma_{\alpha}^2}{2}}, \quad \text{i.e.,} \quad \sigma_{\alpha}^2 = \frac{\sqrt{(1+2\alpha)^2 + 8\alpha} - (1+2\alpha)}{2\alpha}. \qquad (1.13) $$

Specifically, for any $\varepsilon \in \mathbb{R}_+^*$ there exists a sufficiently large $C_\varepsilon \in \mathbb{R}_+^*$ such that

$$ \begin{aligned} \mathcal{D}_{\text{KL}} \left( \frac{F_n}{\|F_n\|_{L^1(\mathbb{R}^d)}} \middle\| \mathbf{F}_\alpha \right) &\le C_\varepsilon ((2 \mathbf{k}_\alpha)^2 + \varepsilon)^n, \\ \left| \frac{\|F_n\|_{L^1(\mathbb{R}^d)}}{\|F_{n-1}\|_{L^1(\mathbb{R}^d)}} - \lambda_\alpha \right| &\le C_\varepsilon (2 \mathbf{k}_\alpha + \varepsilon)^n, \end{aligned} $$

for any $n \in \mathbb{N}$, where $\mathcal{D}_{\text{KL}}$ is the Kullback-Leibler divergence (or relative entropy), i.e.,

$$ \mathcal{D}_{\text{KL}}(P\|Q) := \int_{\mathbb{R}^d} P(x) \log \left( \frac{P(x)}{Q(x)} \right) dx, \qquad (1.14) $$

for any $P, Q \in L_+^1(\mathbb{R}^d) \cap \mathcal{P}(\mathbb{R}^d)$, and the coefficient $\mathbf{k}_\alpha \in (0, \frac{1}{2})$ reads

$$ \mathbf{k}_{\alpha} = \frac{\sigma_{\alpha}^2}{2 + \sigma_{\alpha}^2}, \quad \text{i.e.,} \quad \frac{(3 + 2\alpha) - \sqrt{(1 + 2\alpha)^2 + 8\alpha}}{4}. \qquad (1.15) $$

Note that any solution $(\lambda, F)$ to the eigenproblem (1.11) yields a solution to the time-discrete problem (1.1) through the ansatz (1.10). Then, question (Q1) will readily follow from question (Q2) and, in particular, the unique solution $(\lambda_\alpha, F_\alpha)$ to (1.11) is Gaussian.

**Corollary 1.2.** Assume that $\alpha \in \mathbb{R}_+^*$, then there exists a unique solution $(\lambda_\alpha, F_\alpha)$ to the eigenproblem (1.11), with $F_\alpha$ being a probability measure, namely given by (1.12).

**Remark 1.3** (About the assumption of quadratic selection). By the assumption of a quadratic selection function, we can henceforth push extensively explicit computations of the iterated operator. The latter consists in recursive multiplication and convolution by Gaussian functions. This enables capturing the essence of the relaxation phenomenon, and having a precise description of the behaviour at infinity, which crucially helps to localize the convergence argument. As a by-product, we are able to consider very general initial condition $F_0$. Two of the authors, together with F. Santambrogio, have obtained similar results when the selection function $m$ is more generally assumed to be strongly convex [12]. However, they imposed stringent conditions on the initial datum $F_0$, that is, it should behave at infinity as the equilibrium profile $F_\alpha$ in a very strong sense. It would be of interest to merge the two works, having a general (strongly convex) selection function, and a general initial datum. This is left for future work.

**Remark 1.4** (About the choice of metric). The convergence of the profiles has been quantified in Kullback-Leibler divergence (1.14) in Theorem 1.1, in contrast with the results in [20, 21], where the quadratic Wasserstein distance was used for perturbative regimes of the case without selection. We anticipate that there are several compelling reasons for such a change of metric in our non-perturbative setting:

(i) **(Quadratic Wasserstein distance)** When $\alpha = 0$, the operator $\mathcal{T}$ reduces to $\mathcal{B}$, which is non-expansive with respect to the quadratic Wasserstein distance, and indeed contractive over distributions with common center of mass, cf. Section 2.1. However, when $\alpha > 0$, the multiplicative operator leads to an operator $\mathcal{T}$ which is not even Lipschitz continuous with respect to the quadratic Wasserstein distance, cf Section 2.2. Therefore, the quadratic Wasserstein distance seems to be unadapted to scenarios where reproduction and selection operate together.
levels of the tree. More specifically, a strong convergence of generic solutions {$F_n$} of the time-discrete problem (1.1) towards the steady Gaussian profile $F_\alpha$ solving (1.11) is achieved locally with respect to $x$. Third, we prove an appropriate propagation of quadratic and exponential moments, leading to uniform tightness of the family {$F_n$}$_{n \in \mathbb{N}}$. Finally, we glue all the information together and conclude the final global convergence result in question (Q2) in relative entropy. We refer to Section 2 for a more detailed sketch of our strategy of proof. Specifically, we obtain our main result:

**Theorem 1.1.** Assume that $\alpha \in \mathbb{R}_+^*$ and set any initial non-negative measure $F_0 \in \mathcal{M}_+(\mathbb{R}^d)$. The solution {$F_n$}$_{n \in \mathbb{N}}$ to the time-discrete problem (1.1) verifies that the growth rates $\|F_n\|_{L^1(\mathbb{R}^d)}/\|F_{n-1}\|_{L^1(\mathbb{R}^d)}$ relax towards $\lambda_\alpha$ and the normalized profiles $F_n/\|F_n\|_{L^1(\mathbb{R}^d)}$ relax towards $F_\alpha$ with

$$ \lambda_{\alpha} := \left(1 + \alpha \left(1 + \frac{\sigma_{\alpha}^2}{2}\right)\right)^{-\frac{d}{2}}, \quad F_{\alpha} := G_{0, \sigma_{\alpha}^2}, \qquad (1.12) $$

and the variance $\sigma_\alpha^2 \in \mathbb{R}_+^*$ is the unique positive root of the equation

$$ \frac{1}{\sigma_{\alpha}^2} = \alpha + \frac{1}{1 + \frac{\sigma_{\alpha}^2}{2}}, \quad \text{i.e.,} \quad \sigma_{\alpha}^2 = \frac{\sqrt{(1+2\alpha)^2 + 8\alpha} - (1+2\alpha)}{2\alpha}. \qquad (1.13) $$

Specifically, for any $\varepsilon \in \mathbb{R}_+^*$ there exists a sufficiently large $C_\varepsilon \in \mathbb{R}_+^*$ such that

$$ \begin{aligned} \mathcal{D}_{\text{KL}} \left( \frac{F_n}{\|F_n\|_{L^1(\mathbb{R}^d)}} \middle\| \mathbf{F}_\alpha \right) &\le C_\varepsilon ((2 \mathbf{k}_\alpha)^2 + \varepsilon)^n, \\ \left| \frac{\|F_n\|_{L^1(\mathbb{R}^d)}}{\|F_{n-1}\|_{L^1(\mathbb{R}^d)}} - \lambda_\alpha \right| &\le C_\varepsilon (2 \mathbf{k}_\alpha + \varepsilon)^n, \end{aligned} $$

for any $n \in \mathbb{N}$, where $\mathcal{D}_{\text{KL}}$ is the Kullback-Leibler divergence (or relative entropy), i.e.,

$$ \mathcal{D}_{\text{KL}}(P\|Q) := \int_{\mathbb{R}^d} P(x) \log \left( \frac{P(x)}{Q(x)} \right) dx, \qquad (1.14) $$

for any $P, Q \in L_+^1(\mathbb{R}^d) \cap \mathcal{P}(\mathbb{R}^d)$, and the coefficient $\mathbf{k}_\alpha \in (0, \frac{1}{2})$ reads

$$ \mathbf{k}_{\alpha} = \frac{\sigma_{\alpha}^{2}}{2 + \sigma_{\alpha}^{2}}, \quad \text{i.e.,} \quad \frac{(3 + 2\alpha) - \sqrt{(1 + 2\alpha)^2 + 8\alpha}}{4}. \qquad (1.15) $$

Note that any solution $(\lambda, F)$ to the eigenproblem (1.11) yields a solution to the time-discrete problem (1.1) through the ansatz (1.10). Then, question (Q1) will readily follow from question (Q2) and, in particular, the unique solution $(\lambda_\alpha, F_\alpha)$ to (1.11) is Gaussian.

**Corollary 1.2.** Assume that $\alpha \in \mathbb{R}_+^*$, then there exists a unique solution $(\lambda_\alpha, F_\alpha)$ to the eigenproblem (1.11), with $F_\alpha$ being a probability measure, namely given by (1.12).

**Remark 1.3** (About the assumption of quadratic selection). By the assumption of a quadratic selection function, we can henceforth push extensively explicit computations of the iterated operator. The latter consists in recursive multiplication and convolution by Gaussian functions. This enables capturing the essence of the relaxation phenomenon, and having a precise description of the behaviour at infinity, which crucially helps to localize the convergence argument. As a by-product, we are able to consider very general initial condition $F_0$. Two of the authors, together with F. Santambrogio, have obtained similar results when the selection function $m$ is more generally assumed to be strongly convex [12]. However, they imposed stringent conditions on the initial datum $F_0$, that is, it should behave at infinity as the equilibrium profile $F_\alpha$ in a very strong sense. It would be of interest to merge the two works, having a general (strongly convex) selection function, and a general initial datum. This is left for future work.

**Remark 1.4** (About the choice of metric). The convergence of the profiles has been quantified in Kullback-Leibler divergence (1.14) in Theorem 1.1, in contrast with the results in [20, 21], where the quadratic Wasserstein distance was used for perturbative regimes of the case without selection. We anticipate that there are several compelling reasons for such a change of metric in our non-perturbative setting:

(i) **(Quadratic Wasserstein distance)** When $\alpha = 0$, the operator $\mathcal{T}$ reduces to $\mathcal{B}$, which is non-expansive with respect to the quadratic Wasserstein distance, and indeed contractive over distributions with common center of mass, cf. Section 2.1. However, when $\alpha > 0$, the multiplicative operator leads to an operator $\mathcal{T}$ which is not even Lipschitz continuous with respect to the quadratic Wasserstein distance, cf Section 2.2. Therefore, the quadratic Wasserstein distance seems to be unadapted to scenarios where reproduction and selection operate together.
(ii) **(Log-Lipschitz norm)** As we show in Section 2.3, the operator $\mathcal{T}$ is non-expansive in the log-Lipschits norm $\|\nabla \log \frac{F}{F_\alpha}\|_{L^\infty(\mathbb{R}^d)}$ for all $\alpha \ge 0$, and indeed it is contractive if $\alpha > 0$. However, this contraction only gives actual information when the initial datum $F_0$ has identical tails to the Gaussian density $F_\alpha$ so that initially the log-Lipschitz norm is finite.

(iii) **(Relative entropy distance)** Note that the above log-Lipschitz norm amounts to the natural $L^\infty$ version of the relative Fisher information $\|\nabla \log \frac{F}{F_\alpha}\|_{L^2(\mathbb{R}^d, F)}$. By the log-Soboled inequality, contraction of the log-Lipschits norm readily implies decay of the Kullback-Leibler divergence, which is a more standard metric in relative entropy arguments. However, we emphasize that our use of the Kullback-Leibler divergence is not only aesthetic, but we actually need it in order to go beyond the above structural constraint on the tails of the initial data. Specifically, for generic initial data we need to prove an appropriate shaping of tails over time, which cannot be expressed using uniform norms, but only in an averaged sense (compatible with the relative entropy).

Altogether justifies that the whilst the quadratic Wasserstein distance is useful in perturbative regimes, the use of alternative norms is necessary to quantify contraction in purely non-perturbative settings.

**Remark 1.5** (About the positivity of $\alpha$). Our form of ergodicity, as measured in Theorem 1.1, it breaks down when $\alpha = 0$, simply because the operator is invariant by translation in that case, and it admits a one-parameter family $F_{\alpha=0}(\cdot - \mu)$ with $\mu \in \mathbb{R}^d$ of fixed Gaussian probability densities. Nevertheless, as mentioned in item (i) above, when $\alpha = 0$ and we restrict to centered initial data, there is convergence to the right centered Gaussian probability density $F_{\alpha=0}$ with respect to the Wasserstein distance. It is an open problem to make both approaches meet for $\alpha = 0$, and prove contraction in norms stronger than the Wasserstein distance, but comparable to the log-Lipschitz norm, in this subclass of initial data.

The rest of the paper is organized as follows. In Section 2 we discuss about the generic incompatibility of the Wasserstein distance with a multiplicative operator and we provide a brief outlook of the strategy of our proof. In Section 3 we provide some necessary notation and we introduce the special class of Gaussian solutions of both problems (1.1) and (1.11), which will inspire some parts of the paper. Section 4 is devoted to introduce some main properties of $\mathcal{T}$ regarding the emergence of Gaussian behavior (in the large) from generic initial data, and a suitable propagation of quadratic and exponential moments across generations. In Section 5 we reformulate (1.1) via a high-dimensional integral operator propagating ancestors' information across the different levels of the pedigree chart, which will be the cornerstone to study the long-term dynamics. In Section 6 we prove our main results, namely, Theorem 1.1 and Corollary 1.2. Section 7 contains some numerical experiments that illustrate the results in this paper. In Section 8 we provide some conclusions and perspectives. Finally, Appendix A contains the adimensionalization of the problem, and the relationship between (1.1) and the previous time-continuous analogues in the literature is discussed. A full list of the main notations in the paper is presented in Table 1.

## 2. MOTIVATION AND STRATEGY OF THE PROOF OF THEOREM 1.1

In this section, we discuss the incompatibility of the quadratic Wasserstein distance to quantify directly contractivity under the joint effect of the reproduction operator $\mathcal{B}$ in (1.3) and a generic multiplicative selection $e^{-m}$. Although a small perturbation of the case of flat selection ($m \equiv 0$) could still be considered via a perturbative argument (see discussion above, [21] and also [20]), our novel approach is able to tackle a purely non-perturbative setting. We end the section by briefly discussing the strategy of our proof.

**2.1. Some properties of the sexual reproduction operator.** We start by recalling some of the main properties of the sexual reproduction operator $\mathcal{B}$. Since the fecundity rate has been normalized to 1 (see Appendix A), then $\mathcal{B}$ preserves the mass and center of mass, namely,

$$ \|\mathcal{B}[F]\|_{L^1(\mathbb{R}^d)} = \|F\|_{\mathcal{M}(\mathbb{R}^d)}, \quad \int_{\mathbb{R}^d} x \mathcal{B}[F](x) dx = \int_{\mathbb{R}^d} x F(dx), \qquad (2.1) $$

for any $F \in \mathcal{M}_+(\mathbb{R}^d)$. Furthermore, it is contractive in the space of probability measures with a common center of mass, endowed with the quadratic Wasserstein metric. As discussed above, this property has been used fruitfully by G. RAOUL (cf. [21]) to analyse the long term behavior of the time-continuous problem in the regime of weak (and compactly supported) selection acting on fecundity. For the sake of clarity, we recall this fact and its proof below:

**Lemma 2.1.** Assume that $F_1, F_2 \in \mathcal{P}_2(\mathbb{R}^d)$ and they have the same center of mass. Then,

$$ W_2^2(\mathcal{B}[F_1], \mathcal{B}[F_2]) \le \frac{1}{2} W_2^2(F_1, F_2). $$
(ii) **(Log-Lipschitz norm)** As we show in Section 2.3, the operator $\mathcal{T}$ is non-expansive in the log-Lipschits norm $\|\nabla \log \frac{F}{F_\alpha}\|_{L^\infty(\mathbb{R}^d)}$ for all $\alpha \ge 0$, and indeed it is contractive if $\alpha > 0$. However, this contraction only gives actual information when the initial datum $F_0$ has identical tails to the Gaussian density $F_\alpha$ so that initially the log-Lipschitz norm is finite.

(iii) **(Relative entropy distance)** Note that the above log-Lipschitz norm amounts to the natural $L^\infty$ version of the relative Fisher information $\|\nabla \log \frac{F}{F_\alpha}\|_{L^2(\mathbb{R}^d, F)}$. By the log-Soboled inequality, contraction of the log-Lipschits norm readily implies decay of the Kullback-Leibler divergence, which is a more standard metric in relative entropy arguments. However, we emphasize that our use of the Kullback-Leibler divergence is not only aesthetic, but we actually need it in order to go beyond the above structural constraint on the tails of the initial data. Specifically, for generic initial data we need to prove an appropriate shaping of tails over time, which cannot be expressed using uniform norms, but only in an averaged sense (compatible with the relative entropy).

Altogether justifies that the whilst the quadratic Wasserstein distance is useful in perturbative regimes, the use of alternative norms is necessary to quantify contraction in purely non-perturbative settings.

**Remark 1.5** (About the positivity of $\alpha$). Our form of ergodicity, as measured in Theorem 1.1, it breaks down when $\alpha = 0$, simply because the operator is invariant by translation in that case, and it admits a one-parameter family $F_{\alpha=0}(\cdot - \mu)$ with $\mu \in \mathbb{R}^d$ of fixed Gaussian probability densities. Nevertheless, as mentioned in item (i) above, when $\alpha = 0$ and we restrict to centered initial data, there is convergence to the right centered Gaussian probability density $F_{\alpha=0}$ with respect to the Wasserstein distance. It is an open problem to make both approaches meet for $\alpha = 0$, and prove contraction in norms stronger than the Wasserstein distance, but comparable to the log-Lipschitz norm, in this subclass of initial data.

The rest of the paper is organized as follows. In Section 2 we discuss about the generic incompatibility of the Wasserstein distance with a multiplicative operator and we provide a brief outlook of the strategy of our proof. In Section 3 we provide some necessary notation and we introduce the special class of Gaussian solutions of both problems (1.1) and (1.11), which will inspire some parts of the paper. Section 4 is devoted to introduce some main properties of $\mathcal{T}$ regarding the emergence of Gaussian behavior (in the large) from generic initial data, and a suitable propagation of quadratic and exponential moments across generations. In Section 5 we reformulate (1.1) via a high-dimensional integral operator propagating ancestors' information across the different levels of the pedigree chart, which will be the cornerstone to study the long-term dynamics. In Section 6 we prove our main results, namely, Theorem 1.1 and Corollary 1.2. Section 7 contains some numerical experiments that illustrate the results in this paper. In Section 8 we provide some conclusions and perspectives. Finally, Appendix A contains the adimensionalization of the problem, and the relationship between (1.1) and the previous time-continuous analogues in the literature is discussed. A full list of the main notations in the paper is presented in Table 1.

## 2. MOTIVATION AND STRATEGY OF THE PROOF OF THEOREM 1.1

In this section, we discuss the incompatibility of the quadratic Wasserstein distance to quantify directly contractivity under the joint effect of the reproduction operator $\mathcal{B}$ in (1.3) and a generic multiplicative selection $e^{-m}$. Although a small perturbation of the case of flat selection ($m \equiv 0$) could still be considered via a perturbative argument (see discussion above, [21] and also [20]), our novel approach is able to tackle a purely non-perturbative setting. We end the section by briefly discussing the strategy of our proof.

**2.1. Some properties of the sexual reproduction operator.** We start by recalling some of the main properties of the sexual reproduction operator $\mathcal{B}$. Since the fecundity rate has been normalized to 1 (see Appendix A), then $\mathcal{B}$ preserves the mass and center of mass, namely,

$$ \|\mathcal{B}[F]\|_{L^1(\mathbb{R}^d)} = \|F\|_{\mathcal{M}(\mathbb{R}^d)}, \quad \int_{\mathbb{R}^d} x \mathcal{B}[F](x) dx = \int_{\mathbb{R}^d} x F(dx), \qquad (2.1) $$

for any $F \in \mathcal{M}_+(\mathbb{R}^d)$. Furthermore, it is contractive in the space of probability measures with a common center of mass, endowed with the quadratic Wasserstein metric. As discussed above, this property has been used fruitfully by G. RAOUL (cf. [21]) to analyse the long term behavior of the time-continuous problem in the regime of weak (and compactly supported) selection acting on fecundity. For the sake of clarity, we recall this fact and its proof below:

**Lemma 2.1.** Assume that $F_1, F_2 \in \mathcal{P}_2(\mathbb{R}^d)$ and they have the same center of mass. Then,

$$ W_2^2(\mathcal{B}[F_1], \mathcal{B}[F_2]) \le \frac{1}{2} W_2^2(F_1, F_2). $$
<table><caption>TABLE 1. List of notations</caption><thead><tr><th>Notation</th><th>Meaning</th><th>Page</th></tr></thead><tbody><tr><td><i>M</i>(<b>R</b><sup>d</sup>), <i>M</i><sub>+</sub>(<b>R</b><sup>d</sup>)</td><td>Signed and non-negative finite Radon measures</td><td>4</td></tr><tr><td><i>P</i>(<b>R</b><sup>d</sup>), <i>P</i><sub>2</sub>(<b>R</b><sup>d</sup>)</td><td>Probability measures (with finite 2nd order moment)</td><td>5</td></tr><tr><td><i>P</i> ⊗ <i>Q</i> ∈ <i>M</i>(<b>R</b><sup>2d</sup>)</td><td>Tensor product of the finite measures <i>P</i>, <i>Q</i> ∈ <i>M</i>(<b>R</b><sup>d</sup>)</td><td>43</td></tr><tr><td><i>D</i><sub>KL</sub>(<i>P</i>||<i>Q</i>)</td><td>Kullback-Leibler divergence between <i>P</i>, <i>Q</i> ∈ <i>P</i>(<b>R</b><sup>d</sup>)</td><td>4</td></tr><tr><td><i>W</i><sub>2</sub>(<i>P</i>, <i>Q</i>)</td><td>Quadratic Wasserstein distance between <i>P</i>, <i>Q</i> ∈ <i>P</i><sub>2</sub>(<b>R</b><sup>d</sup>)</td><td>6</td></tr><tr><td><i>G</i><sub>μ,σ<sup>2</sup></sub></td><td>Gaussian with mean <i>μ</i> ∈ <b>R</b><sup>d</sup> and covariance <i>σ</i><sup>2</sup><i>I</i><sub>d</sub> ∈ <b>R</b><sup>d×d</sup></td><td>13</td></tr><tr><td><i>G</i><sub>μ,Σ</sub></td><td>Gaussian with mean <i>μ</i> ∈ <b>R</b><sup>d</sup> and covariance <i>Σ</i> ∈ <b>R</b><sup>d×d</sup></td><td>33</td></tr><tr><td><i>N</i>(<i>μ</i>, <i>Σ</i>)</td><td>Normal with mean <i>μ</i> ∈ <b>R</b><sup>d</sup> and covariance <i>Σ</i> ∈ <b>R</b><sup>d×d</sup></td><td>33</td></tr><tr><td><i>E</i>[<i>X</i>]</td><td>Expectation of a random variable <i>X</i></td><td>34</td></tr><tr><td><i>G</i>(<i>x</i>)</td><td>Gaussian mixing kernel <i>G</i><sub>0,I<sub>d</sub></sub></td><td>2</td></tr><tr><td><i>m</i>(<i>x</i>)</td><td>Quadratic selection function <i>α</i><sub>2</sub>|<i>x</i>|<sup>2</sup></td><td>2</td></tr><tr><td><i>B</i>[<i>F</i>]</td><td>Fisher's infinitesimal operator</td><td>2</td></tr><tr><td><i>T</i>[<i>F</i>]</td><td>Selection-reproduction operator</td><td>1</td></tr><tr><td><i>M</i>[<i>F</i>]</td><td>Normalized multiplicative operator by <i>e</i><sup>-<i>m</i></sup></td><td>8</td></tr><tr><td><i>S</i>[<i>F</i>]</td><td>Scaled selection-reproduction operator</td><td>21</td></tr><tr><td><i>E</i><sub>n</sub>[<i>F</i>]</td><td>High-dimensional integral operators</td><td>35</td></tr><tr><td>{<i>F</i><sub>n</sub>}<sub><i>n</i>∈<i>N</i></sub></td><td>Solution to the time-evolution equation (1.1)</td><td>1</td></tr><tr><td>(<b>λ</b><sub>α</sub>, <b>F</b><sub>α</sub>)</td><td>Gaussian solution to the non-linear eigenproblem (1.11)</td><td>4</td></tr><tr><td><b>σ</b><sub>α</sub><sup>2</sup></td><td>Variance of Gaussian eigenfunction <b>F</b><sub>α</sub></td><td>4</td></tr><tr><td><i>F</i> = <i>e</i><sup><i>m</i></sup><i>F</i> / <i>F</i><sub>α=0</sub></td><td>Normalized profile associated to <i>F</i> ∈ <i>M</i><sub>+</sub>(<b>R</b><sup>d</sup>)</td><td>26</td></tr><tr><td><b>T</b><sup>n</sup>, <b>T</b><sub>*</sub><sup>n</sup>, <b>T</b><sup>n</sup></td><td>Perfect binary tree, rootless tree and leafless tree</td><td>11</td></tr><tr><td><b>L</b><sub>m</sub><sup>n</sup>, <b>L</b><sup>n</sup></td><td>Level <i>m</i> and leaves (level <i>n</i>) of the tree</td><td>11</td></tr><tr><td><i>i</i>1, <i>i</i>2</td><td>Parents of a node <i>i</i> ∈ <b>T</b><sup>n</sup> of the tree</td><td>11</td></tr><tr><td><b>x</b><sub><i>n</i></sub> = (<i>x</i><sub><i>i</i></sub>)<sub><i>i</i>∈<b>T</b><sub>*</sub><sup>n</sup></sub></td><td>Variables indexed by the rootless tree</td><td>13</td></tr><tr><td><b>y</b><sub><i>n</i></sub> = (<i>y</i><sub><i>i</i></sub>)<sub><i>i</i>∈<b>T</b><sub>*</sub><sup>n</sup></sub></td><td>Variables indexed by leafs</td><td>13</td></tr><tr><td><i>z</i><sub><i>n</i></sub> = (<i>z</i><sub><i>j</i></sub>)<sub><i>j</i>∈<b>L</b><sup>n</sup></sub></td><td>Variables indexed by leafs</td><td>13</td></tr><tr><td>||<b>x</b>|| = (<sum><i>x</i><sub><i>i</i></sub></i> ||<i>x</i><sub><i>i</i></sub>||<sup>2</sup>)<sup>1/2</sup></td><td><i>ℓ</i><sub>2</sub> sum of Euclidean norms of <b>x</b> = (<i>x</i><sub>1</sub>, ..., <i>x</i><sub><i>n</i></sub>) ∈ <b>R</b><sup><i>dn</i></sup></td><td>36</td></tr><tr><td><i>Φ</i><sub><i>n</i></sub><sup><i>j</i></sup>(<i>x</i>; <b>y</b><sub><i>n</i></sub>)</td><td>Lineage map from leaf <i>j</i> ∈ <b>L</b><sup>n</sup> to root value <i>x</i> ∈ <b>R</b><sup>d</sup></td><td>30</td></tr><tr><td><i>u</i> ⊗ <i>v</i> ∈ <b>R</b><sup><i>d</i>×<i>d</i></sup></td><td>Kronecker product (<i>u</i><sub><i>i</i></sub><i>v</i><sub><i>j</i></sub>)<sub><i>1</i>≤<i>i</i>,<i>j</i>≤<i>N</i></sub> of vectors <i>u</i>, <i>v</i> ∈ <b>R</b><sup>d</sup></td><td>34</td></tr><tr><td><i>k</i><sub><i>n</i></sub>, <i>κ</i><sub><i>n</i></sub></td><td>Sequences of coefficients in the change of variables</td><td>28</td></tr><tr><td>(2 <b>k</b><sub>α</sub>)<sup>2</sup></td><td>Convergence rate of <i>D</i><sub>KL</sub></td><td>4</td></tr><tr><td>2 <b>k</b><sub>α</sub></td><td>Convergence rate of the log-Lipschitz norm</td><td>8</td></tr><tr><td><b>r</b><sub>α</sub></td><td>Relaxation rates of variances recursion</td><td>15</td></tr></tbody></table>

*Proof.* Recall the following dual characterizations of the quadratic Wasserstein distance (cf. [1, 24]):

$$ 
\begin{aligned}
W_2^2(F_1, F_2) &= \inf \left\{ \int_{\mathbb{R}^{2d}} |x-y|^2 \gamma(dx, dy) : \gamma \in \mathcal{P}(\mathbb{R}^d \times \mathbb{R}^d), \pi_1\#\gamma = F_1, \pi_2\#\gamma = F_2 \right\} \\
&= \sup \left\{ \int_{\mathbb{R}^d} \phi_1 dF_1 + \int_{\mathbb{R}^d} \phi_2 dF_2 : |\phi_1(x) + \phi_2(y)| \le |x-y|^2 \quad \forall x, y \in \mathbb{R}^d \right\},
\end{aligned} 
\quad (2.2) $$

where $\pi_i: \mathbb{R}^d \times \mathbb{R}^d \to \mathbb{R}^d$ is the projection onto the $i$-th component for $i = 1, 2$. Taking any couple $(\phi_1, \phi_2)$ as above and using the specific form of the operator $\mathcal{B}$ we obtain:

$$ 
\begin{aligned}
& \int_{\mathbb{R}^d} \mathcal{B}[F_1](x) \phi_1(x) dx + \int_{\mathbb{R}^d} \mathcal{B}[F_2](y) \phi_2(y) dy \\
&= \int_{\mathbb{R}^{3d}} G\left(x - \frac{x_1 + x_2}{2}\right) \phi_1(x) F_1(dx_1) F_1(dx_2) dx + \int_{\mathbb{R}^{3d}} G\left(y - \frac{y_1 + y_2}{2}\right) \phi_2(y) F_2(dy_1) F_2(dy_2) dy \\
&= \int_{\mathbb{R}^{3d}} G(z) \phi_1\left(z + \frac{x_1 + x_2}{2}\right) F_1(dx_1) F_1(dx_2) dz + \int_{\mathbb{R}^{3d}} G(z) \phi_2\left(z + \frac{y_1 + y_2}{2}\right) F_2(dy_1) F_2(dy_2) dz.
\end{aligned}
$$
<table><caption>TABLE 1. List of notations</caption><thead><tr><th>Notation</th><th>Meaning</th><th>Page</th></tr></thead><tbody><tr><td>M(ℝ<sup>d</sup>), M<sub>+</sub>(ℝ<sup>d</sup>)</td><td>Signed and non-negative finite Radon measures</td><td>4</td></tr><tr><td>P(ℝ<sup>d</sup>), P<sub>2</sub>(ℝ<sup>d</sup>)</td><td>Probability measures (with finite 2nd order moment)</td><td>5</td></tr><tr><td>P ⊗ Q ∈ M(ℝ<sup>2d</sup>)</td><td>Tensor product of the finite measures P, Q ∈ M(ℝ<sup>d</sup>)</td><td>43</td></tr><tr><td>D<sub>KL</sub>(P||Q)</td><td>Kullback-Leibler divergence between P, Q ∈ P(ℝ<sup>d</sup>)</td><td>4</td></tr><tr><td>W<sub>2</sub>(P, Q)</td><td>Quadratic Wasserstein distance between P, Q ∈ P<sub>2</sub>(ℝ<sup>d</sup>)</td><td>6</td></tr><tr><td>G<sub>μ,σ<sup>2</sup></sub></td><td>Gaussian with mean μ ∈ ℝ<sup>d</sup> and covariance σ<sup>2</sup>I<sub>d</sub> ∈ ℝ<sup>d×d</sup></td><td>13</td></tr><tr><td>G<sub>μ,Σ</sub></td><td>Gaussian with mean μ ∈ ℝ<sup>d</sup> and covariance Σ ∈ ℝ<sup>d×d</sup></td><td>33</td></tr><tr><td>N(μ, Σ)</td><td>Normal with mean μ ∈ ℝ<sup>d</sup> and covariance Σ ∈ ℝ<sup>d×d</sup></td><td>33</td></tr><tr><td>E[X]</td><td>Expectation of a random variable X</td><td>34</td></tr><tr><td>G(x)</td><td>Gaussian mixing kernel G<sub>0,I<sub>d</sub></sub></td><td>2</td></tr><tr><td>m(x)</td><td>Quadratic selection function α/2|x|<sup>2</sup></td><td>2</td></tr><tr><td>B[F]</td><td>Fisher's infinitesimal operator</td><td>2</td></tr><tr><td>T[F]</td><td>Selection-reproduction operator</td><td>1</td></tr><tr><td>M[F]</td><td>Normalized multiplicative operator by e<sup>-m</sup></td><td>8</td></tr><tr><td>S[F]</td><td>Scaled selection-reproduction operator</td><td>21</td></tr><tr><td>E<sub>n</sub>[F]</td><td>High-dimensional integral operators</td><td>35</td></tr><tr><td>{F<sub>n</sub>}<sub>n∈N</sub></td><td>Solution to the time-evolution equation (1.1)</td><td>1</td></tr><tr><td>(λ<sub>α</sub>, F<sub>α</sub>)</td><td>Gaussian solution to the non-linear eigenproblem (1.11)</td><td>4</td></tr><tr><td>σ<sub>α</sub><sup>2</sup></td><td>Variance of Gaussian eigenfunction F<sub>α</sub></td><td>4</td></tr><tr><td>F = e<sup>mF</sup> / F<sub>α=0</sub></td><td>Normalized profile associated to F ∈ M<sub>+</sub>(ℝ<sup>d</sup>)</td><td>26</td></tr><tr><td>T<sup>n</sup>, T<sup>n</sup><sub>*</sub>, T<sup>n</sup></td><td>Perfect binary tree, rootless tree and leafless tree</td><td>11</td></tr><tr><td>L<sup>n</sup><sub>m</sub>, L<sup>n</sup></td><td>Level m and leaves (level n) of the tree</td><td>11</td></tr><tr><td>i1, i2</td><td>Parents of a node i ∈ T<sup>n</sup> of the tree</td><td>11</td></tr><tr><td>x<sub>n</sub> = (x<sub>i</sub>)<sub>i∈T<sup>n</sup><sub>*</sub></sub></td><td>Variables indexed by the rootless tree</td><td>13</td></tr><tr><td>y<sub>n</sub> = (y<sub>i</sub>)<sub>i∈T<sup>n</sup><sub>*</sub></sub></td><td>Variables indexed by leafs</td><td>13</td></tr><tr><td>z<sub>n</sub> = (z<sub>j</sub>)<sub>j∈L<sup>n</sup></sub></td><td>Variables indexed by leafs</td><td>13</td></tr><tr><td>|| x || = (Σ<sup>n</sup><sub>i=1</sub> |x<sub>i</sub>|<sup>2</sup>)<sup>1/2</sup></td><td>ℓ<sub>2</sub> sum of Euclidean norms of x = (x<sub>1</sub>,..., x<sub>n</sub>) ∈ ℝ<sup>dn</sup></td><td>36</td></tr><tr><td>Φ<sup>j</sup><sub>n</sub>(x; y<sub>n</sub>)</td><td>Lineage map from leaf j ∈ L<sup>n</sup> to root value x ∈ ℝ<sup>d</sup></td><td>30</td></tr><tr><td>u ⊗ v ∈ ℝ<sup>d×d</sup></td><td>Kronecker product (u<sub>i</sub>v<sub>j</sub>)<sub>1≤i,j≤N</sub> of vectors u, v ∈ ℝ<sup>d</sup></td><td>34</td></tr><tr><td>k<sub>n</sub>, κ<sub>n</sub></td><td>Sequences of coefficients in the change of variables</td><td>28</td></tr><tr><td>(2 k<sub>α</sub>)<sup>2</sup></td><td>Convergence rate of D<sub>KL</sub></td><td>4</td></tr><tr><td>2 k<sub>α</sub></td><td>Convergence rate of the log-Lipschitz norm</td><td>8</td></tr><tr><td>r<sub>α</sub></td><td>Relaxation rates of variances recursion</td><td>15</td></tr></tbody></table>

*Proof.* Recall the following dual characterizations of the quadratic Wasserstein distance (cf. [1, 24]):

$$
\begin{align}
W_2^2(F_1, F_2) &= \inf \left\{ \int_{\mathbb{R}^{2d}} |x-y|^2 \gamma(dx, dy) : \gamma \in \mathcal{P}(\mathbb{R}^d \times \mathbb{R}^d), \pi_{1\#}\gamma = F_1, \pi_{2\#}\gamma = F_2 \right\} \tag{2.2} \\
&= \sup \left\{ \int_{\mathbb{R}^d} \phi_1 dF_1 + \int_{\mathbb{R}^d} \phi_2 dF_2 : |\phi_1(x) + \phi_2(y)| \le |x-y|^2 \quad \forall x, y \in \mathbb{R}^d \right\}, \nonumber
\end{align}
$$

where $\pi_i: \mathbb{R}^d \times \mathbb{R}^d \to \mathbb{R}^d$ is the projection onto the $i$-th component for $i = 1, 2$. Taking any couple $(\phi_1, \phi_2)$ as above and using the specific form of the operator $\mathcal{B}$ we obtain:

$$
\begin{align*}
& \int_{\mathbb{R}^d} \mathcal{B}[F_1](x) \phi_1(x) dx + \int_{\mathbb{R}^d} \mathcal{B}[F_2](y) \phi_2(y) dy \\
&= \int_{\mathbb{R}^{3d}} G\left(x - \frac{x_1 + x_2}{2}\right) \phi_1(x) F_1(dx_1) F_1(dx_2) dx + \int_{\mathbb{R}^{3d}} G\left(y - \frac{y_1 + y_2}{2}\right) \phi_2(y) F_2(dy_1) F_2(dy_2) dy \\
&= \int_{\mathbb{R}^{3d}} G(z) \phi_1\left(z + \frac{x_1 + x_2}{2}\right) F_1(dx_1) F_1(dx_2) dz + \int_{\mathbb{R}^{3d}} G(z) \phi_2\left(z + \frac{y_1 + y_2}{2}\right) F_2(dy_1) F_2(dy_2) dz.
\end{align*}
$$
Consider any transference plan $\gamma \in \mathcal{P}(\mathbb{R}^d \times \mathbb{R}^d)$ between $(F_1, F_2)$ as above. Specifically, we have that $\pi_{1\#}\gamma = F_1$ and $\pi_{2\#}\gamma = F_2$. Then, we can gather both integrals as follows:

$$ \int_{\mathbb{R}^d} \mathcal{B}[F_1](x) \phi_1(x) dx + \int_{\mathbb{R}^d} \mathcal{B}[F_2](y) \phi_2(y) dy \\ = \int_{\mathbb{R}^{5d}} G(z) \left( \phi_1 \left( z + \frac{x_1+x_2}{2} \right) + \phi_2 \left( z + \frac{y_1+y_2}{2} \right) \right) \gamma(dx_1, dy_1) \gamma(dx_2, dy_2) dz. $$

Since the condition $|\phi_1(x) + \phi_2(y)| \le |x - y|^2$ is verified for all $x, y \in \mathbb{R}^d$, then we find

$$ \int_{\mathbb{R}^d} \mathcal{B}[F_1](x) \phi_1(x) dx + \int_{\mathbb{R}^d} \mathcal{B}[F_2](y) \phi_2(y) dy \\ \le \frac{1}{4} \int_{\mathbb{R}^{5d}} G(z) |(x_1+x_2) - (y_1+y_2)|^2 \gamma(dx_1, dy_1) \gamma(dx_2, dy_2) dz \\ = \frac{1}{4} \int_{\mathbb{R}^{4d}} |(x_1+x_2) - (y_1+y_2)|^2 \gamma(dx_1, dy_1) \gamma(dx_2, dy_2) \\ = \frac{1}{4} \int_{\mathbb{R}^{2d}} |x_1 - y_1|^2 \gamma(dx_1, dy_1) + \frac{1}{4} \int_{\mathbb{R}^{2d}} |x_2 - y_2|^2 \gamma(dx_2, dy_2), $$

where in the third line we have used that $G$ is a probability density and in the last line we have used the crucial fact that $F_1$ and $F_2$ share the same center of mass in order to cancel the cross-terms (otherwise the estimate would boil down to a non-expansiveness estimate). Taking supremum over $(\phi_1, \phi_2)$ and infimum over $\gamma$ and using the dual characterizations (2.2) yields the result. $\square$

Note that the fact that both $F_1$ and $F_2$ have the same center of mass has been crucially used to cancel the crossed term. Otherwise, by the Cauchy-Schwarz inequality we would merely obtain non-expansiveness:

$$ W_2(\mathcal{B}[F_1], \mathcal{B}[F_2]) \le W_2(F_1, F_2). \qquad (2.3) $$

Using the contractivity property of Lemma 2.1 along with the conservation of mass and center of mass in (2.1) yields the long term dynamics of (1.1) in the special case $\alpha = 0$.

**Corollary 2.2.** Set any initial datum $F_0 \in \mathcal{M}_+(\mathbb{R}^d)$ such that $\int_{\mathbb{R}^d} |x|^2 F_0(dx) < \infty$ and consider the solution $\{F_n\}_{n \in \mathbb{N}}$ to the time-discrete problem (1.1) with $\alpha = 0$, i.e., $F_n = \mathcal{B}[F_{n-1}]$ for all $n \in \mathbb{N}$. Then,

$$ \frac{\|F_n\|_{L^1(\mathbb{R}^d)}}{\|F_{n-1}\|_{L^1(\mathbb{R}^d)}} = 1, \quad W_2\left(\frac{F_n}{\|F_n\|_{L^1(\mathbb{R}^d)}}, G_{\mu_0,2}\right) \lesssim \frac{1}{2^{n/2}}, $$

for every $n \in \mathbb{N}$, where $\mu_0 := \int_{\mathbb{R}^d} x F_0(dx)$. In particular, the set of stationary distributions under $\mathcal{B}$ is $\{\mathbf{F}_{\alpha=0}(\cdot - \mu) : \mu \in \mathbb{R}^d\}$, where $\mathbf{F}_{\alpha=0} = G_{0,2}$ will denote here on the Gaussian centered at the origin with variance equals 2 in agreement with the notation (1.12) in Theorem 1.1.

The preceding result might be regarded as the counterpart of Theorem 1.1 for $\alpha = 0$. Indeed, by Talagrand transportation inequality for a Gaussian measure [23] we obtain the relation

$$ W_2^2\left(\frac{F_n}{\|F_n\|_{L^1(\mathbb{R}^d)}}, G_{\mu_0,2}\right) \le 4 D_{\text{KL}}\left(\frac{F_n}{\|F_n\|_{L^1(\mathbb{R}^d)}}\| G_{\mu_0,2}\right). $$

However, Theorem 1.1 does not hold when $\alpha = 0$, as mentioned in Remark 1.5, due to two fundamental reasons. First, when $\alpha = 0$ there is translation invariance and therefore, for generic $F_0 \in \mathcal{M}_+(\mathbb{R}^d)$ one cannot expect that the normalized profiles $F_n/\|F_n\|_{L^1(\mathbb{R}^d)}$ always converge to the Gaussian $\mathbf{F}_{\alpha=0}$ centered at the origin (contrarily to what happens when $\alpha > 0$). Otherwise, the centers of mass must get shifted toward the origin, thus breaking the translation invariance. Indeed, as mentioned in Corollary 2.2, the equilibria are not unique when $\alpha = 0$ (contrarily to the case $\alpha > 0$), and the center of mass of the resulting equilibria must stay equal to the initial one. Second, even if we set the initial center of mass at the origin so that we kill the translation invariance, our method of proof of Theorem 1.1 leads to estimates that blow up as $\alpha \to 0$ since we have $2k_{\alpha=0} = 1$.

**2.2. Incompatibility of Wasserstein metric with multiplicative operators.** Note that by definition (1.2), our operator $\mathcal{T}$ is the composition of the sexual resproduction operator $\mathcal{B}$ with the multiplicative operator by the survival probability $e^{-m}$. Then, it might be natural to study perturbations of the previous Lemma 2.1 including the following conservative multiplicative operator.
Consider any transference plan $\gamma \in \mathcal{P}(\mathbb{R}^d \times \mathbb{R}^d)$ between $(F_1, F_2)$ as above. Specifically, we have that $\pi_{1\#}\gamma = F_1$ and $\pi_{2\#}\gamma = F_2$. Then, we can gather both integrals as follows:

$$ \int_{\mathbb{R}^d} \mathcal{B}[F_1](x) \phi_1(x) dx + \int_{\mathbb{R}^d} \mathcal{B}[F_2](y) \phi_2(y) dy \\ = \int_{\mathbb{R}^{5d}} G(z) \left( \phi_1 \left( z + \frac{x_1+x_2}{2} \right) + \phi_2 \left( z + \frac{y_1+y_2}{2} \right) \right) \gamma(dx_1, dy_1) \gamma(dx_2, dy_2) dz. $$

Since the condition $|\phi_1(x) + \phi_2(y)| \le |x - y|^2$ is verified for all $x, y \in \mathbb{R}^d$, then we find

$$ \int_{\mathbb{R}^d} \mathcal{B}[F_1](x) \phi_1(x) dx + \int_{\mathbb{R}^d} \mathcal{B}[F_2](y) \phi_2(y) dy \\ \le \frac{1}{4} \int_{\mathbb{R}^{5d}} G(z) |(x_1+x_2) - (y_1+y_2)|^2 \gamma(dx_1, dy_1) \gamma(dx_2, dy_2) dz \\ = \frac{1}{4} \int_{\mathbb{R}^{4d}} |(x_1+x_2) - (y_1+y_2)|^2 \gamma(dx_1, dy_1) \gamma(dx_2, dy_2) \\ = \frac{1}{4} \int_{\mathbb{R}^{2d}} |x_1 - y_1|^2 \gamma(dx_1, dy_1) + \frac{1}{4} \int_{\mathbb{R}^{2d}} |x_2 - y_2|^2 \gamma(dx_2, dy_2), $$

where in the third line we have used that $G$ is a probability density and in the last line we have used the crucial fact that $F_1$ and $F_2$ share the same center of mass in order to cancel the cross-terms (otherwise the estimate would boil down to a non-expansiveness estimate). Taking supremum over $(\phi_1, \phi_2)$ and infimum over $\gamma$ and using the dual characterizations (2.2) yields the result. $\square$

Note that the fact that both $F_1$ and $F_2$ have the same center of mass has been crucially used to cancel the crossed term. Otherwise, by the Cauchy-Schwarz inequality we would merely obtain non-expansiveness:

$$ W_2(\mathcal{B}[F_1], \mathcal{B}[F_2]) \le W_2(F_1, F_2). \quad (2.3) $$

Using the contractivity property of Lemma 2.1 along with the conservation of mass and center of mass in (2.1) yields the long term dynamics of (1.1) in the special case $\alpha = 0$.

**Corollary 2.2.** Set any initial datum $F_0 \in \mathcal{M}_+(\mathbb{R}^d)$ such that $\int_{\mathbb{R}^d} |x|^2 F_0(dx) < \infty$ and consider the solution $\{F_n\}_{n \in \mathbb{N}}$ to the time-discrete problem (1.1) with $\alpha = 0$, i.e., $F_n = \mathcal{B}[F_{n-1}]$ for all $n \in \mathbb{N}$. Then,

$$ \frac{\|F_n\|_{L^1(\mathbb{R}^d)}}{\|F_{n-1}\|_{L^1(\mathbb{R}^d)}} = 1, \quad W_2\left(\frac{F_n}{\|F_n\|_{L^1(\mathbb{R}^d)}}, G_{\mu_0,2}\right) \lesssim \frac{1}{2^{n/2}}, $$

for every $n \in \mathbb{N}$, where $\mu_0 := \int_{\mathbb{R}^d} x F_0(dx)$. In particular, the set of stationary distributions under $\mathcal{B}$ is $\{\mathbf{F}_{\alpha=0}(\cdot - \mu) : \mu \in \mathbb{R}^d\}$, where $\mathbf{F}_{\alpha=0} = G_{0,2}$ will denote here on the Gaussian centered at the origin with variance equals 2 in agreement with the notation (1.12) in Theorem 1.1.

The preceding result might be regarded as the counterpart of Theorem 1.1 for $\alpha = 0$. Indeed, by Talagrand transportation inequality for a Gaussian measure [23] we obtain the relation

$$ W_2^2 \left( \frac{F_n}{\|F_n\|_{L^1(\mathbb{R}^d)}}, G_{\mu_0,2} \right) \le 4 D_{\text{KL}} \left( \frac{F_n}{\|F_n\|_{L^1(\mathbb{R}^d)}} \| G_{\mu_0,2} \right). $$

However, Theorem 1.1 does not hold when $\alpha = 0$, as mentioned in Remark 1.5, due to two fundamental reasons. First, when $\alpha = 0$ there is translation invariance and therefore, for generic $F_0 \in \mathcal{M}_+(\mathbb{R}^d)$ one cannot expect that the normalized profiles $F_n/\|F_n\|_{L^1(\mathbb{R}^d)}$ always converge to the Gaussian $\mathbf{F}_{\alpha=0}$ centered at the origin (contrarily to what happens when $\alpha > 0$). Otherwise, the centers of mass must get shifted toward the origin, thus breaking the translation invariance. Indeed, as mentioned in Corollary 2.2, the equilibria are not unique when $\alpha = 0$ (contrarily to the case $\alpha > 0$), and the center of mass of the resulting equilibria must stay equal to the initial one. Second, even if we set the initial center of mass at the origin so that we kill the translation invariance, our method of proof of Theorem 1.1 leads to estimates that blow up as $\alpha \to 0$ since we have $2k_{\alpha=0} = 1$.

**2.2. Incompatibility of Wasserstein metric with multiplicative operators.** Note that by definition (1.2), our operator $\mathcal{T}$ is the composition of the sexual resproduction operator $\mathcal{B}$ with the multiplicative operator by the survival probability $e^{-m}$. Then, it might be natural to study perturbations of the previous Lemma 2.1 including the following conservative multiplicative operator.
**Definition 2.3** (Normalization of multiplicative operator).

$$ \mathcal{M}[F] := \frac{e^{-mF}}{\|e^{-mF}\|_{\mathcal{M}(\mathbb{R}^d)}}, \quad F \in \mathcal{M}_+(\mathbb{R}^d) \setminus \{0\}. $$

However, the latter is not Lipschitz continuous with respect to the quadratic Wasserstein metric. Hence the composition of $\mathcal{B}$ and $\mathcal{M}$ is not expected to be contractive, even in the case of weak selection, without any additional restriction. We illustrate such a Lipschitz discontinuity of $\mathcal{M}$ in the following example.

**Example 2.4.** Suppose that the $m \in C_+^1(\mathbb{R})$ is radially symmetric around the origin and set $F_1, F_2 \in P_2(\mathbb{R})$ as the sum of two Dirac masses, $F_1$ being symmetric, and $F_2$ nearly symmetric. More precisely,

$$ F_1 = \frac{1}{2}\delta_{-h} + \frac{1}{2}\delta_h, \quad F_2 = \frac{1}{2}\delta_{-h+\epsilon} + \frac{1}{2}\delta_{h+\epsilon}, $$

where $h > 0$ is fixed so that $m'(h) \neq 0$ and $\epsilon > 0$ is small. Note that

$$ \mathcal{M}[F_1] = F_1, \quad \mathcal{M}[F_2] = (1 - p_{\epsilon})\delta_{-h+\epsilon} + p_{\epsilon}\delta_{h+\epsilon}, $$

where $p_{\epsilon} \in (0, 1)$ is given by

$$ p_{\epsilon} = \frac{e^{-m(h+\epsilon)}}{e^{-m(-h+\epsilon)} + e^{-m(h+\epsilon)}}. $$

On the one hand, we have $W_2(F_1, F_2) = \epsilon$ because $F_2$ is simply deduced from $F_1$ by a translation of size $\epsilon$. On the other hand, assuming $m'(h) > 0$ for simplicity (a similar argument holds if $m'(h) < 0$) and taking $\epsilon > 0$ sufficiently small we obtain $p_{\epsilon} < \frac{1}{2}$ and

$$ W_2(\mathcal{M}[F_1], \mathcal{M}[F_2]) = \sqrt{\epsilon^2 + 4h(h - \epsilon)\left(\frac{1}{2} - p_{\epsilon}\right)} \sim \epsilon + h m'(h)^{1/2} \epsilon^{1/2}, $$

as $\epsilon \to 0$ by symmetry and the mean value theorem. In a sense, the leading order term corresponds to the cost of moving a piece of mass $\frac{1}{2} - p_{\epsilon}$ from $-h + \epsilon$ to $h$ in order to equilibrate the Dirac masses in the transport plan. Since $m'(h) \neq 0$ this leads to the Lipschitz-discontinuity of the operator $\mathcal{M}$.

In [21], this case is ruled out by assuming that the densities are uniformly bounded below on compact intervals. In our case, we will avoid relying on the above quantification in Lemma 2.1 for the contractivity of $\mathcal{B}$ and will propose a different strategy that we discuss in the sequel.

**2.3. Log-Lipschitz contraction estimate.** In this paper, we explore an alternative approach to derive some suitable contraction of the operator $\mathcal{T}$ in appropriate norms in the regime of strong, but quadratic, selection. As anticipated in Remark 1.4, our computations suggest using a log-Lipschitz norm, which measures the uniform deviation of tails of any profile $F$ relative to the Gaussian tails of $F_{\alpha}$. More specifically, we obtain the following result for any $\alpha \ge 0$ (but only useful when $\alpha > 0$).

**Lemma 2.5** (Log-Lipschitz estimate). Let $m$ be the quadratic selection function in (1.5) and consider any value $\alpha \ge 0$. Then, the following estimate holds true

$$ \left\| \nabla \log \frac{T[F]}{F_{\alpha}} \right\|_{L^{\infty}(\mathbb{R}^d)} \le 2 k_{\alpha} \left\| \nabla \log \frac{F}{F_{\alpha}} \right\|_{L^{\infty}(\mathbb{R}^d)}, $$

for any $F \in L_+^1(\mathbb{R}^d) \cap C^1(\mathbb{R}^d)$ such that $\|\nabla \log \frac{F}{F_{\alpha}}\|_{L^{\infty}(\mathbb{R}^d)} < \infty$. The coefficient $k_{\alpha}$ above is given by the numerical value (1.15) in Theorem 1.1.

When applied to $F = F_n$ for any solution $\{F_n\}_{n \in \mathbb{N}}$ of the time-evolution problem (1.1), the above inequality in Lemma 2.5 allows propagating a control of growth/decrease of the log-Lipschitz norm under the flow of the equation. In particular, when $\alpha > 0$, we have $2k_{\alpha} < 1$ (cf. Figure 1), and therefore we obtain an actual contractivity estimate. There is one main drawback though: typically $\nabla \log \frac{F_n}{F_{\alpha}}$ are genuinely unbounded unless $F_0$ and $F_{\alpha}$ have the same Gaussian decay for large $x$ (which is much too restrictive). A major point of our work will precisely be to circumvent this unbounded factors.

However, when $\alpha = 0$ we have $2k_{\alpha} = 1$, and we simply obtain non-expansiveness. We emphasize that in the derivation of such a rough estimate, we do not use carefully the preservation of the center of mass. This results in a non-expansiveness estimate, which is conceptually no better than the previous non-expansiveness property (2.3) in the quadratic Wasserstein metric, obtained as in Lemma 2.1 without assuming that the centers of mass are the same. For this reason, our method of proof will finally not yield satisfying results in the case $\alpha = 0$, but we find this idea illuminating to address the case of a non-trivial effect of selection when $\alpha > 0$. A refinement of Lemma 2.5 leading to real contractivity would require tackling more carefully the center of mass. However, as mentioned in Remark 1.5, we will not address this in the current paper, and we refer to Section 8 for some perspectives and future works in this line.
**Definition 2.3** (Normalization of multiplicative operator).

$$ \mathcal{M}[F] := \frac{e^{-mF}}{\|e^{-mF}\|_{\mathcal{M}(\mathbb{R}^d)}}, \quad F \in \mathcal{M}_+(\mathbb{R}^d) \setminus \{0\}. $$

However, the latter is not Lipschitz continuous with respect to the quadratic Wasserstein metric. Hence the composition of $\mathcal{B}$ and $\mathcal{M}$ is not expected to be contractive, even in the case of weak selection, without any additional restriction. We illustrate such a Lipschitz discontinuity of $\mathcal{M}$ in the following example.

**Example 2.4.** Suppose that the $m \in C_+^1(\mathbb{R})$ is radially symmetric around the origin and set $F_1, F_2 \in P_2(\mathbb{R})$ as the sum of two Dirac masses, $F_1$ being symmetric, and $F_2$ nearly symmetric. More precisely,

$$ F_1 = \frac{1}{2}\delta_{-h} + \frac{1}{2}\delta_h, \quad F_2 = \frac{1}{2}\delta_{-h+\epsilon} + \frac{1}{2}\delta_{h+\epsilon}, $$

where $h > 0$ is fixed so that $m'(h) \neq 0$ and $\epsilon > 0$ is small. Note that

$$ \mathcal{M}[F_1] = F_1, \quad \mathcal{M}[F_2] = (1 - p_{\epsilon})\delta_{-h+\epsilon} + p_{\epsilon}\delta_{h+\epsilon}, $$

where $p_{\epsilon} \in (0, 1)$ is given by

$$ p_{\epsilon} = \frac{e^{-m(h+\epsilon)}}{e^{-m(-h+\epsilon)} + e^{-m(h+\epsilon)}}. $$

On the one hand, we have $W_2(F_1, F_2) = \epsilon$ because $F_2$ is simply deduced from $F_1$ by a translation of size $\epsilon$. On the other hand, assuming $m'(h) > 0$ for simplicity (a similar argument holds if $m'(h) < 0$) and taking $\epsilon > 0$ sufficiently small we obtain $p_{\epsilon} < \frac{1}{2}$ and

$$ W_2(\mathcal{M}[F_1], \mathcal{M}[F_2]) = \sqrt{\epsilon^2 + 4h(h - \epsilon)\left(\frac{1}{2} - p_{\epsilon}\right)} \sim \epsilon + h m'(h)^{1/2} \epsilon^{1/2}, $$

as $\epsilon \to 0$ by symmetry and the mean value theorem. In a sense, the leading order term corresponds to the cost of moving a piece of mass $\frac{1}{2} - p_{\epsilon}$ from $-h + \epsilon$ to $h$ in order to equilibrate the Dirac masses in the transport plan. Since $m'(h) \neq 0$ this leads to the Lipschitz-discontinuity of the operator $\mathcal{M}$.

In [21], this case is ruled out by assuming that the densities are uniformly bounded below on compact intervals. In our case, we will avoid relying on the above quantification in Lemma 2.1 for the contractivity of $\mathcal{B}$ and will propose a different strategy that we discuss in the sequel.

**2.3. Log-Lipschitz contraction estimate.** In this paper, we explore an alternative approach to derive some suitable contraction of the operator $\mathcal{T}$ in appropriate norms in the regime of strong, but quadratic, selection. As anticipated in Remark 1.4, our computations suggest using a log-Lipschitz norm, which measures the uniform deviation of tails of any profile $F$ relative to the Gaussian tails of $F_{\alpha}$. More specifically, we obtain the following result for any $\alpha \ge 0$ (but only useful when $\alpha > 0$).

**Lemma 2.5** (Log-Lipschitz estimate). Let $m$ be the quadratic selection function in (1.5) and consider any value $\alpha \ge 0$. Then, the following estimate holds true

$$ \left\| \nabla \log \frac{T[F]}{F_{\alpha}} \right\|_{L^{\infty}(\mathbb{R}^d)} \le 2 k_{\alpha} \left\| \nabla \log \frac{F}{F_{\alpha}} \right\|_{L^{\infty}(\mathbb{R}^d)}, $$

for any $F \in L_+^1(\mathbb{R}^d) \cap C^1(\mathbb{R}^d)$ such that $\|\nabla \log \frac{F}{F_{\alpha}}\|_{L^{\infty}(\mathbb{R}^d)} < \infty$. The coefficient $k_{\alpha}$ above is given by the numerical value (1.15) in Theorem 1.1.

When applied to $F = F_n$ for any solution $\{F_n\}_{n \in \mathbb{N}}$ of the time-evolution problem (1.1), the above inequality in Lemma 2.5 allows propagating a control of growth/decrease of the log-Lipschitz norm under the flow of the equation. In particular, when $\alpha > 0$, we have $2k_{\alpha} < 1$ (cf. Figure 1), and therefore we obtain an actual contractivity estimate. There is one main drawback though: typically $\nabla \log \frac{F_n}{F_{\alpha}}$ are genuinely unbounded unless $F_0$ and $F_{\alpha}$ have the same Gaussian decay for large $x$ (which is much too restrictive). A major point of our work will precisely be to circumvent this unbounded factors.

However, when $\alpha = 0$ we have $2k_{\alpha} = 1$, and we simply obtain non-expansiveness. We emphasize that in the derivation of such a rough estimate, we do not use carefully the preservation of the center of mass. This results in a non-expansiveness estimate, which is conceptually no better than the previous non-expansiveness property (2.3) in the quadratic Wasserstein metric, obtained as in Lemma 2.1 without assuming that the centers of mass are the same. For this reason, our method of proof will finally not yield satisfying results in the case $\alpha = 0$, but we find this idea illuminating to address the case of a non-trivial effect of selection when $\alpha > 0$. A refinement of Lemma 2.5 leading to real contractivity would require tackling more carefully the center of mass. However, as mentioned in Remark 1.5, we will not address this in the current paper, and we refer to Section 8 for some perspectives and future works in this line.
<figure><img src="image_2.png" /><figcaption>FIGURE 1. Parameter $k_\alpha$ against parameter $\alpha$.</figcaption></figure>

*Proof of Lemma 2.5.* Our starting point resides in the following normalization of the operator $\mathcal{T}$:

$$ \frac{\mathcal{T}[F](x)}{\mathbf{F}_{\alpha}(x)} = \frac{\lambda_{\alpha}}{\int_{\mathbb{R}^d} \frac{F(x')}{\mathbf{F}_{\alpha}(x')} \mathbf{F}_{\alpha}(x') dx'} \iint_{\mathbb{R}^{2d}} \mathbf{P}(x; x_1, x_2) \frac{F(x_1)}{\mathbf{F}_{\alpha}(x_1)} \frac{F(x_2)}{\mathbf{F}_{\alpha}(x_2)} dx_1 dx_2, \quad (2.4) $$

for all $x \in \mathbb{R}^d$, where we have

$$ \mathbf{P}(x; x_1, x_2) = \frac{e^{\frac{1}{2}(\frac{1}{\sigma_{\alpha}^2}-\alpha)|x|^2}}{\lambda_{\alpha}(2\pi\sigma_{\alpha})^d} \exp \left[ -\frac{1}{2} \left| x - \frac{x_1+x_2}{2} \right|^2 - \frac{1}{2\sigma_{\alpha}^2} (|x_1|^2 + |x_2|^2) \right]. \quad (2.5) $$

Above we have exploited the explicit Gaussian shape of $\mathbf{F}_\alpha$ to find $\mathbf{P}$ explicitly. Note that $\mathbf{P}$ consists in a one-step Markov transition kernel representing the probability that parental traits $(x_1, x_2)$ lead to a descendant trait $x$. Indeed, $\iint_{\mathbb{R}^{2d}} \mathbf{P}(x; x_1, x_2) dx_1 dx_2 = 1$ for all $x \in \mathbb{R}^d$ because $(\lambda_\alpha, \mathbf{F}_\alpha)$ solves the non-linear eigenproblem (1.11). We remark that the quadratic form in the exponential in (2.5) reaches its maximum value at $(x_1, x_2) = (kx, kx)$ for $k = \sigma_\alpha^2 / (2 + \sigma_\alpha^2)$ and by definition (1.15) of $k_\alpha$, we infer $k = k_\alpha$. This motivates using the change of variables

$$ x_1 = k_{\alpha} x + y_1, \quad x_2 = k_{\alpha} x + y_2, \quad (2.6) $$

which appropriately centers the quadratic form at the minimum. Specifically, we obtain

$$ \begin{aligned} & \frac{1}{2} \left| x - \frac{x_1 + x_2}{2} \right|^2 - \frac{1}{2\sigma_{\alpha}^2} (|x_1|^2 + |x_2|^2) \\ &= -\frac{1}{2} \left| \frac{y_1 + y_2}{2} \right|^2 - \frac{1}{2\sigma_{\alpha}^2} (|y_1|^2 + |y_2|^2) - \frac{1}{2} \left( (1 - k_{\alpha})^2 + \frac{2k_{\alpha}^2}{\sigma_{\alpha}^2} \right) |x|^2 + \frac{1}{2} \left( (1 - k_{\alpha}) - \frac{2k_{\alpha}}{\sigma_{\alpha}^2} \right) x \cdot (y_1 + y_2), \\ &= -\frac{1}{2} \left| \frac{y_1 + y_2}{2} \right|^2 - \frac{1}{2\sigma_{\alpha}^2} (|y_1|^2 + |y_2|^2) - \frac{1}{2} \left( \frac{1}{\sigma_{\alpha}^2} - \alpha \right) |x|^2, \end{aligned} $$

where in the second line we have noticed that the coefficient of the crossed term $x \cdot (y_1 + y_2)$ vanishes due to the relation $k_\alpha = \frac{\sigma_\alpha^2}{2+\sigma_\alpha^2}$, and the coefficient of the $|x|^2$ factor can be reformulated as

$$ (1 - k_{\alpha})^2 + \frac{2k_{\alpha}^2}{\sigma_{\alpha}^2} = \frac{2}{2 + \sigma_{\alpha}^2} = \frac{1}{1 + \frac{\sigma_{\alpha}^2}{2}} = \frac{1}{\sigma_{\alpha}^2} - \alpha, $$

thanks to the implicit equation (1.13) satisfied by $\sigma_\alpha^2$. Therefore, (2.4)-(2.5) transform into

$$ \frac{\mathcal{T}[F](x)}{\mathbf{F}_{\alpha}(x)} = \frac{\lambda_{\alpha}}{\int_{\mathbb{R}^d} \frac{F(x')}{\mathbf{F}_{\alpha}(x')} \mathbf{F}_{\alpha}(x') dx'} \iint_{\mathbb{R}^{2d}} \tilde{\mathbf{P}}(y_1, y_2) \frac{F(\mathbf{k}_{\alpha} x + y_1)}{\mathbf{F}_{\alpha}(\mathbf{k}_{\alpha} x + y_1)} \frac{F(\mathbf{k}_{\alpha} x + y_2)}{\mathbf{F}_{\alpha}(\mathbf{k}_{\alpha} x + y_2)} dy_1 dy_2, \quad (2.7) $$

$$ \tilde{\mathbf{P}}(y_1, y_2) = \frac{1}{\lambda_{\alpha}(2\pi\sigma_{\alpha})^d} \exp \left[ -\frac{1}{2} \left| \frac{y_1 + y_2}{2} \right|^2 - \frac{1}{2\sigma_{\alpha}^2} (|y_1|^2 + |y_2|^2) \right]. \quad (2.8) $$

We remark that the new Markov transition kernel $\tilde{\mathbf{P}} = \tilde{\mathbf{P}}(y_1, y_2)$ does not depend on $x$ thanks to the explicit cancellation of the $|x|^2$ dependent factors. At this level, we start to observe the ergodicity
<figure><img src="image_2.png" /><figcaption>FIGURE 1. Parameter $k_\alpha$ against parameter $\alpha$.</figcaption></figure>

*Proof of Lemma 2.5.* Our starting point resides in the following normalization of the operator $\mathcal{T}$:

$$ \frac{\mathcal{T}[F](x)}{\mathbf{F}_{\alpha}(x)} = \frac{\lambda_{\alpha}}{\int_{\mathbb{R}^d} \frac{F(x')}{\mathbf{F}_{\alpha}(x')} \mathbf{F}_{\alpha}(x') dx'} \iint_{\mathbb{R}^{2d}} \mathbf{P}(x; x_1, x_2) \frac{F(x_1)}{\mathbf{F}_{\alpha}(x_1)} \frac{F(x_2)}{\mathbf{F}_{\alpha}(x_2)} dx_1 dx_2, \quad (2.4) $$

for all $x \in \mathbb{R}^d$, where we have

$$ \mathbf{P}(x; x_1, x_2) = \frac{e^{\frac{1}{2}(\frac{1}{\sigma_{\alpha}^2} - \alpha)|x|^2}}{\lambda_{\alpha}(2\pi \sigma_{\alpha})^d} \exp \left[ -\frac{1}{2} \left| x - \frac{x_1 + x_2}{2} \right|^2 - \frac{1}{2\sigma_{\alpha}^2} (|x_1|^2 + |x_2|^2) \right]. \quad (2.5) $$

Above we have exploited the explicit Gaussian shape of $\mathbf{F}_{\alpha}$ to find $\mathbf{P}$ explicitly. Note that $\mathbf{P}$ consists in a one-step Markov transition kernel representing the probability that parental traits $(x_1, x_2)$ lead to a descendant trait $x$. Indeed, $\iint_{\mathbb{R}^{2d}} \mathbf{P}(x; x_1, x_2) dx_1 dx_2 = 1$ for all $x \in \mathbb{R}^d$ because $(\lambda_{\alpha}, \mathbf{F}_{\alpha})$ solves the non-linear eigenproblem (1.11). We remark that the quadratic form in the exponential in (2.5) reaches its maximum value at $(x_1, x_2) = (kx, kx)$ for $k = \sigma_{\alpha}^2 / (2 + \sigma_{\alpha}^2)$ and by definition (1.15) of $k_{\alpha}$, we infer $k = k_{\alpha}$. This motivates using the change of variables

$$ x_1 = k_{\alpha} x + y_1, \quad x_2 = k_{\alpha} x + y_2, \quad (2.6) $$

which appropriately centers the quadratic form at the minimum. Specifically, we obtain

$$ \begin{aligned} & \frac{1}{2} \left| x - \frac{x_1 + x_2}{2} \right|^2 - \frac{1}{2\sigma_{\alpha}^2} (|x_1|^2 + |x_2|^2) \\ &= -\frac{1}{2} \left| \frac{y_1 + y_2}{2} \right|^2 - \frac{1}{2\sigma_{\alpha}^2} (|y_1|^2 + |y_2|^2) - \frac{1}{2} \left( (1 - k_{\alpha})^2 + \frac{2k_{\alpha}^2}{\sigma_{\alpha}^2} \right) |x|^2 + \frac{1}{2} \left( (1 - k_{\alpha}) - \frac{2k_{\alpha}}{\sigma_{\alpha}^2} \right) x \cdot (y_1 + y_2), \\ &= -\frac{1}{2} \left| \frac{y_1 + y_2}{2} \right|^2 - \frac{1}{2\sigma_{\alpha}^2} (|y_1|^2 + |y_2|^2) - \frac{1}{2} \left( \frac{1}{\sigma_{\alpha}^2} - \alpha \right) |x|^2, \end{aligned} $$

where in the second line we have noticed that the coefficient of the crossed term $x \cdot (y_1 + y_2)$ vanishes due to the relation $k_{\alpha} = \frac{\sigma_{\alpha}^2}{2+\sigma_{\alpha}^2}$, and the coefficient of the $|x|^2$ factor can be reformulated as

$$ (1 - k_{\alpha})^2 + \frac{2k_{\alpha}^2}{\sigma_{\alpha}^2} = \frac{2}{2 + \sigma_{\alpha}^2} = \frac{1}{1 + \frac{\sigma_{\alpha}^2}{2}} = \frac{1}{\sigma_{\alpha}^2} - \alpha, $$

thanks to the implicit equation (1.13) satisfied by $\sigma_{\alpha}^2$. Therefore, (2.4)-(2.5) transform into

$$ \frac{\mathcal{T}[F](x)}{\mathbf{F}_{\alpha}(x)} = \frac{\lambda_{\alpha}}{\int_{\mathbb{R}^d} \frac{F(x')}{\mathbf{F}_{\alpha}(x')} \mathbf{F}_{\alpha}(x') dx'} \iint_{\mathbb{R}^{2d}} \tilde{\mathbf{P}}(y_1, y_2) \frac{F(\mathbf{k}_{\alpha} x + y_1)}{\mathbf{F}_{\alpha}(\mathbf{k}_{\alpha} x + y_1)} \frac{F(\mathbf{k}_{\alpha} x + y_2)}{\mathbf{F}_{\alpha}(\mathbf{k}_{\alpha} x + y_2)} dy_1 dy_2, \quad (2.7) $$

$$ \tilde{\mathbf{P}}(y_1, y_2) = \frac{1}{\lambda_{\alpha}(2\pi \sigma_{\alpha})^d} \exp \left[ -\frac{1}{2} \left| \frac{y_1 + y_2}{2} \right|^2 - \frac{1}{2\sigma_{\alpha}^2} (|y_1|^2 + |y_2|^2) \right]. \quad (2.8) $$

We remark that the new Markov transition kernel $\tilde{\mathbf{P}} = \tilde{\mathbf{P}}(y_1, y_2)$ does not depend on $x$ thanks to the explicit cancellation of the $|x|^2$ dependent factors. At this level, we start to observe the ergodicity
phenomenon since the dependence of $x$ on the right hand side of (2.7) has shrunk by a factor $\mathbf{k}_\alpha$. A possible strategy to see if there is a quantitative degradation of the dependence on $x$ is to take logarithmic derivatives and try to relate the log-Lipschitz norms of $\mathcal{T}[F]$ and $F$. Specifically, we have

$$ \nabla \log \frac{\mathcal{T}[F](x)}{\mathbf{F}_{\alpha}(x)} = \mathbf{k}_{\alpha} \int_{\mathbb{R}^{2d}} \left( \nabla \log \frac{F(\mathbf{k}_{\alpha} x + y_1)}{\mathbf{F}_{\alpha}(\mathbf{k}_{\alpha} x + y_1)} + \nabla \log \frac{F(\mathbf{k}_{\alpha} x + y_2)}{\mathbf{F}_{\alpha}(\mathbf{k}_{\alpha} x + y_2)} \right) \nu(x; dy_1, dy_2), \quad (2.9) $$

where the $x$-dependent measures $\nu(x; dy_1, dy_2)$ on the variables $(y_1, y_2)$ have the following density with respect to the Lebesgue measure:

$$ \frac{\nu(x; dy_1, dy_2)}{dy_1 dy_2} = \frac{\tilde{\mathbf{P}}(y_1, y_2) \frac{F(\mathbf{k}_\alpha x + y_1)}{\mathbf{F}_\alpha(\mathbf{k}_\alpha x + y_1)} \frac{F(\mathbf{k}_\alpha x + y_2)}{\mathbf{F}_\alpha(\mathbf{k}_\alpha x + y_2)}}{\iint_{\mathbb{R}^{2d}} \tilde{\mathbf{P}}(y'_1, y'_2) \frac{F(\mathbf{k}_\alpha x + y'_1)}{\mathbf{F}_\alpha(\mathbf{k}_\alpha x + y'_1)} \frac{F(\mathbf{k}_\alpha x + y'_2)}{\mathbf{F}_\alpha(\mathbf{k}_\alpha x + y'_2)} dy'_1 dy'_2} $$

Since the integrands of (2.9) are uniformly bounded by our assumptions, we end the proof by taking $L^\infty$ bounds and using that $\nu$ are probability measures on the variables $(y_1, y_2)$. $\square$

## 2.4. Brief description of our strategy.
As advanced before, our strategy is based on a finer understanding of the iterations of $\mathcal{T}$ across generations. Specifically, it relies on a suitable reformulation of the solutions $\{F_n\}_{n \in \mathbb{N}}$ solving the recursion (1.1) for the special quadratic selection function $m$ in (1.5). At this stage we intentionally keep notation simple and intuitive, since our goal is to briefly present the main strategy. However, a rigorous approach with more descriptive notation for the trees of ancestors, which arise from the recursion, is developed in detail in Section 5.

The first step consists in choosing the appropriate normalization extending the previous normalization $F_n/\mathbf{F}_\alpha$ in the proof of Lemma 2.5. In this paper, we have opted for

$$ \bar{F}_n = \frac{e^m F_n}{\mathbf{F}_{\alpha=0}}, $$

(cf. Definition 5.1) but there is some freedom here. In particular, the term $e^m$ is not mandatory, but it is convenient as such for an easier sorting of the various terms. Indeed, we may typically consider any normalization $F_n/G_{0,\sigma^2}$ with $\sigma^2$ larger but arbitrarily close to $\sigma_\alpha^2$ (being $\sigma_\alpha^2$ the variance (1.13) of the expected equilibrium $\mathbf{F}_\alpha$). In other words, we could take any “close-to-optimal” normalization as compared to the “optimal” normalization used in the strategy in Section 2.3. However, for general selection functions $m$ we do not know the expected equilibrium. To account for a more robust viewpoint, we follow a more robust approach so that we do not “optimize” the normalization.

In a nutshell: we iterate the operator $\mathcal{T}$ in the recursion for $F_n$ up to the initial datum and we obtain (modulo a multiplicative constant) a formula of the following type

$$ F_n(x) \propto \int_{\mathbb{R}^{(2^{n+1}-2)d}} \mathbf{P}_n(x; \{x_i\}) \prod_j \bar{F}_0(x_j) d\{x_i\}, $$

for an explicit $n$-step transition kernel $\mathbf{P}_n$ with Gaussian shape:

$$ \mathbf{P}_n(x; \{x_i\}) := \exp(-\mathbf{Q}_n(x; \{x_i\})). $$

Here, the index $i$ spans the $2^{n+1}-2$ members of the genealogical tree with $n$ generations of ancestors from the root (excluded) to the leaves, whilst $j$ only spans the $2^n$ members at the leaves. In addition, $\mathbf{Q}_n = \mathbf{Q}_n(x; \{x_i\})$ is a quadratic form taking as arguments all traits $x$ and $\{x_i\}$ (i.e., all the $2^{n+1}-1$ ancestors in the family chart including the root $x$). As compared to Section 2.3, we have iterated $n$ times, which raises the dimension of the quadratic form to $2^{n+1}-1$. In addition, the change of normalization leads to additional quadratic contributions $\frac{\gamma}{2}|x_i|^2$ (with $\gamma = \frac{\sigma_\alpha^2}{2+\sigma_\alpha^2}$)). Therefore, the minimizers of $\mathbf{Q}_n$ get shifted (cf. Section 5). As a result of successive backwards changes of variables from leaves to the root very much in the spirit of (2.6), the following alternative expression is obtained:

$$ F_n(x) \propto \tilde{\mathbf{F}}_n(x) \int_{\mathbb{R}^{(2^{n+1}-2)d}} \tilde{\mathbf{P}}_n(\{y_i\}) \prod_j \bar{F}_0(\Phi_n^j(x; \{y_i\})) d\{y_i\}. \quad (2.10) $$

Above, for every leaf $j$ the maps $\Phi_n^j$ defined by

$$ \Phi_n^j(x; \{y_i\}) := \kappa_n x + \Lambda^j(\{y_i\}), $$

where $\Lambda_n^j$ are affine transformations with respect to the variables $\{y_i\}$. We shall call them the lineage maps (cf. Definition 5.6), since they contain precise information of the effect of the leaf $j$ of the genealogical
phenomenon since the dependence of $x$ on the right hand side of (2.7) has shrunk by a factor $\mathbf{k}_\alpha$. A possible strategy to see if there is a quantitative degradation of the dependence on $x$ is to take logarithmic derivatives and try to relate the log-Lipschitz norms of $\mathcal{T}[F]$ and $F$. Specifically, we have

$$ \nabla \log \frac{\mathcal{T}[F](x)}{\mathbf{F}_{\alpha}(x)} = \mathbf{k}_{\alpha} \int_{\mathbb{R}^{2d}} \left( \nabla \log \frac{F(\mathbf{k}_{\alpha} x + y_1)}{\mathbf{F}_{\alpha}(\mathbf{k}_{\alpha} x + y_1)} + \nabla \log \frac{F(\mathbf{k}_{\alpha} x + y_2)}{\mathbf{F}_{\alpha}(\mathbf{k}_{\alpha} x + y_2)} \right) \nu(x; dy_1, dy_2), \quad (2.9) $$

where the $x$-dependent measures $\nu(x; dy_1, dy_2)$ on the variables $(y_1, y_2)$ have the following density with respect to the Lebesgue measure:

$$ \frac{\nu(x; dy_1, dy_2)}{dy_1 dy_2} = \frac{\tilde{\mathbf{P}}(y_1, y_2) \frac{F(\mathbf{k}_\alpha x + y_1)}{\mathbf{F}_\alpha(\mathbf{k}_\alpha x + y_1)} \frac{F(\mathbf{k}_\alpha x + y_2)}{\mathbf{F}_\alpha(\mathbf{k}_\alpha x + y_2)}}{\iint_{\mathbb{R}^{2d}} \tilde{\mathbf{P}}(y'_1, y'_2) \frac{F(\mathbf{k}_\alpha x + y'_1)}{\mathbf{F}_\alpha(\mathbf{k}_\alpha x + y'_1)} \frac{F(\mathbf{k}_\alpha x + y'_2)}{\mathbf{F}_\alpha(\mathbf{k}_\alpha x + y'_2)} dy'_1 dy'_2} $$

Since the integrands of (2.9) are uniformly bounded by our assumptions, we end the proof by taking $L^\infty$ bounds and using that $\nu$ are probability measures on the variables $(y_1, y_2)$. $\square$

## 2.4. Brief description of our strategy.
As advanced before, our strategy is based on a finer understanding of the iterations of $\mathcal{T}$ across generations. Specifically, it relies on a suitable reformulation of the solutions $\{F_n\}_{n \in \mathbb{N}}$ solving the recursion (1.1) for the special quadratic selection function $m$ in (1.5). At this stage we intentionally keep notation simple and intuitive, since our goal is to briefly present the main strategy. However, a rigorous approach with more descriptive notation for the trees of ancestors, which arise from the recursion, is developed in detail in Section 5.

The first step consists in choosing the appropriate normalization extending the previous normalization $F_n/\mathbf{F}_\alpha$ in the proof of Lemma 2.5. In this paper, we have opted for

$$ \bar{F}_n = \frac{e^m F_n}{\mathbf{F}_{\alpha=0}}, $$

(cf. Definition 5.1) but there is some freedom here. In particular, the term $e^m$ is not mandatory, but it is convenient as such for an easier sorting of the various terms. Indeed, we may typically consider any normalization $F_n/G_{0,\sigma^2}$ with $\sigma^2$ larger but arbitrarily close to $\sigma_\alpha^2$ (being $\sigma_\alpha^2$ the variance (1.13) of the expected equilibrium $\mathbf{F}_\alpha$). In other words, we could take any “close-to-optimal” normalization as compared to the “optimal” normalization used in the strategy in Section 2.3. However, for general selection functions $m$ we do not know the expected equilibrium. To account for a more robust viewpoint, we follow a more robust approach so that we do not “optimize” the normalization.

In a nutshell: we iterate the operator $\mathcal{T}$ in the recursion for $F_n$ up to the initial datum and we obtain (modulo a multiplicative constant) a formula of the following type

$$ F_n(x) \propto \int_{\mathbb{R}^{(2^{n+1}-2)d}} \mathbf{P}_n(x; \{x_i\}) \prod_j \bar{F}_0(x_j) d\{x_i\}, $$

for an explicit $n$-step transition kernel $\mathbf{P}_n$ with Gaussian shape:

$$ \mathbf{P}_n(x; \{x_i\}) := \exp(-\mathbf{Q}_n(x; \{x_i\})). $$

Here, the index $i$ spans the $2^{n+1}-2$ members of the genealogical tree with $n$ generations of ancestors from the root (excluded) to the leaves, whilst $j$ only spans the $2^n$ members at the leaves. In addition, $\mathbf{Q}_n = \mathbf{Q}_n(x; \{x_i\})$ is a quadratic form taking as arguments all traits $x$ and $\{x_i\}$ (i.e., all the $2^{n+1}-1$ ancestors in the family chart including the root $x$). As compared to Section 2.3, we have iterated $n$ times, which raises the dimension of the quadratic form to $2^{n+1}-1$. In addition, the change of normalization leads to additional quadratic contributions $\frac{\gamma}{2}|x_i|^2$ (with $\gamma = \frac{\sigma_\alpha^2}{2+\sigma_\alpha^2}$)). Therefore, the minimizers of $\mathbf{Q}_n$ get shifted (cf. Section 5). As a result of successive backwards changes of variables from leaves to the root very much in the spirit of (2.6), the following alternative expression is obtained:

$$ F_n(x) \propto \tilde{\mathbf{F}}_n(x) \int_{\mathbb{R}^{(2^{n+1}-2)d}} \tilde{\mathbf{P}}_n(\{y_i\}) \prod_j \bar{F}_0(\Phi_n^j(x; \{y_i\})) d\{y_i\}. \quad (2.10) $$

Above, for every leaf $j$ the maps $\Phi_n^j$ defined by

$$ \Phi_n^j(x; \{y_i\}) := \kappa_n x + \Lambda^j(\{y_i\}), $$

where $\Lambda_n^j$ are affine transformations with respect to the variables $\{y_i\}$. We shall call them the lineage maps (cf. Definition 5.6), since they contain precise information of the effect of the leaf $j$ of the genealogical
tree on the resulting trait $x$. In addition, $P_n$ are $n$-step transition kernel with Gaussian shape:

$$ \tilde{P}_n(\{y_i\}) = \exp(-\tilde{Q}_n(\{y_i\})), $$

for a $x$-independent quadratic form $\tilde{Q}_n$ taking as arguments all traits $\{y_i\}$ (root excluded). In this case, the $|x|^2$ remainders coming from the changes of variables on $P_n$ do not simplify due to our different normalization, and they contribute with an explicit Gaussian factor $\tilde{F}_n(x)$. Then, the dependency upon $x$ is split into two parts: the explicit term $\tilde{F}_n(x)$ outside the integral (which is proven to converge towards $F_\alpha$ by construction) and the contribution $\kappa_n x$ of the lineage maps $\Phi_n^j(x; \{y_i\})$ inside $\bar{F}_0$ for an appropriate sequence $\{\kappa_n\}_{n \in \mathbb{N}} \subseteq \mathbb{R}_+$. The major observation here is that $\kappa_n \ll 2^{-n}$ when $\alpha > 0$ because of the strong shift towards the origin under selection. Otherwise, for $\alpha = 0$ we only have $\kappa_n = 2^{-n}$.

Formula (2.10) for $\alpha > 0$ then suggests a strong form of ergodicity, reminiscent of the contraction of the log-Lipschitz norm in Lemma 2.5, where the dependency on $x$ fades as $n \to \infty$. Indeed, there are still $2^n$ terms in the product indexed by $j$, but the contraction parameter $\kappa_n$ appears to decay fast enough to compensate them. To make this argument quantitative, let us differentiate (2.10) again to obtain

$$ \nabla \log F_n(x) = \nabla \log \tilde{F}_n(x) + \kappa_n \sum_j \int_{\mathbb{R}^{(2^{n+1}-2)d}} \nabla \log \bar{F}_0(\Phi_n^j(x; \{y_i\})) \nu_n(x; d\{y_i\}), \quad (2.11) $$

where the $x$-dependent probability measures $\nu_n(x; d\{y_i\})$ on the variables $\{y_i\}$ have the following density with respect to the Lebesgue measure:

$$ \frac{\nu_n(x; d\{y_i\})}{d\{y_i\}} = \frac{\tilde{P}_n(\{y_i\}) \prod_j \bar{F}_0(\Phi_n^j(x; \{y_i\}))}{\int_{\mathbb{R}^{(2^{n+1}-2)d}} \tilde{P}_n(\{y'_i\}) \prod_j \bar{F}_0(\Phi_n^j(x; \{y'_i\})) d\{y'_i\}}. $$

We remark that a naive repetition of the strategy in Lemma 2.5 under the additional assumption that $\nabla \log \bar{F}_0 \in L^\infty(\mathbb{R}^d)$ would immediately retrieve an exponential decay of the second factor in (2.11) when $\alpha > 0$ (or only a uniform bound if $\alpha = 0$) since the large sum over $j$ would be bounded by $2^n \kappa_n \|\nabla \log \bar{F}_0\|_{L^\infty(\mathbb{R}^d)}$. However, as mentioned before, there is a strong drawback: the factors $\nabla \log \bar{F}_0$ are not bounded, not even after long enough times. In the current formal argument, we have not discussed about the precise shape of the linear components $\Lambda_n^j$ of the lineage maps $\Phi_n^j$, since it was not relevant thus far. However, we anticipate that to overcome the aforementioned complication, we shall require more precise estimates for large values of $\{y_i\}$ in the high-dimensional integral, which are guaranteed by the strong enough decay of $\kappa_n$, and which will be essential in the rigorous proof in Section 6.

More specifically, and this explains why our approach moves from one-step contraction estimates (in the spirit of Section 2.3) to ergodicity results, note that the same reformulation of the iteration as above could be done exactly $n-k$ times up to an advanced enough time step $k$. By doing so, we can grasp on some natural Gaussian shaping under selection, which we expect to lead to unbounded $\nabla \log \bar{F}_k$, but growing sublinearly at infinity. Hence, we require precise compensations in the high-dimensional integral (2.11) which we find by deriving a suitable control of moments of $\nu_n$. By doing so, one has to irremediably move from uniform estimates to averaged estimates, and being able to propagate them for large times. This requires a thorough and highly technical analysis, which becomes the main objective of this paper.

For an easier readability, and to guide the reader along the various steps of the proof, we provide an overall map of it in Figure 2, which allows interconnecting the main fundamental results.

### 3. PRELIMINARIES

In this part, we collect some necessary preliminary tools and results that will be used later on.

3.1. **Perfect binary trees.** In this paper, we shall systematically use indices $i$ that do not range on discrete intervals $\{1, \dots, n\}$ but rather on the vertices of a specific type of trees, which are called perfect binary trees, see [22, 25]. First, we introduce the notion of binary trees, which we present using their universal address system, see [22, Section 11.3.3]. Namely, a binary tree $T$ consists in a finite subset of words $T \subseteq W_2 := \bigcup_{k=0}^\infty \{1, 2\}^k$ with letters in the alphabet $\{1, 2\}$ that verifies:

(i) $\emptyset \in T$.

(ii) If $i1 \in T$ or $i2 \in T$ for some word $i \in W_2$, then $i \in T$.

This implies that the root is the empty word, and $T$ is stable under chopping letters on the right of its words. Given a word $i \in T$, we denote its length (number of letters) or height in the tree by $|i|$. In particular, $|\emptyset| = 0$. A binary tree $T$ is said to be perfect if, in addition, the following properties hold:

(iii) Given a word $i \in W_2$, then $i1 \in T$ if, and only if, $i2 \in T$.
tree on the resulting trait $x$. In addition, $P_n$ are $n$-step transition kernel with Gaussian shape:

$$ \tilde{P}_n(\{y_i\}) = \exp(-\tilde{Q}_n(\{y_i\})), $$

for a $x$-independent quadratic form $\tilde{Q}_n$ taking as arguments all traits $\{y_i\}$ (root excluded). In this case, the $|x|^2$ remainders coming from the changes of variables on $P_n$ do not simplify due to our different normalization, and they contribute with an explicit Gaussian factor $\tilde{F}_n(x)$. Then, the dependency upon $x$ is split into two parts: the explicit term $\tilde{F}_n(x)$ outside the integral (which is proven to converge towards $F_\alpha$ by construction) and the contribution $\kappa_n x$ of the lineage maps $\Phi_n^j(x; \{y_i\})$ inside $\bar{F}_0$ for an appropriate sequence $\{\kappa_n\}_{n \in \mathbb{N}} \subseteq \mathbb{R}_+$. The major observation here is that $\kappa_n \ll 2^{-n}$ when $\alpha > 0$ because of the strong shift towards the origin under selection. Otherwise, for $\alpha = 0$ we only have $\kappa_n = 2^{-n}$.

Formula (2.10) for $\alpha > 0$ then suggests a strong form of ergodicity, reminiscent of the contraction of the log-Lipschitz norm in Lemma 2.5, where the dependency on $x$ fades as $n \to \infty$. Indeed, there are still $2^n$ terms in the product indexed by $j$, but the contraction parameter $\kappa_n$ appears to decay fast enough to compensate them. To make this argument quantitative, let us differentiate (2.10) again to obtain

$$ \nabla \log F_n(x) = \nabla \log \tilde{F}_n(x) + \kappa_n \sum_j \int_{\mathbb{R}^{(2^{n+1}-2)d}} \nabla \log \bar{F}_0(\Phi_n^j(x; \{y_i\})) \nu_n(x; d\{y_i\}), \quad (2.11) $$

where the $x$-dependent probability measures $\nu_n(x; d\{y_i\})$ on the variables $\{y_i\}$ have the following density with respect to the Lebesgue measure:

$$ \frac{\nu_n(x; d\{y_i\})}{d\{y_i\}} = \frac{\tilde{P}_n(\{y_i\}) \prod_j \bar{F}_0(\Phi_n^j(x; \{y_i\}))}{\int_{\mathbb{R}^{(2^{n+1}-2)d}} \tilde{P}_n(\{y'_i\}) \prod_j \bar{F}_0(\Phi_n^j(x; \{y'_i\})) d\{y'_i\}}. $$

We remark that a naive repetition of the strategy in Lemma 2.5 under the additional assumption that $\nabla \log \bar{F}_0 \in L^\infty(\mathbb{R}^d)$ would immediately retrieve an exponential decay of the second factor in (2.11) when $\alpha > 0$ (or only a uniform bound if $\alpha = 0$) since the large sum over $j$ would be bounded by $2^n \kappa_n \|\nabla \log \bar{F}_0\|_{L^\infty(\mathbb{R}^d)}$. However, as mentioned before, there is a strong drawback: the factors $\nabla \log \bar{F}_0$ are not bounded, not even after long enough times. In the current formal argument, we have not discussed about the precise shape of the linear components $\Lambda_n^j$ of the lineage maps $\Phi_n^j$, since it was not relevant thus far. However, we anticipate that to overcome the aforementioned complication, we shall require more precise estimates for large values of $\{y_i\}$ in the high-dimensional integral, which are guaranteed by the strong enough decay of $\kappa_n$, and which will be essential in the rigorous proof in Section 6.

More specifically, and this explains why our approach moves from one-step contraction estimates (in the spirit of Section 2.3) to ergodicity results, note that the same reformulation of the iteration as above could be done exactly $n-k$ times up to an advanced enough time step $k$. By doing so, we can grasp on some natural Gaussian shaping under selection, which we expect to lead to unbounded $\nabla \log \bar{F}_k$, but growing sublinearly at infinity. Hence, we require precise compensations in the high-dimensional integral (2.11) which we find by deriving a suitable control of moments of $\nu_n$. By doing so, one has to irremediably move from uniform estimates to averaged estimates, and being able to propagate them for large times. This requires a thorough and highly technical analysis, which becomes the main objective of this paper.

For an easier readability, and to guide the reader along the various steps of the proof, we provide an overall map of it in Figure 2, which allows interconnecting the main fundamental results.

### 3. PRELIMINARIES

In this part, we collect some necessary preliminary tools and results that will be used later on.

**3.1. Perfect binary trees.** In this paper, we shall systematically use indices $i$ that do not range on discrete intervals $\{1, \dots, n\}$ but rather on the vertices of a specific type of trees, which are called perfect binary trees, see [22, 25]. First, we introduce the notion of binary trees, which we present using their universal address system, see [22, Section 11.3.3]. Namely, a binary tree $T$ consists in a finite subset of words $T \subseteq W_2 := \cup_{k=0}^\infty \{1, 2\}^k$ with letters in the alphabet $\{1, 2\}$ that verifies:

(i) $\emptyset \in T$.

(ii) If $i1 \in T$ or $i2 \in T$ for some word $i \in W_2$, then $i \in T$.

This implies that the root is the empty word, and $T$ is stable under chopping letters on the right of its words. Given a word $i \in T$, we denote its length (number of letters) or height in the tree by $|i|$. In particular, $|\emptyset| = 0$. A binary tree $T$ is said to be perfect if, in addition, the following properties hold:

(iii) Given a word $i \in W_2$, then $i1 \in T$ if, and only if, $i2 \in T$.
<figure><img src="image_3.png" /><figcaption>FIGURE 2. Overall map of the proof of Theorem 1.1</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>FIGURE 3. Perfect binary tree $T^3$</figcaption></figure>

(iv) There exists $n \in \mathbb{N}$ such that $|i| \le n$ for every $i \in T$ and

$$ \#\{i \in T : |i| = n\} = 2^n. $$

This implies that, except for the root, words in $T$ appear in couples $\{i1, i2\}$ and paths in the tree arising from the root achieve the same maximal height $n$. The above four conditions determine a unique tree, which we call the *perfect binary tree with height $n \in \mathbb{N}$* and we denote by $T^n$. See Figure 3 for a graphical representation of the perfect binary tree $T^3$ of height 3. Binary trees are often used to describe the different generations of offsprings after a given individual. In our case, we shall make a reverse use of trees. Namely, a binary tree will represent the *family tree* or *pedigree chart* of an individual, consisting of the different generations of *ancestors* of such an individual. For this purpose, we shall establish the following terminology:

(i) **(Leaves)** We say that a vertex $i \in T^n$ is a leaf of the tree if $|i| = n$. The set of all leaves of $T^n$ will be denoted by

$$ L^n := \{i \in T^n : |i| = n\}. $$

(ii) **(Levels)** We say that a vertex $i \in T^n$ is on the $m$-th level of the tree if $|i| = m$ for some $1 \le m \le n$. The set of all vertices on the $m$-th level will be denoted

$$ L_m^n := \{i \in T^n : |i| = m\}. $$

In particular, note that $L_0^n = \{\emptyset\}$ is the root and $L_n^n = L^n$ are the leaves. For simplicity of notation, we denote the set of vertices of the root-free and leaves-free tree respectively by

$$ T_*^n := T^n \setminus \{\emptyset\}, \quad \hat{T}^n := T^n \setminus L^n. $$
<figure><img src="image_3.png" /><figcaption>FIGURE 2. Overall map of the proof of Theorem 1.1</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>FIGURE 3. Perfect binary tree $T^3$</figcaption></figure>

(iv) There exists $n \in \mathbb{N}$ such that $|i| \le n$ for every $i \in T$ and

$$ \#\{i \in T : |i| = n\} = 2^n. $$

This implies that, except for the root, words in $T$ appear in couples $\{i1, i2\}$ and paths in the tree arising from the root achieve the same maximal height $n$. The above four conditions determine a unique tree, which we call the *perfect binary tree with height $n \in \mathbb{N}$* and we denote by $T^n$. See Figure 3 for a graphical representation of the perfect binary tree $T^3$ of height 3. Binary trees are often used to describe the different generations of offsprings after a given individual. In our case, we shall make a reverse use of trees. Namely, a binary tree will represent the *family tree* or *pedigree chart* of an individual, consisting of the different generations of *ancestors* of such an individual. For this purpose, we shall establish the following terminology:

(i) **(Leaves)** We say that a vertex $i \in T^n$ is a leaf of the tree if $|i| = n$. The set of all leaves of $T^n$ will be denoted by

$$ L^n := \{i \in T^n : |i| = n\}. $$

(ii) **(Levels)** We say that a vertex $i \in T^n$ is on the $m$-th level of the tree if $|i| = m$ for some $1 \le m \le n$. The set of all vertices on the $m$-th level will be denoted

$$ L_m^n := \{i \in T^n : |i| = m\}. $$

In particular, note that $L_0^n = \{\emptyset\}$ is the root and $L_n^n = L^n$ are the leaves. For simplicity of notation, we denote the set of vertices of the root-free and leaves-free tree respectively by

$$ T_*^n := T^n \setminus \{\emptyset\}, \quad \hat{T}^n := T^n \setminus L^n. $$
(iii) (Child) Given a vertex $i \in \mathbb{T}_*^n$, then there exists a unique word $j \in W_2$ such that $i = j1$ or $i = j2$. Such a vertex $j$ is called the child of $i$ and will be denoted by $c(i) = j$.

(iv) (Parents) Given a vertex $\hat{\mathbb{T}}^n$, we define the parents of $i$ as the subset of all vertices that have the same common child $i$, that is,

$$P(i) := \{i1, i2\}.$$

(v) (Mate) Given a vertex $i \in \mathbb{T}_*^n$, we define the mate of $i$ and we denote it $m(i)$ as the only other vertex in $\mathbb{T}^n$ that has the same child as $i$, i.e.,

$$\{i, m(i)\} = \{c(i)1, c(i)2\}.$$

(vi) (Tree order) Given two vertices $i, j \in \mathbb{T}^n$, we say that $i \le j$ if the associated words are so ordered according to the lexicographical order of the set of words $W_2$ with two letters $\{1, 2\}$.

(vii) (Highest common descendant) Given two vertices $i, j \in \hat{\mathbb{T}}^n$, we define the highest common descendant of $i$ and $j$, and we denote it by $i \wedge j$, as

$$i \wedge j := \max\{l \in \mathbb{T}^n : l \le i, l \le j\},$$

where the maximum is considered with respect to above tree (lexicographic) order.

**Remark 3.1** (Tree-indexed variables). Given $n \in \mathbb{N}$, we shall identify $\mathbb{R}^{2(2^n-1)d} \equiv (\mathbb{R}^d)^{\mathbb{T}_*^n}$ and $\mathbb{R}^{2^d} \equiv (\mathbb{R}^d)^{\mathbb{L}^n}$. Specifically, vectors $\mathbf{x}_n, \mathbf{y}_n \in \mathbb{R}^{2(2^n-1)d}$ and $\mathbf{z}_n \in \mathbb{R}^{2^d}$ will be regarded as indexed families

$$\mathbf{x}_n = (x_i)_{i \in \mathbb{T}_*^n}, \quad \mathbf{y}_n = (y_i)_{i \in \mathbb{T}_*^n}, \quad \mathbf{z}_n = (z_j)_{j \in \mathbb{L}^n}, \qquad (3.1)$$

where $x_i, y_i \in \mathbb{R}^d$ and $z_j \in \mathbb{R}^d$ for each $i \in \mathbb{T}_*^n$ and $j \in \mathbb{L}^n$.

**3.2. Gaussian solutions.** In this part, we compute particular solutions of the time-discrete problem (1.1) and the associated eigenproblem (1.11). As we advanced before, we shall exploit the explicit algebraic structure imposed by $G$ (given by the Gaussian mixing kernel (1.4)) and $m$ (given by the quadratic selection function (1.5)). Namely, explicit Gaussian solution will be obtained. We recall first the following stability property of Gaussians under convolutions.

**Lemma 3.2** (Stability of Gaussians). The following relation holds true

$$G_{\mu_1, \sigma_1^2} * G_{\mu_2, \sigma_2^2} = G_{\mu_1 + \mu_2, \sigma_1^2 + \sigma_2^2}, \qquad (3.2)$$

for any couple of means $\mu_1, \mu_2 \in \mathbb{R}^d$ and variances $\sigma_1^2, \sigma_2^2 > 0$.

Using the above result, we can obtain the following explicit evaluation of the operator $\mathcal{T}$ in (1.2) over the class of Gaussian functions.

**Lemma 3.3** (Evaluation on Gaussians). Consider any $\mu \in \mathbb{R}^d$ and $\sigma^2 \in \mathbb{R}_+^*$, then

$$\mathcal{T}[G_{\mu, \sigma^2}] = m_* G_{\mu_*, \sigma_*^2},$$

where the parameters $m_*$, $\mu_*$ and $\sigma_*^2$ are given by:

$$m_* = \frac{e^{-\frac{1}{2} \frac{\alpha |\mu|^2}{1 + \alpha(1 + \frac{\sigma^2}{2})}}}{(1 + \alpha(1 + \frac{\sigma^2}{2}))^{d/2}}, \quad \mu_* = \frac{\mu}{1 + \alpha(1 + \frac{\sigma^2}{2})}, \quad \sigma_*^2 = \frac{1 + \frac{\sigma^2}{2}}{1 + \alpha(1 + \frac{\sigma^2}{2})}. \qquad (3.3)$$

*Proof.* On the one hand, note that by Lemma 3.2

$$\mathcal{B}[G_{\mu, \sigma^2}](x) = \left( G\left(\frac{\cdot}{2}\right) * G_{\mu, \sigma^2} * G_{\mu, \sigma^2} \right) (2x) = G_{\mu, 1 + \frac{\sigma^2}{2}}(x),$$

for each $x \in \mathbb{R}^d$. Therefore, by definition of $\mathcal{T}$ in (1.2) we obtain

$$\mathcal{T}[G_{\mu, \sigma^2}](x) = e^{-\frac{\alpha}{2}|x|^2} G_{\mu, 1 + \frac{\sigma^2}{2}}(x) = \frac{1}{(2\pi(1 + \frac{\sigma^2}{2}))^{d/2}} \exp \left( -\frac{\alpha}{2}|x|^2 - \frac{1}{2} \frac{1}{1 + \frac{\sigma^2}{2}} |x - \mu|^2 \right),$$

for each $x \in \mathbb{R}^d$. By completing the square inside the exponential, we conclude our result.

Consequently, the following explicit Gaussian solution of the eigenproblem (1.11) is found.

**Proposition 3.4** (Gaussian solution of the eigenproblem). Assume that $\alpha \in \mathbb{R}_+^*$, then the eigenproblem (1.11) has a unique Gaussian solution $(\lambda_\alpha, F_\alpha)$, determined by the relation (1.12) in Theorem 1.1.
(iii) (Child) Given a vertex $i \in \mathbb{T}_*^n$, then there exists a unique word $j \in W_2$ such that $i = j1$ or $i = j2$. Such a vertex $j$ is called the child of $i$ and will be denoted by $c(i) = j$.

(iv) (Parents) Given a vertex $\hat{\mathbb{T}}^n$, we define the parents of $i$ as the subset of all vertices that have the same common child $i$, that is,

$$P(i) := \{i1, i2\}.$$

(v) (Mate) Given a vertex $i \in \mathbb{T}_*^n$, we define the mate of $i$ and we denote it $m(i)$ as the only other vertex in $\mathbb{T}^n$ that has the same child as $i$, i.e.,

$$\{i, m(i)\} = \{c(i)1, c(i)2\}.$$

(vi) (Tree order) Given two vertices $i, j \in \mathbb{T}^n$, we say that $i \le j$ if the associated words are so ordered according to the lexicographical order of the set of words $W_2$ with two letters $\{1, 2\}$.

(vii) (Highest common descendant) Given two vertices $i, j \in \hat{\mathbb{T}}^n$, we define the highest common descendant of $i$ and $j$, and we denote it by $i \wedge j$, as

$$i \wedge j := \max\{l \in \mathbb{T}^n : l \le i, l \le j\},$$

where the maximum is considered with respect to above tree (lexicographic) order.

**Remark 3.1** (Tree-indexed variables). Given $n \in \mathbb{N}$, we shall identify $\mathbb{R}^{2(2^n-1)d} \equiv (\mathbb{R}^d)^{\mathbb{T}_*^n}$ and $\mathbb{R}^{2^d} \equiv (\mathbb{R}^d)^{\mathbb{L}^n}$. Specifically, vectors $\mathbf{x}_n, \mathbf{y}_n \in \mathbb{R}^{2(2^n-1)d}$ and $\mathbf{z}_n \in \mathbb{R}^{2^d}$ will be regarded as indexed families

$$\mathbf{x}_n = (x_i)_{i \in \mathbb{T}_*^n}, \quad \mathbf{y}_n = (y_i)_{i \in \mathbb{T}_*^n}, \quad \mathbf{z}_n = (z_j)_{j \in \mathbb{L}^n}, \qquad (3.1)$$

where $x_i, y_i \in \mathbb{R}^d$ and $z_j \in \mathbb{R}^d$ for each $i \in \mathbb{T}_*^n$ and $j \in \mathbb{L}^n$.

**3.2. Gaussian solutions.** In this part, we compute particular solutions of the time-discrete problem (1.1) and the associated eigenproblem (1.11). As we advanced before, we shall exploit the explicit algebraic structure imposed by $G$ (given by the Gaussian mixing kernel (1.4)) and $m$ (given by the quadratic selection function (1.5)). Namely, explicit Gaussian solution will be obtained. We recall first the following stability property of Gaussians under convolutions.

**Lemma 3.2** (Stability of Gaussians). The following relation holds true

$$G_{\mu_1, \sigma_1^2} * G_{\mu_2, \sigma_2^2} = G_{\mu_1 + \mu_2, \sigma_1^2 + \sigma_2^2}, \qquad (3.2)$$

for any couple of means $\mu_1, \mu_2 \in \mathbb{R}^d$ and variances $\sigma_1^2, \sigma_2^2 > 0$.

Using the above result, we can obtain the following explicit evaluation of the operator $\mathcal{T}$ in (1.2) over the class of Gaussian functions.

**Lemma 3.3** (Evaluation on Gaussians). Consider any $\mu \in \mathbb{R}^d$ and $\sigma^2 \in \mathbb{R}_+^*$, then

$$\mathcal{T}[G_{\mu, \sigma^2}] = m_* G_{\mu, \sigma^2},$$

where the parameters $m_*$, $\mu_*$ and $\sigma_*$ are given by:

$$m_* = \frac{e^{-\frac{1}{2} \frac{\alpha |\mu|^2}{1 + \alpha(1 + \frac{\sigma^2}{2})}}}{(1 + \alpha(1 + \frac{\sigma^2}{2}))^{d/2}}, \quad \mu_* = \frac{\mu}{1 + \alpha(1 + \frac{\sigma^2}{2})}, \quad \sigma_*^2 = \frac{1 + \frac{\sigma^2}{2}}{1 + \alpha(1 + \frac{\sigma^2}{2})}. \qquad (3.3)$$

*Proof.* On the one hand, note that by Lemma 3.2

$$\mathcal{B}[G_{\mu, \sigma^2}](x) = \left( G\left(\frac{\cdot}{2}\right) * G_{\mu, \sigma^2} * G_{\mu, \sigma^2} \right) (2x) = G_{\mu, 1 + \frac{\sigma^2}{2}}(x),$$

for each $x \in \mathbb{R}^d$. Therefore, by definition of $\mathcal{T}$ in (1.2) we obtain

$$\mathcal{T}[G_{\mu, \sigma^2}](x) = e^{-\frac{\alpha}{2}|x|^2} G_{\mu, 1 + \frac{\sigma^2}{2}}(x) = \frac{1}{(2\pi(1 + \frac{\sigma^2}{2}))^{d/2}} \exp \left( -\frac{\alpha}{2}|x|^2 - \frac{1}{2} \frac{1}{1 + \frac{\sigma^2}{2}} |x - \mu|^2 \right),$$

for each $x \in \mathbb{R}^d$. By completing the square inside the exponential, we conclude our result.

Consequently, the following explicit Gaussian solution of the eigenproblem (1.11) is found.

**Proposition 3.4** (Gaussian solution of the eigenproblem). Assume that $\alpha \in \mathbb{R}_+^*$, then the eigenproblem (1.11) has a unique Gaussian solution $(\lambda_\alpha, F_\alpha)$, determined by the relation (1.12) in Theorem 1.1.
<figure><img src="image_5.png" /><figcaption>FIGURE 4. Plot of eigenvalue $\lambda_{\alpha}$ and variance at equilibrium $\sigma_{\alpha}^2$ against $\alpha$ for $d = 1$. The case $\alpha = 0$ corresponds to absence of selection, which is conservative ($\lambda_{\alpha=0} = 1$), and the variance at equilibrium $\sigma_{\alpha=0}^2$ is twice the variance of the mixing kernel in $\mathcal{B}$.</figcaption></figure>

*Proof.* We look for $\lambda \in \mathbb{R}_+^*$, $\mu \in \mathbb{R}^d$ and $\sigma^2 \in \mathbb{R}_+^*$ such that $\lambda G_{\mu, \sigma^2} = \mathcal{T}[G_{\mu, \sigma^2}]$. By Lemma 3.3 and bearing in mind parameters $m_*$, $\mu_*$ and $\sigma_*^2$ in (3.3) we obtain that $(\lambda, \mu, \sigma^2)$ must solve the equations:

$$ \lambda = \frac{e^{-\frac{1}{2} \frac{\alpha |\mu|^2}{1+\alpha(1+\frac{\sigma^2}{2})}}}{(1+\alpha(1+\frac{\sigma^2}{2}))^{d/2}}, \quad \mu = \frac{\mu}{1+\alpha(1+\frac{\sigma^2}{2})}, \quad \sigma^2 = \frac{1+\frac{\sigma^2}{2}}{1+\alpha(1+\frac{\sigma^2}{2})}. $$

Hence, the only solution is given by $\mu = 0$ and $\lambda$ and $\sigma^2$ are determined by (1.12) and (1.13). $\square$

**Remark 3.5** (Dependency on selection). The eigenvalue $\lambda_{\alpha} \in \mathbb{R}_+^*$ and the variance $\sigma_{\alpha}^2 \in \mathbb{R}_+^*$ determined by the relations (1.12) and (1.13) for the unique Gaussian solution $(\lambda_{\alpha}, F_{\alpha})$ to the eigenproblem (1.11) are monotonically decreasing with the selection coefficient $\alpha$. Namely, the larger $\alpha$, the smaller $\sigma_{\alpha}^2$ and $\lambda_{\alpha}$. Indeed, we obtain (see also Figure 4)

$$ \alpha \nearrow \infty \implies \sigma_{\alpha}^2 \searrow 0, \quad \lambda_{\alpha} \searrow 0 \\ \alpha \searrow 0 \implies \sigma_{\alpha}^2 \nearrow 2, \quad \lambda_{\alpha} \nearrow 1. $$

Therefore, if selection is strong, then $F_{\alpha}$ is very concentrated around the origin and, if selection ceases, then $F_{\alpha}$ is twice as spread as the mixing kernel $G$ in (1.4), a famous result in quantitative genetics, see e.g. [9]. Note that $\lambda_{\alpha} < 1$ for any $\alpha \in \mathbb{R}_+^*$ and, consequently, the special steady solutions $F_n = \lambda_{\alpha}^n F_{\alpha}$ coming from ansatz (1.10) always get extinct for large $n$. This is an artificial consequence of our parameter reduction in Appendix A. In particular, if we maintain parameter $\beta$ in (A.4) coming from the ratio between birth and mortality rates, then the above special steady solution only extincts if $\beta < (1 + \alpha(1 + \frac{\sigma_{\alpha}^2}{2}))$.

**Remark 3.6** (Eigenproblem with $\alpha = 0$). The same ideas as in Proposition 3.4 yield the Gaussian solutions to the eigenproblem (1.11) in the absence of selection (i.e., $\alpha = 0$). However, there is no longer uniqueness due to the translation invariance of $\mathcal{B}$. Indeed, we obtain the Gaussian solutions:

$$ \lambda = \lambda_{\alpha=0} = 1, \quad F = F_{\alpha=0}(\cdot - \mu) = G_{\mu,2}, $$

for any $\mu \in \mathbb{R}^d$. Indeed, it is a consequence of the contraction property stated in Section 2.1 that these are all possible generic solutions (not only Gaussian). Specifically, by the conservation of mass (2.1) we obtain that the only possible eigenvalue is $\lambda = 1$. Then, the eigenproblem (1.11) reduces to a fixed point equation for $\mathcal{B}$. We can then conclude by Corollary 2.2.

We emphasize that the above presents a crucial difference of behavior between the nonlinear problem (1.11) with sexual reproduction and the analogous linear version with asexual reproduction operator (1.7). Namely, whilst in the former case there exists nontrivial solutions in the absence of selection, in the latter case it can be seen that such a nontrivial solution does not exists (e.g., through Fourier arguments). This suggests that whilst in the linear problem both reproduction and selection stay balanced, in the nonlinear problem reproduction appears to dominates selection structurally.
<figure><img src="image_5.png" /><figcaption>FIGURE 4. Plot of eigenvalue $\lambda_{\alpha}$ and variance at equilibrium $\sigma_{\alpha}^2$ against $\alpha$ for $d = 1$. The case $\alpha = 0$ corresponds to absence of selection, which is conservative ($\lambda_{\alpha=0} = 1$), and the variance at equilibrium $\sigma_{\alpha=0}^2$ is twice the variance of the mixing kernel in $\mathcal{B}$.</figcaption></figure>

*Proof.* We look for $\lambda \in \mathbb{R}_+^*$, $\mu \in \mathbb{R}^d$ and $\sigma^2 \in \mathbb{R}_+^*$ such that $\lambda G_{\mu, \sigma^2} = \mathcal{T}[G_{\mu, \sigma^2}]$. By Lemma 3.3 and bearing in mind parameters $m_*$, $\mu_*$ and $\sigma_*^2$ in (3.3) we obtain that $(\lambda, \mu, \sigma^2)$ must solve the equations:

$$ \lambda = \frac{e^{-\frac{1}{2} \frac{\alpha |\mu|^2}{1+\alpha(1+\frac{\sigma^2}{2})}}}{(1+\alpha(1+\frac{\sigma^2}{2}))^{d/2}}, \quad \mu = \frac{\mu}{1+\alpha(1+\frac{\sigma^2}{2})}, \quad \sigma^2 = \frac{1+\frac{\sigma^2}{2}}{1+\alpha(1+\frac{\sigma^2}{2})}. $$

Hence, the only solution is given by $\mu = 0$ and $\lambda$ and $\sigma^2$ are determined by (1.12) and (1.13). $\square$

**Remark 3.5** (Dependency on selection). The eigenvalue $\lambda_{\alpha} \in \mathbb{R}_+^*$ and the variance $\sigma_{\alpha}^2 \in \mathbb{R}_+^*$ determined by the relations (1.12) and (1.13) for the unique Gaussian solution $(\lambda_{\alpha}, F_{\alpha})$ to the eigenproblem (1.11) are monotonically decreasing with the selection coefficient $\alpha$. Namely, the larger $\alpha$, the smaller $\sigma_{\alpha}^2$ and $\lambda_{\alpha}$. Indeed, we obtain (see also Figure 4)

$$ \begin{aligned} \alpha \nearrow \infty &\Rightarrow \sigma_{\alpha}^2 \searrow 0, & \lambda_{\alpha} \searrow 0 \\ \alpha \searrow 0 &\Rightarrow \sigma_{\alpha}^2 \nearrow 2, & \lambda_{\alpha} \nearrow 1. \end{aligned} $$

Therefore, if selection is strong, then $F_{\alpha}$ is very concentrated around the origin and, if selection ceases, then $F_{\alpha}$ is twice as spread as the mixing kernel $G$ in (1.4), a famous result in quantitative genetics, see e.g. [9]. Note that $\lambda_{\alpha} < 1$ for any $\alpha \in \mathbb{R}_+^*$ and, consequently, the special steady solutions $F_n = \lambda_{\alpha}^n F_{\alpha}$ coming from ansatz (1.10) always get extinct for large $n$. This is an artificial consequence of our parameter reduction in Appendix A. In particular, if we maintain parameter $\beta$ in (A.4) coming from the ratio between birth and mortality rates, then the above special steady solution only extincts if $\beta < (1 + \alpha(1 + \frac{\sigma_{\alpha}^2}{2}))$.

**Remark 3.6** (Eigenproblem with $\alpha = 0$). The same ideas as in Proposition 3.4 yield the Gaussian solutions to the eigenproblem (1.11) in the absence of selection (i.e., $\alpha = 0$). However, there is no longer uniqueness due to the translation invariance of $\mathcal{B}$. Indeed, we obtain the Gaussian solutions:

$$ \lambda = \lambda_{\alpha=0} = 1, \quad F = F_{\alpha=0}(\cdot - \mu) = G_{\mu,2}, $$

for any $\mu \in \mathbb{R}^d$. Indeed, it is a consequence of the contraction property stated in Section 2.1 that these are all possible generic solutions (not only Gaussian). Specifically, by the conservation of mass (2.1) we obtain that the only possible eigenvalue is $\lambda = 1$. Then, the eigenproblem (1.11) reduces to a fixed point equation for $\mathcal{B}$. We can then conclude by Corollary 2.2.

We emphasize that the above presents a crucial difference of behavior between the nonlinear problem (1.11) with sexual reproduction and the analogous linear version with asexual reproduction operator (1.7). Namely, whilst in the former case there exists nontrivial solutions in the absence of selection, in the latter case it can be seen that such a nontrivial solution does not exists (e.g., through Fourier arguments). This suggests that whilst in the linear problem both reproduction and selection stay balanced, in the nonlinear problem reproduction appears to dominates selection structurally.
Finally, we note that by iteration of the operator $\mathcal{T}$ and Lemma 3.3 we can compute the explicit form of solutions of the time-discrete problem (1.1) issued at Gaussian initial data. In fact, the Gaussian structure is preserved along generations, although mass, mean and variance are modified.

**Proposition 3.7** (Gaussian solutions of the time-discrete problem). Consider any $m_0 \in \mathbb{R}_+^*$, $\mu_0 \in \mathbb{R}^d$ and $\sigma_0^2 \in \mathbb{R}_+$, and set the Gaussian initial datum $F_0 = m_0 G_{\mu_0, \sigma_0^2}$. Hence, the solution $\{F_n\}_{n \in \mathbb{N}}$ to the time-discrete problem (1.1) takes the form

$$F_n = m_n G_{\mu_n, \sigma_n^2}, \quad (3.4)$$

for $n \in \mathbb{N}$, where the parameters $m_n \in \mathbb{R}_+^*$, $\mu_n \in \mathbb{R}^d$ and $\sigma_n^2 \in \mathbb{R}_+^*$ are governed by the recursions:

$$m_n = m_{n-1} \frac{e^{-\frac{1}{2} \frac{\alpha |\mu_{n-1}|^2}{1+\alpha \left(1+\frac{\sigma_{n-1}^2}{2}\right)}}}{\left(1+\alpha\left(1+\frac{\sigma_{n-1}^2}{2}\right)\right)^{d/2}}, \quad \mu_n = \frac{\mu_{n-1}}{1+\alpha\left(1+\frac{\sigma_{n-1}^2}{2}\right)}, \quad \frac{1}{\sigma_n^2} = \alpha + \frac{1}{1+\frac{\sigma_{n-1}^2}{2}}. \quad (3.5)$$

**Lemma 3.8** (Relaxation of variance). Assume that $\alpha \in \mathbb{R}_+^*$, consider any $\sigma_0^2 \in \mathbb{R}_+^*$ and define the sequence $\{\sigma_n^2\}_{n \in \mathbb{N}}$ by recursion according to the third recursion in (3.5), i.e.

$$\frac{1}{\sigma_n^2} = \alpha + \frac{1}{1 + \frac{\sigma_{n-1}^2}{2}}, \quad (3.6)$$

for every $n \in \mathbb{N}$. Hence, if $\sigma_0^2 > \sigma_\alpha^2$ then $\sigma_n \searrow \sigma_\alpha^2$ as $n \to \infty$ and, if $\sigma_0^2 < \sigma_\alpha^2$ then $\sigma_n \nearrow \sigma_\alpha^2$ as $n \to \infty$, where $\sigma_\alpha^2$ is given by (1.13). In addition, we obtain the convergence rates

$$|\sigma_n^2 - \sigma_\alpha^2| \le C_v r_\alpha^n, \quad (3.7)$$

for every $n \in \mathbb{N}$, where the constant $C_v \in \mathbb{R}_+$ depends on $\sigma_0^2$ and $\alpha$ (we obtain that $C_v = 0$ if $\sigma_0^2 = \sigma_\alpha^2$), and the ratio $r_\alpha \in (0, \frac{1}{2})$ satisfies $r_\alpha = 2k_\alpha^2$ and is given explicitly by

$$r_\alpha := \frac{8}{\left( (2\alpha + 3) + \sqrt{(2\alpha + 1)^2 + 8\alpha} \right)^2}. \quad (3.8)$$

*Proof.* • STEP 1: Monotone convergence of $\{\sigma_n^2\}_{n \in \mathbb{N}}$.

Consider the function $f: \mathbb{R}_+ \to \mathbb{R}$ given by

$$f(x) = \alpha + \frac{1}{1 + \frac{1}{2x}}, \quad x \in \mathbb{R}_+,$$

and, for simplicity of notation, define

$$x_* := (\sigma_\alpha^2)^{-1}, \quad x_n := (\sigma_n^2)^{-1}, \quad n \in \mathbb{N}. \quad (3.9)$$

Notice that $x_*$ is the unique fixed point of $f$ in $\mathbb{R}_+$ and $\{x_n\}_{n \in \mathbb{N}}$ determines the fixed-point iteration of the map $f$ issued at $x_0 = (\sigma_0^2)^{-1}$, that is, $x_n = f(x_{n-1})$, $n \in \mathbb{N}$. Our goal is to show that $\{x_n\}_{n \in \mathbb{N}}$ converges towards $x_*$ and find convergence rates. By direct computation we obtain $f'(x) = 2/(1+2x)^2$, which is above 1 near $x=0$, and therefore $f$ is not contractive. Hence, we cannot apply the usual Banach contraction principle and a different argument is provided. First, note that $x \in \mathbb{R}_+ \mapsto f(x)$ is strictly increasing and $x \in \mathbb{R}_+^* \mapsto \frac{f(x)}{x}$ is strictly decreasing. Then, we obtain that

$$x < f(x) < x_*, \quad \text{if } 0 \le x < x_*,$$

$$f(x) = x_*, \quad \text{if } x = x_*,$$

$$x_* < f(x) < x, \quad \text{if } x > x_*,$$

which implies the aforementioned monotonicity properties of $\{x_n\}_{n \in \mathbb{N}}$ and, in addition,

$$\min\{x_*, x_0\} \le x_n \le \max\{x_*, x_0\}, \quad (3.10)$$

for any $n \in \mathbb{N}$. In particular, this yields the monotonic convergence of $\{x_n\}_{n \in \mathbb{N}}$ towards $x_*$, or equivalently, the claimed monotonic convergence of $\{\sigma_n^2\}_{n \in \mathbb{N}}$ towards $\sigma_\alpha^2$ as $n \to \infty$.

• STEP 2: Convergence rates.

Second, to compute the convergence rates we shall use the special algebraic structure of function $f$, which is the restriction to $\mathbb{R}_+$ of the Möbius function $M: \mathbb{R} \setminus \{-\frac{1}{2}\} \to \mathbb{R}$ given by

$$M(x) := \frac{2(\alpha + 1)x + \alpha}{2x + 1}, \quad x \in \mathbb{R} \setminus \left\{-\frac{1}{2}\right\}.$$
Finally, we note that by iteration of the operator $\mathcal{T}$ and Lemma 3.3 we can compute the explicit form of solutions of the time-discrete problem (1.1) issued at Gaussian initial data. In fact, the Gaussian structure is preserved along generations, although mass, mean and variance are modified.

**Proposition 3.7** (Gaussian solutions of the time-discrete problem). Consider any $m_0 \in \mathbb{R}_+^*$, $\mu_0 \in \mathbb{R}^d$ and $\sigma_0^2 \in \mathbb{R}_+$, and set the Gaussian initial datum $F_0 = m_0 G_{\mu_0, \sigma_0^2}$. Hence, the solution $\{F_n\}_{n \in \mathbb{N}}$ to the time-discrete problem (1.1) takes the form

$$F_n = m_n G_{\mu_n, \sigma_n^2}, \qquad (3.4)$$

for $n \in \mathbb{N}$, where the parameters $m_n \in \mathbb{R}_+^*$, $\mu_n \in \mathbb{R}^d$ and $\sigma_n^2 \in \mathbb{R}_+^*$ are governed by the recursions:

$$m_n = m_{n-1} \frac{e^{-\frac{1}{2} \frac{\alpha |\mu_{n-1}|^2}{1+\alpha \left(1+\frac{\sigma_{n-1}^2}{2}\right)}}}{\left(1+\alpha \left(1+\frac{\sigma_{n-1}^2}{2}\right)\right)^{d/2}}, \quad \mu_n = \frac{\mu_{n-1}}{1+\alpha \left(1+\frac{\sigma_{n-1}^2}{2}\right)}, \quad \frac{1}{\sigma_n^2} = \alpha + \frac{1}{1+\frac{\sigma_{n-1}^2}{2}}. \qquad (3.5)$$

**Lemma 3.8** (Relaxation of variance). Assume that $\alpha \in \mathbb{R}_+^*$, consider any $\sigma_0^2 \in \mathbb{R}_+^*$ and define the sequence $\{\sigma_n^2\}_{n \in \mathbb{N}}$ by recursion according to the third recursion in (3.5), i.e.

$$\frac{1}{\sigma_n^2} = \alpha + \frac{1}{1 + \frac{\sigma_{n-1}^2}{2}}, \qquad (3.6)$$

for every $n \in \mathbb{N}$. Hence, if $\sigma_0^2 > \sigma_\alpha^2$ then $\sigma_n \searrow \sigma_\alpha^2$ as $n \to \infty$ and, if $\sigma_0^2 < \sigma_\alpha^2$ then $\sigma_n \nearrow \sigma_\alpha^2$ as $n \to \infty$, where $\sigma_\alpha^2$ is given by (1.13). In addition, we obtain the convergence rates

$$|\sigma_n^2 - \sigma_\alpha^2| \le C_v r_\alpha^n, \qquad (3.7)$$

for every $n \in \mathbb{N}$, where the constant $C_v \in \mathbb{R}_+$ depends on $\sigma_0^2$ and $\alpha$ (we obtain that $C_v = 0$ if $\sigma_0^2 = \sigma_\alpha^2$), and the ratio $r_\alpha \in (0, \frac{1}{2})$ satisfies $r_\alpha = 2k_\alpha^2$ and is given explicitly by

$$r_\alpha := \frac{8}{\left( (2\alpha + 3) + \sqrt{(2\alpha + 1)^2 + 8\alpha} \right)^2}. \qquad (3.8)$$

*Proof.* • STEP 1: Monotone convergence of $\{\sigma_n^2\}_{n \in \mathbb{N}}$.

Consider the function $f: \mathbb{R}_+ \to \mathbb{R}$ given by

$$f(x) = \alpha + \frac{1}{1 + \frac{1}{2x}}, \quad x \in \mathbb{R}_+,$$

and, for simplicity of notation, define

$$x_* := (\sigma_\alpha^2)^{-1}, \quad x_n := (\sigma_n^2)^{-1}, \quad n \in \mathbb{N}. \qquad (3.9)$$

Notice that $x_*$ is the unique fixed point of $f$ in $\mathbb{R}_+$ and $\{x_n\}_{n \in \mathbb{N}}$ determines the fixed-point iteration of the map $f$ issued at $x_0 = (\sigma_0^2)^{-1}$, that is, $x_n = f(x_{n-1})$, $n \in \mathbb{N}$. Our goal is to show that $\{x_n\}_{n \in \mathbb{N}}$ converges towards $x_*$ and find convergence rates. By direct computation we obtain $f'(x) = 2/(1+2x)^2$, which is above 1 near $x=0$, and therefore $f$ is not contractive. Hence, we cannot apply the usual Banach contraction principle and a different argument is provided. First, note that $x \in \mathbb{R}_+ \mapsto f(x)$ is strictly increasing and $x \in \mathbb{R}_+^* \mapsto \frac{f(x)}{x}$ is strictly decreasing. Then, we obtain that

$$x < f(x) < x_*, \quad \text{if } 0 \le x < x_*,$$

$$f(x) = x_*, \quad \text{if } x = x_*,$$

$$x_* < f(x) < x, \quad \text{if } x > x_*,$$

which implies the aforementioned monotonicity properties of $\{x_n\}_{n \in \mathbb{N}}$ and, in addition,

$$\min\{x_*, x_0\} \le x_n \le \max\{x_*, x_0\}, \qquad (3.10)$$

for any $n \in \mathbb{N}$. In particular, this yields the monotonic convergence of $\{x_n\}_{n \in \mathbb{N}}$ towards $x_*$, or equivalently, the claimed monotonic convergence of $\{\sigma_n^2\}_{n \in \mathbb{N}}$ towards $\sigma_\alpha^2$ as $n \to \infty$.

• STEP 2: Convergence rates.

Second, to compute the convergence rates we shall use the special algebraic structure of function $f$, which is the restriction to $\mathbb{R}_+$ of the Möbius function $M: \mathbb{R} \setminus \{-\frac{1}{2}\} \to \mathbb{R}$ given by

$$M(x) := \frac{2(\alpha + 1)x + \alpha}{2x + 1}, \quad x \in \mathbb{R} \setminus \left\{-\frac{1}{2}\right\}.$$
<figure><img src="image_6.png" /><figcaption>FIGURE 5. The ratio $r_\alpha$ against parameter $\alpha$.</figcaption></figure>

Note that $M$ has two fixed points:

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

where we note that $x_+ = x_* \in \mathbb{R}_+^*$ and $x_- \in \mathbb{R}^-$. Then, by definition of $\{x_n\}_{n \in \mathbb{N}}$ in (3.9) we obtain

$$x_n - x_+ = \frac{2(\alpha + 1)x_{n-1} + \alpha}{2x_{n-1} + 1} - \frac{2(\alpha + 1)x_+ + \alpha}{2x_+ + 1} = \frac{2}{(2x_{n-1} + 1)(2x_+ + 1)}(x_{n-1} - x_+) \quad (3.11)$$

Since $x_0$ (thus $x_{n-1}$ for initial time steps $n$) can be chosen arbitrary close to 0, then, the best a priori control that we can have on the prefactor in the right hand side of (3.11) is

$$0 \le \frac{2}{(2x_{n-1} + 1)(2x_+ + 1)} \le \frac{2\sigma_{\alpha}^2}{2 + \sigma_{\alpha}^2} = 2k_{\alpha},$$

which gives true contraction in (3.11) when $\alpha > 0$, since $2k_\alpha < 1$. Nevertheless, such a rate is non-optimal, and we show an alternative approach to achieve a sharper one. Specifically, note that

$$x_n - x_- = \frac{2(\alpha + 1)x_{n-1} + \alpha}{2x_{n-1} + 1} - \frac{2(\alpha + 1)x_- + \alpha}{2x_- + 1} = \frac{2}{(2x_{n-1} + 1)(2x_- + 1)}(x_{n-1} - x_-) \quad (3.12)$$

By dividing (3.11) by (3.12) and iterating the identity we obtain

$$\frac{x_n - x_+}{x_n - x_-} = \left( \frac{2x_- + 1}{2x_+ + 1} \right)^n \frac{x_0 - x_+}{x_0 - x_-}, \quad (3.13)$$

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

$$\frac{2x_- + 1}{2x_+ + 1} = \frac{(2\alpha + 3) - \sqrt{(2\alpha + 1)^2 + 8\alpha}}{(2\alpha + 3) + \sqrt{(2\alpha + 1)^2 + 8\alpha}} = \frac{8}{((2\alpha + 3) + \sqrt{(2\alpha + 1)^2 + 8\alpha})^2} = r_{\alpha} < \frac{1}{2}.$$

Therefore, using the uniform control (3.10) and (3.13), we conclude that

$$|x_n - x_*| = \frac{|x_n - x_-|}{|x_0 - x_-|} |x_0 - x_*| \mathbf{r}_\alpha^n \le \frac{\max\{x_*, x_0\} - x_-}{\min\{x_*, x_0\} - x_-} |x_0 - x_*| \mathbf{r}_\alpha^n, \quad (3.14)$$

for every $n \in \mathbb{N}$, thus leading to an improved rate since $\mathbf{r}_\alpha = 2\mathbf{k}_\alpha^2 < 2\mathbf{k}_\alpha$. Finally, notice that

$$|\sigma_n^2 - \sigma_\alpha^2| = \sigma_n^2 \sigma_\alpha^2 |x_n - x_*| \le \max\{\sigma_\alpha^2, \sigma_0^2\}^2 |x_n - x_*|, \quad (3.15)$$

for every $n \in \mathbb{N}$. Hence, joining (3.14), (3.15) along with (3.10) ends the proof. $\square$

In particular, for Gaussian initial data we recover an explicit particular case of the general relaxation result in Theorem 1.1.

**Corollary 3.9 (Relaxation of Gaussian solutions).** Assume that $\alpha \in \mathbb{R}_+^*$, consider any $m_0 \in \mathbb{R}_+$, $\mu_0 \in \mathbb{R}^d$ and $\sigma_0^2 \in \mathbb{R}_+$, and set the Gaussian initial datum $F_0 = m_0 G_{\mu_0, \sigma_0}$. Hence, the Gaussian solution $\{F_n\}_{n \in \mathbb{N}}$ to the time-discrete problem (1.1) verifies that the growth rates $\|F_n\|_{L^1(\mathbb{R}^d)}/\|F_{n-1}\|_{L^1(\mathbb{R}^d)}$ relax towards
<figure><img src="image_6.png" /><figcaption>FIGURE 5. The ratio $r_\alpha$ against parameter $\alpha$.</figcaption></figure>

Note that $M$ has two fixed points:

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

where we note that $x_+ = x_* \in \mathbb{R}_+^*$ and $x_- \in \mathbb{R}^-$. Then, by definition of $\{x_n\}_{n \in \mathbb{N}}$ in (3.9) we obtain

$$x_n - x_+ = \frac{2(\alpha + 1)x_{n-1} + \alpha}{2x_{n-1} + 1} - \frac{2(\alpha + 1)x_+ + \alpha}{2x_+ + 1} = \frac{2}{(2x_{n-1} + 1)(2x_+ + 1)}(x_{n-1} - x_+) \quad (3.11)$$

Since $x_0$ (thus $x_{n-1}$ for initial time steps $n$) can be chosen arbitrary close to 0, then, the best a priori control that we can have on the prefactor in the right hand side of (3.11) is

$$0 \le \frac{2}{(2x_{n-1} + 1)(2x_+ + 1)} \le \frac{2\sigma_{\alpha}^2}{2 + \sigma_{\alpha}^2} = 2k_{\alpha},$$

which gives true contraction in (3.11) when $\alpha > 0$, since $2k_\alpha < 1$. Nevertheless, such a rate is non-optimal, and we show an alternative approach to achieve a sharper one. Specifically, note that

$$x_n - x_- = \frac{2(\alpha + 1)x_{n-1} + \alpha}{2x_{n-1} + 1} - \frac{2(\alpha + 1)x_- + \alpha}{2x_- + 1} = \frac{2}{(2x_{n-1} + 1)(2x_- + 1)}(x_{n-1} - x_-) \quad (3.12)$$

By dividing (3.11) by (3.12) and iterating the identity we obtain

$$\frac{x_n - x_+}{x_n - x_-} = \left( \frac{2x_- + 1}{2x_+ + 1} \right)^n \frac{x_0 - x_+}{x_0 - x_-}, \quad (3.13)$$

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

$$\frac{2x_- + 1}{2x_+ + 1} = \frac{(2\alpha + 3) - \sqrt{(2\alpha + 1)^2 + 8\alpha}}{(2\alpha + 3) + \sqrt{(2\alpha + 1)^2 + 8\alpha}} = \frac{8}{((2\alpha + 3) + \sqrt{(2\alpha + 1)^2 + 8\alpha})^2} = r_\alpha < \frac{1}{2}.$$

Therefore, using the uniform control (3.10) and (3.13), we conclude that

$$|x_n - x_*| = \frac{|x_n - x_-|}{|x_0 - x_-|} |x_0 - x_*| \mathbf{r}_\alpha^n \le \frac{\max\{x_*, x_0\} - x_-}{\min\{x_*, x_0\} - x_-} |x_0 - x_*| \mathbf{r}_\alpha^n, \quad (3.14)$$

for every $n \in \mathbb{N}$, thus leading to an improved rate since $\mathbf{r}_\alpha = 2\mathbf{k}_\alpha^2 < 2\mathbf{k}_\alpha$. Finally, notice that

$$|\sigma_n^2 - \sigma_\alpha^2| = \sigma_n^2 \sigma_\alpha^2 |x_n - x_*| \le \max\{\sigma_\alpha^2, \sigma_0^2\}^2 |x_n - x_*|, \quad (3.15)$$

for every $n \in \mathbb{N}$. Hence, joining (3.14), (3.15) along with (3.10) ends the proof. $\square$

In particular, for Gaussian initial data we recover an explicit particular case of the general relaxation result in Theorem 1.1.

**Corollary 3.9 (Relaxation of Gaussian solutions).** Assume that $\alpha \in \mathbb{R}_+^*$, consider any $m_0 \in \mathbb{R}_+$, $\mu_0 \in \mathbb{R}^d$ and $\sigma_0^2 \in \mathbb{R}_+$, and set the Gaussian initial datum $F_0 = m_0 G_{\mu_0, \sigma_0}$. Hence, the Gaussian solution $\{F_n\}_{n \in \mathbb{N}}$ to the time-discrete problem (1.1) verifies that the growth rates $\|F_n\|_{L^1(\mathbb{R}^d)}/\|F_{n-1}\|_{L^1(\mathbb{R}^d)}$ relax towards
$\lambda_\alpha$ and the normalized profiles $F_n/\|F_n\|_{L^1(\mathbb{R}^d)}$ relax towards $F_\alpha$, where $(\lambda_\alpha, F_\alpha)$ is the unique Gaussian solution (1.12) to the eigenproblem (1.11) as proved in Proposition 3.4. Specifically, for any $\varepsilon \in \mathbb{R}_+^*$

$$ \mathcal{D}_{\text{KL}} \left( \frac{F_n}{\|F_n\|_{L^1(\mathbb{R}^d)}} \| \mathbf{F}_\alpha \right) \le C_{\mu,\varepsilon}(\lambda_\alpha^{4/d} + \varepsilon)^n + C_v \mathbf{r}_\alpha^n, \\ \left| \frac{\|F_n\|_{L^1(\mathbb{R}^d)}}{\|F_{n-1}\|_{L^1(\mathbb{R}^d)}} - \lambda_\alpha \right| \le C_{\mu,\varepsilon}(\lambda_\alpha^{4/d} + \varepsilon)^n + C_v \mathbf{r}_\alpha^n, $$

where $\mathcal{D}_{\text{KL}}$ is the Kullback-Leibler divergence (or relative entropy), $\mathbf{r}_\alpha$ is given by (3.8) in Lemma 3.8, $C_{\mu,\varepsilon} \in \mathbb{R}_+$ depends on $\alpha, \mu_0, \sigma_0^2$ and $\varepsilon$, and $C_v \in \mathbb{R}_+$ depends only on $\alpha$ and $\sigma_0^2$. In fact, $C_{\mu,\varepsilon} = 0$ if $\mu_0 = 0$ and $C_v = 0$ if $\sigma_0^2 = \sigma_\alpha^2$.

Again, in this case the computations of the Kullback-Leibler divergence becomes explicit, and it is based on the general formula below for the divergence between two Gaussian density functions:

$$ \mathcal{D}_{\text{KL}}(G_{\mu_1, \Sigma_1} \| G_{\mu_2, \Sigma_2}) = \frac{1}{2}(\mu_2 - \mu_1)^{\top} \Sigma_2^{-1} (\mu_2 - \mu_1) + \frac{1}{2} \text{trace}(\Sigma_1 \Sigma_2^{-1}) - \frac{d}{2} + \frac{1}{2} \log \det(\Sigma_2 \Sigma_1^{-1}), \quad (3.16) $$

for every $\mu_1, \mu_2 \in \mathbb{R}^d$ and any positive non-singular matrices $\Sigma_1, \Sigma_2 \in \mathbb{R}^{d \times d}$.

*Proof of Corollary 3.9.* First, note that when $\alpha \in \mathbb{R}_+^*$, then the parameters $m_n, \mu_n$ and $\sigma_n^2$ in (3.5) verify

$$ \lim_{n \to \infty} \frac{m_n}{m_{n-1}} = \lambda_\alpha, \quad \lim_{n \to \infty} \mu_n = 0, \quad \lim_{n \to \infty} \sigma_n^2 = \sigma_\alpha^2. $$

In the sequel, we quantify the rates of convergence. Since $\sigma_n^2$ has already been studied in Lemma 3.8 we shall focus only on $\mu_n$ and $m_n$. Let us fix any arbitrarily small $\varepsilon \in \mathbb{R}_+^*$. On the one hand note that

$$ \frac{|\mu_n|}{|\mu_{n-1}|} - \lambda_\alpha^{2/d} = \frac{1}{1 + \alpha \left(1 + \frac{\sigma_{n-1}^2}{2}\right)} - \frac{1}{1 + \alpha \left(1 + \frac{\sigma_\alpha^2}{2}\right)}, $$

for every $n \in \mathbb{N}$. By the mean value theorem and Lemma 3.8 we obtain

$$ \left| \frac{|\mu_n|}{|\mu_{n-1}|} - \lambda_\alpha^{2/d} \right| \le C_v \mathbf{r}_\alpha^n, $$

for an appropriate $C_v \in \mathbb{R}_+^*$ depending on $\alpha$ and $\sigma_0^2$. Since $\lambda_\alpha \in (0, 1)$ then d'Alembert's ratio test readily shows that $\mu_n$ relaxes to zero as a geometric sequence. Indeed, note that

$$ |\mu_n| = \left( \frac{|\mu_n|}{|\mu_{n-1}|} - \lambda_\alpha^{2/d} \right) |\mu_{n-1}| + \lambda_\alpha^{2/d} |\mu_{n-1}| \le (\lambda_\alpha^{2/d} + C_v \mathbf{r}_\alpha^n) |\mu_{n-1}|, $$

for any $n \in \mathbb{N}$. By an inductive argument, this yields

$$ |\mu_n| \le \prod_{k=1}^{n} (\lambda_{\alpha}^{2/d} + C_v \mathbf{r}_{\alpha}^n) |\mu_0| \lesssim C_{\varepsilon}(\lambda_{\alpha}^{2/d} + \varepsilon)^n, \quad (3.17) $$

for sufficiently large $C_{\mu,\varepsilon} \in \mathbb{R}^d$, where we have used that $\mathbf{r}_\alpha \in (0, 1)$ to absorb $C_v \mathbf{r}_\alpha^n$ in an $\varepsilon$-small term. On the other hand, note that

$$ \frac{m_n}{m_{n-1}} - \lambda_\alpha = \frac{e^{-\frac{1}{2} \frac{\alpha |\mu_{n-1}|^2}{1+\alpha \left(1+\frac{\sigma_{n-1}^2}{2}\right)}}}{\left(1+\alpha\left(1+\frac{\sigma_{n-1}^2}{2}\right)\right)^{d/2}} - \frac{1}{\left(1+\alpha\left(1+\frac{\sigma_\alpha^2}{2}\right)\right)^{d/2}} \\ = \frac{e^{-\frac{1}{2} \frac{\alpha |\mu_{n-1}|^2}{1+\alpha \left(1+\frac{\sigma_{n-1}^2}{2}\right)}} - 1}{\left(1+\alpha\left(1+\frac{\sigma_{n-1}^2}{2}\right)\right)^{d/2}} + \left( \frac{1}{\left(1+\alpha\left(1+\frac{\sigma_{n-1}^2}{2}\right)\right)^{d/2}} - \frac{1}{\left(1+\alpha\left(1+\frac{\sigma_\alpha^2}{2}\right)\right)^{d/2}} \right). $$

By the mean value theorem,

$$ \left| \frac{m_n}{m_{n-1}} - \lambda_\alpha \right| \le \frac{\alpha |\mu_{n-1}|^2}{2(1+\alpha)^{d/2+1}} + \frac{d\alpha}{4} \frac{|\sigma_{n-1}^2 - \sigma_\alpha^2|}{(1+\alpha)^{d/2+1}}. $$

Therefore, using contraction of mean (3.17) and contraction of variance (3.7) in Lemma 3.8 entail the rate of convergence for the $m_n/m_{n-1}$. Finally, given that both $F_n/\|F_n\|_{L^1(\mathbb{R}^d)}$ and $F_\alpha$ are Gaussian functions, the relative entropy can be computed explicitly through (3.16). Specifically,

$$ \mathcal{D}_{\text{KL}} \left( \frac{F_n}{\|F_n\|_{L^1(\mathbb{R}^d)}} \| \mathbf{F}_\alpha \right) = \frac{|\mu_n|^2}{2 \sigma_\alpha^2} + \frac{d}{2} \left( \frac{\sigma_n^2}{\sigma_\alpha^2} - 1 \right) - \frac{d}{2} \log \left( \frac{\sigma_n^2}{\sigma_\alpha^2} \right). $$
$\lambda_\alpha$ and the normalized profiles $F_n/\|F_n\|_{L^1(\mathbb{R}^d)}$ relax towards $F_\alpha$, where $(\lambda_\alpha, F_\alpha)$ is the unique Gaussian solution (1.12) to the eigenproblem (1.11) as proved in Proposition 3.4. Specifically, for any $\varepsilon \in \mathbb{R}_+^*$

$$ \mathcal{D}_{\text{KL}} \left( \frac{F_n}{\|F_n\|_{L^1(\mathbb{R}^d)}} \| \mathbf{F}_\alpha \right) \le C_{\mu,\varepsilon}(\lambda_\alpha^{4/d} + \varepsilon)^n + C_v \mathbf{r}_\alpha^n, \\ \left| \frac{\|F_n\|_{L^1(\mathbb{R}^d)}}{\|F_{n-1}\|_{L^1(\mathbb{R}^d)}} - \lambda_\alpha \right| \le C_{\mu,\varepsilon}(\lambda_\alpha^{4/d} + \varepsilon)^n + C_v \mathbf{r}_\alpha^n, $$

where $\mathcal{D}_{\text{KL}}$ is the Kullback-Leibler divergence (or relative entropy), $\mathbf{r}_\alpha$ is given by (3.8) in Lemma 3.8, $C_{\mu,\varepsilon} \in \mathbb{R}_+$ depends on $\alpha, \mu_0, \sigma_0^2$ and $\varepsilon$, and $C_v \in \mathbb{R}_+$ depends only on $\alpha$ and $\sigma_0^2$. In fact, $C_{\mu,\varepsilon} = 0$ if $\mu_0 = 0$ and $C_v = 0$ if $\sigma_0^2 = \sigma_\alpha^2$.

Again, in this case the computations of the Kullback-Leibler divergence becomes explicit, and it is based on the general formula below for the divergence between two Gaussian density functions:

$$ \mathcal{D}_{\text{KL}}(G_{\mu_1, \Sigma_1} \| G_{\mu_2, \Sigma_2}) = \frac{1}{2}(\mu_2 - \mu_1)^{\top} \Sigma_2^{-1} (\mu_2 - \mu_1) + \frac{1}{2} \text{trace}(\Sigma_1 \Sigma_2^{-1}) - \frac{d}{2} + \frac{1}{2} \log \det(\Sigma_2 \Sigma_1^{-1}), \quad (3.16) $$

for every $\mu_1, \mu_2 \in \mathbb{R}^d$ and any positive non-singular matrices $\Sigma_1, \Sigma_2 \in \mathbb{R}^{d \times d}$.

*Proof of Corollary 3.9.* First, note that when $\alpha \in \mathbb{R}_+^*$, then the parameters $m_n, \mu_n$ and $\sigma_n^2$ in (3.5) verify

$$ \lim_{n \to \infty} \frac{m_n}{m_{n-1}} = \lambda_\alpha, \quad \lim_{n \to \infty} \mu_n = 0, \quad \lim_{n \to \infty} \sigma_n^2 = \sigma_\alpha^2. $$

In the sequel, we quantify the rates of convergence. Since $\sigma_n^2$ has already been studied in Lemma 3.8 we shall focus only on $\mu_n$ and $m_n$. Let us fix any arbitrarily small $\varepsilon \in \mathbb{R}_+^*$. On the one hand note that

$$ \frac{|\mu_n|}{|\mu_{n-1}|} - \lambda_\alpha^{2/d} = \frac{1}{1 + \alpha(1 + \frac{\sigma_{n-1}^2}{2})} - \frac{1}{1 + \alpha(1 + \frac{\sigma_\alpha^2}{2})}, $$

for every $n \in \mathbb{N}$. By the mean value theorem and Lemma 3.8 we obtain

$$ \left| \frac{|\mu_n|}{|\mu_{n-1}|} - \lambda_\alpha^{2/d} \right| \le C_v \mathbf{r}_\alpha^n, $$

for an appropriate $C_v \in \mathbb{R}_+^*$ depending on $\alpha$ and $\sigma_0^2$. Since $\lambda_\alpha \in (0, 1)$ then d'Alembert's ratio test readily shows that $\mu_n$ relaxes to zero as a geometric sequence. Indeed, note that

$$ |\mu_n| = \left( \frac{|\mu_n|}{|\mu_{n-1}|} - \lambda_\alpha^{2/d} \right) |\mu_{n-1}| + \lambda_\alpha^{2/d} |\mu_{n-1}| \le (\lambda_\alpha^{2/d} + C_v \mathbf{r}_\alpha^n) |\mu_{n-1}|, $$

for any $n \in \mathbb{N}$. By an inductive argument, this yields

$$ |\mu_n| \le \prod_{k=1}^{n} (\lambda_{\alpha}^{2/d} + C_v \mathbf{r}_{\alpha}^n) |\mu_0| \lesssim C_{\varepsilon}(\lambda_{\alpha}^{2/d} + \varepsilon)^n, \quad (3.17) $$

for sufficiently large $C_{\mu,\varepsilon} \in \mathbb{R}^d$, where we have used that $\mathbf{r}_\alpha \in (0, 1)$ to absorb $C_v \mathbf{r}_\alpha^n$ in an $\varepsilon$-small term. On the other hand, note that

$$ \frac{m_n}{m_{n-1}} - \lambda_\alpha = \frac{e^{-\frac{1}{2} \frac{\alpha |\mu_{n-1}|^2}{1+\alpha \left(1+\frac{\sigma_{n-1}^2}{2}\right)}}}{\left(1+\alpha\left(1+\frac{\sigma_{n-1}^2}{2}\right)\right)^{d/2}} - \frac{1}{\left(1+\alpha\left(1+\frac{\sigma_\alpha^2}{2}\right)\right)^{d/2}} \\ = \frac{e^{-\frac{1}{2} \frac{\alpha |\mu_{n-1}|^2}{1+\alpha \left(1+\frac{\sigma_{n-1}^2}{2}\right)}} - 1}{\left(1+\alpha\left(1+\frac{\sigma_{n-1}^2}{2}\right)\right)^{d/2}} + \left( \frac{1}{\left(1+\alpha\left(1+\frac{\sigma_{n-1}^2}{2}\right)\right)^{d/2}} - \frac{1}{\left(1+\alpha\left(1+\frac{\sigma_\alpha^2}{2}\right)\right)^{d/2}} \right). $$

By the mean value theorem,

$$ \left| \frac{m_n}{m_{n-1}} - \lambda_\alpha \right| \le \frac{\alpha |\mu_{n-1}|^2}{2(1+\alpha)^{d/2+1}} + \frac{d\alpha}{4} \frac{|\sigma_{n-1}^2 - \sigma_\alpha^2|}{(1+\alpha)^{d/2+1}}. $$

Therefore, using contraction of mean (3.17) and contraction of variance (3.7) in Lemma 3.8 entail the rate of convergence for the $m_n/m_{n-1}$. Finally, given that both $F_n/\|F_n\|_{L^1(\mathbb{R}^d)}$ and $F_\alpha$ are Gaussian functions, the relative entropy can be computed explicitly through (3.16). Specifically,

$$ \mathcal{D}_{\text{KL}} \left( \frac{F_n}{\|F_n\|_{L^1(\mathbb{R}^d)}} \| \mathbf{F}_\alpha \right) = \frac{|\mu_n|^2}{2 \sigma_\alpha^2} + \frac{d}{2} \left( \frac{\sigma_n^2}{\sigma_\alpha^2} - 1 \right) - \frac{d}{2} \log \left( \frac{\sigma_n^2}{\sigma_\alpha^2} \right). $$
<figure><img src="image_7.png" /><figcaption>FIGURE 6. Comparison of convergence rates in Theorem 1.1 and Corollary 3.9</figcaption></figure>

Hence, using again (3.17) and (3.7) concludes our result.

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

The proof of Theorem 1.1 for generic non-Gaussian initial data $F_0 \in \mathcal{M}_+(\mathbb{R}^d)$ will be the core of this paper and we postpone it to Section 6.

**Remark 3.10** (Optimality of convergence rates). As illustrated in Figure 6, the convergence rate of the normalized profiles obtained in our main Theorem 1.1 (blue line) is sharp, compared to the explicit one found in Corollary 3.9 (red line) for the special class of Gaussian solutions. In fact, one can easily check that the identity $\lambda_{\alpha}^{4/d} = (2\mathbf{k}_{\alpha})^2$ holds. However, there is a mismatch between the convergence rate of the rate of growth of mass in Theorem 1.1 (orange line) and the sharp rate for Gaussian solutions found in Corollary 3.9 (again red line). On the one hand, it stands to reason that the relative entropy has a certain quadratic structure, whilst the rate of growth of mass is related to $L^1$ norms, and then it is not quadratic. We will see that more clearly in the proof of Theorem 1.1 in Section 6, where we use explicitly the following relation

$$ \left| \frac{\|F_n\|_{L^1(\mathbb{R}^d)}}{\|F_{n-1}\|_{L^1(\mathbb{R}^d)}} - \lambda_{\alpha} \right| \lesssim \sqrt{D_{\mathrm{KL}} \left( \frac{F_n}{\|F_n\|_{L^1(\mathbb{R}^d)}} \|\mathbf{F}_{\alpha}\| \right)}, $$

based on Pinsker's inequality. However, a certain quadratic structure is still present in the rate of convergence of $m_n/m_{n-1}$ in Corollary 3.9, which may become explicit for generic solutions to (1.1) if one finds the hidden cancellations. However, for simplicity we do not address this technical detail here.

## 4. SOME PROPERTIES OF THE OPERATOR $\mathcal{T}$

First, note that for any $F \in \mathcal{M}_+(\mathbb{R}^d)$ we obtain that $\mathcal{T}[F] \in W^{k,1}(\mathbb{R}^d) \cap W^{k,\infty}(\mathbb{R}^d)$ for each $k \in \mathbb{N}$ thanks to the fact that $F$ has finite mass and the Gaussian mixing kernel $G$ in (1.4) is smooth and has bounded and integrable derivatives of any order. In particular, solutions $\{F_n\}_{n \in \mathbb{N}}$ to the time-discrete problem (1.1) become instantaneously smooth after the first generation for any generic initial datum $F_0 \in \mathcal{M}_+(\mathbb{R}^d)$. In the sequel, we will derive a quantitative control on the emergence of Gaussian tails for $\mathcal{T}[F]$. In addition, we will quantify the propagation of quadratic and exponential moments for the normalized profile $\mathcal{T}[F]/\|\mathcal{T}[F]\|_{L^1(\mathbb{R}^d)}$. Both a priori estimates will be required later in Section 6.

### 4.1. Emergence of Gaussian tails.

**Lemma 4.1** (Emergence of tails). Assume that $\alpha \in \mathbb{R}_+^*$, consider $F \in \mathcal{M}_+(\mathbb{R}^d)$ and set $\bar{\sigma}^2, \underline{\sigma}^2 \in \mathbb{R}_+^*$ by

$$ \bar{\sigma}^2 := \frac{1}{\alpha}, \quad \underline{\sigma}^2 \in \left(0, \frac{1}{1+\alpha}\right). $$

Then, the following properties are fulfilled:

(i) **(Upper control of tails)** There exists $C = C(\alpha, F) > 0$ such that

$$ \mathcal{T}[F](x) \leq C G_{0,\bar{\sigma}^2}(x), \quad (4.1) $$

for every $x \in \mathbb{R}^d$.
<figure><img src="image_7.png" /><figcaption>FIGURE 6. Comparison of convergence rates in Theorem 1.1 and Corollary 3.9</figcaption></figure>

Hence, using again (3.17) and (3.7) concludes our result.

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

The proof of Theorem 1.1 for generic non-Gaussian initial data $F_0 \in \mathcal{M}_+(\mathbb{R}^d)$ will be the core of this paper and we postpone it to Section 6.

**Remark 3.10 (Optimality of convergence rates).** As illustrated in Figure 6, the convergence rate of the normalized profiles obtained in our main Theorem 1.1 (blue line) is sharp, compared to the explicit one found in Corollary 3.9 (red line) for the special class of Gaussian solutions. In fact, one can easily check that the identity $\lambda_{\alpha}^{4/d} = (2\mathbf{k}_{\alpha})^2$ holds. However, there is a mismatch between the convergence rate of the rate of growth of mass in Theorem 1.1 (orange line) and the sharp rate for Gaussian solutions found in Corollary 3.9 (again red line). On the one hand, it stands to reason that the relative entropy has a certain quadratic structure, whilst the rate of growth of mass is related to $L^1$ norms, and then it is not quadratic. We will see that more clearly in the proof of Theorem 1.1 in Section 6, where we use explicitly the following relation

$$ \left| \frac{\|F_n\|_{L^1(\mathbb{R}^d)}}{\|F_{n-1}\|_{L^1(\mathbb{R}^d)}} - \lambda_{\alpha} \right| \lesssim \sqrt{D_{\text{KL}} \left( \frac{F_n}{\|F_n\|_{L^1(\mathbb{R}^d)}} \|\mathbf{F}_{\alpha}\| \right)}, $$

based on Pinsker's inequality. However, a certain quadratic structure is still present in the rate of convergence of $m_n/m_{n-1}$ in Corollary 3.9, which may become explicit for generic solutions to (1.1) if one finds the hidden cancellations. However, for simplicity we do not address this technical detail here.

## 4. SOME PROPERTIES OF THE OPERATOR $\mathcal{T}$

First, note that for any $F \in \mathcal{M}_+(\mathbb{R}^d)$ we obtain that $\mathcal{T}[F] \in W^{k,1}(\mathbb{R}^d) \cap W^{k,\infty}(\mathbb{R}^d)$ for each $k \in \mathbb{N}$ thanks to the fact that $F$ has finite mass and the Gaussian mixing kernel $G$ in (1.4) is smooth and has bounded and integrable derivatives of any order. In particular, solutions $\{F_n\}_{n \in \mathbb{N}}$ to the time-discrete problem (1.1) become instantaneously smooth after the first generation for any generic initial datum $F_0 \in \mathcal{M}_+(\mathbb{R}^d)$. In the sequel, we will derive a quantitative control on the emergence of Gaussian tails for $\mathcal{T}[F]$. In addition, we will quantify the propagation of quadratic and exponential moments for the normalized profile $\mathcal{T}[F]/\|\mathcal{T}[F]\|_{L^1(\mathbb{R}^d)}$. Both a priori estimates will be required later in Section 6.

### 4.1. Emergence of Gaussian tails.

**Lemma 4.1 (Emergence of tails).** Assume that $\alpha \in \mathbb{R}_+^*$, consider $F \in \mathcal{M}_+(\mathbb{R}^d)$ and set $\bar{\sigma}^2, \underline{\sigma}^2 \in \mathbb{R}_+^*$ by

$$ \bar{\sigma}^2 := \frac{1}{\alpha}, \quad \underline{\sigma}^2 \in \left(0, \frac{1}{1+\alpha}\right). $$

Then, the following properties are fulfilled:

(i) **(Upper control of tails)** There exists $C = C(\alpha, F) > 0$ such that

$$ \mathcal{T}[F](x) \leq C G_{0,\bar{\sigma}^2}(x), \quad (4.1) $$

for every $x \in \mathbb{R}^d$.
(ii) (Lower control of tails) There exists $c = c(\alpha, \underline{\sigma}^2, F) > 0$ such that

$$ \mathcal{T}[F](x) \geq c G_{0, \underline{\sigma}^2}(x), \qquad (4.2) $$

for every $x \in \mathbb{R}^d$.

(iii) (Control of log-derivative) Assume that $F$ is absolutely continuous with respect to the Lebesgue measure and it has Gaussian tail, i.e., there exists $\sigma^2 \in \mathbb{R}_+^*$ and $C' \in \mathbb{R}_+^*$ such that $F(x) \leq C'G_{0,\sigma^2}$ for every $x \in \mathbb{R}^d$. Then, there exists $C'' = C''(\alpha, \underline{\sigma}^2, \sigma^2, C', F) > 0$ such that

$$ |\nabla \log \mathcal{T}[F](x)| \leq C''(1 + |x|), \qquad (4.3) $$

for every $x \in \mathbb{R}^d$.

*Proof.* First, notice that $F \in \mathcal{M}_+(\mathbb{R}^d)$ is any generic distribution but $G \in L^\infty(\mathbb{R}^d)$. Then, $\mathcal{B}[F] \in L^\infty(\mathbb{R}^d)$ and, by definition of $\mathcal{T}$ in (1.2), we obtain

$$ \mathcal{T}[F](x) \leq \|\mathcal{B}[F]\|_{L^\infty(\mathbb{R}^d)} e^{-\frac{\alpha}{2}|x|^2} = \|\mathcal{B}[F]\|_{L^\infty(\mathbb{R}^d)} (2\pi\alpha^{-1})^{d/2} G_{0,\bar{\sigma}^2}(x), $$

for every $x \in \mathbb{R}^d$. Hence, (4.1) holds for appropriate $C > 0$. Second, note that

$$ \begin{aligned} \mathcal{T}[F](x) &= \frac{e^{-\frac{\alpha}{2}|x|^2}}{(2\pi)^{d/2} \|F\|_{\mathcal{M}_+(\mathbb{R}^d)}} \int_{\mathbb{R}^{2d}} \exp\left(-\frac{1}{2}\left|x - \frac{x_1+x_2}{2}\right|^2\right) F(dx_1)F(dx_2) \\ &= \frac{e^{-\frac{\alpha}{2}|x|^2}}{(2\pi)^{d/2} \|F\|_{\mathcal{M}_+(\mathbb{R}^d)}} \int_{\mathbb{R}^{2d}} \exp\left(-\frac{1}{2}|x|^2 - \frac{1}{8}|x_1+x_2|^2 + \frac{1}{2}x \cdot (x_1+x_2)\right) F(dx_1)F(dx_2) \\ &\geq \frac{e^{-\frac{\alpha}{2}|x|^2}}{(2\pi)^{d/2} \|F\|_{\mathcal{M}_+(\mathbb{R}^d)}} \int_{\mathbb{R}^{2d}} \exp\left(-\frac{2+\varepsilon}{4}|x|^2 - \frac{2+\varepsilon}{8\varepsilon}|x_1+x_2|^2\right) F(dx_1)F(dx_2) \\ &\geq \frac{e^{-\frac{\alpha}{2}|x|^2}}{(2\pi)^{d/2} \|F\|_{\mathcal{M}_+(\mathbb{R}^d)}} \int_{\mathbb{R}^{2d}} \exp\left(-\frac{2+\varepsilon}{4}|x|^2 - \frac{2+\varepsilon}{4\varepsilon}|x_1|^2 - \frac{2+\varepsilon}{4\varepsilon}|x_2|^2\right) F(dx_1)F(dx_2), \end{aligned} $$

for every $\varepsilon > 0$, where in the third and last lines we have used Cauchy–Schwarz’s and Young’s inequalities. Consequently, we obtain that

$$ \mathcal{T}[F](x) \geq (2\pi)^{d/2} \left(\frac{2\varepsilon}{2+\varepsilon}\right)^d \frac{\|G_{0, \frac{2\varepsilon}{2+\varepsilon}} F\|_{\mathcal{M}_+(\mathbb{R}^d)}^2}{\|F\|_{\mathcal{M}_+(\mathbb{R}^d)}} e^{-\frac{\alpha}{2}|x|^2} e^{-\frac{2+\varepsilon}{4}|x|^2}, $$

for every $x \in \mathbb{R}^d$. Taking $\varepsilon > 0$ small enough and an appropriate constant $c > 0$ yields (4.2). Finally, taking derivatives on (1.2) we obtain

$$ \begin{aligned} \nabla \mathcal{T}[F](x) &= -\alpha x \mathcal{T}[F](x) - \frac{e^{-\frac{\alpha}{2}|x|^2}}{\|F\|_{L^1(\mathbb{R}^d)}} \int_{\mathbb{R}^{2d}} \left(x - \frac{x_1+x_2}{2}\right) G\left(x - \frac{x_1+x_2}{2}\right) F(dx_1) F(dx_2) \\ &= -(1+\alpha) x \mathcal{T}[F](x) - \frac{e^{-\frac{\alpha}{2}|x|^2}}{\|F\|_{L^1(\mathbb{R}^d)}} \int_{\mathbb{R}^{2d}} \left(\frac{x_1+x_2}{2}\right) G\left(x - \frac{x_1+x_2}{2}\right) F(dx_1) F(dx_2). \end{aligned} $$

Therefore, dividing by $\mathcal{T}[F](x)$ and taking norms yield the following estimate

$$ |\nabla \log \mathcal{T}[F](x)| \leq (1+\alpha)|x| + \frac{1}{\sqrt{2}\|F\|_{L^1(\mathbb{R}^d)}} \frac{e^{-\frac{\alpha}{2}|x|^2}}{\mathcal{T}[F](x)} \int_{\mathbb{R}^{2d}} |(x_1, x_2)| G\left(x - \frac{x_1+x_2}{2}\right) F(dx_1)F(dx_2), \qquad (4.4) $$

for every $x \in \mathbb{R}^d$. Consider any $R > 0$ and split the integral in the right hand side of (4.4) into the subsets $\{(x_1, x_2) \in \mathbb{R}^{2d} : |(x_1, x_2)| \leq \sqrt{2}R|x|\}$ and $\{(x_1, x_2) \in \mathbb{R}^{2d} : |(x_1, x_2)| > \sqrt{2}R|x|\}$. For the first subset, we use the definition (1.2) of the operator $\mathcal{T}$. For the second subset, we use the uniform bound
(ii) (Lower control of tails) There exists $c = c(\alpha, \underline{\sigma}^2, F) > 0$ such that

$$ \mathcal{T}[F](x) \geq c G_{0, \underline{\sigma}^2}(x), \qquad (4.2) $$

for every $x \in \mathbb{R}^d$.

(iii) (Control of log-derivative) Assume that $F$ is absolutely continuous with respect to the Lebesgue measure and it has Gaussian tail, i.e., there exists $\sigma^2 \in \mathbb{R}_+^*$ and $C' \in \mathbb{R}_+^*$ such that $F(x) \le C' G_{0, \sigma^2}$ for every $x \in \mathbb{R}^d$. Then, there exists $C'' = C''(\alpha, \underline{\sigma}^2, \sigma^2, C', F) > 0$ such that

$$ |\nabla \log \mathcal{T}[F](x)| \le C''(1 + |x|), \qquad (4.3) $$

for every $x \in \mathbb{R}^d$.

*Proof.* First, notice that $F \in \mathcal{M}_+(\mathbb{R}^d)$ is any generic distribution but $G \in L^\infty(\mathbb{R}^d)$. Then, $\mathcal{B}[F] \in L^\infty(\mathbb{R}^d)$ and, by definition of $\mathcal{T}$ in (1.2), we obtain

$$ \mathcal{T}[F](x) \le \|\mathcal{B}[F]\|_{L^\infty(\mathbb{R}^d)} e^{-\frac{\alpha}{2}|x|^2} = \|\mathcal{B}[F]\|_{L^\infty(\mathbb{R}^d)} (2\pi\alpha^{-1})^{d/2} G_{0, \underline{\sigma}^2}(x), $$

for every $x \in \mathbb{R}^d$. Hence, (4.1) holds for appropriate $C > 0$. Second, note that

$$ \begin{aligned} \mathcal{T}[F](x) &= \frac{e^{-\frac{\alpha}{2}|x|^2}}{(2\pi)^{d/2} \|F\|_{\mathcal{M}_+(\mathbb{R}^d)}} \int_{\mathbb{R}^{2d}} \exp\left(-\frac{1}{2}\left|x - \frac{x_1+x_2}{2}\right|^2\right) F(dx_1)F(dx_2) \\ &= \frac{e^{-\frac{\alpha}{2}|x|^2}}{(2\pi)^{d/2} \|F\|_{\mathcal{M}_+(\mathbb{R}^d)}} \int_{\mathbb{R}^{2d}} \exp\left(-\frac{1}{2}|x|^2 - \frac{1}{8}|x_1+x_2|^2 + \frac{1}{2}x \cdot (x_1+x_2)\right) F(dx_1)F(dx_2) \\ &\ge \frac{e^{-\frac{\alpha}{2}|x|^2}}{(2\pi)^{d/2} \|F\|_{\mathcal{M}_+(\mathbb{R}^d)}} \int_{\mathbb{R}^{2d}} \exp\left(-\frac{2+\varepsilon}{4}|x|^2 - \frac{2+\varepsilon}{8\varepsilon}|x_1+x_2|^2\right) F(dx_1)F(dx_2) \\ &\ge \frac{e^{-\frac{\alpha}{2}|x|^2}}{(2\pi)^{d/2} \|F\|_{\mathcal{M}_+(\mathbb{R}^d)}} \int_{\mathbb{R}^{2d}} \exp\left(-\frac{2+\varepsilon}{4}|x|^2 - \frac{2+\varepsilon}{4\varepsilon}|x_1|^2 - \frac{2+\varepsilon}{4\varepsilon}|x_2|^2\right) F(dx_1)F(dx_2), \end{aligned} $$

for every $\varepsilon > 0$, where in the third and last lines we have used Cauchy–Schwarz’s and Young’s inequalities. Consequently, we obtain that

$$ \mathcal{T}[F](x) \ge (2\pi)^{d/2} \left(\frac{2\varepsilon}{2+\varepsilon}\right)^d \frac{\|G_{0, \frac{2\varepsilon}{2+\varepsilon}} F\|_{\mathcal{M}_+(\mathbb{R}^d)}^2}{\|F\|_{\mathcal{M}_+(\mathbb{R}^d)}} e^{-\frac{\alpha}{2}|x|^2} e^{-\frac{2+\varepsilon}{4}|x|^2}, $$

for every $x \in \mathbb{R}^d$. Taking $\varepsilon > 0$ small enough and an appropriate constant $c > 0$ yields (4.2). Finally, taking derivatives on (1.2) we obtain

$$ \begin{aligned} \nabla \mathcal{T}[F](x) &= -\alpha x \mathcal{T}[F](x) - \frac{e^{-\frac{\alpha}{2}|x|^2}}{\|F\|_{L^1(\mathbb{R}^d)}} \int_{\mathbb{R}^{2d}} \left(x - \frac{x_1+x_2}{2}\right) G\left(x - \frac{x_1+x_2}{2}\right) F(dx_1) F(dx_2) \\ &= -(1+\alpha) x \mathcal{T}[F](x) - \frac{e^{-\frac{\alpha}{2}|x|^2}}{\|F\|_{L^1(\mathbb{R}^d)}} \int_{\mathbb{R}^{2d}} \left(\frac{x_1+x_2}{2}\right) G\left(x - \frac{x_1+x_2}{2}\right) F(dx_1) F(dx_2). \end{aligned} $$

Therefore, dividing by $\mathcal{T}[F](x)$ and taking norms yield the following estimate

$$ |\nabla \log \mathcal{T}[F](x)| \le (1+\alpha)|x| + \frac{1}{\sqrt{2}\|F\|_{L^1(\mathbb{R}^d)}} \frac{e^{-\frac{\alpha}{2}|x|^2}}{\mathcal{T}[F](x)} \int_{\mathbb{R}^{2d}} |(x_1, x_2)| G\left(x - \frac{x_1+x_2}{2}\right) F(dx_1)F(dx_2), \qquad (4.4) $$

for every $x \in \mathbb{R}^d$. Consider any $R > 0$ and split the integral in the right hand side of (4.4) into the subsets $\{(x_1, x_2) \in \mathbb{R}^{2d} : |(x_1, x_2)| \le \sqrt{2}R|x|\}$ and $\{(x_1, x_2) \in \mathbb{R}^{2d} : |(x_1, x_2)| > \sqrt{2}R|x|\}$. For the first subset, we use the definition (1.2) of the operator $\mathcal{T}$. For the second subset, we use the uniform bound
of $G$ along with the assumptions on $F$ and the preceding lower bound (4.2) of $\mathcal{T}[F]$. Then, we obtain

$$
\begin{align*}
|\nabla \log \mathcal{T}[F](x)| &\le (1+\alpha+R)|x| + \frac{C^2 \sigma^d e^{-\frac{\alpha}{2}|x|^2} e^{\frac{1}{2\sigma^2}|x|^2}}{\sqrt{2} c \|F\|_{L^1(\mathbb{R}^d)}} \int_{|(x_1,x_2)|>\sqrt{2}R|x|} |(x_1,x_2)| e^{-\frac{1}{2\sigma^2}|(x_1,x_2)|^2} dx_1 dx_2 \\
&= (1+\alpha+R)|x| + \frac{C^2 \sigma^d e^{-\frac{\alpha}{2}|x|^2} e^{\frac{1}{2\sigma^2}|x|^2}}{\sqrt{2} c \|F\|_{L^1(\mathbb{R}^d)}} \int_{\sqrt{2}R|x|}^{+\infty} r^{2d} e^{-\frac{1}{2\sigma^2}r^2} dr \\
&\le (1+\alpha+R)|x| + \frac{C^2 \sigma^d e^{-\frac{\alpha}{2}|x|^2} e^{\frac{1}{2\sigma^2}|x|^2}}{\sqrt{2} c \|F\|_{L^1(\mathbb{R}^d)}} \left(\frac{8d-4}{e}\right)^{d-\frac{1}{2}} \sigma^{2d-1} \int_{\sqrt{2}R|x|}^{+\infty} r e^{-\frac{1}{4\sigma^2}r^2} dr \\
&= (1+\alpha+R)|x| + \frac{\sqrt{2} C^2 \sigma^d \sigma^{2d+1}}{c \|F\|_{L^1(\mathbb{R}^d)}} \left(\frac{8d-4}{e}\right)^{d-\frac{1}{2}} e^{-\frac{1}{2}\left(\alpha+\frac{R^2}{\sigma^2}-\frac{1}{\sigma^2}\right)|x|^2},
\end{align*}
$$

for every $x \in \mathbb{R}^d$, where in the third line we have used the inequality

$$
r^{2d-1} \le \left( \frac{2d-1}{e} \right)^{d-\frac{1}{2}} \left( \frac{\sigma}{\varepsilon} \right)^{2d-1} \exp \left( \frac{\varepsilon}{2\sigma^2} r^2 \right),
$$

for every $r > 0$ and $\varepsilon > 0$ in the particular case of $\varepsilon = 1/2$. Taking $R > 0$ large enough ends the proof. $\square$

The previous result can be iterated to obtain similar properties for solutions {$F_n$}$_{n \in \mathbb{N}}$ to the time-
discrete problem (1.1) issued at a generic initial datum $F_0 \in \mathcal{M}_+(\mathbb{R}^d)$. In particular, we note that
Gaussian tails emerge instantaneously after the first generation $n = 1$ and linear growth of the log-
derivative is guaranteed after the second generation $n = 2$.

**Corollary 4.2** (Propagation of tails). Assume that $\alpha \in \mathbb{R}_+^*$, 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)$ and set $\bar{\sigma}_1^2, \underline{\sigma}_1^2 \in \mathbb{R}_+^*$ by

$$
\bar{\sigma}_1^2 := \frac{1}{\alpha}, \quad \underline{\sigma}_1^2 \in \left(0, \frac{1}{1+\alpha}\right).
$$

Define the associated sequences of variances {$\bar{\sigma}_n^2$}$_{n \in \mathbb{N}}$ and {$\underline{\sigma}_n^2$}$_{n \in \mathbb{N}}$ by the following recursive relations

$$
\frac{1}{\bar{\sigma}_{n+1}^2} := \alpha + \frac{1}{1 + \frac{\bar{\sigma}_n^2}{2}}, \quad \frac{1}{\underline{\sigma}_{n+1}^2} := \alpha + \frac{1}{1 + \frac{\underline{\sigma}_n^2}{2}}, \qquad (4.5)
$$

for every $n \in \mathbb{N}$. Then, the following properties are fulfilled:

(i) (Upper control of tails of $F_n$) There exist $C_n = C_n(\alpha, F_0) > 0$ such that

$$
F_n(x) \le C_n G_{0, \bar{\sigma}_n^2}(x), \qquad (4.6)
$$

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

(ii) (Lower control of tails of $F_n$) There exist $c_n = c_n(\alpha, \underline{\sigma}_1^2, F_0) > 0$ such that

$$
F_n(x) \ge c_n G_{0, \underline{\sigma}_n^2}(x), \qquad (4.7)
$$

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

(iii) (Control of log-derivative of $F_n$) There exists $D_n = D_n(\alpha, \underline{\sigma}_1^2, F_0)$ such that

$$
|\nabla \log F_n(x)| \le D_n(1 + |x|), \qquad (4.8)
$$

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

Proof. Note that once (4.6) and (4.7) are proved, we can readily apply (4.3) in Lemma 4.1 to $F = F_{n-1}$ with $n \ge 2$ (which satisfies the required upper control by a Gaussian function) and we recover the estimate (4.8) for the log-derivative of $F_n = \mathcal{T}[F_{n-1}]$. Then, we just focus on proving estimates (4.6) and (4.7) through an inductive argument. First, note that Lemma 4.1 applied to $F = F_0$ yields (4.6) and (4.7) with $n = 1$. Let us assume that (4.6) and (4.7) hold for some $n \in \mathbb{N}$ and let us prove it for $n+1$. Specifically, using the induction hypothesis we obtain that

$$
\begin{align*}
F_{n+1}(x) &= e^{-\frac{\alpha}{2}|x|^2} \mathcal{B}[F_n](x) \\
&\le \frac{C_n^2}{\|F_n\|_{L^1(\mathbb{R}^d)}} e^{-\frac{\alpha}{2}|x|^2} \left(G\left(\frac{\cdot}{2}\right) * G_{0, \bar{\sigma}_n^2} * G_{0, \underline{\sigma}_n^2}\right)(2x) \\
&= \frac{C_n^2}{\|F_n\|_{L^1(\mathbb{R}^d)}} e^{-\frac{\alpha}{2}|x|^2} G_{0, 1+\frac{\bar{\sigma}_n^2}{2}}(x) \le C_{n+1} G_{0, \bar{\sigma}_{n+1}^2}(x),
\end{align*}
$$
of $G$ along with the assumptions on $F$ and the preceding lower bound (4.2) of $\mathcal{T}[F]$. Then, we obtain

$$
\begin{align*}
|\nabla \log \mathcal{T}[F](x)| &\le (1+\alpha+R)|x| + \frac{C^2 \sigma^d e^{-\frac{\alpha}{2}|x|^2} e^{\frac{1}{2\sigma^2}|x|^2}}{\sqrt{2} c \|F\|_{L^1(\mathbb{R}^d)}} \int_{|(x_1,x_2)|>\sqrt{2}R|x|} |(x_1,x_2)| e^{-\frac{1}{2\sigma^2}|(x_1,x_2)|^2} dx_1 dx_2 \\
&= (1+\alpha+R)|x| + \frac{C^2 \sigma^d e^{-\frac{\alpha}{2}|x|^2} e^{\frac{1}{2\sigma^2}|x|^2}}{\sqrt{2} c \|F\|_{L^1(\mathbb{R}^d)}} \int_{\sqrt{2}R|x|}^{+\infty} r^{2d} e^{-\frac{1}{2\sigma^2}r^2} dr \\
&\le (1+\alpha+R)|x| + \frac{C^2 \sigma^d e^{-\frac{\alpha}{2}|x|^2} e^{\frac{1}{2\sigma^2}|x|^2}}{\sqrt{2} c \|F\|_{L^1(\mathbb{R}^d)}} \left(\frac{8d-4}{e}\right)^{d-\frac{1}{2}} \sigma^{2d-1} \int_{\sqrt{2}R|x|}^{+\infty} r e^{-\frac{1}{4\sigma^2}r^2} dr \\
&= (1+\alpha+R)|x| + \frac{\sqrt{2} C^2 \sigma^d \sigma^{2d+1}}{c \|F\|_{L^1(\mathbb{R}^d)}} \left(\frac{8d-4}{e}\right)^{d-\frac{1}{2}} e^{-\frac{1}{2}\left(\alpha+\frac{R^2}{\sigma^2}-\frac{1}{\sigma^2}\right)|x|^2},
\end{align*}
$$

for every $x \in \mathbb{R}^d$, where in the third line we have used the inequality

$$
r^{2d-1} \le \left( \frac{2d-1}{e} \right)^{d-\frac{1}{2}} \left( \frac{\sigma}{\varepsilon} \right)^{2d-1} \exp \left( \frac{\varepsilon}{2\sigma^2} r^2 \right),
$$

for every $r > 0$ and $\varepsilon > 0$ in the particular case of $\varepsilon = 1/2$. Taking $R > 0$ large enough ends the proof. $\square$

The previous result can be iterated to obtain similar properties for solutions {$F_n$}$_{n \in \mathbb{N}}$ to the time-
discrete problem (1.1) issued at a generic initial datum $F_0 \in \mathcal{M}_+(\mathbb{R}^d)$. In particular, we note that
Gaussian tails emerge instantaneously after the first generation $n = 1$ and linear growth of the log-
derivative is guaranteed after the second generation $n = 2$.

**Corollary 4.2** (Propagation of tails). Assume that $\alpha \in \mathbb{R}_+^*$, 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)$ and set $\bar{\sigma}_1^2, \underline{\sigma}_1^2 \in \mathbb{R}_+^*$ by

$$
\bar{\sigma}_1^2 := \frac{1}{\alpha}, \quad \underline{\sigma}_1^2 \in \left(0, \frac{1}{1+\alpha}\right).
$$

Define the associated sequences of variances {$\bar{\sigma}_n^2$}$_{n \in \mathbb{N}}$ and {$\underline{\sigma}_n^2$}$_{n \in \mathbb{N}}$ by the following recursive relations

$$
\frac{1}{\bar{\sigma}_{n+1}^2} := \alpha + \frac{1}{1 + \frac{\bar{\sigma}_n^2}{2}}, \quad \frac{1}{\underline{\sigma}_{n+1}^2} := \alpha + \frac{1}{1 + \frac{\underline{\sigma}_n^2}{2}}, \qquad (4.5)
$$

for every $n \in \mathbb{N}$. Then, the following properties are fulfilled:

(i) (Upper control of tails of $F_n$) There exist $C_n = C_n(\alpha, F_0) > 0$ such that

$$
F_n(x) \le C_n G_{0, \bar{\sigma}_n^2}(x), \qquad (4.6)
$$

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

(ii) (Lower control of tails of $F_n$) There exist $c_n = c_n(\alpha, \underline{\sigma}_1^2, F_0) > 0$ such that

$$
F_n(x) \ge c_n G_{0, \underline{\sigma}_n^2}(x), \qquad (4.7)
$$

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

(iii) (Control of log-derivative of $F_n$) There exists $D_n = D_n(\alpha, \underline{\sigma}_1^2, F_0)$ such that

$$
|\nabla \log F_n(x)| \le D_n(1 + |x|), \qquad (4.8)
$$

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

Proof. Note that once (4.6) and (4.7) are proved, we can readily apply (4.3) in Lemma 4.1 to $F = F_{n-1}$ with $n \ge 2$ (which satisfies the required upper control by a Gaussian function) and we recover the estimate (4.8) for the log-derivative of $F_n = \mathcal{T}[F_{n-1}]$. Then, we just focus on proving estimates (4.6) and (4.7) through an inductive argument. First, note that Lemma 4.1 applied to $F = F_0$ yields (4.6) and (4.7) with $n = 1$. Let us assume that (4.6) and (4.7) hold for some $n \in \mathbb{N}$ and let us prove it for $n+1$. Specifically, using the induction hypothesis we obtain that

$$
\begin{align*}
F_{n+1}(x) &= e^{-\frac{\alpha}{2}|x|^2} \mathcal{B}[F_n](x) \\
&\le \frac{C_n^2}{\|F_n\|_{L^1(\mathbb{R}^d)}} e^{-\frac{\alpha}{2}|x|^2} \left( G\left(\frac{\cdot}{2}\right) * G_{0, \bar{\sigma}_n^2} * G_{0, \underline{\sigma}_n^2} \right) (2x) \\
&= \frac{C_n^2}{\|F_n\|_{L^1(\mathbb{R}^d)}} e^{-\frac{\alpha}{2}|x|^2} G_{0, 1+\frac{\bar{\sigma}_n^2}{2}}(x) \le C_{n+1} G_{0, \bar{\sigma}_{n+1}^2}(x),
\end{align*}
$$
for every $x \in \mathbb{R}^d$ and appropriate $C_{n+1}$, where we have used Lemma 3.2 for the stability of Gaussian under convolutions and the definition (4.5) of $\bar{\sigma}_{n+1}^2$. This proves (4.6) for $n+1$ and a similar argument yields the lower estimate (4.7). $\square$

By Lemma 3.8 we note that the sequence of variances $\{\bar{\sigma}_n^2\}_{n \in \mathbb{N}}$ and $\{\underline{\sigma}_n^2\}_{n \in \mathbb{N}}$ relax towards the asymptotic value $\sigma_\alpha^2$ with a geometric convergence rate. This is consistent with our main result in Theorem 1.1 and Corollary 1.2. In fact, this suggests that any solution $\{F_n\}_{n \in \mathbb{N}}$ to the time-discrete problem (1.1) must relax towards the asymptotic profile $F_\alpha$. In addition, recall that for any solution $(\lambda, F)$ of the eigenproblem (1.11) we recover a particular solution of the time-discrete problem (1.1) via the ansatz (1.10), i.e., $F_n = \lambda^n F$. Hence, we expect that $(\lambda_\alpha, F_\alpha)$ must indeed be the unique solution to the eigenvalue problem (1.11). However, we are still far from characterizing the full long-term dynamics for the solution $\{F_n\}_{n \in \mathbb{N}}$ in Theorem 1.1 and the uniqueness result in Corollary 1.2. Namely, the above coefficients $C_n$ and $c_n$ contains crucial information about the balance of mass and their values in Corollary 4.2 are not necessarily optimal. Indeed, note that they are given by explicit recursive formulas

$$C_{n+1} = (1 - \alpha \bar{\sigma}_n^2)^{d/2} \frac{C_n^2}{\|F_n\|_{L^1(\mathbb{R}^d)}}, \quad c_{n+1} = (1 - \alpha \bar{\sigma}_n^2)^{d/2} \frac{c_n^2}{\|F_n\|_{L^1(\mathbb{R}^d)}},$$

but they require further knowledge about the behavior of $\|F_n\|_{L^1(\mathbb{R}^d)}$. In the next paragraph, we overcome this important issue of scaling factor by focusing on renormalized profiles, for which we prove uniform propagation of moments.

**4.2. Propagation of quadratic and exponential moments.** In this section we analyze the propagation of quadratic and exponential moments along the normalized profiles $F_n/\|F_n\|_{L^1(\mathbb{R}^d)}$, where $\{F_n\}_{n \in \mathbb{N}}$ is a solution to (1.1). As we anticipated in Section 2.4, we need to move from uniform estimates of the high-dimensional integral (2.11) to estimates in an average sense. This will be precisely the point in proof of the main Theorem 1.1 in Section 6.2 where a suitable control of moments will be crucial, as sketched in the overall map in Figure 2. Before addressing the case of generic initial data $F_0 \in \mathcal{M}_+(\mathbb{R}^d)$, we illustrate the particular explicit case of Gaussian initial data $F_0 = m_0 G_{\mu_0, \sigma_0^2}$.

**Remark 4.3 (Gaussian case).** Consider any $m_0 \in \mathbb{R}_+$, $\mu_0 \in \mathbb{R}^d$ and $\sigma_0^2 \in \mathbb{R}_+$, set $\theta \in \mathbb{R}_+^*$ with $\theta < \frac{1}{2\max\{\sigma_0^2, \sigma_\alpha^2\}}$ and consider the solution $\{F_n\}_{n \in \mathbb{N}}$ to the time-discrete problem (1.1) with initial datum $F_0 = m_0 G_{\mu_0, \sigma_0^2}$. Then, by direct calculation we obtain

$$\int_{\mathbb{R}^d} |x|^2 \frac{F_n(x)}{\|F_n\|_{L^1(\mathbb{R}^d)}} dx = |\mu_n|^2 + \sigma_n^2, \\ \int_{\mathbb{R}^d} e^{\theta|x|^2} \frac{F_n(x)}{\|F_n\|_{L^1(\mathbb{R}^d)}} dx = \frac{e^{\frac{\theta}{1-2\theta\sigma_n^2} |\mu_n|^2}}{(1-2\theta\sigma_n^2)^{d/2}},$$

Here, $\{\mu_n\}_{n \in \mathbb{N}}$ and $\{\sigma_n^2\}_{n \in \mathbb{N}}$ are the mean and variance of $\{F_n\}_{n \in \mathbb{N}}$ according to formula (3.5) in Proposition 3.7. As studied in the proof of Corollary 3.9 we have the uniform bounds

$$|\mu_n| \le |\mu_0| \quad \text{and} \quad \min\{\sigma_0^2, \sigma_\alpha^2\} \le \sigma_n^2 \le \max\{\sigma_0^2, \sigma_\alpha^2\},$$

for every $n \in \mathbb{N}$. Then, we obtain that both the quadratic and exponential moments are uniformly bounded.

In the sequel, we explore the case of generic initial data $F_0 \in \mathcal{M}_+(\mathbb{R})$. For simplicity of notation, we define the following conservative operator associated with $\mathcal{T}$.

**Definition 4.4 (Normalization of $\mathcal{T}$).**

$$S[F] := \frac{\mathcal{T}[F]}{\|\mathcal{T}[F]\|_{L^1(\mathbb{R}^d)}}, \quad F \in \mathcal{M}_+(\mathbb{R}^d) \setminus \{0\}.$$

Therefore, we obtain the following control of quadratic moments under $\mathcal{S}$.

**Lemma 4.5 (Control of quadratic moments).** Assume that $\alpha \in \mathbb{R}_+^*$ and set any parameter $\eta \in (\frac{1}{2(1+\alpha)^2}, 1)$. Then, there exists a constant $M = M(\alpha, \eta) \in \mathbb{R}_+^*$ such that

$$\int_{\mathbb{R}^d} |x|^2 S[F](x) dx \le M + \eta \int_{\mathbb{R}^d} |x|^2 \frac{F(dx)}{\|F\|_{\mathcal{M}(\mathbb{R}^d)}},$$

for any measure $F \in \mathcal{M}_+(\mathbb{R}^d)$.
for every $x \in \mathbb{R}^d$ and appropriate $C_{n+1}$, where we have used Lemma 3.2 for the stability of Gaussian under convolutions and the definition (4.5) of $\bar{\sigma}_{n+1}^2$. This proves (4.6) for $n+1$ and a similar argument yields the lower estimate (4.7). $\square$

By Lemma 3.8 we note that the sequence of variances $\{\bar{\sigma}_n^2\}_{n \in \mathbb{N}}$ and $\{\underline{\sigma}_n^2\}_{n \in \mathbb{N}}$ relax towards the asymptotic value $\sigma_\alpha^2$ with a geometric convergence rate. This is consistent with our main result in Theorem 1.1 and Corollary 1.2. In fact, this suggests that any solution $\{F_n\}_{n \in \mathbb{N}}$ to the time-discrete problem (1.1) must relax towards the asymptotic profile $F_\alpha$. In addition, recall that for any solution $(\lambda, F)$ of the eigenproblem (1.11) we recover a particular solution of the time-discrete problem (1.1) via the ansatz (1.10), i.e., $F_n = \lambda^n F$. Hence, we expect that $(\lambda_\alpha, F_\alpha)$ must indeed be the unique solution to the eigenvalue problem (1.11). However, we are still far from characterizing the full long-term dynamics for the solution $\{F_n\}_{n \in \mathbb{N}}$ in Theorem 1.1 and the uniqueness result in Corollary 1.2. Namely, the above coefficients $C_n$ and $c_n$ contains crucial information about the balance of mass and their values in Corollary 4.2 are not necessarily optimal. Indeed, note that they are given by explicit recursive formulas

$$C_{n+1} = (1 - \alpha \bar{\sigma}_n^2)^{d/2} \frac{C_n^2}{\|F_n\|_{L^1(\mathbb{R}^d)}}, \quad c_{n+1} = (1 - \alpha \bar{\sigma}_n^2)^{d/2} \frac{c_n^2}{\|F_n\|_{L^1(\mathbb{R}^d)}},$$

but they require further knowledge about the behavior of $\|F_n\|_{L^1(\mathbb{R}^d)}$. In the next paragraph, we overcome this important issue of scaling factor by focusing on renormalized profiles, for which we prove uniform propagation of moments.

**4.2. Propagation of quadratic and exponential moments.** In this section we analyze the propagation of quadratic and exponential moments along the normalized profiles $F_n/\|F_n\|_{L^1(\mathbb{R}^d)}$, where $\{F_n\}_{n \in \mathbb{N}}$ is a solution to (1.1). As we anticipated in Section 2.4, we need to move from uniform estimates of the high-dimensional integral (2.11) to estimates in an average sense. This will be precisely the point in proof of the main Theorem 1.1 in Section 6.2 where a suitable control of moments will be crucial, as sketched in the overall map in Figure 2. Before addressing the case of generic initial data $F_0 \in \mathcal{M}_+(\mathbb{R}^d)$, we illustrate the particular explicit case of Gaussian initial data $F_0 = m_0 G_{\mu_0, \sigma_0^2}$.

**Remark 4.3 (Gaussian case).** Consider any $m_0 \in \mathbb{R}_+$, $\mu_0 \in \mathbb{R}^d$ and $\sigma_0^2 \in \mathbb{R}_+$, set $\theta \in \mathbb{R}_+^*$ with $\theta < \frac{1}{2\max\{\sigma_0^2, \sigma_\alpha^2\}}$ and consider the solution $\{F_n\}_{n \in \mathbb{N}}$ to the time-discrete problem (1.1) with initial datum $F_0 = m_0 G_{\mu_0, \sigma_0^2}$. Then, by direct calculation we obtain

$$\int_{\mathbb{R}^d} |x|^2 \frac{F_n(x)}{\|F_n\|_{L^1(\mathbb{R}^d)}} dx = |\mu_n|^2 + \sigma_n^2, \\ \int_{\mathbb{R}^d} e^{\theta|x|^2} \frac{F_n(x)}{\|F_n\|_{L^1(\mathbb{R}^d)}} dx = \frac{e^{\frac{\theta}{1-2\theta\sigma_n^2} |\mu_n|^2}}{(1-2\theta\sigma_n^2)^{d/2}},$$

Here, $\{\mu_n\}_{n \in \mathbb{N}}$ and $\{\sigma_n^2\}_{n \in \mathbb{N}}$ are the mean and variance of $\{F_n\}_{n \in \mathbb{N}}$ according to formula (3.5) in Proposition 3.7. As studied in the proof of Corollary 3.9 we have the uniform bounds

$$|\mu_n| \le |\mu_0| \quad \text{and} \quad \min\{\sigma_0^2, \sigma_\alpha^2\} \le \sigma_n^2 \le \max\{\sigma_0^2, \sigma_\alpha^2\},$$

for every $n \in \mathbb{N}$. Then, we obtain that both the quadratic and exponential moments are uniformly bounded.

In the sequel, we explore the case of generic initial data $F_0 \in \mathcal{M}_+(\mathbb{R})$. For simplicity of notation, we define the following conservative operator associated with $\mathcal{T}$.

**Definition 4.4 (Normalization of $\mathcal{T}$).**

$$S[F] := \frac{\mathcal{T}[F]}{\|\mathcal{T}[F]\|_{L^1(\mathbb{R}^d)}}, \quad F \in \mathcal{M}_+(\mathbb{R}^d) \setminus \{0\}.$$

Therefore, we obtain the following control of quadratic moments under $\mathcal{S}$.

**Lemma 4.5 (Control of quadratic moments).** Assume that $\alpha \in \mathbb{R}_+^*$ and set any parameter $\eta \in (\frac{1}{2(1+\alpha)^2}, 1)$. Then, there exists a constant $M = M(\alpha, \eta) \in \mathbb{R}_+^*$ such that

$$\int_{\mathbb{R}^d} |x|^2 S[F](x) dx \le M + \eta \int_{\mathbb{R}^d} |x|^2 \frac{F(dx)}{\|F\|_{\mathcal{M}(\mathbb{R}^d)}},$$

for any measure $F \in \mathcal{M}_+(\mathbb{R}^d)$.
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). \quad (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. \quad (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), \quad (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). \quad (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)}}}. \quad (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.