Vincent Calvez
email: vincent.calvez@math.cnrs.fr

Thomas Lepoutre
email: thomas.lepoutre@inria.fr

David Poyato
email: davidpoyato@ugr.es

ERGODICITY OF THE FISHER INFINITESIMAL MODEL WITH QUADRATIC SELECTION

Keywords: . 2020 Mathematics Subject Classification. 35B40, 35P30, 35Q92, 47G20, 92D15 Integro-differential equations, asymptotic behavior, nonlinear spectral theory, quantitative genetics

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.

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 [START_REF] Fisher | The correlation between relatives on the supposition of mendelian inheritance[END_REF]. 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 reconcilling 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 [START_REF] Barton | The infinitesimal model: Definition, derivation, and implications[END_REF] 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 ∈ R d in a population undergoing sexual reproduction and the effect of natural selection. Specifically, starting at any initial configuration F 0 ∈ M + (R d ) of trait distribution, we analyze the long term dynamics of trait distributions {F n } n∈N across successive generations n ∈ N, which solves the following recursion

F n = T [F n-1 ], (1.1)
for any n ∈ N. The operator T encodes the balanced effect of sexual reproduction in the population (under the infinitesimal model) and natural selection. Specifically, T is defined by

T [F ] := e -m B[F ], (1.2)
for any trait distribution F ∈ M + (R d ). On the one hand, m = m(x) ≥ 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 B is chosen to be Fisher's infinitesimal operator and it takes the form

B[F ](x) := R 2d G x - x 1 + x 2 2 F (x 1 )F (x 2 ) R d F (x ′ ) dx ′ dx 1 dx 2 , x ∈ R d , (1.3)
for any trait distribution F ∈ M + (R d ). Here G = G(x) is a probability density (the mixing kernel) and the factor G(x -x1+x2

2

) represents the transition probability that two given individuals with trait values x 1 , x 2 ∈ R d will mate and yield a descendant with trait value x ∈ R d . In other words the resulting trait distributes around the mean value x1+x2 2 of parents' trait with law G. By definition, B[F ] quantifies the number of births after all possible matches of any couple of individuals according to the trait distribution F . Altogether, 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 B has recently pulled the attention of both the applied and more theoretical communities, cf. [START_REF] Barton | The infinitesimal model: Definition, derivation, and implications[END_REF][START_REF] Calvez | Asymptotic analysis of a quantitative genetics model with nonlinear integral operator[END_REF][START_REF] Mirrahimi | Dynamics of sexual populations structured by a space variable and a phenotypical trait[END_REF][START_REF] Patout | The cauchy problem for the infinitesimal model in the regime of small variance[END_REF][START_REF] Raoul | Macroscopic limit from a structured population model to the Kirkpatrick-Barton model[END_REF][START_REF]Exponential convergence to a steady-state for a population genetics model with sexual reproduction and selection[END_REF]. In this paper we shall restrict to Gaussian mixing kernel and quadratic selection function, i.e.,

G(x) := 1 (2π) d/2 e -|x| 2 2 , x ∈ R d , (1.4
)

m(x) := α 2 |x| 2 , x ∈ R d , (1.5)
where α ∈ 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 [START_REF] Barton | The infinitesimal model: Definition, derivation, and implications[END_REF]. 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):

∂ t f = -m(x)f + R[f ], t ≥ 0, x ∈ R d , f (0, x) = f 0 (x), x ∈ R d , (1.6)
where, m = m(x) is again the mortality rate and R = R[f ] is the reproduction operator. Hence, R[f ](x) determines the amount of births with trait value x ∈ 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 R across generations. Many studies consider a linear reproduction operator R, associated with a probability density K = K(x) characterizing the mutational effects at birth, of the form

R[f ](x) := R d K(x -y)f (t, y) dy, x ∈ R d .
(1.7)

The factor K(x -y) determines the probability that an individual with trait value y ∈ R d produces a descendant with trait value x ∈ 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 [START_REF] Barles | Concentrations and constrained Hamilton-Jacobi equations arising in adaptive dynamics[END_REF][START_REF] Barles | Concentration in Lotka-Volterra parabolic or integral equations: a general convergence result[END_REF][START_REF] Diekmann | The dynamics of adaptation: an ilumination example and a Hamilton-Jacobi approach[END_REF], including an additional density-dependent competition term that makes the analysis non standard (see e.g. [START_REF] Calvez | Uniqueness of the viscosity solution of a constrained Hamilton-Jacobi equation[END_REF] 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 B in (1.3) as reproduction operator R = B. Several asymptotic regimes have been addressed: large reproduction rate [START_REF] Mirrahimi | Dynamics of sexual populations structured by a space variable and a phenotypical trait[END_REF][START_REF] Raoul | Macroscopic limit from a structured population model to the Kirkpatrick-Barton model[END_REF][START_REF]Exponential convergence to a steady-state for a population genetics model with sexual reproduction and selection[END_REF], small variance asymptotics [START_REF] Calvez | Asymptotic analysis of a quantitative genetics model with nonlinear integral operator[END_REF][START_REF] Patout | The cauchy problem for the infinitesimal model in the regime of small variance[END_REF]. 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 α in the present work.

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

f (t, x) = e λt F (x), (t, x) ∈ R + × R d , (1.8)
The parameter λ ∈ R represents the rate at which the number of individuals grows (if λ ≥ 0) or decreases (if λ ≤ 0), and F = F (x) ≥ 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 (λ, F ):

λF (x) + m(x)F (x) = B[F ](x), x ∈ R d , R d F (x ′ ) dx ′ = 1.
(1.9)

Note that the operator B is genuinely non-linear so that methods based on the Krein-Rutman theory or maximum principles cannot be applied straightforwardly, see [START_REF] Berestycki | Persistence criteria for populations with non-local dispersion[END_REF][START_REF] Coville | On eigenvalue problems arising from nonlocal diffusion models[END_REF] and the references therein for the linear case. Further, usual extensions of the Krein-Rutman theory to 1-homogeneous operators [START_REF] Mahadevan | A note on a non-linear Krein-Rutman theorem[END_REF] cannot be applied neither because 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 → (λ + m) -1 B[F ] over an appropriate cone of L 1 (R d ) that is conserved by the nonlinear operator, see [START_REF] Bourgeron | Existence of recombination-selection equilibria for sexual populations[END_REF]. However, uniqueness cannot be achieved by this method. Moreover, it has been proven in [START_REF] Calvez | Asymptotic analysis of a quantitative genetics model with nonlinear integral operator[END_REF]Corollary 1.5] that several equilibrium states (λ, 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 R = B and a trait dependent fecundity rate [START_REF]Exponential convergence to a steady-state for a population genetics model with sexual reproduction and selection[END_REF]. 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 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 [START_REF]Exponential convergence to a steady-state for a population genetics model with sexual reproduction and selection[END_REF]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 [START_REF] Bolley | Convergence to equilibrium in Wasserstein distance for Fokker-Planck equations[END_REF]. 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 ≡ 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) = λ n F (x), (n, x) ∈ N × R d .
(1.10)

Again, λ ∈ R * + is the rate of growth (if λ ≥ 1) or decrease (if λ ≤ 1) of individuals, and F = F (x) is an unknown probability density. This yields the following generalized eigenproblem for the couple (λ, F ):

λF (x) = T [F ](x), x ∈ R d , R d F (x ′ ) dx ′ = 1.
( 1.11) In this paper we aim to address the following questions:

(Q1) Does the eigenproblem (1.11) have a unique solution (λ α , F α ), with λ α ∈ R and F α being a probability measure, for each α ∈ R * + ? (Q2) Consider any generic initial datum F 0 ∈ M + (R d ) and its associated solution {F n } n∈N of the time-discrete problem (1.1). Do the renormalized profiles F n /∥F n ∥ L 1 (R d ) converge to the unique steady profile F α solving (1.11) when n → ∞? We shall prove that the answer to both questions is affirmative. It stands to reason that similar existence and local uniqueness results like in [START_REF] Bourgeron | Existence of recombination-selection equilibria for sexual populations[END_REF][START_REF] Calvez | Asymptotic analysis of a quantitative genetics model with nonlinear integral operator[END_REF] 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 T , leading to quantitative estimates for the relaxation of profiles {F n } n∈N towards F α . First, we prove that an explicit (Gaussian) solution to the eigenproblem (1.11) exists. Second, by computing n iterations of the operator T , that is F n = 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 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∈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 α 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∈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 α ∈ R * + and set any initial non-negative measure

F 0 ∈ M + (R d ).
The solution {F n } n∈N to the time-discrete problem (1.1) verifies that the growth rates

∥F n ∥ L 1 (R d ) /∥F n-1 ∥ L 1 (R d )
relax towards λ α and the normalized profiles

F n /∥F n ∥ L 1 (R d ) relax towards F α with λ α := 1 + α 1 + σ 2 α 2 -d 2 , F α := G 0,σ 2 α , (1.12)
and the variance σ 2 α ∈ R * + is the unique positive root of the equation

1 σ 2 α = α + 1 1 + σ 2 α 2 , i.e., σ 2 α = (1 + 2α) 2 + 8α -(1 + 2α) 2α . (1.13)
Specifically, for any ε ∈ R * + there exists a sufficiently large

C ε ∈ R * + such that D KL F n ∥F n ∥ L 1 (R d ) F α ≤ C ε ((2 k α ) 2 + ε) n , ∥F n ∥ L 1 (R d ) ∥F n-1 ∥ L 1 (R d ) -λ α ≤ C ε (2 k α +ε) n ,
for any n ∈ N, where D KL is the Kullback-Leibler divergence (or relative entropy), i.e.,

D KL (P ∥Q) := R d P (x) log P (x) Q(x) dx, (1.14)
for any P, Q ∈ L 1 + (R d ) ∩ P(R d ), and the coefficient k α ∈ (0, 1 2 ) reads k α = σ 2 α 2 + σ 2 α , i.e., (3 + 2α) -(1 + 2α) 2 + 8α 4 . (1.15)
Note that any solution (λ, 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 (λ α , F α ) to (1.11) is Gaussian.

Corollary 1.2. Assume that α ∈ R * + , then there exists a unique solution (λ α , F α ) to the eigenproblem (1.11), with F α 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 [START_REF] Calvez | Uniform contractivity of the Fisher infinitesimal model with strongly convex selection[END_REF]. However, they imposed stringent conditions on the initial datum F 0 , that is, it should behave at infinity as the equilibrium profile F α 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 [START_REF] Raoul | Macroscopic limit from a structured population model to the Kirkpatrick-Barton model[END_REF][START_REF]Exponential convergence to a steady-state for a population genetics model with sexual reproduction and selection[END_REF], 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 α = 0, the operator T reduces to B, which is nonexpansive with respect to the quadratic Wasserstein distance, and indeed contractive over distributions with common center of mass, cf. Section 2.1. However, when α > 0, the multiplicative operator leads to an operator 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 T is non-expansive in the log-Lipschits norm ∥∇ log F F α ∥ L ∞ (R d ) for all α ≥ 0, and indeed it is contractive if α > 0. However, this contraction only gives actual information when the initial datum F 0 has identical tails to the Gaussian density F α so that initially the log-Lipschitz norm is finite. (iii) (Relative entropy distance) Note that the above log-Lipschit norm amounts to the natural L ∞ version of the relative Fisher information

∥∇ log F F α ∥ L 2 (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 α). Our form of ergodicity, as measured in Theorem 1.1, it breaks down when α = 0, simply because the operator is invariant by translation in that case, and it admits a one-parameter family F α=0 (• -µ) with µ ∈ R d of fixed Gaussian probability densities. Nevertheless, as mentioned in item (i) above, when α = 0 and we restrict to centered initial data, there is convergence to the right centered Gaussian probability density F α=0 with respect to the Wasserstein distance. It is an open problem to make both approaches meet for α = 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 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.

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 B in (1.3) and a generic multiplicative selection e -m . Although a small perturbation of the case of flat selection (m ≡ 0) could still be considered via a perturbative argument (see discussion above, [START_REF]Exponential convergence to a steady-state for a population genetics model with sexual reproduction and selection[END_REF] and also [START_REF] Raoul | Macroscopic limit from a structured population model to the Kirkpatrick-Barton model[END_REF]), 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 B. Since the fecundity rate has been normalized to 1 (see Appendix A), then B preserves the mass and center of mass, namely,

∥B[F ]∥ L 1 (R d ) = ∥F ∥ M(R d ) , R d x B[F ](x) dx = R d x F (dx), (2.1)
for any F ∈ M + (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. [START_REF]Exponential convergence to a steady-state for a population genetics model with sexual reproduction and selection[END_REF]) 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 ∈ P 2 (R d ) and they have the same center of mass. Then,

W 2 2 (B[F 1 ], B[F 2 ]) ≤ 1 2 W 2 2 (F 1 , F 2 ) . Notation Meaning Page M(R d ), M + (R d )
Signed and non-negative finite Radon measures

P(R d ), P 2 (R d )
Probability measures (with finite 2nd order moment)

P ⊗ Q ∈ M(R 2d )
Tensor product of the finite measures

P, Q ∈ M(R d ) D KL (P ∥Q) Kullback-Leibler divergence between P, Q ∈ P(R d ) W 2 (P, Q) Quadratic Wasserstein distance between P, Q ∈ P 2 (R d ) G µ,σ 2
Gaussian with mean µ ∈ R d and covariance

σ 2 I d ∈ R d×d G µ,Σ Gaussian with mean µ ∈ R d and covariance Σ ∈ R d×d N (µ, Σ) Normal with mean µ ∈ R d and covariance Σ ∈ R d×d E[X] Expectation of a random variable X G(x) Gaussian mixing kernel G 0,I d m(x) Quadratic selection function α 2 |x| 2 B[F ] Fisher's infinitesimal operator T [F ] Selection-reproduction operator M[F ] Normalized multiplicative operator by e -m S[F ] Scaled selection-reproduction operator E n [F ] High-dimensional integral operators {F n } n∈N Solution to the time-evolution equation (1.1) (λ α , F α )
Gaussian solution to the non-linear eigenproblem (1.11)

σ 2 α Variance of Gaussian eigenfunction F α F = e m F F α=0 Normalized profile associated to F ∈ M + (R d ) T n , T n * , T n
Perfect binary tree, rootless tree and leafless tree L n m , L n Level m and leaves (level n) of the tree i1, i2

Parents of a node i ∈ T n of the tree

x n = (x i ) i∈T n * y n = (y i ) i∈T n *
Variables indexed by the rootless tree

z n = (z j ) j∈L n
Variables indexed by leafs

∥ x ∥ = n i=1 |x i | 2 1/2 ℓ 2 sum of Euclidean norms of x = (x 1 , . . . , x n ) ∈ R dn Φ j n (x; y n ) Lineage map from leaf j ∈ L n to root value x ∈ R d u ⊗ v ∈ R d×d Kronecker product (u i v j ) 1≤i,j≤N of vectors u, v ∈ R d k n , κ n
Sequences of coefficients in the change of variables (2 k α ) 2 Convergence rate of D KL 2 k α Convergence rate of the log-Lipschitz norm r α Relaxation rates of variances recursion Table 1. List of notations Proof. Recall the following dual characterizations of the quadratic Wasserstein distance (cf. [START_REF] Ambrosio | Gradient flows in metric spaces and in the space of probability measures[END_REF][START_REF] Villani | Optimal transport: old and new[END_REF]):

W 2 2 (F 1 , F 2 ) = inf R 2d |x -y| 2 γ(dx, dy) : γ ∈ P(R d × R d ), π 1# γ = F 1 , π 2# γ = F 2 = sup R d ϕ 1 dF 1 + R d ϕ 2 dF 2 : |ϕ 1 (x) + ϕ 2 (y)| ≤ |x -y| 2 ∀ x, y ∈ R d , (2.2)
where π i : R d × R d → R d is the projection onto the i-th component for i = 1, 2. Taking any couple (ϕ 1 , ϕ 2 ) as above and using the specific form of the operator B we obtain:

R d B[F 1 ](x) ϕ 1 (x) dx + R d B[F 2 ](y) ϕ 2 (y) dy = R 3d G x - x 1 + x 2 2 ϕ 1 (x) F 1 (dx 1 ) F 1 (dx 2 ) dx + R 3d G y - y 1 + y 2 2 ϕ 2 (y) F 2 (dy 1 ) F 2 (dy 2 ) dy = R 3d G(z) ϕ 1 z + x 1 + x 2 2 F 1 (dx 1 ) F 1 (dx 2 ) dz + R 3d G(z) ϕ 2 z + y 1 + y 2 2 F 2 (dy 1 ) F 2 (dy 2 ) dz.
Consider any transference plan γ ∈ P(R d × R d ) between (F 1 , F 2 ) as above. Specifically, we have that π 1# γ = F 1 and π 2# γ = F 2 . Then, we can gather both integrals as follows:

R d B[F 1 ](x) ϕ 1 (x) dx + R d B[F 2 ](y) ϕ 2 (y) dy = R 5d G(z) ϕ 1 z + x 1 + x 2 2 + ϕ 2 z + y 1 + y 2 2 γ(dx 1 , dy 1 ) γ(dx 2 , dy 2 ) dz.
Since the condition |ϕ 1 (x) + ϕ 2 (y)| ≤ |x -y| 2 is verified for all x, y ∈ R d , then we find

R d B[F 1 ](x) ϕ 1 (x) dx + R d B[F 2 ](y) ϕ 2 (y) dy ≤ 1 4 R 5d G(z) |(x 1 + x 2 ) -(y 1 + y 2 )| 2 γ(dx 1 , dy 1 ) γ(dx 2 , dy 2 ) dz = 1 4 R 4d |(x 1 + x 2 ) -(y 1 + y 2 )| 2 γ(dx 1 , dy 1 ) γ(dx 2 , dy 2 ) = 1 4 R 2d |x 1 -y 1 | 2 γ(dx 1 , dy 1 ) + 1 4 R 2d |x 2 -y 2 | 2 γ(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 (ϕ 1 , ϕ 2 ) and infimum over γ and using the dual characterizations (2.2) yields the result. □

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 nonexpansiveness:

W 2 (B[F 1 ], B[F 2 ]) ≤ W 2 (F 1 , F 2 ). (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 α = 0.

Corollary 2.2. Set any initial datum F 0 ∈ M + (R d ) such that R d |x| 2 F 0 (dx) < ∞ and consider the solution {F n } n∈N to the time-discrete problem (1.1) with α = 0, i.e., F n = B[F n-1 ] for all n ∈ N. Then,

∥F n ∥ L 1 (R d ) ∥F n-1 ∥ L 1 (R d ) = 1, W 2 F n ∥F n ∥ L 1 (R d ) , G µ0,2 ≲ 1 2 n/2 ,
for every n ∈ N, where µ 0 := R d x F 0 (dx). In particular, the set of stationary distributions under B is {F α=0 (• -µ) : µ ∈ R d }, where F α=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 α = 0. Indeed, by Talagrand transportation inequality for a Gaussian measure [START_REF] Talagrand | Transportation cost for Gaussian and other measures[END_REF] we obtain the relation

W 2 2 F n ∥F n ∥ L 1 (R d ) , G µ0,2 ≤ 4 D KL F n ∥F n ∥ L 1 (R d ) G µ0,2 .
However, Theorem 1.1 does not hold when α = 0, as mentioned in Remark 1.5, due to two fundamental reasons. First, when α = 0 there is translation invariance and therefore, for generic F 0 ∈ M + (R d ) one cannot expect that the normalized profiles F n /∥F n ∥ L 1 (R d ) always converge to the Gaussian F α=0 centered at the origin (contrarily to what happens when α > 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 α = 0 (contrarily to the case α > 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 α → 0 since we have 2 k α=0 = 1.

2.2.

Incompatibility of Wasserstein metric with multiplicative operators. Note that by definition (1.2), our operator T is the composition of the sexual resproduction operator 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)

.

M[F ] := e -m F ∥e -m F ∥ M(R d ) , F ∈ M + (R d ) \ {0} .
However, the latter is not Lipschitz continuous with respect to the quadratic Wasserstein metric. Hence the composition of B and 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 M in the following example.

Example 2.4. Suppose that the m ∈ C 1 + (R) is radially symmetric around the origin and set F 1 , F 2 ∈ P 2 (R) as the sum of two Dirac masses, F 1 being symmetric, and F 2 nearly symmetric. More precisely,

F 1 = 1 2 δ -h + 1 2 δ h , F 2 = 1 2 δ -h+ε + 1 2 δ h+ε ,
where h > 0 is fixed so that m ′ (h) ̸ = 0 and ε > 0 is small. Note that

M[F 1 ] = F 1 , M[F 2 ] = (1 -p ε )δ -h+ε + p ε δ h+ε ,
where p ε ∈ (0, 1) is given by p ε = e -m(h+ε) e -m(-h+ε) + e -m(h+ε) . On the one hand, we have W 2 (F 1 , F 2 ) = ε because F 2 is simply deduced from F 1 by a translation of size ε. On the other hand, assuming m ′ (h) > 0 for simplicity (a similar argument holds if m ′ (h) < 0) and taking ε > 0 sufficiently small we obtain p ε < 1 2 and

W 2 (M[F 1 ], M[F 2 ]) = ε 2 + 4h(h -ε) 1 2 -p ε ∼ ε + h m ′ (h) 1/2 ε 1/2 ,
as ε → 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 1 2 -p ε from -h + ε to h in order to equilibrate the Dirac masses in the transport plan. Since m ′ (h) ̸ = 0 this leads to the Lipschitz-discontinuity of the operator M.

In [START_REF]Exponential convergence to a steady-state for a population genetics model with sexual reproduction and selection[END_REF], 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 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 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 α . More specifically, we obtain the following result for any α ≥ 0 (but only useful when α > 0). Lemma 2.5 (Log-Lipschitz estimate). Let m be the quadratic selection function in (1.5) and consider any value α ≥ 0. Then, the following estimate holds true

∇ log T [F ] F α L ∞ (R d ) ≤ 2 k α ∇ log F F α L ∞ (R d ) ,
for any

F ∈ L 1 + (R d ) ∩ C 1 (R d ) such that ∥∇ log F F α ∥ L ∞ (R d ) < ∞.
The coefficient k α 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∈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 α > 0, we have 2 k α < 1 (cf. Figure 1), and therefore we obtain an actual contractivity estimate. There is one main drawback though: typically ∇ log Fn F α are genuinely unbounded unless F 0 and F α 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 α = 0 we have 2 k α = 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 α = 0, but we find this idea illuminating to address the case of a non-trivial effect of selection when α > 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. Proof of Lemma 2.5. Our starting point resides in the following normalization of the operator T :

T [F ](x) F α (x) = λ α R d F (x ′ ) F α(x ′ ) F α (x ′ ) dx ′ R 2d P(x; x 1 , x 2 ) F (x 1 ) F α (x 1 ) F (x 2 ) F α (x 2 ) dx 1 dx 2 , (2.4)
for all x ∈ R d , where we have

P(x; x 1 , x 2 ) = e 1 2 ( 1 σ 2 α -α)|x| 2 λ α (2π σ α ) d exp - 1 2 x - x 1 + x 2 2 2 - 1 2 σ 2 α (|x 1 | 2 + |x 2 | 2 ) . (2.5)
Above we have exploited the explicit Gaussian shape of F α to find P explicitly. Note that 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, R 2d P(x; x 1 , x 2 ) dx 1 dx 2 = 1 for all x ∈ R d because (λ α , F α ) 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 = σ 2 α /(2 + σ 2 α ) and by definition (1.15) of k α , we infer k = k α . This motivates using the change of variables

x 1 = k α x + y 1 , x 2 = k α x + y 2 , (2.6)
which appropriately centers the quadratic form at the minimum. Specifically, we obtain

1 2 x - x 1 + x 2 2 2 - 1 2 σ 2 α (|x 1 | 2 + |x 2 | 2 ) = - 1 2
y 1 + y 2 2 2 - 1 2 σ 2 α (|y 1 | 2 + |y 2 | 2 ) - 1 2 (1 -k α ) 2 + 2 k 2 α σ 2 α |x| 2 + 1 2 (1 -k α ) - 2 k α σ 2 α x • (y 1 + y 2 ), = - 1 2
y 1 + y 2 2 2 - 1 2 σ 2 α (|y 1 | 2 + |y 2 | 2 ) - 1 2 1 σ 2 α -α |x| 2 ,
where in the second line we have noticed that the coefficient of the crossed term x • (y 1 + y 2 ) vanishes due to the relation

k α = σ 2 α 2+σ 2 α
, and the coefficient of the |x| 2 factor can be reformulated as

(1 -k α ) 2 + 2 k 2 α σ 2 α = 2 2 + σ 2 α = 1 1 + σ 2 α 2 = 1 σ 2 α -α,
thanks to the implicit equation (1.13) satisfied by σ 2 α . Therefore, (2.4)-(2.5) transform into

T [F ](x) F α (x) = λ α R d F (x ′ ) F α (x ′ ) F α (x ′ ) dx ′ R 2d P(y 1 , y 2 ) F (k α x + y 1 ) F α (k α x + y 1 ) F (k α x + y 2 ) F α (k α x + y 2 ) dy 1 dy 2 , (2.7)
P(y 1 , y 2 ) = 1 λ α (2π σ α ) d exp - 1 2 y 1 + y 2 2 2 - 1 2 σ 2 α (|y 1 | 2 + |y 2 | 2 ) .
(2.8)

We remark that the new Markov transition kernel P = 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 k α . 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 T [F ] and F . Specifically, we have

∇ log T [F ](x) F α (x) = k α R 2d ∇ log F (k α x + y 1 ) F α (k α x + y 1 ) + ∇ log F (k α x + y 2 ) F α (k α x + y 2 ) ν(x; dy 1 , dy 2 ), (2.9)
where the x-dependent measures ν(x; dy 1 , dy 2 ) on the variables (y 1 , y 2 ) have the following density with respect to the Lebesgue measure:

ν(x; dy 1 , dy 2 ) dy 1 dy 2 = P(y 1 , y 2 ) F (k α x + y 1 ) F α (k α x + y 1 ) F (k α x + y 2 ) F α (k α x + y 2 ) R 2d P(y ′ 1 , y ′ 2 ) F (k α x + y ′ 1 ) F α (k α x + y ′ 1 ) F (k α x + y ′ 2 ) F α (k α 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 ∞ bounds and using that ν are probability measures on the variables (y 1 , y 2 ). □ 2.4. Brief description of our strategy. As advanced before, our strategy is based on a finer understanding of the iterations of T across generations. Specifically, it relies on a suitable reformulation of the solutions {F n } n∈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 / F α in the proof of Lemma 2.5. In this paper, we have opted for Fn = e m F n F α=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,σ 2 with σ 2 larger but arbitrarily close to σ 2 α (being σ 2 α the variance (1.13) of the expected equilibrium F α ). 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 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) ∝ R (2 n+1 -2)d P n (x; {x i }) j F0 (x j ) d{x i },
for an explicit n-step transition kernel P n with Gaussian shape:

P n (x; {x i }) := exp(-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,

Q n = 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 γ

2 |x i | 2 (with γ = σ 2 α 2+σ 2 α
)). Therefore, the minimizers of 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) ∝ F n (x) R (2 n+1 -2)d P n ({y i }) j F0 (Φ j n (x; {y i })) d{y i } . (2.10)
Above, for every leaf j the maps Φ j n defined by Φ j n (x; {y i }) := κ n x + Λ j ({y i }), where Λ j n 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:

P n ({y i }) = exp(-Q n ({y i })),
for a x-independent quadratic form 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 F n (x). Then, the dependency upon x is split into two parts: the explicit term F n (x) outside the integral (which is proven to converge towards F α by construction) and the contribution κ n x of the lineage maps Φ j n (x; {y i }) inside F0 for an appropriate sequence {κ n } n∈N ⊆ R + . The major observation here is that κ n ≪ 2 -n when α > 0 because of the strong shift towards the origin under selection. Otherwise, for α = 0 we only have κ n = 2 -n . Formula (2.10) for α > 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 → ∞. Indeed, there are still 2 n terms in the product indexed by j, but the contraction parameter κ n appears to decay fast enough to compensate them. To make this argument quantitative, let us differentiate (2.10) again to obtain

∇ log F n (x) = ∇ log F n (x) + κ n j R (2 n+1 -2)d ∇ log F0 (Φ j n (x; {y i })) ν n (x; d{y i }) , (2.11)
where the x-dependent probability measures ν n (x; d{y i }) on the variables {y i } have the following density with respect to the Lebesgue measure:

ν n (x; d{y i }) d{y i } = P n ({y i }) j F0 (Φ j n (x; {y i })) R (2 n+1 -2)d P n ({y ′ i }) j F0 (Φ j n (x; {y ′ i })) d{y ′ i } .
We remark that a naive repetition of the strategy in Lemma 2.5 under the additional assumption that ∇ log F0 ∈ L ∞ (R d ) would immediately retrieve an exponential decay of the second factor in (2.11) when α > 0 (or only a uniform bound if α = 0) since the large sum over j would be bounded by

2 n κ n ∥∇ log F0 ∥ L ∞ (R d )
. However, as mentioned before, there is a strong drawback: the factors ∇ log F0 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 Λ j n of the lineage maps Φ j n , 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 κ 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 ∇ log Fk , 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 ν 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.

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, . . . , n} but rather on the vertices of a specific type of trees, which are called perfect binary trees, see [START_REF] Rosen | Discrete Mathematics and Its Applications[END_REF][START_REF] Williamson | Combinatorics for computer science[END_REF]. First, we introduce the notion of binary trees, which we present using their universal address system, see [START_REF] Rosen | Discrete Mathematics and Its Applications[END_REF]Section 11.3.3]. Namely, a binary tree T consists in a finite subset of words T ⊆ W 2 := ∪ ∞ k=0 {1, 2} k with letters in the alphabet {1, 2} that verifies:

(i) ∅ ∈ T. (ii) If i1 ∈ T or i2 ∈ T for some word i ∈ W 2 , then i ∈ 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 ∈ T, we denote its length (number of letters) or height in the tree by |i|. In particular, |∅| = 0. A binary tree T is said to be perfect if, in addition, the following properties hold:

(iii) Given a word i ∈ W 2 , then i1 ∈ T if, and only if, i2 ∈ T.

Asymptotics of the high-dimensional integral (Lemma 6. 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 ∈ 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 ∈ 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 ∈ T n : |i| = n}. (ii) (Levels) We say that a vertex i ∈ T n is on the m-th level of the tree if |i| = m for some 1 ≤ m ≤ n.
The set of all vertices on the m-th level will be denoted

L n m := {i ∈ T n : |i| = m}.
In particular, note that L n 0 = {∅} 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 \{∅}, T n := T n \ L n .
(iii) (Child) Given a vertex i ∈ T n * , then there exists a unique word j ∈ 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 ∈ 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 ∈ T n
* , we define the mate of i and we denote it m(i) as the only other vertex in 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 ∈ T n , we say that i ≤ 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 ∈ T n , we define the highest common descendant of i and j, and we denote it by i ∧ j, as

i ∧ j := max{l ∈ T n : l ≤ i, l ≤ j},
where the maximum is considered with respect to above tree (lexicographic) order.

Remark 3.1 (Tree-indexed variables). Given n ∈ N, we shall identify R 2(2 n -1)d ≡ (R d ) T n * and R 2 n d ≡ (R d ) L n . Specifically, vectors x n , y n ∈ R 2(2 n -1)d and z n ∈ R 2 n d
will be regarded as indexed families

x n = (x i ) i∈T n * , y n = (y i ) i∈T n * , z n = (z j ) j∈L n , (3.1)
where x i , y i ∈ R d and z j ∈ R d for each i ∈ T n * and j ∈ 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 µ1,σ 2 1 * G µ2,σ 2 2 = G µ1+µ2,σ 2 1 +σ 2 2 , (3.2)
for any couple of means µ 1 , µ 2 ∈ R d and variances σ 2 1 , σ 2 2 > 0. Using the above result, we can obtain the following explicit evaluation of the operator T in (1.2) over the class of Gaussian functions.

Lemma 3.3 (Evaluation on Gaussians). Consider any

µ ∈ R d and σ 2 ∈ R * + , then T [G µ,σ 2 ] = m * G µ * ,σ 2
* , where the parameters m * , µ * and σ 2 * are given by:

m * = e -1 2 α|µ| 2 1+α(1+ σ 2 2 ) (1 + α(1 + σ 2 2 )) d/2 , µ * = µ 1 + α(1 + σ 2 2 ) , σ 2 * = 1 + σ 2 2 1 + α(1 + σ 2 2 ) . (3.3)
Proof. On the one hand, note that by Lemma 3.2

B[G µ,σ 2 ](x) = G • 2 * G µ,σ 2 * G µ,σ 2 (2x) = G µ,1+ σ 2 2 (x), for each x ∈ R d . Therefore, by definition of T in (1.2) we obtain T [G µ,σ 2 ](x) = e -α 2 |x| 2 G µ,1+ σ 2 2 (x) = 1 (2π(1 + σ 2 2 )) d/2 exp - α 2 |x| 2 - 1 2 1 1 + σ 2 2 |x -µ| 2 ,
for each x ∈ 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. Proof. We look for

λ ∈ R * + , µ ∈ R d and σ 2 ∈ R * + such that λG µ,σ 2 = T [G µ,σ 2 ]
. By Lemma 3.3 and bearing in mind parameters m * , µ * and σ 2 * in (3.3) we obtain that (λ, µ, σ 2 ) must solve the equations:

λ = e -1 2 α|µ| 2 1+α(1+ σ 2 2 ) (1 + α(1 + σ 2 2 )) d/2 , µ = µ 1 + α(1 + σ 2 2 ) , σ 2 = 1 + σ 2 2 1 + α(1 + σ 2 2 )
.

Hence, the only solution is given by µ = 0 and λ and σ 2 are determined by (1.12) and (1.13). □ Remark 3.5 (Dependency on selection). The eigenvalue λ α ∈ R * + and the variance σ 2 α ∈ R * + determined by the relations (1.12) and (1.13) for the the unique Gaussian solution (λ α , F α ) to the eigenproblem (1.11) are monotonically decreasing with the selection coefficient α. Namely, the larger α, the smaller σ 2 α and λ α . Indeed, we obtain (see also Figure 4)

α ↗ ∞ =⇒ σ 2 α ↘ 0, λ α ↘ 0 α ↘ 0 =⇒ σ 2 α ↗ 2, λ α ↗ 1.
Therefore, if selection is strong, then F α is very concentrated around the origin and, if selection ceases, then F α is twice as spread as the mixing kernel G in (1.4), a famous result in quantitative genetics, see e.g. [START_REF] Bulmer | The mathematical theory of quantitative genetics[END_REF]. Note that λ α < 1 for any α ∈ R * + and, consequently, the special steady solutions F n = λ n α F α 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 β in (A.4) coming from the ratio between birth and mortality rates, then the above special steady solution only extincts if

β < (1 + α(1 + σ 2 α 2 )).
Remark 3.6 (Eigenproblem with α = 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., α = 0). However, there is no longer uniqueness due to the translation invariance of B. Indeed, we obtain the Gaussian solutions:

λ = λ α=0 = 1, F = F α=0 (• -µ) = G µ,2 ,
for any µ ∈ 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 λ = 1. Then, the eigenproblem (1.11) reduces to a fixed point equation for 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 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 ∈ R * + , µ 0 ∈ R d and σ 2 0 ∈ R + , and set the Gaussian initial datum F 0 = m 0 G µ0,σ 2 0 . Hence, the solution {F n } n∈N to the time-discrete problem (1.1) takes the form

F n = m n G µn,σ 2 n , (3.4)
for n ∈ N, where the parameters m n ∈ R * + , µ n ∈ R d and σ 2 n ∈ R * + are governed by the recursions:

m n = m n-1 e -1 2 α|µ n-1 | 2 1+α 1+ σ 2 n-1 2 (1 + α(1 + σ 2 n-1 2 )) d/2 , µ n = µ n-1 1 + α(1 + σ 2 n-1 2 ) , 1 σ 2 n = α + 1 1 + σ 2 n-1 2 . ( 3
= α + 1 1 + σ 2 n-1 2 , (3.6
)

for every n ∈ N. Hence, if σ 2 0 > σ 2 α then σ n ↘ σ 2 α as n → ∞ and, if σ 2 0 < σ 2 α then σ n ↗ σ 2 α as n → ∞, where σ 2
α is given by (1.13). In addition, we obtain the convergence rates

|σ 2 n -σ 2 α | ≤ C v r n α , (3.7)
for every n ∈ N, where the constant C v ∈ R + depends on σ 2 0 and α (we obtain that

C v = 0 if σ 2 0 = σ 2 α
), and the ratio r α ∈ (0, 1 2 ) satisfies r α = 2 k 2 α and is given explicitly by

r α := 8 (2α + 3) + (2α + 1) 2 + 8α 2 . (3.8)
Proof.

• Step 1: Monotone convergence of {σ 2 n } n∈N . Consider the function f : R + -→ R given by

f (x) = α + 1 1 + 1 2x , x ∈ R + ,
and, for simplicity of notation, define

x * := (σ 2 α ) -1 , x n := (σ 2 n ) -1 , n ∈ N. (3.9)
Notice that x * is the unique fixed point of f in R + and {x n } n∈N determines the fixed-point iteration of the map f issued at

x 0 = (σ 2 0 ) -1 , that is, x n = f (x n-1 ), n ∈ N.
Our goal is to show that {x n } n∈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 ∈ R + → f (x) is strictly increasing and x ∈ R * + → f (x)
x is strictly decreasing. Then, we obtain that

x < f (x) < x * , if 0 ≤ x < x * , f (x) = x * , if x = x * , x * < f (x) < x, if x > x * ,
which implies the aforementioned monotonicity properties of {x n } n∈N and, in addition,

min{x * , x 0 } ≤ x n ≤ max{x * , x 0 }, (3.10)
for any n ∈ N. In particular, this yields the monotonic convergence of {x n } n∈N towards x * , or equivalently, the claimed monotonic convergence of {σ 2 n } n∈N towards σ 2 α as n → ∞. • 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 R + of the Möbius function M : R \ {-1 2 } -→ R given by Note that M has two fixed points:

M (x) := 2(α + 1)x + α 2x + 1 , x ∈ R \ - 1 2 .
x ± := 2α + 1 ± (2α + 1) 2 + 8α 4 ,
where we note that x + = x * ∈ R * + and x -∈ R -. Then, by definition of {x n } n∈N in (3.9) we obtain

x n -x + = 2(α + 1)x n-1 + α 2x n-1 + 1 - 2(α + 1)x + + α 2x + + 1 = 2 (2x n-1 + 1)(2x + + 1) (x n-1 -x + ). (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 ≤ 2 (2x n-1 + 1)(2x + + 1) ≤ 2 σ 2 α 2 + σ 2 α = 2 k α ,
which gives true contraction in (3.11) when α > 0, since 2 k α < 1. Nevertheless, such a rate is nonoptimal, and we show an alternative approach to achieve a sharper one. Specifically, note that

x n -x -= 2(α + 1)x n-1 + α 2x n-1 + 1 - 2(α + 1)x -+ α 2x -+ 1 = 2 (2x n-1 + 1)(2x -+ 1) (x n-1 -x -). (3.12)
By dividing (3.11) by (3.12) and iterating the identity we obtain

x n -x + x n -x - = 2x -+ 1 2x + + 1 n x 0 -x + x 0 -x - , (3.13)
for every n ∈ N. In fact, the basis can be restated as follows

2x -+ 1 2x + + 1 = (2α + 3) -(2α + 1) 2 + 8α (2α + 3) + (2α + 1) 2 + 8α = 8 (2α + 3) + (2α + 1) 2 + 8α 2 = r α < 1 2 .
Therefore, using the uniform control (3.10) and (3.13), we conclude that

|x n -x * | = |x n -x -| |x 0 -x -| |x 0 -x * | r n α ≤ max{x * , x 0 } -x - min{x * , x 0 } -x - |x 0 -x * | r n α , (3.14)
for every n ∈ N, thus leading to an improved rate since

r α = 2 k 2 α < 2 k α . Finally, notice that |σ 2 n -σ 2 α | = σ 2 n σ 2 α |x n -x * | ≤ max{σ 2 α , σ 2 0 } 2 |x n -x * | , (3.15)
for every n ∈ N. Hence, joining (3.14), (3.15) along with (3.10) ends the proof. □

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 α ∈ R * + , consider any m 0 ∈ R + , µ 0 ∈ R d and σ 2 0 ∈ R + , and set the Gaussian initial datum F 0 = m 0 G µ0,σ0 . Hence, the Gaussian solution {F n } n∈N to the time-discrete problem (1.1) verifies that the growth rates

∥F n ∥ L 1 (R d ) /∥F n-1 ∥ L 1 (R d ) relax towards λ α and the normalized profiles F n /∥F n ∥ L 1 (R d ) relax towards F α ,
where (λ α , F α ) is the unique Gaussian solution (1.12) to the eigenproblem (1.11) as proved in Proposition 3.4. Specifically, for any ε

∈ R * + D KL F n ∥F n ∥ L 1 (R d ) F α ≤ C µ,ε (λ 4/d α +ε) n + C v r n α , ∥F n ∥ L 1 (R d ) ∥F n-1 ∥ L 1 (R d ) -λ α ≤ C µ,ε (λ 4/d α +ε) n + C v r n α ,
where D KL is the Kullback-Leibler divergence (or relative entropy), r α is given by (3.8) in Lemma 3.8, C µ,ε ∈ R + depends on α, µ 0 , σ 2 0 and ε, and

C v ∈ R + depends only α and σ 2 0 . In fact, C µ,ε = 0 if µ 0 = 0 and C v = 0 if σ 2 0 = σ 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:

D KL (G µ1,Σ1 ∥G µ2,Σ2 ) = 1 2 (µ 2 -µ 1 ) ⊤ Σ -1 2 (µ 2 -µ 1 ) + 1 2 trace(Σ 1 Σ -1 2 ) - d 2 + 1 2 log det(Σ 2 Σ -1 1 ), (3.16)
for every µ 1 , µ 2 ∈ R d and any positive non-singular matrices Σ 1 , Σ 2 ∈ R d×d .

Proof of Corollary 3.9. First, note that when α ∈ R * + , then the parameters m n , µ n and σ 2 n in (3.5) verify

lim n→∞ m n m n-1 = λ α , lim n→∞ µ n = 0, lim n→∞ σ 2 n = σ 2 α .
In the sequel, we quantify the rates of convergence. Since σ 2 n has already been studied in Lemma 3.8 we shall focus only on µ n and m n . Let us fix any arbitrarily small ε ∈ R * + . On the one hand note that

|µ n | |µ n-1 | -λ 2/d α = 1 1 + α 1 + σ 2 n-1 2 - 1 1 + α 1 + σ 2 α 2 ,
for every n ∈ N. By the mean value theorem and Lemma 3.8 we obtain

|µ n | |µ n-1 | -λ 2/d α ≤ C v r n α ,
for an appropriate C v ∈ R * + depending on α and σ 2 0 . Since λ α ∈ (0, 1) then d'Alembert's ratio test readily shows that µ n relaxes to zero as a geometric sequence. Indeed, note that

|µ n | = |µ n | |µ n-1 | -λ 2/d α |µ n-1 | + λ 2/d α |µ n-1 | ≤ (λ 2/d α +C v r n α )|µ n-1 |,
for any n ∈ N. By an inductive argument, this yields

|µ n | ≤ n k=1 (λ 2/d α +C v r n α )|µ 0 | ≲ C ε (λ 2/d α +ε) n , (3.17)
for sufficiently large C µ,ε ∈ R d , where we have used that r α ∈ (0, 1) to absorb C v r n α in an ε-small term. On the other hand, note that

m n m n-1 -λ α = e -1 2 α|µ n-1 | 2 1+α 1+ σ 2 n-1 2 (1 + α(1 + σ 2 n-1 2 )) d/2 - 1 (1 + α(1 + σ 2 α 2 )) d/2 = e -1 2 α|µ n-1 | 2 1+α 1+ σ 2 n-1 2 -1 (1 + α(1 + σ 2 n-1 2 )) d/2 +   1 (1 + α(1 + σ 2 n-1 2 )) d/2 - 1 (1 + α(1 + σ 2 α 2 )) d/2   .
By the mean value theorem,

m n m n-1 -λ α ≤ α|µ n-1 | 2 2(1 + α) d 2 +1
+ dα 4

|σ 2 n-1 -σ 2 α | (1 + α) 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 (R d ) and F α are Gaussian functions, the relative entropy can be computed explicitly through (3.16). Specifically, Hence, using again (3.17) and (3.7) concludes our result. □

D KL F n ∥F n ∥ L 1 (R d ) F α = |µ n | 2 2 σ 2 α + d 2 σ 2 n σ 2 α -1 - d 2 log σ 2 n σ 2 α . 5
The proof of Theorem 1.1 for generic non-Gaussian initial data F 0 ∈ M + (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 λ 4/d α = (2 k α ) 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

∥F n ∥ L 1 (R d ) ∥F n-1 ∥ L 1 (R d ) -λ α ≲ D KL F n ∥F n ∥ L 1 (R d ) F α ,
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.

Some properties of the operator T

First, note that for any

F ∈ M + (R d ) we obtain that T [F ] ∈ W k,1 (R d ) ∩ W k,∞ (R d )
for each k ∈ 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∈N to the time-discrete problem (1.1) become instantaneously smooth after the first generation for any generic initial datum

F 0 ∈ M + (R d ).
In the sequel, we will derive a quantitative control on the emergence of Gaussian tails for T [F ]. In addition, we will quantify the propagation of quadratic and exponential moments for the normalized profile

T [F ]/∥T [F ]∥ L 1 (R d )
. Both a priori estimates will be required later in Section 6.

σ 2 := 1 α , σ 2 ∈ 0, 1 1 + α .
Then, the following properties are fulfilled: (i) (Upper control of tails) There exists C = C(α, F ) > 0 such that

T [F ](x) ≤ C G 0,σ 2 (x), (4.1)
for every x ∈ R d .

(ii) (Lower control of tails) There exists c = c(α, σ 2 , F ) > 0 such that

T [F ](x) ≥ c G 0,σ 2 (x), (4.2)
for every x ∈ 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

σ 2 ∈ R * + and C ′ ∈ R * + such that F (x) ≤ C ′ G 0,σ 2 for every x ∈ R d . Then, there exists C ′′ = C ′′ (α, σ 2 , σ 2 , C ′ , F ) > 0 such that |∇ log T [F ](x)| ≤ C ′′ (1 + |x|), (4.3)
for every x ∈ R d .

Proof. First, notice that

F ∈ M + (R d ) is any generic distribution but G ∈ L ∞ (R d ). Then, B[F ] ∈ L ∞ (R d )
and, by definition of T in (1.2), we obtain

T [F ](x) ≤ ∥B[F ]∥ L ∞ (R d ) e -α 2 |x| 2 = ∥B[F ]∥ L ∞ (R d ) 2πα -1 d/2 G 0,σ 2 (x),
for every x ∈ R d . Hence, (4.1) holds for appropriate C > 0. Second, note that

T [F ](x) = e -α 2 |x| 2 (2π) d/2 ∥F ∥ M+(R d ) R 2d exp - 1 2 x - x 1 + x 2 2 2 F (dx 1 )F (dx 2 ) = e -α 2 |x| 2 (2π) d/2 ∥F ∥ M+(R d ) R 2d exp - 1 2 |x| 2 - 1 8 |x 1 + x 2 | 2 + 1 2 x • (x 1 + x 2 ) F (dx 1 )F (dx 2 ) ≥ e -α 2 |x| 2 (2π) d/2 ∥F ∥ M+(R d ) R 2d exp - 2 + ε 4 |x| 2 - 2 + ε 8ε |x 1 + x 2 | 2 F (dx 1 )F (dx 2 ) ≥ e -α 2 |x| 2 (2π) d/2 ∥F ∥ M+(R d ) R 2d exp - 2 + ε 4 |x| 2 - 2 + ε 4ε |x 1 | 2 - 2 + ε 4ε |x 2 | 2 F (dx 1 )F (dx 2 ),
for every ε > 0, where in the third and last lines we have used Cauchy-Schwarz's and Young's inequalities. Consequently, we obtain that

T [F ](x) ≥ (2π) d/2 2ε 2 + ε d ∥G 0, 2ε 2+ε F ∥ 2 M+(R d ) ∥F ∥ M+(R d ) e -α 2 |x| 2 e -2+ε 4 |x| 2 ,
for every x ∈ R d . Taking ε > 0 small enough and an appropriate constant c > 0 yields (4.2). Finally, taking derivatives on (1.2) we obtain

∇T [F ](x) = -α x T [F ](x) - e -α 2 |x| 2 ∥F ∥ L 1 (R d ) R 2d x - x 1 + x 2 2 G x - x 1 + x 2 2 F (dx 1 ) F (dx 2 ) = -(1 + α) x T [F ](x) - e -α 2 |x| 2 ∥F ∥ L 1 (R d ) R 2d x 1 + x 2 2 G x - x 1 + x 2 2 F (dx 1 ) F (dx 2 ).
Therefore, dividing by T [F ](x) and taking norms yield the following estimate

|∇ log T [F ](x)| ≤ (1 + α)|x| + 1 √ 2∥F ∥ L 1 (R d ) e -α 2 |x| 2 T [F ](x) R 2d |(x 1 , x 2 )| G x - x 1 + x 2 2 F (dx 1 )F (dx 2 ), (4.4)
for every x ∈ 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 ) ∈ R 2d : |(x 1 , x 2 )| ≤ √ 2R|x|} and {(x 1 , x 2 ) ∈ R 2d : |(x 1 , x 2 )| > √ 2R|x|}.
For the first subset, we use the definition (1.2) of the operator 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 T [F ]. Then, we obtain

|∇ log T [F ](x)| ≤ (1 + α + R)|x| + C 2 σ d e -α 2 |x| 2 e 1 2σ 2 |x| 2 √ 2c∥F ∥ L 1 (R d ) |(x1,x2)|> √ 2R|x| |(x 1 , x 2 )|e -1 2σ 2 |(x1,x2)| 2 dx 1 dx 2 = (1 + α + R)|x| + C 2 σ d e -α 2 |x| 2 e 1 2σ 2 |x| 2 √ 2c∥F ∥ L 1 (R d ) +∞ √ 2R|x| r 2d e -1 2σ 2 r 2 dr ≤ (1 + α + R)|x| + C 2 σ d e -α 2 |x| 2 e 1 2σ 2 |x| 2 √ 2c∥F ∥ L 1 (R d ) 8d -4 e d-1 2 σ 2d-1 +∞ √ 2R|x| re -1 4σ 2 r 2 dr = (1 + α + R)|x| + √ 2C 2 σ d σ 2d+1 c∥F ∥ L 1 (R d ) 8d -4 e d-1 2 e -1 2 α+ R 2 σ 2 -1 σ 2 |x| 2 ,
for every x ∈ R d , where in the third line we have used the inequality

r 2d-1 ≤ 2d -1 e d-1 2 σ ε 2d-1 exp ε 2σ 2 r 2 ,
for every r > 0 and ε > 0 in the particular case of ε = 1/2. Taking R > 0 large enough ends the proof. □

The previous result can be iterated to obtain similar properties for solutions {F n } n∈N to the timediscrete problem (1.1) issued at a generic initial datum F 0 ∈ M + (R d ). In particular, we note that Gaussian tails emerge instantaneously after the first generation n = 1 and linear growth of the logderivative is guaranteed after the second generation n = 2.

Corollary 4.2 (Propagation of tails). Assume that

α ∈ R * + , consider the solution {F n } n∈N of (1.1) issued at a generic initial datum F 0 ∈ M + (R d ) and set σ 2 1 , σ 2 1 ∈ R * + by σ 2 1 := 1 α , σ 2 1 ∈ 0, 1 1 + α .
Define the associated sequences of variances {σ 2 n } n∈N and {σ 2 n } n∈N by the following recursive relations

1 σ 2 n+1 := α + 1 1 + σ 2 n 2 , 1 σ 2 n+1 := α + 1 1 + σ 2 n 2 , (4.5)
for every n ∈ N. Then, the following properties are fulfilled: (i) (Upper control of tails of F n ) There exist

C n = C n (α, F 0 ) > 0 such that F n (x) ≤ C n G 0,σ 2 n (x), (4.6)
for each x ∈ R d and every n ∈ N.

(ii) (Lower control of tails of F n ) There exist c n = c n (α, σ 2 1 , F 0 ) > 0 such that F n (x) ≥ c n G 0,σ 2 n (x), (4.7)
for each x ∈ R d and every n ∈ N.

(iii) (Control of log-derivative of F n ) There exists D n = D n (α, σ 2 1 , F 0 ) such that |∇ log F n (x)| ≤ D n (1 + |x|), (4.8)
for each x ∈ R d and every n ≥ 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 ≥ 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 = 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 ∈ N and let us prove it for n + 1. Specifically, using the induction hypothesis we obtain that

F n+1 (x) = e -α 2 |x| 2 B[F n ](x) ≤ C 2 n ∥F n ∥ L 1 (R d ) e -α 2 |x| 2 G • 2 * G 0,σ 2 n * G 0,σ 2 n (2x) = C 2 n ∥F n ∥ L 1 (R d ) e -α 2 |x| 2 G 0,1+ σ 2 n 2 (x) ≤ C n+1 G 0,σ 2 n+1 (x),
for every x ∈ 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 σ 2 n+1 . This proves (4.6) for n + 1 and a similar argument yields the lower estimate (4.7). □ By Lemma 3.8 we note that the sequence of variances {σ 2 n } n∈N and {σ 2 n } n∈N relax towards the asymptotic value σ 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∈N to the time-discrete problem (1.1) must relax towards the asymptotic profile F α . In addition, recall that for any solution (λ, 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 = λ n F . Hence, we expect that (λ α , F α ) 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∈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 -ασ 2 n ) d/2 C 2 n ∥F n ∥ L 1 (R d ) , c n+1 = (1 -ασ 2 n ) d/2 c 2 n ∥F n ∥ L 1 (R d ) ,
but they require further knowledge about the behavior of

∥F n ∥ L 1 (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 (R d )
, where {F n } n∈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 ∈ M + (R d ), we illustrate the particular explicit case of Gaussian initial data

F 0 = m 0 G µ0,σ 2 0 . Remark 4.3 (Gaussian case). Consider any m 0 ∈ R + , µ 0 ∈ R d and σ 2 0 ∈ R + , set θ ∈ R * + with θ < 1 2 max{σ 2 0 ,σ 2
α } and consider the solution {F n } n∈N to the time-discrete problem (1.1) with initial datum F 0 = m 0 G µ0,σ 2 0 . Then, by direct calculation we obtain

R d |x| 2 F n (x) ∥F n ∥ L 1 (R d ) dx = |µ n | 2 + σ 2 n , R d e θ|x| 2 F n (x) ∥F n ∥ L 1 (R d ) dx = e θ 1-2θσ 2 n |µn| 2 (1 -2θσ 2 n ) d/2
, Here, {µ n } n∈N and {σ 2 n } n∈N are the mean and variance of {F n } n∈N according to formula (3.5) in Proposition 3.7. As studied in the proof of Corollary 3.9 we have the uniform bounds

|µ n | ≤ |µ 0 | and min{σ 2 0 , σ 2 α } ≤ σ 2 n ≤ max{σ 2 0 , σ 2
α }, for every n ∈ 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 ∈ M + (R). For simplicity of notation, we define the following conservative operator associated with T .

Definition 4.4 (Normalization of T ).

S[F

] := T [F ] ∥T [F ]∥ L 1 (R d ) , F ∈ M + (R d ) \ {0} .
Therefore, we obtain the following control of quadratic moments under S.

= M (α, η) ∈ R * + such that R d |x| 2 S[F ](x) dx ≤ M + η R d |x| 2 F (dx) ∥F ∥ M(R d ) ,
for any measure

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

R d |x| 2 M[F ](x) dx ≤ R d |x| 2 F (dx) . (4.9)
This inequality is intuitively clear: by applying the selection operator 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:

R d |x| 2 M[F ](x) dx - R d |x| 2 F (dx) = 1 2 R 2d |x| 2 -|y| 2 e -m(x) -e -m(y) F (dx)F (dy) 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 nondecreasing as it is the case for our quadratic choice (1.5). A similar estimate can also be tested on the operator B. Specifically, by a simple change of variable the following identity holds true:

R d |x| 2 B[F ](x) dx = Var G + 1 2 R d |x| 2 F (dx) + 1 2 R d x F (dx) 2 . (4.10)
Then, it might seem natural to combine these two relationships in order to control the quadratic moments of the composition S = M • 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 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 M, as the reproduction component B is invariant by translation. However, it cannot be derived solely from 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 B, so that the combination of M and 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 ∈ M + (R d ). First, note that R d |x| 2 S[F ](x) dx = R 3d |x| 2 e -m(x) G x - x 1 + x 2 2 F (dx 1 ) ∥F ∥ M(R d ) F (dx 2 ) ∥F ∥ M(R d ) dx R 3d e -m(x) G x - x 1 + x 2 2 F (dx 1 ) ∥F ∥ M(R d ) F (dx 2 ) ∥F ∥ M(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 ∈ R d we define x := x1+x2 2 and we obtain

R d |x| 2 e -m(x) G(x -x) dx = 1 (1 + α) d/2 exp - 1 2 α 1 + α |x| 2 1 1 + α + |x| 2 (1 + α) 2 , (4.11) R d e -m(x) G(x -x) dx = 1 (1 + α) d/2 exp - 1 2 α 1 + α |x| 2 . (4.12)
Therefore, we can write alternatively

R d |x| 2 S[F ](x) dx = 1 1 + α + 1 (1 + α) 2 R 2d |x| 2 exp - 1 2 α 1 + α |x| 2 F (dx 1 ) ∥F ∥ M(R d ) F (dx 2 ) ∥F ∥ M(R d ) R 2d exp - 1 2 α 1 + α |x| 2 F (dx 1 ) ∥F ∥ M(R d ) F (dx 2 ) ∥F ∥ M(R d )
.

(4.13)

Our ultimate goal is to bound the second term in (4.13). To this end, we provide a lower bound for the denominator and an upper bound for the numerator in two separate steps.

• Step 1: Consider an arbitrary R 1 ∈ R * + to be determined later. Then, we obtain

R 2d exp - 1 2 α 1 + α |x| 2 F (dx 1 ) ∥F ∥ M(R d ) F (dx 2 ) ∥F ∥ M(R d ) ≥ |x|≤R1 exp - 1 2 α 1 + α |x| 2 F (x 1 ) ∥F ∥ L 1 (R d ) F (x 2 ) ∥F ∥ L 1 (R d ) dx 1 dx 2 ≥ exp - 1 2 α 1 + α R 2 1 1 - |x|>R1 F (dx 1 ) ∥F ∥ M(R d ) F (dx 2 ) ∥F ∥ M(R d ) ≥ exp - 1 2 α 1 + α R 2 1 1 - 1 2R 2 1 R d |x| 2 F (dx) ∥F ∥ M(R d ) , (4.14)
where in last line we have used Young's and Chebyshev's inequalities.

•

Step 2: Consider an arbitrary R 2 ∈ R * + to be determined later. Then, we obtain

R 2d |x| 2 exp - 1 2 α 1 + α |x| 2 F (dx 1 ) ∥F ∥ M(R d ) F (dx 2 ) ∥F ∥ M(R d ) = I 1 + I 2 ,
where each term reads

I 1 := |x|≤R2 |x| 2 exp - 1 2 α 1 + α |x| 2 F (dx 1 ) ∥F ∥ M(R d ) F (dx 2 ) ∥F ∥ M(R d ) , I 2 := |x|>R2 |x| 2 exp - 1 2 α 1 + α |x| 2 F (dx 1 ) ∥F ∥ M(R d ) F (dx 2 ) ∥F ∥ M(R d )
.

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

I 1 ≤ R 2 2 R 2d exp - 1 2 α 1 + α |x| 2 F (dx 1 ) ∥F ∥ M(R d ) F (dx 2 ) ∥F ∥ M(R d ) . (4.15)
On the other hand, the second term can be estimated by

I 2 ≤ 2 eγ 1 + α α |x|>R2 exp - 1 -γ 2 α 1 + α |x| 2 F (dx 1 ) ∥F ∥ M(R d ) F (dx 2 ) ∥F ∥ M(R d ) ≤ 2 eγ 1 + α α exp - 1 -γ 2 α 1 + α R 2 2 , (4.16)
for any γ ∈ (0, 1), where in the first line we have used that e θs ≥ θes with s = |x| 2 and θ = γ 2 α 1+α , and in the second line we have exploited that F/∥F ∥ M(R d ) are probability measures, thus normalized.

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

R d |x| 2 S[F ](x) dx ≤ 1 1 + α + R 2 2 (1 + α) 2 + 2 eγ 1 α(1 + α) exp 1 2 α 1+α R 2 1 exp -1-γ 2 α 1+α R 2 2 1 -1 2R 2 1 R d |x| 2 F (dx) ∥F ∥ M(R d ) ,
for any R 1 , R 2 ∈ R * + and γ ∈ (0, 1). We set the values

R 2 1 = 1 2(1 -δ) R d |x| 2 F (dx) ∥F ∥ M(R d ) , R 2 2 = R 2 1 1 -γ ,
with γ, δ ∈ (0, 1). Then, we obtain

R d |x| 2 S[F ](x) dx ≤ 1 1 + α + 2 eγδ 1 α(1 + α) + 1 2(1 -δ)(1 -γ)(1 + α) 2 R d |x| 2 F (dx) ∥F ∥ M(R d ) ,
thus ending the proof by the arbitrariness of γ, δ ∈ (0, 1). □

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

Corollary 4.6 (Propagation of quadratic moments). Assume that α ∈ R * + , set any parameter η ∈ 1 2(1+α) 2 , 1 and consider the solution {F n } n∈N of (1.1) issued at a generic initial datum

F 0 ∈ M + (R d ) with R d |x| 2 F 0 (dx) < ∞. Then, there exists M = M (α, η) such that R d |x| 2 F n (x) ∥F n ∥ L 1 (R d ) dx ≤ M 1 -η + η n R d |x| 2 F 0 (dx) ∥F 0 ∥ M(R d ) ,
for any n ∈ N.

Let us note that, on the one hand the factor η n makes the dependence on the variance of the initial datum F 0 negligible as the amount of generations increases. On the other hand, the constant M does not depend on the chosen initial profile F 0 . This is compatible with our ergodicity property in Theorem 1.1. Lemma 4.7 (Control of exponential moments). Assume that α ∈ R * + , set θ ∈ R * + so that θ < α 2 and any parameter χ > α 1+α θ 2(1+α-2θ)(α-2θ) . Then, there exists a constant C = C(α, θ, χ) such that

R d e θ|x| 2 S[F ](x) dx ≤ C 1 + exp χ R d |x| 2 F (dx) ∥F ∥ M(R d ) ,
for any measure

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

R d e θ|x| 2 S[F ](x) dx = 1 + α 1 + α -2θ d/2 R 2d exp - 1 2 α -2θ 1 + α -2θ |x| 2 F (dx 1 ) ∥F ∥ M(R d ) F (dx 2 ) ∥F ∥ M(R d ) R 2d exp - 1 2 α 1 + α |x| 2 F (dx 1 ) ∥F ∥ M(R d ) F (dx 2 ) ∥F ∥ M(R d )
, (4.17) where x = x1+x2 2 . Again, we provide bounds for the numerator and denominator in (4.17). Regarding the denominator, we recover estimate (4.14) for any R 1 ∈ R * + . Thus, we focus on the bound of the numerator

R 2d exp - 1 2 α -2θ 1 + α -2θ |x| 2 F (dx 1 ) ∥F ∥ M(R d ) F (dx 2 ) ∥F ∥ M(R d ) = I 1 + I 2 ,
where each term reads

I 1 := |x|≤R2 exp - 1 2 α -2θ 1 + α -2θ |x| 2 F (dx 1 ) ∥F ∥ M(R d ) F (dx 2 ) ∥F ∥ M(R d ) , I 2 := |x|>R2 exp - 1 2 α -2θ 1 + α -2θ |x| 2 F (dx 1 ) ∥F ∥ M(R d ) F (dx 2 ) ∥F ∥ M(R d )
.

On the one hand, note that

I 1 ≤ exp θR 2 2 (1 + α -2θ) 2 R 2d exp - 1 2 α 1 + α |x| 2 F (dx 1 ) ∥F ∥ M(R d ) F (dx 2 ) ∥F ∥ M(R d ) , (4.18)
where we have used the mean value theorem applied to the function r ∈ R * + → r 1+r to derive the relation

α 1 + α - α -2θ 1 + α -2θ ≤ 2θ (1 + α -2θ) 2 .
On the other hand,

I 2 ≤ exp - 1 2 α -2θ 1 + α -2θ R 2 2 . ( 4
R d e θ|x| 2 S[F ](x) dx ≤ 1 + α 1 + α -2θ d/2    exp θR 2 2 (1 + α -2θ) 2 + exp -1 2 α-2θ 1+α-2θ R 2 2 exp 1 2 α 1+α R 2 1 1 -1 2R 2 1 R d |x| 2 F (dx) ∥F ∥ M(R d )    , for any R 1 , R 2 ∈ R * + .
As in the proof of Lemma 4.5 we choose them as follows

R 2 1 = 1 2(1 -δ) R d |x| 2 F (dx) ∥F ∥ M(R d ) , R 2 2 = α 1 + α 1 + α -2θ α -2θ R 2 1 ,
where δ ∈ (0, 1). Then, we obtain

R d e θ|x| 2 S[F ](x) dx ≤ 1 + α 1 + α -2θ d/2 1 δ + exp α 1 + α θ (1 + α -2θ)(α -2θ) 1 2(1 -δ) R d |x| 2 F (dx) ∥F ∥ M(R d ) ,
and we end the proof by arbitrariness of δ ∈ (0, 1). □ Therefore, small exponential moments can be controlled by quadratic moments. We can then couple Lemma 4.7 and Corollary 4.6 to obtain the propagation of exponential moments in the time-discrete problem (1.1).

Corollary 4.8 (Propagation of exponential moments). Assume that α ∈ R *

+ , set θ ∈ R * + with θ < α 2 and any parameter η ∈ 1 2(1+α) 2 , 1 . Consider the solution {F n } n∈N of (1.1) issued at a generic initial datum

F 0 ∈ M + (R d ) with R d |x| 2 F 0 (dx) < ∞.
Then, there exist C = C(α, θ, η) and

C ′ = C ′ (α, θ, η) such that R d e θ|x| 2 F n (x) ∥F n ∥ L 1 (R d ) dx ≤ C 1 + exp C ′ 1 + η n-1 R d |x| 2 F 0 (dx) ∥F 0 ∥ M(R d ) ,
for any n ∈ N.

Reformulating the recursion

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

⋄ Iteration 1:
By the explicit definition of T in (1.2) we readily obtain

F 1 (x) = 1 ∥F 0 ∥ M(R d ) R 2d e -m(x) G x - x 1 + x 2 2 F 0 (x 1 )F 0 (x 2 ) dx 1 , (5.1)
where we denote

x 1 = (x 1 , x 2 ) ∈ R 2d .
⋄ Iteration 2: By (5.1) with F 0 replaced by F 1 and using F 1 = T [F 0 ] in the right hand side we get

F 2 (x) = 1 ∥F 1 ∥ L 1 (R d ) ∥F 0 ∥ 2 M(R d ) R 6d e -(m(x)+m(x1)+m(x2)) G x - x 1 + x 2 2 × G x 1 - x 11 + x 12 2 G x 2 - x 21 + x 22 2 F 0 (x 11 )F 0 (x 12 )F 0 (x 21 )F 0 (x 22 ) dx 2 , (5.2)
where we denote x 2 = (x 1 , x 2 , x 11 , x 12 , x 21 , x 22 ) ∈ R 6d .

⋄ Iteration n ∈ N: By a clear inductive process, we can iterate the operator T as many times n ∈ N as needed to recover an explicit dependence of F n on the initial datum F 0 . Of course, this generates a high-dimensional integral involving 2(2 n -1) variables that we can label along the vertices of the perfect binary tree T n according to the universal address system notation in Section 3.1. Namely, we obtain

F n (x) = 1 n-1 m=0 ∥F m ∥ 2 n-1-m L 1 (R d ) × R 2(2 n -1)d exp   - i∈ T n m(x i )   i∈ T n G x i - x i1 + x i2 2 i∈L n F 0 (x i ) dx n , (5.3)
where we have used the tree-indexed notation x n = (x i ) i∈T n * ∈ R 2(2 n -1)d in Remark 3.1 and x ∅ := x. The ultimate goal of this section is to find a suitable change of variables as discussed in Section 2.4. We shall see that the effect of such a change of variable on the above high-dimensional integral in (5.3) will reveal that actual dependence of {F n } n∈N on the shape of F 0 . We will see that such a dependence is indeed weak and the initial datum is rapidly "forgotten" across the different levels of the tree T n . This ergodicity property will be crucial in our analysis and will be exploited later in Section 6 to derive our main result in Theorem 1.1.

5.1.

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

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

F (dx) := e m(x) F (dx) F α=0 (x) , x ∈ R d .
Here, we remind that F α=0 = G 0,2 is the Gaussian eigenfunction corresponding to (1.12) with α = 0. We also recall that there is some freedom in the normalization, and in particular the term e m is not mandatory but it is convenient for an easier sorting of the various terms, as we anticipated in Section 2.4. By substituting (1.4) and (1.5) into (5.1) and writting the integral in the right hand side in terms of the rescaled function F0 in Definition 5.1, we obtain that

(e m F 2 )(x) = 1 (2π) 3d/2 1 (4π) 2d 1 ∥F 1 ∥ L 1 (R d ) ∥F 0 ∥ 2 M(R d ) × R 6d exp - α 2 |x 1 | 2 + |x 2 | 2 - 1 2 x - x 1 + x 2 2 2 × exp - α 0 2 |x 11 | 2 + |x 12 | 2 + |x 21 | 2 + |x 22 | 2 - 1 2 x 1 - x 11 + x 12 2 2 - 1 2 x 2 - x 21 + x 22 2 2 × F0 (x 11 ) F0 (x 12 ) F0 (x 21 ) F0 (x 22 ) dx 2 .
(5.4)

Here, the parameter α 0 ∈ R * + has been defined by

α 0 := 1 2 + α,
so that the quadratic terms in the fist term of the third line of (5.4) correct the rescaling F0 of F 0 . The new formulation will be obtained by an appropriate change of variables, so that the quadratic forms inside the above exponential get appropriately centered. To do so, notice that this formula involves as many variables x i as indices i in the perfect binary tree T 2 . Indeed, we have sorted the different terms in (5.4) in such a way that the second line only involves variables x 1 and x 2 in the first level L 2 1 of the tree, whilst the third line contains the terms involving the variables x 11 , x 12 , x 21 and x 22 in the second level of the tree L 2 2 (i.e., the leaves L 2 ). We will divide the method into two steps. First, we will address the change of the variables indexed by the leaves L 2 . Second, we will perform the change of variables indexed by the first level of the tree L 2 1 .

• Step 1: Change of variables for x 11 , x 12 , x 21 and x 22 . Consider some coefficient k 1 > 0 to be determined later and define the change of variables x 11 → y 11 , x 12 → y 12 , x 21 → y 21 and x 22 → y 22 given by

x 11 = k 1 x 1 + y 11 , x 12 = k 1 x 1 + y 12 , x 21 = k 1 x 2 + y 21 , x 22 = k 1 x 2 + y 12 .
On the one hand, using the change of variables in the terms inside the exponential of (5.4) which involve x 11 and x 12 , we obtain

α 0 2 |x 11 | 2 + α 0 2 |x 12 | 2 + 1 2 x 1 - x 11 + x 12 2 2 = α 0 2 |k 1 x 1 + y 11 | 2 + α 0 2 |k 1 x 1 + y 12 | 2 + 1 2 (1 -k 1 )x 1 - y 11 + y 12 2 2 = 1 2 2α 0 k 2 1 + (1 -k 1 ) 2 |x 1 | 2 + α 0 2 |y 11 | 2 + α 0 2 |y 12 | 2 + 1 8 |y 11 + y 12 | 2 + α 0 k 1 - 1 -k 1 2 x 1 • y 11 + α 0 k 1 - 1 -k 1 2 x 1 • y 12 .
On the other hand, using the change of variables in the terms which involve x 21 and x 22 , we get

α 0 2 |x 21 | 2 + α 0 2 |x 22 | 2 + 1 2 x 2 - x 21 + x 22 2 2 = α 0 2 |k 1 x 2 + y 21 | 2 + α 0 2 |k 1 x 2 + y 22 | 2 + 1 2 (1 -k 1 )x 2 - y 21 + y 22 2 2 = 1 2 2α 0 k 2 1 + (1 -k 1 ) 2 |x 2 | 2 + α 0 2 |y 21 | 2 + α 0 2 |y 22 | 2 + 1 8 |y 21 + y 22 | 2 + α 0 k 1 - 1 -k 1 2 x 2 • y 21 + α 0 k 1 - 1 -k 1 2 x 2 • y 22 .
Notice that one can eliminate the crossed terms in both expressions by choosing

k 1 := 1 1 + 2α 0 = 1 2(1 + α) .
In that case, adding both terms we obtain that

α 0 2 |x 11 | 2 + α 0 2 |x 12 | 2 + 1 2 x 1 - x 11 + x 12 2 2 + α 0 2 |x 21 | 2 + α 0 2 |x 22 | 2 + 1 2 x 2 - x 21 + x 22 2 2 = 1 2 (1 -k 1 )|x 1 | 2 + 1 4k 1 |y 11 | 2 + 1 4k 1 |y 12 | 2 - 1 8 |y 11 -y 12 | 2 + 1 2 (1 -k 1 )|x 2 | 2 + 1 4k 1 |y 21 | 2 + 1 4k 1 |y 22 | 2 - 1 8 |y 21 -y 22 | 2 .
Putting everything together into (5.4) yields where the parameter α 1 ∈ R * + has been defined by

(e m F 2 )(x) = 1 (2π) 3d/2 1 (4π) 2d 1 ∥F 1 ∥ L 1 (R d ) ∥F 0 ∥ 2 M(R d ) × R 6d exp - α 1 2 |x 1 | 2 + |x 2 | 2 - 1 2 x - x 1 + x 2 2 2 × exp - 1
α 1 := 1 -k 1 + α,
in order to recombine the initial terms in the second line of (5.4) involving variables indexed by the first level L 2 1 with the new remainders that have appeared from the previous step.

• Step 2: Change of variables for x 1 and x 2 . Consider some coefficient k 2 > 0 to be determined later and define the change of variables x 1 → y 1 and x 2 → y 2 given by

x 1 = k 2 x + y 1 , x 2 = k 2 x + y 2 .
This time, using the change of variables in the terms of the exponential in the second line of (5.5) yields

α 1 2 |x 1 | 2 + α 1 2 |x 2 | 2 + 1 2 x - x 1 + x 2 2 2 = α 1 2 |k 2 x + y 1 | 2 + α 1 2 |k 2 x + y 2 | 2 + 1 2 (1 -k 2 )x - y 1 + y 2 2 2 = 1 2 2α 1 k 2 2 + (1 -k 2 ) 2 |x| 2 + α 1 2 |y 1 | 2 + α 1 2 |y 2 | 2 + 1 8 |y 1 + y 2 | 2 + α 1 k 2 - 1 -k 2 2 x • y 1 + α 1 k 2 - 1 -k 2 2 x • y 2 .
Again, we can cancel the crossed term by choosing

k 2 := 1 1 + 2α 1 = 1 3 + 2α -2k 1 .
Namely, we obtain that

α 1 2 |x 1 | 2 + α 1 2 |x 2 | 2 + 1 2 x - x 1 + x 2 2 2 = 1 2 (1 -k 2 )|x| 2 + 1 4k 2 |y 1 | 2 + 1 4k 2 |y 2 | 2 - 1 8 |y 1 -y 2 | 2 .
Then, putting everything together into (5.5) yields

F 2 (x) = 1 (2π) 3d/2 1 (4π) 2d 1 ∥F 1 ∥ L 1 (R d ) ∥F 0 ∥ 2 M(R d ) e -1+α-k 2 2 |x| 2 × R 6d exp - 1 4k 2 |y 1 | 2 + |y 2 | 2 + 1 8 |y 1 -y 2 | 2 × exp - 1 4k 1 |y 11 | 2 + |y 12 | 2 + |y 21 | 2 + |y 22 | 2 + 1 8 |y 11 -y 12 | 2 + |y 21 -y 22 | 2 × F0 (k 1 k 2 x + k 1 y 1 + y 11 ) F0 (k 1 k 2 x + k 1 y 1 + y 12 ) × F0 (k 1 k 2 x + k 1 y 2 + y 21 ) F0 (k 1 k 2 x + k 1 y 2 + y 22 ) dy 2 , (5.6)
where we denote again y 2 = (y 1 , y 2 , y 11 , y 12 , y 21 , y 22 ) ∈ R 6d .

Before stating the main result for general n ∈ N, we collect some natural notation according to the preceding computations, that will be useful here on. First, we define the following sequences of coefficients. Definition 5.2 (Coefficients). Consider any α ∈ R + .

(1) The coefficients {k n } n∈N are defined by the recursive formula

k 1 := 1 2(1 + α) , k n := 1 3 + 2α -2k n-1
, for n ≥ 2.

(5.7)

(2) The coefficients {κ n } n∈N are defined by the recursive formula

κ 0 := 1, κ n := k 1 • • • k n , for n ≥ 1. (5.8)
We note that for n = 2 the above coefficients k 1 , k 2 reduce to those appearing in the previous reformulation (5.6) of the recursion. Since it will be used later, we study the asymptotic behavior of such sequences of coefficients as n → ∞. (5.9)

Specifically, k α is given explicitly by formula (1.15) and it is related to the coefficient r α in formula (3.8) of Lemma 3.8 by r α = 2 k 2 α . In addition,

k n -k α ≤ C r n-1 α , (5.10)
for every n ∈ N, where r α is and C ∈ R * + depends only on α. (ii) (Coefficients κ n ) The coefficients {κ n } n∈N are positive numbers and decay geometrically with maximal rate k α . Specifically, for every ε ∈ R * + there exists

C ε ∈ R * + such that κ n ≤ C ε (k α +ε) n , (5.11)
for every n ∈ N.

Proof. On the one hand, the properties of {κ n } n∈N readily follow from those of {k n } n∈N and the relation (5.8) in Definition 5.2. Specifically, for any given q ∈ N and any n ≥ q we get the decomposition

κ n = q m=0 k m n m=q+1 k m ≤ k q 1 k n-q q+1 = k q 1 k q q+1 k n q+1 ≤ k q 1 k q q+1 (C r q+1 α + k α ) n .
Therefore, taking q sufficiently large so that C r q+1 α ≤ ε we conclude (5.11). We then focus on the properties of {k n } n∈N and we use a similar strategy like in Lemma 3.8.

• Step 1: Monotone convergence of {k n } n∈N . Define the following function

f (x) := 1 3 + 2α -2x , x ∈ -∞, 3 2 + α .
Then, by (5.7) {k n } n∈N obeys the following recursive relation

k n = f (k n-1 ), n > 1.
Let us consider the fixed points x -< x + of f , i.e,

x ± = (3 + 2α) ± (1 + 2α) 2 + 8α 4 .
By inspection, it is clear that 0

< k α = x -< x + < 3 2 + α and x -< f (x) < x, if x ∈ (x -, x + ),
where k α is given in (1.15). Since k 1 ∈ (x -, x + ), then we conclude that {k n } n∈N is a well defined, positive and monotonically decreasing sequence contained in the compact interval [x -, x + ]. Therefore, it must converge towards some limit ℓ, which is a solution of f (x) = x i.e., ℓ ∈ {x -, x + }. Since the full sequence {k n } n∈N is below x + and decreasing, we then conclude that ℓ = x -= k α . In particular, we obtain

k α < k n ≤ k 1 , (5.12)
for any n ∈ N.

•

Step 2: Convergence rates. As for Lemma 3.8, we shall use the special algebraic structure of the function f , which is the restriction to -∞, 3 2 + α of the Möbius function M : R \ { 3 2 + α} -→ R given by

M (x) := 1 3 + 2α -2x , x ∈ R \ 3 2 + α .
By definition of {k n } n∈N in (5.7) we obtain

k n -x + = 1 3 + 2α -2k n-1 - 1 3 + 2α -2x + = 2 (3 + 2α -2k n-1 )(3 + 2α -2x + ) (k n-1 -x + ).
Similarly, we obtain

k n -x -= 1 3 + 2α -2k n-1 - 1 3 + 2α -2x - = 2 (3 + 2α -2k n-1 )(3 + 2α -2x -) (k n-1 -x -).
By dividing both expressions and iterating such an identity we get for every n ∈ N. In fact, the basis can be restated as follows

k n -x - k n -x + = 3 + 2α -2x + 3 + 2α -2x - n-1 k 1 -x - k 1 -x + , (5.13)
3 + 2α -2x + 3 + 2α -2x - = 3 + 2α -(1 + 2α) 2 + 8α 3 + 2α + (1 + 2α) 2 + 8α = 8 3 + 2α + (1 + 2α) 2 + 8α 2 ≡ r α < 1,
where r α is determined by (3.8) in Lemma 3.8. Therefore, using (5.12) and (5.13) we conclude that

k n -k α = k 1 -x - x + -k 1 (x + -k n ) r n-1 α ≤ k 1 -x - x + -k 1 (x + -x -) r n-1 α ,
for any n ∈ N, thus ending the proof. □ Definition 5.4 (Quadratic forms). Consider the sequence {k n } n∈N in Definition 5.2. Then, we define the quadratic form

Q n = Q n (y n ) by Q n (y n ) := n-1 m=0 i∈L n m 1 4k n-m (|y i1 | 2 + |y i2 | 2 ) - 1 8 |y i1 -y i2 | 2 , y n ∈ R 2(2 n -1)d , (5.14)
where we are using the tree-indexed notation y n = (y i ) i∈T n * ∈ R 2(2 n -1)d in Remark 3.1. Again, note that when n = 2 the above quadratic form Q 2 (y 2 ) reduces to the one inside the exponential of (5.6). We now show the following uniform control of the quadratic forms Q n as n → ∞.

β min := 1 + 2α 4 , β max := 1 4 k α ,
where k α is given by formula (1.15). Then, we obtain that

β min ∥ y n ∥ 2 ≤ Q n (y n ) ≤ β max ∥ y n ∥ 2 ,
where we are using the tree-indexed notation

y n = (y i ) i∈T n * ∈ R 2(2 n -1)d in Remark 3.1. Proof.
By virtue of Young's inequality, we obtain the following lower and upper bound for

Q n n-1 m=0 i∈L n m 1 4 1 k n-m -1 (|y i1 | 2 + |y i2 | 2 ) ≤ Q n (y n ) ≤ n-1 m=0 i∈L n m 1 4 1 k n-m (|y i1 | 2 + |y i2 | 2 ),
for every

y n = (y i ) i∈T n * ∈ R 2(2 n -1)d .
Then, the result follows from Lemma 5.3, which guarantees the uniform control k α < k m-n ≤ k 1 for every m = 0, . . . , n -1. □ Definition 5.6 (Lineage maps). We define the lineage maps Φ j n = Φ j n (x; y n ) associated to the leave j ∈ L n (see Figure 7) as follows

Φ j n (x; y n ) := κ n x + n-1 m=0 κ m y c m (j) , y n ∈ R 2(2 n -1)d , (5.15)
where x ∈ R d represents the root value, We are now ready to state the main result of this section extending formula (5.6) to any n ∈ N. The starting point is again formula (5.3). In the particular case of Gaussian mixing kernel (1.4) and quadratic selection function (1.5), it takes the explicit form

y n = (y i ) i∈T n * ∈ R 2(2 n -1)
(e m F n )(x) = 1 (2π) (2 n -1)d/2 1 (4π) 2 n-1 d 1 n-1 m=0 ∥F m ∥ 2 n-1-m L 1 (R d ) × R 2(2 n -1)d exp   - n-2 m=0 i∈L n m α 2 (|x i1 | 2 + |x i2 | 2 ) + 1 2 x i - x i1 + x i2 2 2   × exp   - i∈L n n-1 1 2 (α + 1 2 )(|x i1 | 2 + |x i2 | 2 ) + 1 2 x i - x i1 + x i2 2 2   × j∈L n F0 (x j ) dx n , (5.16)
for any x ∈ R d and n ∈ N. Above, we have used the tree-indexed notation x n = (x i ) i∈T n * ∈ R 2(2 n -1)d in Remark 3.1, we have set x ∅ := x and we have considered the rescaled distribution F0 according to Definition 5.1. Again, notice that the exponential terms in the fist term of the third line of (5.16) have been introduced to appropriately correct the rescaled distribution F0 of F 0 . As a consequence, the factor (4π) -2 n-1 d arises from the normalization by F α=0 (x j ) whilst the factor (2π) -(2 n -1)d/2 comes from the repeated products of the Gaussian mixing kernel G. The main result then reads as follows Proposition 5.7 (Reformulation of the iterations I). Assume that α ∈ R + and set any initial datum

F 0 ∈ M + (R d ).
Hence, the solution {F n } n∈N to the time-discrete problem (1.1) admits the following form

F n (x) = e -1+α-kn 2 |x| 2 (4π) 2 n-1 d (2π) (2 n -1)d/2 1 n-1 m=0 ∥F m ∥ 2 n-1-m L 1 (R d ) R 2(2 n -1)d e -Qn(y n ) j∈L n F0 (Φ j n (x; y n )) dy n , (5.17)
for every x ∈ R d and n ∈ N, where we are using the tree-indexed notation y n = (y i ) i∈T n * ∈ R 2(2 n -1)d in Remark 3.1 along with further notation from Definitions 5.1, 5.2, 5.4 and 5.6.

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

x i1 = k n-m x i + y i1 , x i2 = k n-m x i + y i2 ,
for every i ∈ L n m , which we will apply in a backwards way starting at indices i ∈ L n n-1 and ending at i = ∅. Again, the objective of such a change of variables is to appropriately move the reference frame so that the above quadratic forms 1 2

x i -xi1+xi2 2
2 inside the exponential of (5.16) get appropriately centered and no crossed terms remain when the selection parts α 2 (|x i1 | 2 + |x i2 | 2 ) are taken into account. For simplicity, along the proof we shall restrict to tracking how the various terms inside the exponentials in (5.16) get modified after the change of variables.

Proof of Proposition 5.7.

• Step 1: Change of variables x i1 → y i1 and x i2 → y i2 for i ∈ L n n-1 . We consider the change of variables given by

x i1 = k 1 x i + y i1 , x i2 = k 1 x i + y i2 ,
where k 1 is given by Definition 5.2. The collection of all the terms for i ∈ L n n-1 in (5.16) then read

α 0 2 |x i1 | 2 + α 0 2 |x i2 | 2 + 1 2 x i - x i1 + x i2 2 2 = α 0 2 |k 1 x i + y i1 | 2 + α 0 2 |k 1 x i + y i2 | 2 + 1 2 (1 -k 1 )x i - y i1 + y i2 2 2 = 1 2 2α 0 k 2 1 + (1 -k 1 ) 2 |x i | 2 + α 0 2 |y i1 | 2 + α 0 2 |y i2 | 2 + 1 8 |y i1 + y i2 | 2 + α 0 k 1 - 1 -k 1 2 x i • y i1 + α 0 k 1 - 1 -k 1 2 x i • y i2 = 1 2 2α 0 k 2 1 + (1 -k 1 ) 2 |x i | 2 + α 0 2 |y i1 | 2 + α 0 2 |y i2 | 2 + 1 8 |y i1 + y i2 | 2 = 1 2 (1 -k 1 )|x i | 2 + 1 4k 1 |y i1 | 2 + 1 4k 1 |y i2 | 2 - 1 8 |y i1 -y i2 | 2 ,
where we have defined the coefficient α 0 ∈ R * + by

α 0 := α + 1 2 ,
to recombine terms in the second and third lines of (5.16), and we have used the following relation between k 1 and α 0 in order to cancel the crossed terms, i.e.,

k 1 = 1 1 + 2α 0 .
• Step 2: Change of variables x i1 → y i1 and x i2 → y i2 for i ∈ L n n-2 . Again, we consider the change of variables given by

x i1 = k 2 x i + y i1 , x i2 = k 2 x i + y i2 ,
where k 2 is given by Definition 5.2. Putting together the terms for i ∈ L n n-2 in (5.16) and the above x i -dependent remainder in the above expression in Step 1 yields the following term under the above change of variables

α 1 2 |x i1 | 2 + α 1 2 |x i2 | 2 + 1 2 x i - x i1 + x i2 2 2 = α 1 2 |k 2 x i + y i1 | 2 + α 1 2 |k 2 x i + y i2 | 2 + 1 2 (1 -k 2 )x i - y i1 + y i2 2 2 = 1 2 2α 1 k 2 2 + (1 -k 2 ) 2 |x i | 2 + α 1 2 |y i1 | 2 + α 1 2 |y i2 | 2 + 1 8 |y i1 + y i2 | 2 + α 1 k 2 - 1 -k 2 2 x i • y i1 + α 1 k 2 - 1 -k 2 2 x i • y i2 = 1 2 2α 1 k 2 2 + (1 -k 2 ) 2 |x i | 2 + α 1 2 |y i1 | 2 + α 1 2 |y i2 | 2 + 1 8 |y i1 + y i2 | 2 = 1 2 (1 -k 2 )|x i | 2 + 1 4k 2 |y i1 | 2 + 1 4k 2 |y i2 | 2 - 1 8 |y i1 -y i2 | 2 ,
where we have defined again coefficient the

α 1 ∈ R * + by α 1 := 1 + α -k 1 ,
in order to absorb the above-mentioned x i -dependent remainder. In addition, note that we have used again the following relation between k 2 and α 1 to cancel the crossed term, i.e.,

k 2 = 1 1 + 2α 1 .
• Step 3: Change of variables x i1 → y i1 and x i2 → y i2 for i ∈ L n n-3 . We consider the change of variables given by

x i1 = k 3 x i + y i1 , x i2 = k 3 x i + y i2 ,
where k 3 is given by Definition 5.2. Putting together the terms for i ∈ L n n-3 in (5.16) and the above x i -dependent remainder in the above expression in Step 2 yields the following term under the above change of variables

α 2 2 |x i1 | 2 + α 2 2 |x i2 | 2 + 1 2 x i - x i1 + x i2 2 2 = α 2 2 |k 3 x i + y i1 | 2 + α 2 2 |k 3 x i + y i2 | 2 + 1 2 (1 -k 3 )x i - y i1 + y i2 2 2 = 1 2 2α 2 k 2 3 + (1 -k 3 ) 2 |x i | 2 + α 2 2 |y i1 | 2 + α 2 2 |y i2 | 2 + 1 8 |y i1 + y i2 | 2 + α 2 k 3 - 1 -k 3 2 x i • y i1 + α 2 k 3 - 1 -k 3 2 x i • y i2 = 1 2 2α 2 k 2 3 + (1 -k 3 ) 2 |x i | 2 + α 2 2 |y i1 | 2 + α 2 2 |y i2 | 2 + 1 8 |y i1 + y i2 | 2 = 1 2 (1 -k 3 )|x i | 2 + 1 4k 3 |y i1 | 2 + 1 4k 3 |y i2 | 2 - 1 8 |y i1 -y i2 | 2 ,
where we have defined again coefficient the

α 2 ∈ R * + α 2 := 1 + α -k 2 ,
in order to absorb the above-mentioned x i -dependent remainder. In addition, note that we have used again the following relation between k 3 and α 2 to cancel the crossed term, i.e.,

k 3 = 1 1 + 2α 2 .
Following a similar recursive process, we readily identify the quadratic form Q n in Definition 5.4 in the exponential of the integrand in the final expression (5.17).

• Step 4: Identifying the lineage maps Φ j n and the exponential factor. On the one hand, the hierarchy of changes of variables immediately leads to x j = Φ j n (x; y n ), for any leaf j ∈ L n , thanks to Definition 5.2 for the coefficients {κ n } n∈N and Definition 5.6 for the lineage maps Φ j n . On the other hand, since the Jacobian determinant of each change of variables is 1, no further factors appear during the iterative process except for the exponential x-dependent remainder e -1+α-kn 2 |x| 2 of the last step of the recurrence, which is not absorbed in the quadratic form Q n . □ 5.2. Probabilistic reinterpretation. For general purposes, in this part we reformulate the result in Proposition 5.7 in appropriate probabilistic terms. More precisely, this reformulation will not only provide shorter formulas, but it will also allow identifying a problematic key point when studying the asymptotics of the high-dimensional integral (5.17) in Section 6, namely, the presence of non-negligible correlations between the various factors indexed with indexes over leaves j ∈ L n of the tree, which do not dissipate even for long time n → ∞. We shall use the following notation.

Definition 5.8 (High-dimensional normal distributions). Consider any integer n ∈ N and

Σ m n :=     (2 -k n-m )k n-m 1 -k n-m I d - k 2 n-m 1 -k n-m I d - k 2 n-m 1 -k n-m I d (2 -k n-m )k n-m 1 -k n-m I d     , (5.18)
for every m = 0, . . . , n -1. Then, we will define the random vector Y n = (Y i ) i∈T n * distributed according to the following multivariate normal distribution

G n (y n ) := n-1 m=0 i∈L n m 1 (2π) d det(Σ m n ) exp - 1 2 (y T i1 , y T i2 )(Σ m n ) -1 y i1 y i2 , y n ∈ R 2(2 n -1)d , (5.19)
where again we are using the tree-indexed notation y n = (y i ) i∈T n * ∈ R 2(2 n -1)d in Remark 3.1. Notice that for every index i ∈ L n m with m = 0, . . . , n -1 the pair (Y i1 , Y i2 ) is normally distributed according to N (0, Σ m n ) and Y i1 and Y i2 are certainly correlated. Nevertheless, the vector (Y i1 , Y i2 ) is independent on any other random vector Y j for indices j ̸ = i1 and j ̸ = i2. We are now ready to introduce the following probabilistic reformulation of Proposition 5.7. Proposition 5.9 (Reformulation of the iterations II). Assume that α ∈ R * + and set any initial datum F 0 ∈ M + (R d ). Hence, the solution {F n } n∈N to the time-discrete problem (1.1) admits the following form

F n (x) = e -1+α-kn 2 |x| 2 (2π) d/2 2 2 n-1 d     n-1 m=0 4k 2 n-m 1-kn-m 2 m-1 d ∥F m ∥ 2 n-1-m L 1 (R d )     E   j∈L n F0 (Φ j n (x; Y n ))   , (5.20)
for every x ∈ R d and n ∈ N, where we are using the high-dimensional random vector Y n = (Y i ) i∈T n * in Definition 5.8 along with further notation from Definitions 5.1, 5.2, 5.4 and 5.6.

Proof. By Proposition 5.7 we obtain

F n (x) = e -1+α-kn 2 |x| 2 (4π) 2 n-1 d (2π) (2 n -1)d/2 1 n-1 m=0 ∥F m ∥ 2 n-1-m L 1 (R d ) R 2(2 n -1)d e -Qn(y n ) j∈L n F0 (Φ j n (x; y n )) dy n = (2π) (2 n -1)d/2 (4π) 2 n-1 d e -1+α-kn 2 |x| 2     n-1 m=0 4k 2 n-m 1-kn-m 2 m-1 d ∥F m ∥ 2 n-1-m L 1 (R d )     R 2(2 n -1)d G n (y n ) j∈L n F0 (Φ j n (x; y n )) dy n ,
for every x ∈ R d and each n ∈ N, where in the last line we have used that the inverse and determinant of the covariance matrices Σ m n in (5.18) take the form

(Σ m n ) -1 :=    2 -k n-m 4k n-m I d 1 4 I d 1 4 I d 2 -k n-m 4k n-m I d    , det(Σ m n ) = 4k 2 n-m 1 -k n-m d .
Then, the result follows by applying the law of the unconscious statistician (LOTUS) to relate the preceding integral to the expectation of the random variable j∈L n F0 (Φ j n (x; Y n )). □

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

Remark 5.10 (Non-negligible correlations). Since Y n in Definition 5.8 is normally distributed and the lineage map Φ j n (x; y n ) in Definition 5.6 are affine transformations of y n , then all (Φ j n (x; Y n )) j∈L n are also normally distributed. Unfortunately, the components of (Φ j n (x; Y n )) j∈L n are correlated since, in particular, all the components (i1, i2) are correlated in view of the structure (5.18) of their covariance matrices. Specifically, a straightforward computation shows that the covariance matrix 1 for the notation ⊗ of the Kronecker product) takes the following explicit form

K ij = E (Φ i n (x; Y n ) -E Φ i n (x; Y n )) ⊗ (Φ j n (x; Y n ) -E Φ j n (x; Y n )) , between any couple (i, j) ∈ L n × L n of components (see Table
K ij =                      n-1 q=0 κ 2 q (2 -k q+1 )k q+1 1 -k q+1 I d , if |i ∧ j| = n,   -κ 2 n-|i∧j|-1 k 2 n-|i∧j| 1 -k n-|i∧j| + n-1 q=n-|i∧j| κ 2 q (2 -k q+1 )k q+1 1 -k q+1   I d , if 0 < |i ∧ j| < n, -κ 2 n-1 k 2 n 1 -k n I d , if |i ∧ j| = 0.
Here, |i ∧ j| represents the level of the highest common descendant i ∧ j of i and j (cf. Section 3.1). Note that correlations only disappears for leaves (i, j) such that n -|i ∧ j| → ∞, that is, when the length of the path from i (or j) to the highest common descendant i ∧ j diverges as the height n of the tree increases. Unfortunately, correlations are non-negligible in all other cases even for n → ∞.

Uniqueness of equilibria and quantitative ergodicity property

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

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

E n [F ](x) := E   j∈L n F (Φ j n (x; Y n ))   , x ∈ R d , (6.1)
for any n ∈ N and any non-negative measurable function F : R d → R.

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

Γ(a + 1, x) ≤ 2e -x/2 , (6.3)
for every x, a ∈ R * + verifying the constraints

x > e, 2a < x log(x) . (6.4)

Proof. Let us restate the incomplete Gamma function as follows

Γ(a + 1, x) = +∞ x exp(-ψ a (s)s) ds,
where ψ a (s) is the function

ψ a (s) := 1 -a log(s) s , s > 0.
Notice that ψ a is monotonically increasing for s ∈ [e, ∞). Hence, we obtain that

ψ a (s) ≥ ψ a (x) ≥ 1 2 ,
for every s ≥ x, where we have used (6.4)

F ∈ C 1 (R d ) such that 0 < F (x) ≤ C, |∇ log F (x)| ≤ D(1 + |x|), (6.5)
for every x ∈ R d and appropriate C, D ∈ R * + . Then, for every ε ∈ R * + , there exists a large enough constant C ε ∈ R * + such that the following property holds true

|∇ log E n [F ](x)| ≤ C ε (2 k α +ε) n + (r α +ε) n |x| + e -(2+ε) n ε e Cε(2 rα +ε) n |x| 2 , (6.6)
for every x ∈ R d and any n ∈ N, where E n is the operator in Definition 6.1.

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

F (x) = log F (0) + 1 0 ∇(log F )(θx) • x dθ,
for any x ∈ R d . Using the linear growth assumption in (6.5) for the log-derivative ∇ log F we obtain that log F (x) ≥ log F (0) -D|x|

1 0 (1 + θ|x|) dθ ≥ log F (0) - D 2 -D|x| 2 ,
for any x ∈ R d , where in the last inequality we have integrated with respect to θ and we have used Young's inequality to interpolate |x| by |x| 2 . Therefore, we obtain the following Gaussian lower bound

F (x) ≥ ce -β|x| 2 , (6.7)
for each x ∈ R d and some c, β ∈ R * + depending only on F (0) and D. For simplicity of notation we define

E n := E n [F ],
for any n ∈ N. In the sequel, we shall study the behavior of ∇ log E n as n → ∞.

• in the tree contributes to the amount of 2 m leaves. Then, we obtain

j∈L n Φ j n (x; y n ) ≤ 2 n κ n |x| +   n-1 m=0 i∈L n n-m (2 m κ m ) 2   1/2   n-1 m=0 i∈L n n-m |y i | 2   1/2 = 2 n κ n |x| + n-1 m=0 2 n-m (2 m κ m ) 2 1/2 ∥ y n ∥ ≤ 2 n κ n |x| + √ 2 2 n/2 ∥ y n ∥, (6.8)
for any y n ∈ R 2(2 n -1)d , where in the first line we have used the Cauchy-Schwarz inequality on the joint sum over indices m and i, in the second line we have used the definition ∥ y n ∥ = ( i∈T n * |y i | 2 ) 1/2 of the ℓ 2 norm (cf. Table 1), and in the last line we have used that

n-1 m=0 2 n-m (2 m κ m ) 2 ≤ 2 n n-1 m=0 2 m k 2m 1 ≤ 2 n ∞ m=0 2 m k 2m 1 = 2 n 1 -2k 2 1 ≤ 2 n+1 .
Above, we have used that k n ≤ k 1 for all n ∈ N by virtue of Lemma 5.3(i) and also that

1 1 -2k 2 1 = 1 + 4α + 2α 2 2(1 + α) 2 = 2(1 + α 2 ) (1 + α) 2 + 2α + α 2 < 2.
• Step 2: Derivative of log E n . Taking derivatives on (6.1), using the second property in (6.5) along with (6.8) entails

|∇ log E n (x)| ≤ κ n E j∈L n |∇ log F (Φ j n (x; Y n ))| j∈L n F (Φ j n (x; Y n )) E j∈L n F (Φ j n (x; Y n )) ≲ κ n E j∈L n (1 + |Φ j n (x; Y n )|) j∈L n F (Φ j n (x; Y n )) E j∈L n F (Φ j n (x; Y n )) ≲ 2 n κ n + 2 n κ 2 n |x| + 2 n/2 κ n R n (x), (6.9)
for every x ∈ R d , where the remainder R n reads

R n (x) := E ∥ Y n ∥ j∈L n F (Φ j n (x; Y n )) E j∈L n F (Φ j n (x; Y n )) , x ∈ R d .
This estimation provides a first insight about the ergodicity property. Indeed the correlations, as measured by κ n , decreases faster than 2 -n (cf. Lemma 5.3), meaning that the contribution 2 n κ n decays fast to zero. It is also the case of the other contributions, locally in x, but the last one, namely 2 n/2 κ n R n (x), requires more work. Consider any √ 2 < q < 1 √ 2 kα , which exists because 2 k α < 1 for α ∈ R * + , and split

R n (x) = R 1 n (x) + R 2 n (x) in terms of the functions R 1 n (x) := 1 Z n (x) ∥ y n ∥≤q n ∥ y n ∥ j∈L n F (Φ j n (x; y n ))e -Qn(y n ) dy n , (6.10)
R 2 n (x) := 1 Z n (x) ∥ y n ∥>q n ∥ y n ∥ j∈L n F (Φ j n (x; y n ))e -Qn(y n ) dy n , (6.11)
and the normalization factor

Z n (x) := R 2(2 n -1)d j∈L n
F (Φ j n (x; y n ))e -Qn(y n ) dy n . (6.12)

• Step 3: Lower bound for Z n . Using the above lower bound by a Gaussian in (6.7) along with the upper bound for Q n in Lemma 5.5 implies the following estimate of Z n in (6.12):

Z n (x) ≥ c 2 n R 2(2 n -1)d exp   -β j∈L n |Φ j n (x; y n )| 2   exp(-β max ∥ y n ∥ 2 ) dy n ≥ c 2 n e -2 2n+1 κ 2 n β∥x∥ 2 R 2(2 n -1)d exp -(2 n+2 β + β max )∥ y n ∥ 2 dy n = c 2 n e -2 2n+1 κ 2 n β|x| 2 S 2(2 n -1)d-1 +∞ 0 r 2(2 n -1)d-1 exp -(2 n+2 β + β max )r 2 dr = c 2 n e -2 2n+1 κ 2 n β|x| 2 2(2 n+2 β + β max ) (2 n -1)d S 2(2 n -1)d-1 +∞ 0 s (2 n -1)d-1 e -s ds = c 2 n e -2 2n+1 κ 2 n β|x| 2 2(2 n+2 β + β max ) (2 n -1)d S 2(2 n -1)d-1 Γ((2 n -1)d),
where in the last step we have used the Gamma function, that is defined by Γ(a) := +∞ 0 s a-1 e -s ds, a > 0.

Recall that area of the hypersphere is

|S d-1 | = 2π d/2 Γ( d 2 )
. Then, we obtain

Z n (x) ≥ c 2 n π (2 n -1)d (2 n+2 β + β max ) (2 n -1)d e -2 2n+1 κ 2 n β|x| 2 , ( 6.13)
for every x ∈ R d .

• Step 4: Upper bound for R 1 n and R 2 n . On the one hand, notice that R 1 n (x) ≤ q n , (6.14) for all x ∈ R d by definition of R 1 n in (6.10). On the other hand, using the lower bound of Q n in Lemma 5.5 we can estimate R 2 n in (6.11) by

R 2 n (x) ≤ C 2 n q n Z n (x) ∥ y n ∥≥q n ∥y n ∥ 2 e -βmin∥ y n ∥ 2 dy n = C 2 n q n Z n (x) S 2(2 n -1)d-1 +∞ q n r 2(2 n -1)d+1 e -βminr 2 dr = C 2 n q n Z n (x) S 2(2 n -1)d-1 1 2β (2 n -1)d+1 min +∞ βminq 2n s (2 n -1)d e -s ds = C 2 n π (2 n -1)d q n β (2 n -1)d+1 min Z n (x) Γ (2 n -1)d + 1, β min q 2n Γ((2 n -1)d) , (6.15)
for all x ∈ R d , where we have used again the incomplete Gamma function (6.2) in Lemma 6.2. Putting the lower estimate (6.13) for Z n into (6.15) yields

R 2 n (x) ≲ e 2 2n+1 κ 2 n β|x| 2 C 2 n (2 n+2 β + β max ) (2 n -1)d q n c 2 n β (2 n -1)d min Γ (2 n -1)d + 1, β min q 2n Γ((2 n -1)d) , (6.16)
for all x ∈ R d . Now, take a = (2 n -1)d and x = β min q 2n and notice that the constraints (6.4) in Lemma 6.2 hold true for large enough n ∈ N as long as we choose q > √ 2. Hence, we obtain the estimate

Γ (2 n -1)d + 1, β min q 2n ≤ 2e -β min 2 q 2n
, for all n ≥ n * = n * (β min , d, q). Using Stirling's formula for the remaining Gamma function in (6.16), we conclude the estimate

R 2 n (x) ≲ e 2 2n+1 κ 2 n β|x| 2 C 2 n (2 n+2 β + β max ) (2 n -1)d q n c 2 n β (2 n -1)d min e -β min 2 q 2n (2 n -1) 1/2 (2 n -1)d e (2 n -1)d , (6.17)
for every x ∈ R d and each n ≥ n * .

• Step 5: Final conclusion. Let us put the preceding estimates (6.14) and (6.17) into (6.9) to achieve

|∇ log E n (x)| ≲ 2 n κ n + 2 n/2 q n κ n + 2 n κ 2 n |x| + M 2 n e β min 2 q 2n e 2 2n+1 κ 2 n β|x| 2 , ( 6.18)
for every x ∈ R d and each n ≥ n * , where M = M (C, c, β, β min , β max , d, q) is a universal constant that we have introduced to absorb all the double exponential factors of the form a 2 n with a > 0 in the fourth term along with any further lower order factor. On the one hand, choosing q > √ 2 sufficiently close to √ 2 and using (5.11) in Lemma 5.3 we obtain the following asymptotics for the factors in the first three terms of (6.18):

2 n κ n ≤ C ε (2 k α +ε) n , 2 n/2 q n κ n ≤ C ε (2 k α +ε) n , 2 n κ 2 n ≤ C ε (2(k α +ε) 2 ) n ≤ C ε (r α +ε) n , for arbitrarily small ε ∈ R * + and sufficiently large C ε ∈ R * + .
Note that in the last identity we have used the relation r α = 2 k 2 α . On the other hand, since q has been taken larger than √ 2 and arbitrarily close to √ 2, then the double exponential in the denominator of the last term in (6.18) kills the double exponential on the numerator independently on the explicit value of M . Therefore, we obtain

M 2 n e β min 2 q 2n ≤ C ε e -(2+ε) n ε ,
by appropriately increasing C ε if necessary. Similarly, we infer the following control for the remaining x-dependent exponential factor in the fourth term of (6.18)

e 2 2n+1 κ 2 n β|x| 2 ≤ e Cε(4(kα +ε) 2 ) n |x| 2 ≤ e Cε(2 rα +ε) n |x| 2 ,
which yields the announced estimate (6.6). □ Remark 6.4. Under the assumptions in Lemma 6.3, notice that the above result implies in particular that for every ε ∈ R * + and R ∈ R * + , there exists a large enough constant

C ε,R ∈ R * + such that sup |x|≤R |∇ log E n [F ](x)| ≤ C ε,R (2 k α ) n ,
for any n ∈ N. In particular, ∇ log E n [F ] converge to zero as n → ∞ uniformly on compact sets.

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

∈ R * + and R ∈ R * + , there exists a large enough constant C ε,R ∈ R * + such that sup |x|≤R F n (x) ∥F n ∥ L 1 (R d ) -F α (x) ≤ C ε,R (2 k α ) n , (6.19)
for any n ∈ N, where F α is the Gaussian profile (1.12). In particular, F n /∥F n ∥ L 1 (R d ) converge to F α as n → ∞ uniformly on compact sets.

In the proof of Corollary 6.5 we will require a careful control of the asymptotic behavior of the ratios F n (0)/∥F n ∥ L 1 (R d ) as n → ∞, which we provide below. First, note that if {F n } n∈N is the solution (3.4) of the time-discrete problem (1.1) issued at a Gaussian initial datum F 0 like in Proposition 3.7, we have

F n (0) ∥F n ∥ L 1 (R d ) = G µn,σ 2 n (0) → F α (0) = 1 (2π σ 2 α ) d/2
, as n → ∞ explicitly, where µ n and σ 2 n are determined by the recursive relations (3.5) and we have used that µ n → 0 and σ 2 n → σ 2 α according to Lemma 3.8. In the following lemma, we provide a similar quantitative result with convergence rates for generic initial data F 0 ∈ M + (R d ). Lemma 6.6. Assume that α ∈ R * + , set any initial datum F 0 ∈ M + (R d ) and consider the solution {F n } n∈N to the time-discrete problem (1.1). For any ε ∈ R * + there exists C ε ∈ R * + sufficiently large with

F n (0) ∥F n ∥ L 1 (R d ) -F α (0) ≤ C ε (2 k α +ε) n , (6.20)
for any n ∈ N.

Proof. Let us fix any ε ∈ R * + arbitrarily small so that 2 k α +ε < 1. Again, this is possible because 2 k α < 1 for α ∈ R * + . For simplicity of notation, we shall define the following sequence of coefficients:

c n := F n (0) ∥F n ∥ L 1 (R d )
, n ∈ N.

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

1 c n = R d e -1+α-kn 2 |x| 2 E n [ F0 ](x) dx E n [ F0 ](0) = |x|≤Rn e -1+α-kn 2 |x| 2 E n [ F0 ](x) E n [ F0 ](0) dx + 1 c n |x|>Rn F n (x) ∥F n ∥ L 1 (R d ) dx,
where we set the sequence of radii {R n } n∈N as follows

R n := 1 (2 k α ) n 2m+1 , n ∈ N, (6.21)
for a fixed value m ∈ N, which we take large enough so that

(2 k α ) 2m 2m+1 ≤ 2 k α + ε 2 . (6.22)
Whilst our choice of R n is not justified at first glance, we claim that it has been taken as to minimize the decay rate on the error terms E n and E n below. By solving the implicit equation on c n we infer

c n = 1 -Ẽn (2π σ 2 α ) d/2 + E n , (6.23)
for every n ∈ N, where the error terms Ẽn and E n take the form

Ẽn := |x|>Rn F n (x) ∥F n ∥ L 1 (R d ) dx, E n := |x|≤Rn e -1+α-kn 2 |x| 2 E n [ F0 ](x) E n [ F0 ](0) dx -(2π σ 2 α ) d/2 , (6.24)
Additionally, we can split the second error term as

E n = E n,1 -E n,2 + E n,3 with E n,1 := 2π 1 + α -k n d/2 - 2π 1 + α -k α d/2 , E n,2 := |x|>Rn e -1+α-kn 2 |x| 2 dx, E n,3 := |x|≤Rn e -1+α-kn 2 |x| 2 E n [ F0 ](x) E n [ F0 ](0) -1 dx, (6.25)
where we have used the relationship σ 2 α (1 + α -k α ) = 1, which results from (1.13) and (1.15), and therefore the decomposition of E n becomes clear since E n,1 can be reformulated as

E n,1 = R d e -1+α-kn 2 |x| 2 dx -(2π σ 2 α ) d/2 .
On the one hand, given any fixed value θ ∈ R + with θ < α 2 (for instance θ = α 4 ) the propagation of exponential moments in Corollary 4.8 implies that

E θ := sup n∈N R d e θ|x| 2 F n (x) ∥F n ∥ L 1 (R d ) dx < ∞.
By expanding the exponential in power series, one obtains in particular uniformly bounded moments of any order, and in particular, of order 2m, namely

M m := sup n∈N R d |x| 2m F n (x) ∥F n ∥ L 1 (R d ) dx ≤ E θ m! θ m < ∞.
Hence, we infer the following control on the first error term in (6.24)

Ẽn ≤ 1 R 2m n R d |x| 2m F n (x) ∥F n ∥ L 1 (R d ) dx ≲ 1 R 2m n , (6.26)
for any n ∈ N. On the other hand, by Lemma 5.3 we have that E n,1 in (6.25) can be bounded by

|E n,1 | ≲ k n -k α ≲ r n α . (6.27)
By direct calculation we also obtain

E n,2 = |S d-1 | +∞ Rn r d-1 e -1+α-kn 2 r 2 dr ≲ Γ d 2 , 1 + α -k n 2 R 2 n ≲ exp - 1 + α -k n 4(2 k α ) 2n 2m+1 , (6.28)
where Γ(a, x) is the incomplete Gamma function (6.2) and we have used estimate (6.3) in Lemma 6.2. Note that the constraint (6.4) is trivially satisfied since R n → ∞ by our choice (6.21). Finally, we control the error term E n,3 in (6.25) under the addition assumption that F0 satisfies the hypothesis (6.5) in Lemma 6.3. Whilst this is not always true, we show at the end of the proof that we can always assume so without loss of generality by replacing the argument on F0 by an advance enough time step Fm so that selection has properly shaped the Gaussian tails. Under this condition, Lemma 6.3 implies that given any ε

′ ∈ R * + there is C ε ′ ∈ R * + so that, for |x| ≤ R n , we have 1 0 (∇ log E n )[ F0 ](θx) • x dθ ≲ (2 k α +ε ′ ) n R n + (r α +ε ′ ) n R 2 n + e -(2+ε ′ ) n ε ′ e C ε ′ (2 rα +ε ′ ) n R 2 n R n ≲ 2 k α (2 k α ) 1 2m+1 + ε 2 n + r α (2 k α ) 2 2m+1 + ε 2 n + e -(2+ε ′ ) n ε ′ exp C ε ′ 2 r α (2 k α ) 2 2m+1 + ε 2 n 1 (2 k α ) n 2m+1 ≲ (2 k α ) 2m 2m+1 + ε 2 n ≲ (2 k α +ε) n .
In the second inequality above we have taken ε ′ small enough compared with ε. In the third inequality we have used the relation r α = 2 k 2 α , which guarantees that all the contributions in the third line are dominated by the first term ((2 k α ) 2m 2m+1 + ε 2 ) n . Finally, in the last inequality we have used our choice of m ∈ N large enough so that (6.22) holds. Therefore, the mean value theorem implies

E n [ F0 ](x) E n [ F0 ](0) -1 = exp 1 0 (∇ log E n [ F0 ])(θx) • x dθ -1 ≲ (2 k α +ε) n .
Hence, E n,3 in (6.25) can be controlled by Thereby, using the mean value theorem on (6.23) along with the bounds (6.26) and (6.30) concludes that

|E n,3 | ≤ |x|≤Rn e -1+α-kn 2 |x| 2 E n [ F0 ](x) E n [ F0 ](0) -1 dx ≲ (2 k α +ε) n . ( 6
c n - 1 (2π σ 2 α ) d/2 ≲ Ẽn + |E n | ≲ (2 k α +ε) n .
To end the proof, we show that we can always assume that the hypothesis (6.5) in Lemma 6.3 is satisfied without loss of generality. Indeed, note that after two iterations F2 already satisfies (6.5) 2 thanks to Corollary 4.2, but only two iterations are not enough to guarantee (6.5) 1 , in general. Nevertheless, applying Corollary 4.2 leads to the following upper bound on the normalized profiles:

Fm (x) ≤ C m e 1 2 α+ 1 2 -1 σ 2 m |x| 2
, for appropriate C m ∈ R * + , with variances {σ 2 m } m∈N defined by the recurrence (4.5), i.e.,

1 σ 2 m+1 = α + 1 1 + σ 2 m 2 , m ∈ N,
and with initial datum σ 2 1 = 1/α. The boundedness condition (6.5) 2 is then satisfied by Fm if we can show that the prefactor α+ 1 2 -1 σ 2 m in the above exponential bound becomes non-positive for large enough m. This is where the precise choice of normalization in Definition 5.1 plays a role. Specifically, recall that so defined σ 2 m → σ 2 α and we have precise relaxation rates by Lemma 3.8. Hence,

α + 1 2 - 1 σ 2 m = α + 1 2 - 1 σ 2 α + 1 σ 2 α - 1 σ 2 m ≤ α + 1 2 - 1 σ 2 α + C v r m α ,
for any m ∈ N. Since we chose a normalization Fn of the profiles by a Gaussian G 0,σ 2 with variance σ 2 strictly larger that the variance σ 2 α of the equilibrium F α (more specifically 1 σ 2 = 1 σ 2 = α + 1 2 ), we can guarantee that the right hand side above is non-positive if m ≥ n α for a sufficiently large n α ≥ 2 (depending only on α). This justifies that Fnα satisfies both conditions in (6.5). □

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

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

∇ log F n (x) = -(1 + α -k n )x + ∇ log E n [ F0 ](x),
for any x ∈ R d and every n ∈ N. Using the mean value theorem we then achieve log

F n (x) F n (0) -log F α (x) F α (0) = 1 0 (∇ log F n )(θx) • x dθ + 1 2 σ 2 α |x| 2 = - 1 + α -k n 2 |x| + 1 0 ∇ log E n [ F0 ](θx) • x dθ + 1 2 σ 2 α |x| 2 = k n -k α 2 |x| 2 + 1 0 ∇ log E n [ F0 ](θx) • x dθ, (6.31)
where we have used the relation 1/ σ 2 α = 1 + α -k α . On the one hand, the first term in the right hand side of (6.31) converges to zero thanks to Lemma 5.3. On the other hand, for the second term we shall apply Lemma 6.3. To do so we need to guarantee again that F0 verifies the hypothesis (6.5) of such a lemma. Recall that we can always assume the condition (6.5) without loss of generality (cf. last step in the proof of Lemma 6.6). Therefore, we can apply Remark 6.4 with F = F0 and obtain that ∇ log E n [ F0 ] converges to zero uniformly over compact sets with explicit convergence rates. Putting it into (6.31) we obtain more explicitly

sup |x|≤R log F n (x) F n (0) -log F α (x) F α (0) ≤ C ε,R (2 k α ) n ,
for any R ∈ R * + , any n ∈ N and a large enough C ε,R ∈ R * + . This, together with the above control on the asymptotics of ∥F n ∥ L 1 (R d ) /F n (0) in Lemma 6.6, allow proving (6.19). □

In particular, note that the above Corollary 6.5 is enough to prove the uniqueness of solutions to the eigenproblem (1.11), as stated in Corollary 1.2. Whilst the above local convergence result is enough to identify asymptotically the profile F α , a global result with quantitative convergence rates is still missing. In particular, note that the constants C ε,R above blow up when R → ∞ as we see explicitly in Lemma 6.3. A second drawback of this we are unable to characterize the long-time behavior of the mass ∥F n ∥ L 1 (R n ) in terms of the eigenvalue λ α via this method. In the following section, we give an answer to both questions by better exploiting the previous fundamental Lemma 6.3 and using the propagation of quadratic and exponential moments in Section 4.2. 6.2. Global convergence result. We are now ready to prove our main result. Let us emphasize that our final convergence result in Theorem 1.1 is presented using the Kullback-Leibler divergence, which is a very different metric from the log-Lipschitz type metrics used in the previous Section 6.1 for the local convergence results. As anticipated in Remark 1.4, this decision does not obey aesthetic reasons only, but we actually need to move from uniform norms (like the log-Lipschitz norm) to averaged norms (like the Kullback-Leibler divergence) in order to address the deficiency encountered in Lemma 6.3. Specifically, recall that for generic initial data F 0 ∈ M + (R d ), the log-Lipschitz norms of the high-dimensional integral E n [ F0 ] cannot be controlled uniformly due to additional exponentially growing terms.

Proof of Theorem 1.1.

• Step 1: Convergence of the profiles F n /∥F n ∥ L 1 (R d ) . Since F α is Gaussian (thus strongly log-concave), then the logarithmic-Sobolev inequality holds true and therefore we obtain the following control of the relative entropy by the relative Fisher information

D KL F n ∥F n ∥ F α ≤ σ 2 α 2 R d ∇ log F n (x) F α (x) 2 F n (x) ∥F n ∥ L 1 (R d )
dx, see Corollary 5.7.2 and Section 9.3.1 in [START_REF] Bakry | Analysis and Geometry of Markov Diffusion Operators[END_REF] for details. By the reformulation of F n in formula (5.20) of Proposition 5.9 and using the notation in Definition 6.1 we have that

∇ log F n (x) F α (x) = (k n -k α ) x + ∇ log E n [ F0 ](x).
Therefore, we obtain the following upper bound

D KL F n ∥F n ∥ F α ≲ D n,1 + D n,2 , (6.32)
where each factor takes the form

D n,1 := |k n -k α | 2 R d |x| 2 F n (x) ∥F n ∥ L 1 (R d ) dx, D n,2 := R d ∇ log E n [ F0 ](x) 2 F n (x) ∥F n ∥ L 1 (R d )
dx.

On the one hand, for D n,1 we can use the propagation of quadratic moments in Corollary 4.6 together with the explicit relaxation rates of {k n } n∈N to k α in Lemma 5.3 to show that D n,1 ≲ r 2n α . (6.33)

On the other hand, for D n,2 we shall refine the argument in the proof of Corollary 6.5. Specifically, recall again that we can assume that F0 satisfies the hypothesis (6.5) without loss of generality (cf. proof of Corollary 6.5). Then, applying Lemma 6.3 in order to control ∇ log E n [ F0 ] with fixed ε ∈ R * + yields

D n,2 ≲ (2 k α ) 2 + ε n R d 1 + |x| 2 + e 2 Cε(2 rα +ε) n |x| 2 F n (x) ∥F n ∥ L 1 (R d ) dx,
for some sufficiently large C ε ∈ R * + . Above we have used that (2 k α ) n is the decay rate that controls all the others in formula (6.6). Take n ≥ n ε sufficiently large such that

θ ε := sup n≥nε 2 C ε (2 r α +ε) n < α 2 .
Then, Corollaries 4.6 and 4.8 imply that the above quadratic and exponential moments are uniformly bounded with respect to n for n ≥ n ε . Therefore, we have

D n,2 ≲ ((2 k α ) 2 + ε) n . (6.34)
Finally, putting the estimates (6.33) and (6.34) into the split (6.32), and noticing that r 2 α < (2 k α ) 2 thanks to the relation r α = 2 k 2 α , ends this part of the proof. • Step 2: Convergence of the growth rates

∥F n ∥ L 1 (R d ) /∥F n-1 ∥ L 1 (R d ) .
For simplicity of notation we shall define the following sequence of coefficients:

λ n := ∥F n ∥ L 1 (R d ) ∥F n-1 ∥ L 1 (R d )
, n ∈ N.

Our goal then reduces to studying the asymptotic behavior of {λ n } n∈N and obtaining quantitative convergence rates. By definition of operator T in (1.2) we obtain that

λ n = R d e -m(x) R 2d G x - x 1 + x 2 2 F n-1 (x 1 ) ∥F n-1 ∥ L 1 (R d ) F n-1 (x 2 ) ∥F n-1 ∥ L 1 (R d ) dx 1 dx 2 dx.
By direct computation of the integral with respect to x as it was done in (4.11) we find that In addition, recalling the explicit form of the eigen-pair (λ α , F α ) in (1.12) implies

λ n = R 2d H(x 1 , x 2 ) F n-1 (x 1 ) ∥F n-1 ∥ L 1 (R d ) F n-1 (x 2 ) ∥F n-1 ∥ L 1 (R d ) dx 1 dx 2 , ( 6
λ α = R 2d
H(x 1 , x 2 ) F α (x 1 ) F α (x 2 ) dx 1 dx 2 . (6.37) Therefore, taking the difference of (6.35) and (6.37) and noticing that H in (6.36) is bounded yields for all φ ∈ C c (R d ) (cf. Table 1). We do not have a direct convergence result of F n /∥F n ∥ L 1 (R d ) in L 1 norms, but we do have convergence of the Kullback-Leibler divergence thanks to Step 1. By Pinsker's inequality, the latter metric controls the former, namely

|λ n -λ α | ≲ F n ∥F n ∥ L 1 (R d ) ⊗ F n ∥F n ∥ L 1 (R d ) -F α ⊗ F α L 1 (R 2d ) , ( 6
F n ∥F n ∥ L 1 (R d ) ⊗ F n ∥F n ∥ L 1 (R d ) -F α ⊗ F α L 1 (R 2d ) ≤ 1 √ 2 D KL F n ∥F n ∥ L 1 (R d ) ⊗ F n ∥F n ∥ L 1 (R d ) F α ⊗ F α = 1 √ 2 2 D KL F n ∥F n ∥ L 1 (R d )
⊗ F α , (6.39) the expected value for Gaussian solutions, namely, r α ≈ 0.472. These numerical results illustrate the two-step process arising in the relaxation dynamics: convergence towards a Gaussian profile occurs much faster than relaxation of the mean towards the origin due to (weak) selection. Alternatively speaking, the equilibrium variance builds up much faster than the center of the distribution gets to the origin.

Strong selection.

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

λ α ≈ 0.7944, σ 2 α ≈ 0.9221. The latter substantially differs from the value at linkage equilibrium (i.e., σ 2 α=0 = 2). In Figure 11a we note that the solution resembles a Gaussian distribution after a couple of iterations. After 15 iterations we obtain that the normalized profile F n /∥F n ∥ L 1 (R d ) has approximate mean and variance respectively given by ≈ 0.0017 and ≈ 0.9221. In Figure 11b we observe the convergence of the growth rates ∥F n ∥ L 1 (R d ) /∥F n-1 ∥ L 1 (R d ) towards the eigenvalue λ α . After 15 iterations the growth rate becomes ≈ 0.7944, which again agrees with the exact value λ α above up to 5 digits.

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

Conclusions and perspectives

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

Figure 1 .

 1 Figure 1. Parameter k α against parameter α.

Figure 2 .Figure 3 .

 23 Figure 2. Overall map of the proof of Theorem 1.1

Proposition 3 . 4 (

 34 Gaussian solution of the eigenproblem). Assume that α ∈ R * + , then the eigenproblem (1.11) has a unique Gaussian solution (λ α , F α ), determined by the relation (1.12) in Theorem 1.1.

  The malthusian exponent λα against parameter α for d = 1. The equilibrium variance σ 2 α against parameter α.

Figure 4 .

 4 Figure 4. Plot of eigenvalue λ α and variance at equilibrium σ 2 α against α for d = 1. The case α = 0 corresponds to absence of selection, which is conservative (λ α=0 = 1), and the variance at equilibrium σ 2 α=0 is twice the variance of the mixing kernel in B.

Figure 5 .

 5 Figure 5. The ratio r α against parameter α.

Figure 6 .

 6 Figure 6. Comparison of convergence rates in Theorem 1.1 and Corollary 3.9

4. 1 .

 1 Emergence of Gaussian tails.

Lemma 4 . 1 (

 41 Emergence of tails). Assume that α ∈ R * + , consider F ∈ M + (R d ) and set σ 2 , σ 2 ∈ R * + by

Lemma 4 . 5 ( 1 2

 451 Control of quadratic moments). Assume that α ∈ R * + and set any parameter η ∈ (1+α) 2 , 1 . Then, there exists a constant M

Definition 5 . 1 (

 51 Rescaled distributions). Consider anyF ∈ M + (R d ).Then, we define

Lemma 5 . 3 (

 53 Asymptotic behavior of the coefficients). Consider the sequences {k n } n∈N and {κ n } n∈N in Definition 5.2. Then, the following properties hold true: (i) (Coefficients k n ) The coefficients {k n } n∈N are positive numbers and k n ↘ k α as n → ∞, where k α ∈ R * + is the smallest root of the equation 1 3 + 2α -2 k α = k α .

Figure 7 .

 7 Figure 7. Unique path joining the leaf j = 121 to the root in the perfect binary tree T 3 . The corresponding lineage map takes the form Φ 121 n (x; y 3 ) = κ 3 x + κ 2 y 1 + κ 1 y 12 + y 121 .

Lemma 5 . 5 (

 55 Uniform positive definite quadratic forms). Consider the quadratic form Q n = Q n (y n ) in Definition 5.4 and define the couple of coefficients β min , β max ∈ R * + by

  d is represented according to the treeindexed notation in Remark 3.1, {κ n } n∈N is given in Definition 5.2, and c m = c • • • • • c is the m times iterated map c : T n * → T n which, to any vertex i ∈ T n * , it associates its child c(i) ∈ T n (cf. Section 3.1).

Lemma 6 . 2 .

 62 Consider the (upper) incomplete Gamma function with parameter a ∈ R * + , i.e., Γ(a, x) := +∞ x s a-1 e -s ds, x ∈ R * + . (6.2)

Step 1 :κ n x + n- 1 m=0κ

 11 Control of j∈L n |Φ j n (x; y n )|. By definition (5.15) of the lineage map Φ j n , we obtain that j∈L n |Φ j n (x; y n )| = j∈L n m y c m (j) ≤ 2 n κ n |x| + where we are using the tree-indexed notation y n = (y i ) i∈T n * ∈ R 2(2 n -1)d in Remark 3.1 and c : T n * -→ T n is defined in Section 3.1. Note that in the last inequality we have used that each index of the level L n n-m

Corollary 6 . 5 (

 65 Local convergence result). Assume that α ∈ R * + , set any initial datum F 0 ∈ M + (R d ) and consider the solution {F n } n∈N to the time-discrete problem (1.1). Then, for every ε

Proof of Corollary 1 . 2 .

 12 Let (λ, F ) be any solution of the eigenproblem(1.11). Then, the ansatz (1.10) defines a solution F n = λ n F of the time-discrete problem (1.1). By Corollary 6.5 we obtain thatlim n→∞ sup |x|≤R λ n F (x) λ n -F α (x) = 0,for any R ∈ R * + . Hence, F ≡ F α , and thus λ = λ α . □

1 + x 2 2 2 ,

 12 .35) where the function H(x 1 , x 2 ) takes the formH(x 1 , x 2 ) := 1 (1 + α) d/2 exp -α (1 + α) x (x 1 , x 2 ) ∈ R d .(6.36)

  .38) for any n ∈ N, where P ⊗ Q denotes the Kronecker product of two measuresP, Q ∈ M(R d ), namely R 2d φ(x, y) (P ⊗ Q)(dx, dy) = R d R dφ(x, y) P (dx) Q(dy),

Figure 8 .

 8 Figure 8. Relaxation of the normalized profiles F n /∥F n ∥ L 1 (R d ) towards the eigenfunction F α along the time iterations n = 0, 1, 2, 3, 4, 7 and 150 for a step function (7.1) as initial datum and weak selection parameter α = 0.015. The vertical dotted line represents the location of the mean of the equilibrium profile F α .

Figure 9 .

 9 Figure 9. Relaxation of the growth rates ∥F n ∥ L 1 (R d ) /∥F n-1 ∥ L 1 (R d ) towards the eigenvalue λ α along time iterations 0 ≤ n ≤ 150 for a step function (7.1) as initial datum and weak selection parameter α = 0.015.

Figure 10 .

 10 Figure 10. Numerical computation of the rates of convergence of the normalized profiles and the growth rates: (a) semi-log plot of the errors ε prof n := D KL ( F n /∥F n ∥ L 1 ∥ F α ), and (b) semi-log plot of the errors ε massn

  4k 1 |y 11 | 2 + |y 12 | 2 + |y 21 | 2 + |y 22 | 2 + 1 8 |y 11 -y 12 | 2 + |y 21 -y 22 | 2 × F0 (k 1 x 1 + y 11 ) F0 (k 1 x 1 + y 12 ) F0 (k 1 x 2 + y 21 ) F0 (k 1 x 2 + y 22 ) dx 1 dx 2 dy 11 dy 12 dy 21 dy 22 ,

	(5.5)

  Therefore, we obtain the following weak dependency of E n [F ] on the variable x as n → ∞. Lemma 6.3 (Asymptotics of E n ). Assume that α ∈ R * + and consider any

	on each of the above inequalities. Consequently, we get
	Γ(a + 1, x) ≤	+∞	e -s 2 ds,
		x
	and this ends the proof.		□

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 founding from the European Union's Horizon Europe research and innovation program under the Marie Sk lodowska-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.

where in last step we have used the tensorization property of the Kullback-Leilber divergence. The result then follows from (6.38)-(6.39) and the above explicit convergence rates of the normalized profiles F n /∥F n ∥ L 1 (R d ) in Step 1. □

Numerical experiments

For our numerical simulation, we have restricted to one-dimensional traits (i.e., d = 1) and we have considered a step function of the following form as initial datum: ), (7.1) where Z ∈ R * + is a normalizing factor so that F 0 ∈ P(R). Let {F n } n∈N be the solution to the timediscrete problem (1.1) starting at F 0 . Our numerical simulation is performed with Python on the finite computation domain [-15, 60] for the trait variable x, which contains the support of the previous initial datum F 0 and all the mass of each F n (except for a negligible Gaussian tail). In our simulation, we set ∆x = 0.001 as the step for the discretization of our computational domain. In particular, with such a step we compute numerical integrals with respect to x according to the left rectangle rule.

In the following, we explore numerically two different scenarios: weak selection and strong selection. Specifically, we obtain numerical approximations for the profiles F n in both regimes and we illustrate numerically the convergence of the growth rates

towards the eigenvalue λ α , and the relaxation of the normalized profiles F n /∥F n ∥ L 1 (R d ) towards the eigenfunction F α . As a consequence, we derive numerical approximations of the convergence rates to be compared with the theoretical results in this paper. More specifically, we note that the theoretical convergence rates in Theorem 1.1 are sharp, except a mismatch for the rates of growth of mass, which was discussed previously in Remark 3.10 (see also Figure 6). Indeed, we actually attain numerically the same convergence rates as in Corollary 3.9, which were sharp for Gaussian initial data. 7.1. Weak selection. In this part, we consider a small value α = 0.015. This leads to the following numerical values of the features of the equilibrium:

8897. characterizing the eigenpair (λ α , F α ), according to (1.12) and (1.13). Note that since the selection parameter has been taken very small, then the eigenvalue and the variance of the eigenfunction are close to those at linkage equilibrium, namely, λ α=0 = 1 and σ 2 α=0 = 2 (see Remark 3.5 and Figure 4). In Figure 8 we observe the relaxation of the normalized profiles F n /∥F n ∥ L 1 (R d ) towards the eigenfunction F α along the time iterations n = 0, 1, 2, 3, 4, 7, 150. We remark on the strong contraction of the variance during the first few iterations leading to a well-identified Gaussian-shaped profile at time n = 3. At time n = 7 the variance of the profile is already close to that of the equilibrium F α , but its mean is still substantially shifted to the right. Thereafter, we observe the gradual motion of the profiles towards the left at a much lower speed. After 150 iterations, the approximate mean and variance of F n /∥F n ∥ L 1 (R d ) are given by ≈ 0.0474 and ≈ 1.8897 respectively, where the latter agrees with the above exact value σ 2 α up to 5 digits.

In Figure 9 we represent the convergence of the growth rates

towards the eigenvalue λ α . After 150 iterations we obtain that the growth rate takes the approximate value ≈ 0.9857, which again agrees with the exact value λ α above up to 5 digits.

In Figure 10 we have represented the errors

in semi-log plots so that the horizontal axis appears in the natural scale and represents each time iteration, and the vertical axis contains the logarithm of the errors. As observed, both plots reduce to essentially straight lines, suggesting exponential relaxation. Interestingly, the numerical rate of convergence in the Kullback-Leibler divergence measured in Figure 10a is approximately 0.9441. It coincides numerically with the rate of relaxation among the subclass of Gaussian solutions (see Corollary 3.9), namely, λ 4 α ≈ 0.9441, which in turns is identical to the theoretical rate (2 k α ) 2 obtained in Theorem 1.1. These results are in perfect agreement with the other convergence results in Figure 10b: the numerical rate of convergence in the rate of growth of mass is approximately 0.9442, which is close to the one among the class of Gaussian solutions λ 4 α ≈ 0.9441, in contrast with the theoretical upper bound obtained in Theorem 1.1, namely 2 k α ≈ 0.9717. Moreover, the numerical rate of convergence of the variance of the normalized profiles

the survival probability of individuals. Our main result provides quantitative convergence rates of the renormalized distributions towards a unique stationary profile. Indeed, rates are exponential, which can be interpreted loosely as a spectral gap in this non-linear context. It is noticeable that the sexual reproduction operator is contractive under the Wasserstein distance (see Lemma 2.1). However, the generic incompatibility of multiplicative selection with transport distances becomes an apparent obstruction to the use of a direct perturbative approach in the regime of weak selection (see Example 2.4). This obstacle was circumvented by G. Raoul [START_REF]Exponential convergence to a steady-state for a population genetics model with sexual reproduction and selection[END_REF], assuming that selection is restricted to a compact support. We follow a different route in a purely non-perturbative setting. To this end, we restrict to the specific choice of a quadratic selection function of the trait, for the sake of simplicity. Our alternative approach relies on an appropriate study of the propagation of information along large binary trees of ancestors, by a suitable reformulation of high-dimensional integrals, inspired by the changes of variables performed in [START_REF] Calvez | Asymptotic analysis of a quantitative genetics model with nonlinear integral operator[END_REF] in the regime of small variance. This reveals an ergodicity property where the exact shape of the initial distribution is quickly forgotten across generations.

We remark that the above heuristic arguments indicate that the quadratic Wasserstein metric is not fully appropriate for dealing with the problem at hand. However, we could not identify yet a proper metric that extends the contraction property of the neutral case to the quadratic selection case.

Several perspectives are envisaged. First, there is an apparent price to pay with our method, in which we fully exploit the Gaussian structure induced by quadratic selection in order to perform tractable computations within the high-dimensional integrals. We believe though that the restriction to quadratic selection might be overcome in future works to allow for more general selection functions. Second, as studied in [START_REF] Calvez | Asymptotic analysis of a quantitative genetics model with nonlinear integral operator[END_REF], the case of multiple minima on the selection function leads to non-uniqueness of stable equilibria. Then, uncovering the hidden metastability and quantifying the relaxation towards a specific equilibrium is of great interest for its applications in quantitative genetics. Finally, a major problem is to transcend non-overlapping generations and tackle the full time-continuous model as presented in previous literature.

Appendix A. Nondimensionalization and derivation of the time-discrete version

In this section, we shall nondimensionalize the time-continuous model (1.6)-(1.3) with Gaussian mixing kernel G and quadratic selection function m. By time discretization on the Duhamel formulation, we will derive the time-discrete version (1.1) which has been central in this paper. Indeed, we shall show that we can reduce parameters into only one, namely, parameter α ∈ R + in the quadratic selection function (1.5). Bearing in mind all the biological parameters and dimensions, our evolution problem reads:

Here, r m , r b ∈ R + represent the mortality and birth rates and have frequency units, whilst the dimensionless corrections of m and B take the form

Here, G 0,σ 2 represents the Gaussian centered at 0 with variance σ 2 , and σ m and σ have the same units as the quantitative trait. On the one hand, σ m can be regarded as a characteristic unit quantifying the effective range of selection. Namely, if |x| ≥ σ m then m(x) ≥ 1. On the other hand, σ 2 is the genetic variance and 2σ 2 represents the variance at linkage equilibrium often denoted by σ 2 LE := 2σ 2 (cf. [START_REF] Bulmer | The mathematical theory of quantitative genetics[END_REF]). We nondimensionalize our system by appropriately scaling our variables as follows

Here, T and L are characteristic units for time and trait. For simplicity, we set them as follows in terms of the parameters of the system:

In other words, we scale time according to the mortality rate r m and the trait variable according to the genetic variance σ 2 . Dropping hats for simplicity yields the following dimensionless form of the equation

where B = B 1 involves again unit genetic variance thanks to our scaling assumption. This reduces the amount of free parameters to only two of them α, β ∈ R + defined as follows:

In the sequel, we justify our time-discrete version (1.11). Specifically, let us integrate equation (A.2) 1 for t ∈ [t n , t n+1 ] where {t n } n∈N ⊆ R + is any increasing sequence. Then, we recover Duhamel's formula

for any n ∈ N and each x ∈ R d . Using the rectangle rule for the integral in the right hand side leads to the approximations F n (x) ≃ f (t n , x) of the the above time-continuous problem for f = f (t, x):

for any n ∈ N and each x ∈ R d . Here, we have set unitary time steps t n -t n-1 = 1 for simplicity and t 0 = 0. We emphasize that the right hand side of (A.3) consists of two different terms. On the one hand, the first term describes the amount of old individuals from the previously generation n -1 having resisted the effect of selection in generation n. On the other hand, the second term computes the amount of offspring conceived by individuals in the previous generation n -1, where an eventual decline due to selection has also been taken into account. Let us emphasize that in (A.3) individuals from generation n -1 can merge (after breeding) with their offspring to form next generation n. This is typical in most biological populations, which have overlapping generations (e.g., some multivoltine flies like Drosophila melanogaster). However, some other systems have non-overlapping generations so that the full generation n -1 gets extinct after breeding and only their offspring survive at generation n (e.g., some univoltine insects like the Dawson's burrowing bee). We refer to [START_REF] Masaki | Geographic variation of diapause in insects[END_REF] for a discussion on the various types of voltinism and the role of diapause. In that case, by canceling the first term (A.3) reduces to

for any n ∈ N and each x ∈ R. We remark that two different parameters α and β are still present in (A.4). However, given that F → e -α 2 |x| 2 B[F ] is a 1-homogeneous operator we can kill the parameter β by rescaling the trait distributions F n . Specifically, define F n = β -n F n and note that (A.4) reduces to our time-discrete problem (1.1) with only one parameter α.