Ergodicity of the Fisher infinitesimal model with quadratic selection

Nonlinear Analysis 238 (2024) 113392 Contents lists available at ScienceDirect Nonlinear Analysis www.elsevier.com/locate/na Ergodicity of the Fisher infinitesimal model with quadratic selection Vincent Calvez a, Thomas Lepoutre b, David Poyato c,∗ aInstitut Camille Jordan (ICJ), UMR 5208 CNRS & Université Claude Bernard Lyon 1, 69100 Villeurbanne, France bUniv Lyon, Inria, Université Claude Bernard Lyon 1, CNRS UMR5208, Institut Camille Jordan, F-69603 Villeurbanne, France cDepartamento de Matemática Aplicada and Research Unit “Modeling Nature” (MNat), Facultad de Ciencias, Universidad de Granada, 18071 Granada, Spain a r t i c l e i n f o Article history: Received 24 February 2023 Accepted 13 September 2023 Communicated by Enrico Valdinoci MSC: 35B40 35P30 35Q92 47G20 92D15 Keywords: Integro-differential equations Asymptotic behavior Nonlinear spectral theory Quantitative genetics abstract We study the convergence towards a unique equilibrium distribution of the solutions to a time-discrete model with non-overlapping generations arising in quantitative genetics. The model describes the dynamics of a phenotypic distribution with respect to a multi-dimensional trait, which is shaped by selection and Fisher’s infinitesimal model of sexual reproduction. We extend some previous works devoted to the time-continuous analogs, 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. ©2023 Elsevier Ltd. All rights reserved. 1. Introduction Fisher’s infinitesimal model (also known as the polygenic model) is a widely used statistical model in quantitative genetics initially proposed by R. Fisher [1]. 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 ∗Corresponding author. E-mail addresses: vincent.calv[email protected] (V. Calvez), thomas.lep[email protected] (T. Lepoutre), davidpoy[email protected] (D. Poyato). https://doi.org/10.1016/j.na.2023.113392 0362-546X/©2023 Elsevier Ltd. All rights reserved. V. Calvez, T. Lepoutre and D. Poyato Nonlinear Analysis 238 (2024) 113392 value of parents’ traits, with a constant (genetic) variance across generations. This model allowed reconciling Mendelian inheritance and the continuous trait variations documented by F. Galton via a Central Limit Theorem. More specifically, by taking limits when the number of underlying loci tends to infinity on a model with Mendelian inheritance, N. Barton,A. Etheridge and A. V´ eber [2] 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∈Rd in a population undergoing sexual reproduction and the effect of natural selection. Specifically, starting at any initial configuration F0∈ M+(Rd) of trait distribution, we analyze the long term dynamics of trait distributions {Fn}n∈Nacross successive generations n∈N, which solves the following recursion Fn=T[Fn−1],(1.1) for any n∈N. The operator Tencodes the balanced effect of sexual reproduction in the population (under the infinitesimal model) and natural selection. Specifically, Tis defined by T[F] := e−mB[F],(1.2) for any trait distribution F∈ M+(Rd). 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 Bis chosen to be Fisher’s infinitesimal operator and it takes the form B[F](x) := ∫R2d G(x−x1+x2 2)F(x1)F(x2) ∫RdF(x′)dx′dx1dx2, x ∈Rd,(1.3) for any trait distribution F∈ M+(Rd). 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 x1, x2∈Rdwill mate and yield a descendant with trait value x∈Rd. In other words the resulting trait distributes around the mean value x1+x2 2of 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 Fhaving resisted the effect of selection. The above sexual reproduction operator Bhas recently pulled the attention of both the applied and more theoretical communities, cf. [2–7]. In this paper we shall restrict to Gaussian mixing kernel and quadratic selection function, i.e., G(x) := 1 (2π)d/2e−|x|2 2, x ∈Rd,(1.4) m(x) := α 2|x|2, x ∈Rd,(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 [2]. For simplicity, we set a Gaussian Gwith unit variance, but any value of the genetic variance could also be considered (see nondimensionalization in Appendix). Before introducing our results, we shall relate the previous time-discrete problem (1.1) to analogous timecontinuous 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): {∂tf=−m(x)f+R[f], t ≥0, x ∈Rd, f(0, x) = f0(x), x ∈Rd,(1.6) 2 V. Calvez, T. Lepoutre and D. Poyato Nonlinear Analysis 238 (2024) 113392 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∈Rdper 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 Racross 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) := ∫Rd K(x−y)f(t, y)dy, x ∈Rd.(1.7) The factor K(x−y) determines the probability that an individual with trait value y∈Rdproduces a descendant with trait value x∈Rd(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 [8–10], including an additional density-dependent competition term that makes the analysis non standard (see e.g. [11] and references therein for the well-posedness of the constrained Hamilton–Jacobi equation derived in the limit). Recently, inspired by the infinitesimal model, the case of sexual reproduction has been addressed by invoking the preceding nonlinear version Bin (1.3) as reproduction operator R=B. Several asymptotic regimes have been addressed: large reproduction rate [4,6,7], small variance asymptotics [3,5]. 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λtF(x),(t, x)∈R+×Rd,(1.8) The parameter λ∈Rrepresents 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 ∈Rd, ∫RdF(x′)dx′= 1.(1.9) Note that the operator Bis genuinely non-linear so that methods based on the Krein–Rutman theory or maximum principles cannot be applied straightforwardly, see [12,13] and the references therein for the linear case. Further, usual extensions of the Krein–Rutman theory to 1-homogeneous operators [14] cannot be applied neither because Bis 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)−1B[F] over an appropriate cone of L1(Rd) that is conserved by the nonlinear operator, see [15]. However, uniqueness cannot be achieved by this method. Moreover, it has been proven in [3, Corollary 1.5] that several equilibrium states (λ, F) can co-exist in the presence of multiple local minima of m(provided the variance of Gis 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=Band a trait dependent fecundity rate [7]. 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 Bin the space of probability measures sharing the same center of mass. Unfortunately, the Wasserstein metric is 3 V. Calvez, T. Lepoutre and D. Poyato Nonlinear Analysis 238 (2024) 113392 not fully compatible with multiplicative operators, such as trait-dependent fecundity (see the discussion in [7, Section 3.4] and Section 2below). This can be circumvented under the additional assumption that the trait density is locally uniformly bounded below, following [16]. 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 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 Fn(x) = λnF(x),(n, x)∈N×Rd.(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 ∈Rd, ∫RdF(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 λα∈Rand Fαbeing a probability measure, for each α∈R∗ +? (Q2) Consider any generic initial datum F0∈ M+(Rd) and its associated solution {Fn}n∈Nof the timediscrete problem (1.1). Do the renormalized profiles Fn/∥Fn∥L1(Rd)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 [3,15] 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 {Fn}n∈Ntowards Fα. First, we prove that an explicit (Gaussian) solution to the eigenproblem (1.11) exists. Second, by computing niterations of the operator T, that is Fn=Tn[F0], we notice that information of Fnat the trait value xis propagated from the initial datum F0across 2nancestors over a binary tree with height nand rooted at x(the pedigree chart). Interestingly, an appropriate reformulation of Tnin the case of Gaussian mixing (1.4) and quadratic selection (1.5) shows that the dependence of the solution {Fn}n∈Non the initial datum F0is rapidly lost across the different levels of the tree. More specifically, a strong convergence of generic solutions {Fn}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 {Fn}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 2for 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 F0∈ M+(Rd). The solution {Fn}n∈Nto the time-discrete problem (1.1) verifies that the growth rates ∥Fn∥L1(Rd)/∥Fn−1∥L1(Rd)relax 4 V. Calvez, T. Lepoutre and D. Poyato Nonlinear Analysis 238 (2024) 113392 towards λαand the normalized profiles Fn/∥Fn∥L1(Rd)relax towards Fαwith λα:= (1 + α(1 + σ2 α 2))−d 2 ,Fα:= G0,σ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 DKL (Fn ∥Fn∥L1(Rd)Fα)≤Cε((2kα)2+ε)n, ⏐⏐⏐⏐⏐∥Fn∥L1(Rd) ∥Fn−1∥L1(Rd)−λα⏐⏐⏐⏐⏐≤Cε(2kα+ε)n, for any n∈N, where DKL is the Kullback–Leibler divergence (or relative entropy), i.e., DKL(P∥Q) := ∫Rd P(x) log (P(x) Q(x))dx, (1.14) for any P, Q ∈L1 +(Rd)∩P(Rd), 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 behavior at infinity, which crucially helps to localize the convergence argument. As a by-product, we are able to consider very general initial condition F0. Two of the authors, together with F. Santambrogio, have obtained similar results when the selection function mis more generally assumed to be strongly convex [17]. However, they imposed stringent conditions on the initial datum F0, 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 [6,7], 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: 5 V. Calvez, T. Lepoutre and D. Poyato Nonlinear Analysis 238 (2024) 113392 (i) (Quadratic Wasserstein distance) When α= 0, the operator Treduces 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 Twhich 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 Tis non-expansive in the log-Lipschits norm ∥∇log F Fα∥L∞(Rd)for all α≥0, and indeed it is contractive if α > 0. However, this contraction only gives actual information when the initial datum F0has 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α∥L2(Rd,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 µ∈Rdof 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 2we 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 3we 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 4is devoted to introduce some main properties of Tregarding 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 5we 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 6we prove our main results, namely, Theorem 1.1 and Corollary 1.2. Section 7contains some numerical experiments that illustrate the results in this paper. In Section 8we provide some conclusions and perspectives. Finally, Appendix contains the adimensionalization of the problem, and the relationship between (1.1) and the previous time-continuous analogs in the literature is discussed. A full list of the main notations in the paper is presented in Table 1. 2. Motivation and strategy of the proof of Theorem 1.1 In this section, we discuss the incompatibility of the quadratic Wasserstein distance to quantify directly contractivity under the joint effect of the reproduction operator Bin (1.3) and a generic multiplicative 6 V. Calvez, T. Lepoutre and D. Poyato Nonlinear Analysis 238 (2024) 113392 Table 1 List of notations. Notation Meaning Reference M(Rd),M+(Rd)Signed and non-negative finite Radon measures Theorem 1.1 P(Rd),P2(Rd)Probability measures (with finite 2nd order moment) Lemma 2.1 P⊗Q∈ M(R2d)Tensor product of the finite measures P, Q ∈ M(Rd)Eq. (6.39) DKL(P|Q)Kullback–Leibler divergence between P, Q ∈ P(Rd)Eq. (1.14) W2(P, Q)Quadratic Wasserstein distance between P, Q ∈ P2(Rd)Eq. (2.2) Gµ, σ2Gaussian with mean µ∈Rdand covariance σ2Id∈Rd×dSection 3.2 Gµ, ΣGaussian with mean µ∈Rdand covariance Σ∈Rd×dDefinition 5.8 N(µ, Σ)Normal with mean µ∈Rdand covariance Σ∈Rd×dDefinition 5.8 E[X]Expectation of a random variable XEq. (5.20) G(x)Gaussian mixing kernel G0, IdEq. (1.4) m(x)Quadratic selection function α 2|x|2Eq. (1.5) B[F]Fisher’s infinitesimal operator Eq. (1.3) T[F]Selection-reproduction operator Eq. (1.2) M[F]Normalized multiplicative operator by e−mDefinition 2.3 S[F]Scaled selection-reproduction operator Definition 4.4 En[F]High-dimensional integral operators Definition 6.1 {Fn}n∈NSolution to the time-evolution equation (1.1) Eq. (1.1) (λα,Fα)Gaussian solution to the non-linear eigenproblem (1.11) Eq. (1.12) σ2 αVariance of Gaussian eigenfunction FαEq. (1.13) ¯ F=emF Fα=0 Normalized profile associated to F∈ M+(Rd)Definition 5.1 Tn,Tn ∗,ˆTnPerfect binary tree, rootless tree and leafless tree Section 3.1 Ln m,LnLevel mand leaves (level n) of the tree Section 3.1 i1, i2Parents of a node i∈ˆTnof the tree Section 3.1 xn= (xi)i∈Tn ∗ yn= (yi)i∈Tn ∗ Variables indexed by the rootless tree Remark 3.1 zn= (zj)j∈LnVariables indexed by leafs Remark 3.1 ∥x∥=(∑n i=1 |xi|2)1/2ℓ2sum of Euclidean norms of x= (x1, ...,xn)∈Rdn Eq. (6.8) Φj n(x;yn)Lineage map from leaf j∈Lnto root value x∈RdDefinition 5.6 u⊗v∈Rd×dKronecker product (uivj)1≤i, j≤Nof vectors u, v ∈RdRemark 5.10 kn,κnSequences of coefficients in the change of variables Definition 5.2 (2kα)2Convergence rate of DKL Eq. (1.15) 2kαConvergence rate of the log-Lipschitz norm Lemma 2.5 rαRelaxation rates of variances recursion Lemma 3.8 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, [7] and also [6]), 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), then Bpreserves the mass and center of mass, namely, ∥B[F]∥L1(Rd)=∥F∥M(Rd),∫Rd xB[F](x)dx =∫Rd x F(dx),(2.1) for any F∈ M+(Rd). 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. [7]) to analyze 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: 7 V. Calvez, T. Lepoutre and D. Poyato Nonlinear Analysis 238 (2024) 113392 Lemma 2.1. Assume that F1, F2∈ P2(Rd)and they have the same center of mass. Then, W2 2(B[F1],B[F2]) ≤1 2W2 2(F1, F2). Proof. Recall the following dual characterizations of the quadratic Wasserstein distance (cf. [18,19]): W2 2(F1, F2) = inf {∫R2d|x−y|2γ(dx, dy) : γ∈ P(Rd×Rd), π1#γ=F1, π2#γ=F2} = sup {∫Rd ϕ1dF1+∫Rd ϕ2dF2:|ϕ1(x) + ϕ2(y)| ≤ |x−y|2∀x, y ∈Rd}, (2.2) where πi:Rd×Rd→Rdis the projection onto the ith component for i= 1,2. Taking any couple (ϕ1, ϕ2) as above and using the specific form of the operator Bwe obtain: ∫RdB[F1](x)ϕ1(x)dx +∫RdB[F2](y)ϕ2(y)dy =∫R3d G(x−x1+x2 2)ϕ1(x)F1(dx1)F1(dx2)dx +∫R3d G(y−y1+y2 2)ϕ2(y)F2(dy1)F2(dy2)dy =∫R3d G(z)ϕ1(z+x1+x2 2)F1(dx1)F1(dx2)dz +∫R3d G(z)ϕ2(z+y1+y2 2)F2(dy1)F2(dy2)dz. Consider any transference plan γ∈ P(Rd×Rd) between (F1, F2) as above. Specifically, we have that π1#γ=F1and π2#γ=F2. Then, we can gather both integrals as follows: ∫RdB[F1](x)ϕ1(x)dx +∫RdB[F2](y)ϕ2(y)dy =∫R5d G(z)(ϕ1(z+x1+x2 2)+ϕ2(z+y1+y2 2))γ(dx1, dy1)γ(dx2, dy2)dz. Since the condition |ϕ1(x) + ϕ2(y)|≤|x−y|2is verified for all x, y ∈Rd, then we find ∫RdB[F1](x)ϕ1(x)dx +∫RdB[F2](y)ϕ2(y)dy ≤1 4∫R5d G(z)|(x1+x2)−(y1+y2)|2γ(dx1, dy1)γ(dx2, dy2)dz =1 4∫R4d|(x1+x2)−(y1+y2)|2γ(dx1, dy1)γ(dx2, dy2) =1 4∫R2d|x1−y1|2γ(dx1, dy1) + 1 4∫R2d|x2−y2|2γ(dx2, dy2), where in the third line we have used that Gis a probability density and in the last line we have used the crucial fact that F1and F2share 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 F1and F2have the same center of mass has been crucially used to cancel the crossed term. Otherwise, by the Cauchy–Schwarz inequality we would merely obtain non-expansiveness: W2(B[F1],B[F2]) ≤ W2(F1, F2).(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. 8 V. Calvez, T. Lepoutre and D. Poyato Nonlinear Analysis 238 (2024) 113392 Corollary 2.2. Set any initial datum F0∈ M+(Rd)such that ∫Rd|x|2F0(dx)<∞and consider the solution {Fn}n∈Nto the time-discrete problem (1.1) with α= 0, i.e., Fn=B[Fn−1]for all n∈N. Then, ∥Fn∥L1(Rd) ∥Fn−1∥L1(Rd) = 1,W2(Fn ∥Fn∥L1(Rd) , Gµ0,2)≲1 2n/2, for every n∈N, where µ0:= ∫Rdx F0(dx). In particular, the set of stationary distributions under Bis {Fα=0(· −µ) : µ∈Rd}, where Fα=0 =G0,2will denote here on the Gaussian centered at the origin with variance equals 2in 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 [20] we obtain the relation W2 2(Fn ∥Fn∥L1(Rd) , Gµ0,2)≤4DKL (Fn ∥Fn∥L1(Rd)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 F0∈ M+(Rd) one cannot expect that the normalized profiles Fn/∥Fn∥L1(Rd)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 towards 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 2kα=0 = 1. 2.2. Incompatibility of Wasserstein metric with multiplicative operators Note that by definition (1.2), our operator Tis the composition of the sexual reproduction operator Bwith 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−mF ∥e−mF∥M(Rd) , F ∈ M+(Rd)\{0}. However, the latter is not Lipschitz continuous with respect to the quadratic Wasserstein metric. Hence the composition of Band Mis not expected to be contractive, even in the case of weak selection, without any additional restriction. We illustrate such a Lipschitz discontinuity of Min the following example. Example 2.4. Suppose that the m∈C1 +(R) is radially symmetric around the origin and set F1, F2∈ P2(R) as the sum of two Dirac masses, F1being symmetric, and F2nearly symmetric. More precisely, F1=1 2δ−h+1 2δh, F2=1 2δ−h+ε+1 2δh+ε, where h > 0 is fixed so that m′(h)= 0 and ε > 0 is small. Note that M[F1] = F1,M[F2] = (1 −pε)δ−h+ε+pεδh+ε, 9 V. Calvez, T. Lepoutre and D. Poyato Nonlinear Analysis 238 (2024) 113392 (v) (Mate) Given a vertex i∈Tn ∗, we define the mate of iand we denote it m(i) as the only other vertex in Tnthat 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 ∈Tn, we say that i≤jif the associated words are so ordered according to the lexicographical order of the set of words W2with two letters {1,2}. (vii) (Highest common descendant) Given two vertices i, j ∈ˆ Tn, we define the highest common descendant of iand j, and we denote it by i∧j, as i∧j:= max{l∈Tn: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 R2(2n−1)d≡(Rd)Tn ∗and R2nd≡ (Rd)Ln. Specifically, vectors xn,yn∈R2(2n−1)dand zn∈R2ndwill be regarded as indexed families xn= (xi)i∈Tn ∗,yn= (yi)i∈Tn ∗,zn= (zj)j∈Ln,(3.1) where xi, yi∈Rdand zj∈Rdfor each i∈Tn ∗and j∈Ln. 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∈Rdand variances σ2 1, σ2 2>0. Using the above result, we can obtain the following explicit evaluation of the operator Tin (1.2) over the class of Gaussian functions. Lemma 3.3 (Evaluation on Gaussians).Consider any µ∈Rdand σ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), 16 V. Calvez, T. Lepoutre and D. Poyato Nonlinear Analysis 238 (2024) 113392 for each x∈Rd. Therefore, by definition of Tin (1.2) we obtain T[Gµ,σ2](x) = e−α 2|x|2Gµ,1+ σ2 2 (x) = 1 (2π(1 + σ2 2))d/2exp (−α 2|x|2−1 2 1 1 + σ2 2|x−µ|2), for each x∈Rd. By completing the square inside the exponential, we conclude our result. □ Consequently, the following explicit Gaussian solution of the eigenproblem (1.11) is found. Proposition 3.4 (Gaussian Solution of the Eigenproblem).Assume that α∈R∗ +, then the eigenproblem (1.11) has a unique Gaussian solution (λα,Fα), determined by the relation (1.12) in Theorem 1.1. Proof. We look for λ∈R∗ +,µ∈Rdand σ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 σ2are 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 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 Fig. 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 Gin (1.4), a famous result in quantitative genetics, see e.g. [23]. Note that λα<1 for any α∈R∗ +and, consequently, the special steady solutions Fn=λn αFαcoming from ansatz (1.10) always get extinct for large n. This is an artificial consequence of our parameter reduction in Appendix. 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 µ∈Rd. 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 17 V. Calvez, T. Lepoutre and D. Poyato Nonlinear Analysis 238 (2024) 113392 Fig. 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. 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 Tand 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 m0∈R∗ +,µ0∈Rd and σ2 0∈R+, and set the Gaussian initial datum F0=m0Gµ0,σ2 0. Hence, the solution {Fn}n∈Nto the time-discrete problem (1.1) takes the form Fn=mnGµn,σ2 n,(3.4) for n∈N, where the parameters mn∈R∗ +,µn∈Rdand σ2 n∈R∗ +are governed by the recursions: mn=mn−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.5) Lemma 3.8 (Relaxation of Variance).Assume that α∈R∗ +, consider any σ2 0∈R∗ +and define the sequence {σ2 n}n∈Nby recursion according to the third recursion in (3.5), i.e. 1 σ2 n =α+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 α| ≤ Cvrn α,(3.7) for every n∈N, where the constant Cv∈R+depends on σ2 0and α(we obtain that Cv= 0 if σ2 0=σ2 α), and the ratio rα∈(0,1 2)(see Fig. 5) satisfies rα= 2k2 αand is given explicitly by rα:= 8 ((2α+ 3) + √(2α+ 1)2+ 8α)2.(3.8) 18 V. Calvez, T. Lepoutre and D. Poyato Nonlinear Analysis 238 (2024) 113392 Proof. •Step 1: Monotone convergence of {σ2 n}n∈N. Consider the function f:R+−→ Rgiven by f(x) = α+1 1 + 1 2x , x ∈R+, and, for simplicity of notation, define x∗:= (σ2 α)−1, xn:= (σ2 n)−1, n ∈N.(3.9) Notice that x∗is the unique fixed point of fin R+and {xn}n∈Ndetermines the fixed-point iteration of the map fissued at x0= (σ2 0)−1, that is, xn=f(xn−1), n ∈N. Our goal is to show that {xn}n∈Nconverges 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 fis 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) xis 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 {xn}n∈Nand, in addition, min{x∗, x0} ≤ xn≤max{x∗, x0},(3.10) for any n∈N. In particular, this yields the monotonic convergence of {xn}n∈Ntowards x∗, or equivalently, the claimed monotonic convergence of {σ2 n}n∈Ntowards σ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¨obius function M:R\{−1 2} −→ Rgiven by M(x) := 2(α+ 1)x+α 2x+ 1 , x ∈R\{−1 2}. Note that Mhas two fixed points: x±:= 2α+ 1 ±√(2α+ 1)2+ 8α 4, where we note that x+=x∗∈R∗ +and x−∈R−. Then, by definition of {xn}n∈Nin (3.9) we obtain xn−x+=2(α+ 1)xn−1+α 2xn−1+ 1 −2(α+ 1)x++α 2x++ 1 =2 (2xn−1+ 1)(2x++ 1)(xn−1−x+).(3.11) Since x0(thus xn−1for 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 (2xn−1+ 1)(2x++ 1) ≤2σ2 α 2 + σ2 α = 2kα, which gives true contraction in (3.11) when α > 0, since 2kα<1. Nevertheless, such a rate is non-optimal, and we show an alternative approach to achieve a sharper one. Specifically, note that xn−x−=2(α+ 1)xn−1+α 2xn−1+ 1 −2(α+ 1)x−+α 2x−+ 1 =2 (2xn−1+ 1)(2x−+ 1)(xn−1−x−).(3.12) 19 V. Calvez, T. Lepoutre and D. Poyato Nonlinear Analysis 238 (2024) 113392 Fig. 5. The ratio rαagainst parameter α. By dividing (3.11) by (3.12) and iterating the identity we obtain xn−x+ xn−x− =(2x−+ 1 2x++ 1)nx0−x+ x0−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 |xn−x∗|=|xn−x−| |x0−x−||x0−x∗|rn α≤max{x∗, x0}−x− min{x∗, x0}−x−|x0−x∗|rn α,(3.14) for every n∈N, thus leading to an improved rate since rα= 2k2 α<2kα. Finally, notice that |σ2 n−σ2 α|=σ2 nσ2 α|xn−x∗| ≤ max{σ2 α, σ2 0}2|xn−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 m0∈R+,µ0∈Rd and σ2 0∈R+, and set the Gaussian initial datum F0=m0Gµ0,σ0. Hence, the Gaussian solution {Fn}n∈N to the time-discrete problem (1.1) verifies that the growth rates ∥Fn∥L1(Rd)/∥Fn−1∥L1(Rd)relax towards λα and the normalized profiles Fn/∥Fn∥L1(Rd)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∗ + DKL (Fn ∥Fn∥L1(Rd)Fα)≤Cµ,ε(λ4/d α+ε)n+Cvrn α, ⏐⏐⏐⏐⏐∥Fn∥L1(Rd) ∥Fn−1∥L1(Rd)−λα⏐⏐⏐⏐⏐≤Cµ,ε(λ4/d α+ε)n+Cvrn α, 20 V. Calvez, T. Lepoutre and D. Poyato Nonlinear Analysis 238 (2024) 113392 where DKL is the Kullback–Leibler divergence (or relative entropy), rαis given by (3.8) in Lemma 3.8, Cµ,ε ∈R+depends on α,µ0,σ2 0and ε, and Cv∈R+depends only αand σ2 0. In fact, Cµ,ε = 0 if µ0= 0 and Cv= 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: DKL(Gµ1,Σ1∥Gµ2,Σ2) = 1 2(µ2−µ1)⊤Σ−1 2(µ2−µ1) + 1 2trace(Σ1Σ−1 2)−d 2+1 2log det(Σ2Σ−1 1),(3.16) for every µ1, µ2∈Rdand any positive non-singular matrices Σ1,Σ2∈Rd×d. Proof of Corollary 3.9.First, note that when α∈R∗ +, then the parameters mn,µnand σ2 nin (3.5) verify lim n→∞ mn mn−1 =λα,lim n→∞µn= 0,lim n→∞σ2 n=σ2 α. In the sequel, we quantify the rates of convergence. Since σ2 nhas already been studied in Lemma 3.8 we shall focus only on µnand mn. 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 α⏐⏐⏐⏐≤Cvrn α, for an appropriate Cv∈R∗ +depending on αand σ2 0. Since λα∈(0,1) then d’Alembert’s ratio test readily shows that µnrelaxes to zero as a geometric sequence. Indeed, note that |µn|=(|µn| |µn−1|−λ2/d α)|µn−1|+λ2/d α|µn−1| ≤ (λ2/d α+Cvrn α)|µn−1|, for any n∈N. By an inductive argument, this yields |µn| ≤ n ∏ k=1 (λ2/d α+Cvrn α)|µ0|≲Cε(λ2/d α+ε)n,(3.17) for sufficiently large Cµ,ε ∈Rd, where we have used that rα∈(0,1) to absorb Cvrn αin an ε-small term. On the other hand, note that mn mn−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, ⏐⏐⏐⏐ mn mn−1−λα⏐⏐⏐⏐≤α|µn−1|2 2(1 + α)d 2+1 +dα 4|σ2 n−1−σ2 α| (1 + α)d 2+1 . 21 V. Calvez, T. Lepoutre and D. Poyato Nonlinear Analysis 238 (2024) 113392 Fig. 6. Comparison of convergence rates in Theorem 1.1 and Corollary 3.9. Therefore, using contraction of mean (3.17) and contraction of variance (3.7) in Lemma 3.8 entail the rate of convergence for the mn/mn−1. Finally, given that both Fn/∥Fn∥L1(Rd)and Fαare Gaussian functions, the relative entropy can be computed explicitly through (3.16). Specifically, DKL (Fn ∥Fn∥L1(Rd)Fα)=|µn|2 2σ2 α +d 2(σ2 n σ2 α−1)−d 2log (σ2 n σ2 α). Hence, using again (3.17) and (3.7) concludes our result. □ The proof of Theorem 1.1 for generic non-Gaussian initial data F0∈ M+(Rd) will be the core of this paper and we postpone it to Section 6. Remark 3.10 (Optimality of Convergence Rates).As illustrated in Fig. 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 α= (2kα)2holds. 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 L1norms, 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 ⏐⏐⏐⏐⏐∥Fn∥L1(Rd) ∥Fn−1∥L1(Rd)−λα⏐⏐⏐⏐⏐≲   √DKL (Fn ∥Fn∥L1(Rd)Fα), based on Pinsker’s inequality. However, a certain quadratic structure is still present in the rate of convergence of mn/mn−1in Corollary 3.9, which may become explicit for generic solutions to (1.1) if one finds the hidden cancellations. However, for simplicity we do not address this technical detail here. 4. Some properties of the operator T First, note that for any F∈ M+(Rd) we obtain that T[F]∈Wk,1(Rd)∩Wk,∞(Rd) for each k∈N thanks to the fact that Fhas finite mass and the Gaussian mixing kernel Gin (1.4) is smooth and has 22 V. Calvez, T. Lepoutre and D. Poyato Nonlinear Analysis 238 (2024) 113392 bounded and integrable derivatives of any order. In particular, solutions {Fn}n∈Nto the time-discrete problem (1.1) become instantaneously smooth after the first generation for any generic initial datum F0∈ M+(Rd). 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]∥L1(Rd). Both a priori estimates will be required later in Section 6. 4.1. Emergence of Gaussian tails Lemma 4.1 (Emergence of Tails).Assume that α∈R∗ +, consider F∈ M+(Rd)and set σ2, σ2∈R∗ +by σ2:= 1 α, σ2∈(0,1 1 + α). Then, the following properties are fulfilled: (i) (Upper control of tails) There exists C=C(α, F)>0such that T[F](x)≤C G0,σ2(x),(4.1) for every x∈Rd. (ii) (Lower control of tails) There exists c=c(α, σ2, F)>0such that T[F](x)≥c G0,σ2(x),(4.2) for every x∈Rd. (iii) (Control of log-derivative) Assume that Fis 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′G0,σ2for every x∈Rd. Then, there exists C′′ =C′′(α, σ2, σ2, C′, F)>0such that |∇log T[F](x)| ≤ C′′(1 + |x|),(4.3) for every x∈Rd. Proof. First, notice that F∈ M+(Rd) is any generic distribution but G∈L∞(Rd). Then, B[F]∈L∞(Rd) and, by definition of Tin (1.2), we obtain T[F](x)≤ ∥B[F]∥L∞(Rd)e−α 2|x|2=∥B[F]∥L∞(Rd)(2πα−1)d/2G0,σ2(x), for every x∈Rd. Hence, (4.1) holds for appropriate C > 0. Second, note that T[F](x) = e−α 2|x|2 (2π)d/2∥F∥M+(Rd)∫R2d exp (−1 2⏐⏐⏐⏐x−x1+x2 2⏐⏐⏐⏐ 2)F(dx1)F(dx2) =e−α 2|x|2 (2π)d/2∥F∥M+(Rd)∫R2d exp (−1 2|x|2−1 8|x1+x2|2+1 2x·(x1+x2))F(dx1)F(dx2) ≥e−α 2|x|2 (2π)d/2∥F∥M+(Rd)∫R2d exp (−2 + ε 4|x|2−2 + ε 8ε|x1+x2|2)F(dx1)F(dx2) ≥e−α 2|x|2 (2π)d/2∥F∥M+(Rd)∫R2d exp (−2 + ε 4|x|2−2 + ε 4ε|x1|2−2 + ε 4ε|x2|2)F(dx1)F(dx2), 23 V. Calvez, T. Lepoutre and D. Poyato Nonlinear Analysis 238 (2024) 113392 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∥G0,2ε 2+εF∥2 M+(Rd) ∥F∥M+(Rd) e−α 2|x|2e−2+ε 4|x|2, for every x∈Rd. 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∥L1(Rd)∫R2d(x−x1+x2 2)G(x−x1+x2 2)F(dx1)F(dx2) =−(1 + α)xT[F](x)−e−α 2|x|2 ∥F∥L1(Rd)∫R2d(x1+x2 2)G(x−x1+x2 2)F(dx1)F(dx2). Therefore, dividing by T[F](x) and taking norms yield the following estimate |∇log T[F](x)| ≤ (1 + α)|x|+1 √2∥F∥L1(Rd) e−α 2|x|2 T[F](x)∫R2d|(x1, x2)|G(x−x1+x2 2)F(dx1)F(dx2),(4.4) for every x∈Rd. Consider any R > 0 and split the integral in the right hand side of (4.4) into the subsets {(x1, x2)∈R2d:|(x1, x2)| ≤ √2R|x|} and {(x1, x2)∈R2d:|(x1, x2)|>√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 Galong with the assumptions on Fand the preceding lower bound (4.2) of T[F]. Then, we obtain |∇log T[F](x)| ≤ (1 + α+R)|x|+C2σde−α 2|x|2e 1 2σ2|x|2 √2c∥F∥L1(Rd)∫|(x1,x2)|>√2R|x||(x1, x2)|e−1 2σ2|(x1,x2)|2dx1dx2 = (1 + α+R)|x|+C2σde−α 2|x|2e 1 2σ2|x|2 √2c∥F∥L1(Rd)∫+∞ √2R|x| r2de−1 2σ2r2dr ≤(1 + α+R)|x|+C2σde−α 2|x|2e 1 2σ2|x|2 √2c∥F∥L1(Rd)(8d−4 e)d−1 2 σ2d−1∫+∞ √2R|x| re−1 4σ2r2dr = (1 + α+R)|x|+√2C2σdσ2d+1 c∥F∥L1(Rd)(8d−4 e)d−1 2 e−1 2(α+R2 σ2−1 σ2)|x|2 , for every x∈Rd, where in the third line we have used the inequality r2d−1≤(2d−1 e)d−1 2(σ ε)2d−1exp (ε 2σ2r2), 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 {Fn}n∈Nto the time-discrete problem (1.1) issued at a generic initial datum F0∈ M+(Rd). In particular, we note that Gaussian tails emerge instantaneously after the first generation n= 1 and linear growth of the log-derivative is guaranteed after the second generation n= 2. Corollary 4.2 (Propagation of Tails).Assume that α∈R∗ +, consider the solution {Fn}n∈Nof (1.1) issued at a generic initial datum F0∈ M+(Rd)and set σ2 1, σ2 1∈R∗ +by σ2 1:= 1 α, σ2 1∈(0,1 1 + α). 24 V. Calvez, T. Lepoutre and D. Poyato Nonlinear Analysis 238 (2024) 113392 Define the associated sequences of variances {σ2 n}n∈Nand {σ2 n}n∈Nby 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 Fn) There exist Cn=Cn(α, F0)>0such that Fn(x)≤CnG0,σ2 n(x),(4.6) for each x∈Rdand every n∈N. (ii) (Lower control of tails of Fn) There exist cn=cn(α, σ2 1, F0)>0such that Fn(x)≥cnG0,σ2 n(x),(4.7) for each x∈Rdand every n∈N. (iii) (Control of log-derivative of Fn) There exists Dn=Dn(α, σ2 1, F0)such that |∇log Fn(x)| ≤ Dn(1 + |x|),(4.8) for each x∈Rdand 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=Fn−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 Fn=T[Fn−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=F0yields (4.6) and (4.7) with n= 1. Let us assume that (4.6) and (4.7) hold for some n∈Nand let us prove it for n+ 1. Specifically, using the induction hypothesis we obtain that Fn+1(x) = e−α 2|x|2B[Fn](x) ≤C2 n ∥Fn∥L1(Rd) e−α 2|x|2(G(· 2)∗G0,σ2 n∗G0,σ2 n)(2x) =C2 n ∥Fn∥L1(Rd) e−α 2|x|2G0,1+ σ2 n 2 (x)≤Cn+1G0,σ2 n+1 (x), for every x∈Rdand appropriate Cn+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∈Nand {σ2 n}n∈Nrelax 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 {Fn}n∈Nto 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.,Fn=λnF. 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 {Fn}n∈Nin Theorem 1.1 and the uniqueness result in Corollary 1.2. Namely, the above coefficients Cnand cncontains 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 Cn+1 = (1 −ασ2 n)d/2C2 n ∥Fn∥L1(Rd) , cn+1 = (1 −ασ2 n)d/2c2 n ∥Fn∥L1(Rd) , but they require further knowledge about the behavior of ∥Fn∥L1(Rd). 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. 25 V. Calvez, T. Lepoutre and D. Poyato Nonlinear Analysis 238 (2024) 113392 Here, we remind that Fα=0 =G0,2is 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 emis 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 ¯ F0in Definition 5.1, we obtain that (emF2)(x) = 1 (2π)3d/2 1 (4π)2d 1 ∥F1∥L1(Rd)∥F0∥2 M(Rd) ×∫R6d exp (−α 2(|x1|2+|x2|2)−1 2⏐⏐⏐⏐x−x1+x2 2⏐⏐⏐⏐ 2) ×exp (−α0 2(|x11|2+|x12|2+|x21|2+|x22|2)−1 2⏐⏐⏐⏐x1−x11 +x12 2⏐⏐⏐⏐ 2 −1 2⏐⏐⏐⏐x2−x21 +x22 2⏐⏐⏐⏐ 2) ׯ F0(x11)¯ F0(x12)¯ F0(x21)¯ F0(x22)dx2. (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 ¯ F0of F0. 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 xias indices iin the perfect binary tree T2. Indeed, we have sorted the different terms in (5.4) in such a way that the second line only involves variables x1and x2in the first level L2 1of the tree, whilst the third line contains the terms involving the variables x11,x12,x21 and x22 in the second level of the tree L2 2(i.e., the leaves L2). We will divide the method into two steps. First, we will address the change of the variables indexed by the leaves L2. Second, we will perform the change of variables indexed by the first level of the tree L2 1. •Step 1: Change of variables for x11,x12,x21 and x22. Consider some coefficient k1>0 to be determined later and define the change of variables x11 →y11, x12 →y12,x21 →y21 and x22 →y22 given by x11 =k1x1+y11, x12 =k1x1+y12, x21 =k1x2+y21, x22 =k1x2+y12. On the one hand, using the change of variables in the terms inside the exponential of (5.4) which involve x11 and x12, we obtain α0 2|x11|2+α0 2|x12|2+1 2⏐⏐⏐⏐x1−x11 +x12 2⏐⏐⏐⏐ 2 =α0 2|k1x1+y11|2+α0 2|k1x1+y12|2+1 2⏐⏐⏐⏐(1 −k1)x1−y11 +y12 2⏐⏐⏐⏐ 2 =1 2(2α0k2 1+ (1 −k1)2)|x1|2+α0 2|y11|2+α0 2|y12|2+1 8|y11 +y12|2 +(α0k1−1−k1 2)x1·y11 +(α0k1−1−k1 2)x1·y12. 32 V. Calvez, T. Lepoutre and D. Poyato Nonlinear Analysis 238 (2024) 113392 On the other hand, using the change of variables in the terms which involve x21 and x22, we get α0 2|x21|2+α0 2|x22|2+1 2⏐⏐⏐⏐x2−x21 +x22 2⏐⏐⏐⏐ 2 =α0 2|k1x2+y21|2+α0 2|k1x2+y22|2+1 2⏐⏐⏐⏐(1 −k1)x2−y21 +y22 2⏐⏐⏐⏐ 2 =1 2(2α0k2 1+ (1 −k1)2)|x2|2+α0 2|y21|2+α0 2|y22|2+1 8|y21 +y22|2 +(α0k1−1−k1 2)x2·y21 +(α0k1−1−k1 2)x2·y22. Notice that one can eliminate the crossed terms in both expressions by choosing k1:= 1 1+2α0 =1 2(1 + α). In that case, adding both terms we obtain that α0 2|x11|2+α0 2|x12|2+1 2⏐⏐⏐⏐x1−x11 +x12 2⏐⏐⏐⏐ 2 +α0 2|x21|2+α0 2|x22|2+1 2⏐⏐⏐⏐x2−x21 +x22 2⏐⏐⏐⏐ 2 =1 2(1 −k1)|x1|2+1 4k1|y11|2+1 4k1|y12|2−1 8|y11 −y12|2 +1 2(1 −k1)|x2|2+1 4k1|y21|2+1 4k1|y22|2−1 8|y21 −y22|2. Putting everything together into (5.4) yields (emF2)(x) = 1 (2π)3d/2 1 (4π)2d 1 ∥F1∥L1(Rd)∥F0∥2 M(Rd) ×∫R6d exp (−α1 2(|x1|2+|x2|2)−1 2⏐⏐⏐⏐x−x1+x2 2⏐⏐⏐⏐ 2) ×exp (−1 4k1(|y11|2+|y12|2+|y21|2+|y22|2)+1 8(|y11 −y12|2+|y21 −y22|2)) ׯ F0(k1x1+y11)¯ F0(k1x1+y12)¯ F0(k1x2+y21)¯ F0(k1x2+y22)dx1dx2dy11 dy12 dy21 dy22, (5.5) where the parameter α1∈R∗ +has been defined by α1:= 1 −k1+α, in order to recombine the initial terms in the second line of (5.4) involving variables indexed by the first level L2 1with the new remainders that have appeared from the previous step. •Step 2: Change of variables for x1and x2. Consider some coefficient k2>0 to be determined later and define the change of variables x1→y1and x2→y2given by x1=k2x+y1, x2=k2x+y2. This time, using the change of variables in the terms of the exponential in the second line of (5.5) yields α1 2|x1|2+α1 2|x2|2+1 2⏐⏐⏐⏐x−x1+x2 2⏐⏐⏐⏐ 2 =α1 2|k2x+y1|2+α1 2|k2x+y2|2+1 2⏐⏐⏐⏐(1 −k2)x−y1+y2 2⏐⏐⏐⏐ 2 33 V. Calvez, T. Lepoutre and D. Poyato Nonlinear Analysis 238 (2024) 113392 =1 2(2α1k2 2+ (1 −k2)2)|x|2+α1 2|y1|2+α1 2|y2|2+1 8|y1+y2|2 +(α1k2−1−k2 2)x·y1+(α1k2−1−k2 2)x·y2. Again, we can cancel the crossed term by choosing k2:= 1 1+2α1 =1 3+2α−2k1 . Namely, we obtain that α1 2|x1|2+α1 2|x2|2+1 2⏐⏐⏐⏐x−x1+x2 2⏐⏐⏐⏐ 2 =1 2(1 −k2)|x|2+1 4k2|y1|2+1 4k2|y2|2−1 8|y1−y2|2. Then, putting everything together into (5.5) yields F2(x) = 1 (2π)3d/2 1 (4π)2d 1 ∥F1∥L1(Rd)∥F0∥2 M(Rd) e−1+α−k2 2|x|2 ×∫R6d exp (−1 4k2(|y1|2+|y2|2)+1 8|y1−y2|2) ×exp (−1 4k1(|y11|2+|y12|2+|y21|2+|y22|2)+1 8(|y11 −y12|2+|y21 −y22|2)) ׯ F0(k1k2x+k1y1+y11)¯ F0(k1k2x+k1y1+y12) ׯ F0(k1k2x+k1y2+y21)¯ F0(k1k2x+k1y2+y22)dy2, (5.6) where we denote again y2= (y1, y2, y11, y12, y21, y22)∈R6d. 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 {kn}n∈Nare defined by the recursive formula k1:= 1 2(1 + α), kn:= 1 3+2α−2kn−1 ,for n≥2. (5.7) (2) The coefficients {κn}n∈Nare defined by the recursive formula κ0:= 1, κn:= k1···kn,for n≥1.(5.8) We note that for n= 2 the above coefficients k1, k2reduce 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→ ∞. Lemma 5.3 (Asymptotic Behavior of the Coefficients).Consider the sequences {kn}n∈Nand {κn}n∈Nin Definition 5.2. Then, the following properties hold true: 34 V. Calvez, T. Lepoutre and D. Poyato Nonlinear Analysis 238 (2024) 113392 (i) (Coefficients kn) The coefficients {kn}n∈Nare positive numbers and kn↘kαas n→ ∞, where kα∈R∗ +is the smallest root of the equation 1 3+2α−2kα =kα.(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α= 2k2 α. In addition, kn−kα≤Crn−1 α,(5.10) for every n∈N, where rαis and C∈R∗ +depends only on α. (ii) (Coefficients κn) The coefficients {κn}n∈Nare 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∈Nreadily follow from those of {kn}n∈Nand the relation (5.8) in Definition 5.2. Specifically, for any given q∈Nand any n≥qwe get the decomposition κn=(q ∏ m=0 km)( n ∏ m=q+1 km)≤kq 1kn−q q+1 =kq 1 kq q+1 kn q+1 ≤kq 1 kq q+1 (Crq+1 α+kα)n. Therefore, taking qsufficiently large so that Crq+1 α≤εwe conclude (5.11). We then focus on the properties of {kn}n∈Nand we use a similar strategy like in Lemma 3.8. •Step 1: Monotone convergence of {kn}n∈N. Define the following function f(x) := 1 3+2α−2x, x ∈(−∞,3 2+α). Then, by (5.7) {kn}n∈Nobeys the following recursive relation kn=f(kn−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 k1∈(x−, x+), then we conclude that {kn}n∈Nis 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) = xi.e.,ℓ∈ {x−, x+}. Since the full sequence {kn}n∈Nis below x+and decreasing, we then conclude that ℓ=x−=kα. In particular, we obtain kα< kn≤k1,(5.12) for any n∈N. 35 V. Calvez, T. Lepoutre and D. Poyato Nonlinear Analysis 238 (2024) 113392 •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¨obius function M:R\{3 2+α} −→ Rgiven by M(x) := 1 3+2α−2x, x ∈R\{3 2+α}. By definition of {kn}n∈Nin (5.7) we obtain kn−x+=1 3+2α−2kn−1−1 3+2α−2x+ =2 (3 + 2α−2kn−1)(3 + 2α−2x+)(kn−1−x+). Similarly, we obtain kn−x−=1 3+2α−2kn−1−1 3+2α−2x− =2 (3 + 2α−2kn−1)(3 + 2α−2x−)(kn−1−x−). By dividing both expressions and iterating such an identity we get kn−x− kn−x+ =(3+2α−2x+ 3+2α−2x−)n−1k1−x− k1−x+ ,(5.13) for every n∈N. In fact, the basis can be restated as follows 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 kn−kα=k1−x− x+−k1 (x+−kn)rn−1 α≤k1−x− x+−k1 (x+−x−)rn−1 α, for any n∈N, thus ending the proof. □ Definition 5.4 (Quadratic Forms).Consider the sequence {kn}n∈Nin Definition 5.2. Then, we define the quadratic form Qn=Qn(yn) by Qn(yn) := n−1 ∑ m=0 ∑ i∈Ln m(1 4kn−m (|yi1|2+|yi2|2)−1 8|yi1−yi2|2),yn∈R2(2n−1)d,(5.14) where we are using the tree-indexed notation yn= (yi)i∈Tn ∗∈R2(2n−1)din Remark 3.1. Again, note that when n= 2 the above quadratic form Q2(y2) reduces to the one inside the exponential of (5.6). We now show the following uniform control of the quadratic forms Qnas n→ ∞. Lemma 5.5 (Uniform Positive Definite Quadratic Forms).Consider the quadratic form Qn=Qn(yn)in Definition 5.4 and define the couple of coefficients βmin, βmax ∈R∗ +by βmin := 1+2α 4, βmax := 1 4kα , where kαis given by formula (1.15). Then, we obtain that βmin∥yn∥2≤Qn(yn)≤βmax∥yn∥2, where we are using the tree-indexed notation yn= (yi)i∈Tn ∗∈R2(2n−1)din Remark 3.1. 36 V. Calvez, T. Lepoutre and D. Poyato Nonlinear Analysis 238 (2024) 113392 Fig. 7. Unique path joining the leaf j= 121 to the root in the perfect binary tree T3. The corresponding lineage map takes the form Φ121 n(x;y3) = κ3x+κ2y1+κ1y12 +y121. Proof. By virtue of Young’s inequality, we obtain the following lower and upper bound for Qn n−1 ∑ m=0 ∑ i∈Ln m 1 4(1 kn−m−1)(|yi1|2+|yi2|2)≤Qn(yn)≤ n−1 ∑ m=0 ∑ i∈Ln m 1 4 1 kn−m (|yi1|2+|yi2|2), for every yn= (yi)i∈Tn ∗∈R2(2n−1)d. Then, the result follows from Lemma 5.3, which guarantees the uniform control kα< km−n≤k1for every m= 0, . . . , n −1. □ Definition 5.6 (Lineage Maps).We define the lineage maps Φj n=Φj n(x;yn) associated to the leave j∈Ln (see Fig. 7) as follows Φj n(x;yn) := κnx+ n−1 ∑ m=0 κmycm(j),yn∈R2(2n−1)d,(5.15) where x∈Rdrepresents the root value, yn= (yi)i∈Tn ∗∈R2(2n−1)dis represented according to the treeindexed notation in Remark 3.1,{κn}n∈Nis given in Definition 5.2, and cm=c◦ ··· ◦ cis the mtimes iterated map c:Tn ∗→ˆ Tnwhich, to any vertex i∈Tn ∗, it associates its child c(i)∈ˆ Tn(cf. Section 3.1). 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 (emFn)(x) = 1 (2π)(2n−1)d/2 1 (4π)2n−1d 1 ∏n−1 m=0 ∥Fm∥2n−1−m L1(Rd) ×∫R2(2n−1)dexp ⎛ ⎝− n−2 ∑ m=0 ∑ i∈Ln m[α 2(|xi1|2+|xi2|2) + 1 2⏐⏐⏐⏐xi−xi1+xi2 2⏐⏐⏐⏐ 2]⎞ ⎠ ×exp ⎛ ⎝−∑ i∈Ln n−1[1 2(α+1 2)(|xi1|2+|xi2|2) + 1 2⏐⏐⏐⏐xi−xi1+xi2 2⏐⏐⏐⏐ 2]⎞ ⎠ ×∏ j∈Ln ¯ F0(xj)dxn, (5.16) for any x∈Rdand n∈N. Above, we have used the tree-indexed notation xn= (xi)i∈Tn ∗∈R2(2n−1)d in Remark 3.1, we have set x∅:= xand we have considered the rescaled distribution ¯ F0according 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 ¯ F0of F0. As a consequence, the factor (4π)−2n−1darises from the normalization by Fα=0(xj) whilst the factor (2π)−(2n−1)d/2comes from the repeated products of the Gaussian mixing kernel G. The main result then reads as follows 37 V. Calvez, T. Lepoutre and D. Poyato Nonlinear Analysis 238 (2024) 113392 Proposition 5.7 (Reformulation of the Iterations I).Assume that α∈R+and set any initial datum F0∈ M+(Rd). Hence, the solution {Fn}n∈Nto the time-discrete problem (1.1) admits the following form Fn(x) = e−1+α−kn 2|x|2 (4π)2n−1d(2π)(2n−1)d/2 1 ∏n−1 m=0 ∥Fm∥2n−1−m L1(Rd)∫R2(2n−1)de−Qn(yn)∏ j∈Ln ¯ F0(Φj n(x;yn)) dyn,(5.17) for every x∈Rdand n∈N, where we are using the tree-indexed notation yn= (yi)i∈Tn ∗∈R2(2n−1)din 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 xi1=kn−mxi+yi1, xi2=kn−mxi+yi2, for every i∈Ln m, which we will apply in a backwards way starting at indices i∈Ln n−1and 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⏐⏐xi−xi1+xi2 2⏐⏐2inside the exponential of (5.16) get appropriately centered and no crossed terms remain when the selection parts α 2(|xi1|2+|xi2|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 xi1→yi1and xi2→yi2for i∈Ln n−1. We consider the change of variables given by xi1=k1xi+yi1, xi2=k1xi+yi2, where k1is given by Definition 5.2. The collection of all the terms for i∈Ln n−1in (5.16) then read α0 2|xi1|2+α0 2|xi2|2+1 2⏐⏐⏐⏐xi−xi1+xi2 2⏐⏐⏐⏐ 2 =α0 2|k1xi+yi1|2+α0 2|k1xi+yi2|2+1 2⏐⏐⏐⏐(1 −k1)xi−yi1+yi2 2⏐⏐⏐⏐ 2 =1 2(2α0k2 1+ (1 −k1)2)|xi|2+α0 2|yi1|2+α0 2|yi2|2+1 8|yi1+yi2|2 +(α0k1−1−k1 2)xi·yi1+(α0k1−1−k1 2)xi·yi2 =1 2(2α0k2 1+ (1 −k1)2)|xi|2+α0 2|yi1|2+α0 2|yi2|2+1 8|yi1+yi2|2 =1 2(1 −k1)|xi|2+1 4k1|yi1|2+1 4k1|yi2|2−1 8|yi1−yi2|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 k1and α0in order to cancel the crossed terms, i.e., k1=1 1+2α0 . 38 V. Calvez, T. Lepoutre and D. Poyato Nonlinear Analysis 238 (2024) 113392 •Step 2: Change of variables xi1→yi1and xi2→yi2for i∈Ln n−2. Again, we consider the change of variables given by xi1=k2xi+yi1, xi2=k2xi+yi2, where k2is given by Definition 5.2. Putting together the terms for i∈Ln n−2in (5.16) and the above xidependent remainder in the above expression in Step 1 yields the following term under the above change of variables α1 2|xi1|2+α1 2|xi2|2+1 2⏐⏐⏐⏐xi−xi1+xi2 2⏐⏐⏐⏐ 2 =α1 2|k2xi+yi1|2+α1 2|k2xi+yi2|2+1 2⏐⏐⏐⏐(1 −k2)xi−yi1+yi2 2⏐⏐⏐⏐ 2 =1 2(2α1k2 2+ (1 −k2)2)|xi|2+α1 2|yi1|2+α1 2|yi2|2+1 8|yi1+yi2|2 +(α1k2−1−k2 2)xi·yi1+(α1k2−1−k2 2)xi·yi2 =1 2(2α1k2 2+ (1 −k2)2)|xi|2+α1 2|yi1|2+α1 2|yi2|2+1 8|yi1+yi2|2 =1 2(1 −k2)|xi|2+1 4k2|yi1|2+1 4k2|yi2|2−1 8|yi1−yi2|2, where we have defined again coefficient the α1∈R∗ +by α1:= 1 + α−k1, in order to absorb the above-mentioned xi-dependent remainder. In addition, note that we have used again the following relation between k2and α1to cancel the crossed term, i.e., k2=1 1+2α1 . •Step 3: Change of variables xi1→yi1and xi2→yi2for i∈Ln n−3. We consider the change of variables given by xi1=k3xi+yi1, xi2=k3xi+yi2, where k3is given by Definition 5.2. Putting together the terms for i∈Ln n−3in (5.16) and the above xidependent remainder in the above expression in Step 2 yields the following term under the above change of variables α2 2|xi1|2+α2 2|xi2|2+1 2⏐⏐⏐⏐xi−xi1+xi2 2⏐⏐⏐⏐ 2 =α2 2|k3xi+yi1|2+α2 2|k3xi+yi2|2+1 2⏐⏐⏐⏐(1 −k3)xi−yi1+yi2 2⏐⏐⏐⏐ 2 =1 2(2α2k2 3+ (1 −k3)2)|xi|2+α2 2|yi1|2+α2 2|yi2|2+1 8|yi1+yi2|2 +(α2k3−1−k3 2)xi·yi1+(α2k3−1−k3 2)xi·yi2 =1 2(2α2k2 3+ (1 −k3)2)|xi|2+α2 2|yi1|2+α2 2|yi2|2+1 8|yi1+yi2|2 =1 2(1 −k3)|xi|2+1 4k3|yi1|2+1 4k3|yi2|2−1 8|yi1−yi2|2, 39 V. Calvez, T. Lepoutre and D. Poyato Nonlinear Analysis 238 (2024) 113392 where we have defined again coefficient the α2∈R∗ + α2:= 1 + α−k2, in order to absorb the above-mentioned xi-dependent remainder. In addition, note that we have used again the following relation between k3and α2to cancel the crossed term, i.e., k3=1 1+2α2 . Following a similar recursive process, we readily identify the quadratic form Qnin Definition 5.4 in the exponential of the integrand in the final expression (5.17). •Step 4: Identifying the lineage maps Φj nand the exponential factor. On the one hand, the hierarchy of changes of variables immediately leads to xj=Φj n(x;yn), for any leaf j∈Ln, thanks to Definition 5.2 for the coefficients {κn}n∈Nand 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 Qn.□ 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∈Lnof 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∈Nand Σm n:= ⎛ ⎜ ⎜ ⎝ (2 −kn−m)kn−m 1−kn−m Id−k2 n−m 1−kn−m Id −k2 n−m 1−kn−m Id (2 −kn−m)kn−m 1−kn−m Id ⎞ ⎟ ⎟ ⎠,(5.18) for every m= 0, . . . , n −1. Then, we will define the random vector Yn= (Yi)i∈Tn ∗distributed according to the following multivariate normal distribution Gn(yn) := n−1 ∏ m=0 ∏ i∈Ln m 1 (2π)d√det(Σm n)exp (−1 2(yT i1, yT i2)(Σm n)−1(yi1 yi2)),yn∈R2(2n−1)d,(5.19) where again we are using the tree-indexed notation yn= (yi)i∈Tn ∗∈R2(2n−1)din Remark 3.1. Notice that for every index i∈Ln mwith m= 0, . . . , n −1 the pair (Yi1, Yi2) is normally distributed according to N(0,Σm n) and Yi1and Yi2are certainly correlated. Nevertheless, the vector (Yi1, Yi2) is independent on any other random vector Yjfor indices j=i1 and j=i2. We are now ready to introduce the following probabilistic reformulation of Proposition 5.7. 40 V. Calvez, T. Lepoutre and D. Poyato Nonlinear Analysis 238 (2024) 113392 Proposition 5.9 (Reformulation of the Iterations II).Assume that α∈R∗ +and set any initial datum F0∈ M+(Rd). Hence, the solution {Fn}n∈Nto the time-discrete problem (1.1) admits the following form Fn(x) = e−1+α−kn 2|x|2 (2π)d/222n−1d ⎛ ⎜ ⎜ ⎜ ⎜ ⎝ n−1 ∏ m=0 (4k2 n−m 1−kn−m)2m−1d ∥Fm∥2n−1−m L1(Rd) ⎞ ⎟ ⎟ ⎟ ⎟ ⎠ E⎡ ⎣∏ j∈Ln ¯ F0(Φj n(x;Yn))⎤ ⎦,(5.20) for every x∈Rdand n∈N, where we are using the high-dimensional random vector Yn= (Yi)i∈Tn ∗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 Fn(x) = e−1+α−kn 2|x|2 (4π)2n−1d(2π)(2n−1)d/2 1 ∏n−1 m=0 ∥Fm∥2n−1−m L1(Rd)∫R2(2n−1)de−Qn(yn)∏ j∈Ln ¯ F0(Φj n(x;yn)) dyn =(2π)(2n−1)d/2 (4π)2n−1de−1+α−kn 2|x|2⎛ ⎜ ⎜ ⎜ ⎜ ⎝ n−1 ∏ m=0 (4k2 n−m 1−kn−m)2m−1d ∥Fm∥2n−1−m L1(Rd) ⎞ ⎟ ⎟ ⎟ ⎟ ⎠∫R2(2n−1)dGn(yn)∏ j∈Ln ¯ F0(Φj n(x;yn)) dyn, for every x∈Rdand each n∈N, where in the last line we have used that the inverse and determinant of the covariance matrices Σm nin (5.18) take the form (Σm n)−1:= ⎛ ⎜ ⎝ 2−kn−m 4kn−m Id 1 4Id 1 4Id 2−kn−m 4kn−m Id ⎞ ⎟ ⎠,det(Σm n) = (4k2 n−m 1−kn−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∈Ln¯ F0(Φj n(x;Yn)). □ As illustrated in the overall map in Fig. 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 {Fn}n∈Nto the timediscrete 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 powerful ideas, which we introduce later in Lemma 6.3. Remark 5.10 (Non-negligible Correlations).Since Ynin Definition 5.8 is normally distributed and the lineage map Φj n(x;yn)inDefinition 5.6 are affine transformations of yn, then all (Φj n(x;Yn))j∈Lnare also normally distributed. Unfortunately, the components of (Φj n(x;Yn))j∈Lnare 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 Kij =E[(Φi n(x;Yn)−EΦi n(x;Yn)) ⊗(Φj n(x;Yn)−EΦj n(x;Yn))], 41 V. Calvez, T. Lepoutre and D. Poyato Nonlinear Analysis 238 (2024) 113392 Proof. Let us fix any ε∈R∗ +arbitrarily small so that 2kα+ε < 1. Again, this is possible because 2kα<1 for α∈R∗ +. For simplicity of notation, we shall define the following sequence of coefficients: cn:= Fn(0) ∥Fn∥L1(Rd) , n ∈N. Our goal then reduces to studying the asymptotic behavior of {cn}n∈Nand obtaining quantitative convergence rates. Note that by formula (5.20) in Proposition 5.9 for the recursion, we obtain 1 cn =∫Rde−1+α−kn 2|x|2En[¯ F0](x)dx En[¯ F0](0) =∫|x|≤Rn e−1+α−kn 2|x|2En[¯ F0](x) En[¯ F0](0) dx +1 cn∫|x|>Rn Fn(x) ∥Fn∥L1(Rd) dx, where we set the sequence of radii {Rn}n∈Nas follows Rn:= 1 (2kα)n 2m+1 , n ∈N,(6.21) for a fixed value m∈N, which we take large enough so that (2kα)2m 2m+1 ≤2kα+ε 2.(6.22) Whilst our choice of Rnis not justified at first glance, we claim that it has been taken as to minimize the decay rate on the error terms ˜ Enand Enbelow. By solving the implicit equation on cnwe infer cn=1−˜ En (2πσ2 α)d/2+En ,(6.23) for every n∈N, where the error terms ˜ Enand Entake the form ˜ En:= ∫|x|>Rn Fn(x) ∥Fn∥L1(Rd) dx, En:= ∫|x|≤Rn e−1+α−kn 2|x|2En[¯ F0](x) En[¯ F0](0) dx −(2πσ2 α)d/2, (6.24) Additionally, we can split the second error term as En=En,1−En,2+En,3with En,1:= (2π 1 + α−kn)d/2 −(2π 1 + α−kα)d/2 , En,2:= ∫|x|>Rn e−1+α−kn 2|x|2dx, En,3:= ∫|x|≤Rn e−1+α−kn 2|x|2(En[¯ F0](x) En[¯ 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 Enbecomes clear since En,1can be reformulated as En,1=∫Rd e−1+α−kn 2|x|2dx −(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∫Rd eθ|x|2Fn(x) ∥Fn∥L1(Rd) dx < ∞. 48 V. Calvez, T. Lepoutre and D. Poyato Nonlinear Analysis 238 (2024) 113392 By expanding the exponential in power series, one obtains in particular uniformly bounded moments of any order, and in particular, of order 2m, namely Mm:= sup n∈N∫Rd|x|2mFn(x) ∥Fn∥L1(Rd) dx ≤Eθm! θm<∞. Hence, we infer the following control on the first error term in (6.24) ˜ En≤1 R2m n∫Rd|x|2mFn(x) ∥Fn∥L1(Rd) dx ≲1 R2m n ,(6.26) for any n∈N. On the other hand, by Lemma 5.3 we have that En,1in (6.25) can be bounded by |En,1|≲kn−kα≲rn α.(6.27) By direct calculation we also obtain En,2=|Sd−1|∫+∞ Rn rd−1e−1+α−kn 2r2dr ≲Γ(d 2,1 + α−kn 2R2 n)≲exp (−1 + α−kn 4(2kα)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 Rn→ ∞ by our choice (6.21). Finally, we control the error term En,3in (6.25) under the addition assumption that ¯ F0satisfies 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 ¯ F0by an advance enough time step ¯ Fmso 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| ≤ Rn, we have ⏐⏐⏐⏐∫1 0 (∇log En)[ ¯ F0](θx)·x dθ⏐⏐⏐⏐ ≲(2kα+ε′)nRn+ (rα+ε′)nR2 n+e−(2+ε′)n ε′eCε′(2rα+ε′)nR2 nRn ≲(2kα (2kα)1 2m+1 +ε 2)n +(rα (2kα)2 2m+1 +ε 2)n +e−(2+ε′)n ε′exp [Cε′(2rα (2kα)2 2m+1 +ε 2)n]1 (2kα)n 2m+1 ≲((2kα)2m 2m+1 +ε 2)n≲(2kα+ε)n. In the second inequality above we have taken ε′small enough compared with ε. In the third inequality we have used the relation rα= 2k2 α, which guarantees that all the contributions in the third line are dominated by the first term ((2kα)2m 2m+1 +ε 2)n. Finally, in the last inequality we have used our choice of m∈Nlarge enough so that (6.22) holds. Therefore, the mean value theorem implies ⏐⏐⏐⏐En[¯ F0](x) En[¯ F0](0) −1⏐⏐⏐⏐=⏐⏐⏐⏐exp (∫1 0 (∇log En[¯ F0])(θx)·x dθ)−1⏐⏐⏐⏐≲(2kα+ε)n. Hence, En,3in (6.25) can be controlled by |En,3| ≤ ∫|x|≤Rn e−1+α−kn 2|x|2⏐⏐⏐⏐En[¯ F0](x) En[¯ F0](0) −1⏐⏐⏐⏐dx ≲(2kα+ε)n.(6.29) Putting (6.27),(6.28) and (6.29) together yields |En|≲(2kα+ε)n.(6.30) 49 V. Calvez, T. Lepoutre and D. Poyato Nonlinear Analysis 238 (2024) 113392 Thereby, using the mean value theorem on (6.23) along with the bounds (6.26) and (6.30) concludes that ⏐⏐⏐⏐cn−1 (2πσ2 α)d/2⏐⏐⏐⏐≲˜ En+|En|≲(2kα+ε)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 ¯ F2already satisfies (6.5)2thanks 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)≤Cme 1 2(α+1 2−1 σ2 m)|x|2 , for appropriate Cm∈R∗ +, with variances {σ2 m}m∈Ndefined 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)2is then satisfied by ¯ Fmif we can show that the prefactor α+1 2−1 σ2 min 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 α)+Cvrm α, for any m∈N. Since we chose a normalization ¯ Fnof the profiles by a Gaussian G0,σ2with 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 Fn(x) = −(1 + α−kn)x+∇log En[¯ F0](x), for any x∈Rdand every n∈N. Using the mean value theorem we then achieve log Fn(x) Fn(0) −log Fα(x) Fα(0) =∫1 0 (∇log Fn)(θx)·x dθ +1 2σ2 α|x|2 =−1 + α−kn 2|x|+∫1 0∇log En[¯ F0](θx)·x dθ +1 2σ2 α|x|2 =kn−kα 2|x|2+∫1 0∇log En[¯ 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 ¯ F0verifies 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=¯ F0and obtain that ∇log En[¯ F0] converges 50 V. Calvez, T. Lepoutre and D. Poyato Nonlinear Analysis 238 (2024) 113392 to zero uniformly over compact sets with explicit convergence rates. Putting it into (6.31) we obtain more explicitly sup |x|≤R⏐⏐⏐⏐log Fn(x) Fn(0) −log Fα(x) Fα(0) ⏐⏐⏐⏐≤Cε,R (2kα)n, for any R∈R∗ +, any n∈Nand a large enough Cε,R ∈R∗ +. This, together with the above control on the asymptotics of ∥Fn∥L1(Rd)/Fn(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. Proof of Corollary 1.2.Let (λ, F) be any solution of the eigenproblem (1.11). Then, the ansatz (1.10) defines a solution Fn=λnFof the time-discrete problem (1.1).ByCorollary 6.5 we obtain that lim n→∞ sup |x|≤R⏐⏐⏐⏐ λnF(x) λn−Fα(x)⏐⏐⏐⏐= 0, for any R∈R∗ +. Hence, F≡Fα, and thus λ=λα.□ 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 ∥Fn∥L1(Rn)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 F0∈ M+(Rd), the log-Lipschitz norms of the high-dimensional integral En[¯ F0] cannot be controlled uniformly due to additional exponentially growing terms. Proof of Theorem 1.1. •Step 1: Convergence of the profiles Fn/∥Fn∥L1(Rd). 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 DKL (Fn ∥Fn∥Fα)≤σ2 α 2∫Rd⏐⏐⏐⏐∇log (Fn(x) Fα(x))⏐⏐⏐⏐ 2Fn(x) ∥Fn∥L1(Rd) dx, see Corollary 5.7.2 and Section 9.3.1 in [24] for details. By the reformulation of Fnin formula (5.20) of Proposition 5.9 and using the notation in Definition 6.1 we have that ∇log (Fn(x) Fα(x))= (kn−kα)x+∇log En[¯ F0](x). 51 V. Calvez, T. Lepoutre and D. Poyato Nonlinear Analysis 238 (2024) 113392 Therefore, we obtain the following upper bound DKL (Fn ∥Fn∥Fα)≲Dn,1+Dn,2,(6.32) where each factor takes the form Dn,1:= |kn−kα|2∫Rd|x|2Fn(x) ∥Fn∥L1(Rd) dx, Dn,2:= ∫Rd⏐⏐∇log En[¯ F0](x)⏐⏐2Fn(x) ∥Fn∥L1(Rd) dx. On the one hand, for Dn,1we can use the propagation of quadratic moments in Corollary 4.6 together with the explicit relaxation rates of {kn}n∈Nto kαin Lemma 5.3 to show that Dn,1≲r2n α.(6.33) On the other hand, for Dn,2we shall refine the argument in the proof of Corollary 6.5. Specifically, recall again that we can assume that ¯ F0satisfies the hypothesis (6.5) without loss of generality (cf. proof of Corollary 6.5). Then, applying Lemma 6.3 in order to control ∇log En[¯ F0] with fixed ε∈R∗ +yields Dn,2≲((2kα)2+ε)n∫Rd(1 + |x|2+e2Cε(2rα+ε)n|x|2)Fn(x) ∥Fn∥L1(Rd) dx, for some sufficiently large Cε∈R∗ +. Above we have used that (2kα)nis the decay rate that controls all the others in formula (6.6). Take n≥nεsufficiently large such that θε:= sup n≥nε 2Cε(2rα+ε)n<α 2. Then, Corollaries 4.6 and 4.8 imply that the above quadratic and exponential moments are uniformly bounded with respect to nfor n≥nε. Therefore, we have Dn,2≲((2kα)2+ε)n.(6.34) Finally, putting the estimates (6.33) and (6.34) into the split (6.32), and noticing that r2 α<(2kα)2thanks to the relation rα= 2k2 α, ends this part of the proof. •Step 2: Convergence of the growth rates ∥Fn∥L1(Rd)/∥Fn−1∥L1(Rd). For simplicity of notation we shall define the following sequence of coefficients: λn:= ∥Fn∥L1(Rd) ∥Fn−1∥L1(Rd) , n ∈N. Our goal then reduces to studying the asymptotic behavior of {λn}n∈Nand obtaining quantitative convergence rates. By definition of operator Tin (1.2) we obtain that λn=∫Rd e−m(x)∫R2d G(x−x1+x2 2)Fn−1(x1) ∥Fn−1∥L1(Rd) Fn−1(x2) ∥Fn−1∥L1(Rd) dx1dx2dx. By direct computation of the integral with respect to xas it was done in (4.11) we find that λn=∫R2d H(x1, x2)Fn−1(x1) ∥Fn−1∥L1(Rd) Fn−1(x2) ∥Fn−1∥L1(Rd) dx1dx2,(6.35) where the function H(x1, x2) takes the form H(x1, x2) := 1 (1 + α)d/2exp (−α (1 + α)⏐⏐⏐⏐ x1+x2 2⏐⏐⏐⏐ 2),(x1, x2)∈Rd.(6.36) 52 V. Calvez, T. Lepoutre and D. Poyato Nonlinear Analysis 238 (2024) 113392 In addition, recalling the explicit form of the eigen-pair (λα,Fα) in (1.12) implies λα=∫R2d H(x1, x2)Fα(x1)Fα(x2)dx1dx2.(6.37) Therefore, taking the difference of (6.35) and (6.37) and noticing that Hin (6.36) is bounded yields |λn−λα|≲ Fn ∥Fn∥L1(Rd)⊗Fn ∥Fn∥L1(Rd)−Fα⊗FαL1(R2d) ,(6.38) for any n∈N, where P⊗Qdenotes the Kronecker product of two measures P, Q ∈ M(Rd), namely ∫R2d φ(x, y) (P⊗Q)(dx, dy) = ∫Rd(∫Rd φ(x, y)P(dx))Q(dy),(6.39) for all φ∈Cc(Rd) (cf. Table 1). We do not have a direct convergence result of Fn/∥Fn∥L1(Rd)in L1norms, 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  Fn ∥Fn∥L1(Rd)⊗Fn ∥Fn∥L1(Rd)−Fα⊗FαL1(R2d) ≤1 √2   √DKL (Fn ∥Fn∥L1(Rd)⊗Fn ∥Fn∥L1(Rd)Fα⊗Fα) =1 √2   √2DKL (Fn ∥Fn∥L1(Rd)⊗Fα), (6.40) where in last step we have used the tensorization property of the Kullback–Leibler divergence. The result then follows from (6.38)–(6.40) and the above explicit convergence rates of the normalized profiles Fn/∥Fn∥L1(Rd)in Step 1.□ 7. Numerical experiments For our numerical simulation, we have restricted to one-dimensional traits (i.e.,d= 1) and we have considered a step function of the following form as initial datum: F0=1 Z(30 1[−7,−3] + 20 1[7.5,12.5] + 50 1[30,40] + 30 1[52,5,57,5]),(7.1) where Z∈R∗ +is a normalizing factor so that F0∈ P(R). Let {Fn}n∈Nbe the solution to the time-discrete problem (1.1) starting at F0. 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 F0and all the mass of each Fn(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 xaccording 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 Fnin both regimes and we illustrate numerically the convergence of the growth rates ∥Fn∥L1(Rd)/∥Fn−1∥L1(Rd)towards the eigenvalue λα, and the relaxation of the normalized profiles Fn/∥Fn∥L1(Rd)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 Fig. 6). Indeed, we actually attain numerically the same convergence rates as in Corollary 3.9, which were sharp for Gaussian initial data. 53 V. Calvez, T. Lepoutre and D. Poyato Nonlinear Analysis 238 (2024) 113392 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: λα≈0.9857,σ2 α≈1.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 Fig. 4). In Fig. 8 we observe the relaxation of the normalized profiles Fn/∥Fn∥L1(Rd)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 Fn/∥Fn∥L1(Rd)are given by ≈0.0474 and ≈1.8897 respectively, where the latter agrees with the above exact value σ2 αup to 5 digits. In Fig. 9 we represent the convergence of the growth rates ∥Fn∥L1(Rd)/∥Fn−1∥L1(Rd)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 Fig. 10 we have represented the errors εprof n:= DKL (Fn ∥Fn∥L1(Rd)Fα), εmass n:= ⏐⏐⏐⏐⏐∥Fn∥L1(Rd) ∥Fn−1∥L1(Rd)−λα⏐⏐⏐⏐⏐, 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 Fig. 10(a) 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 (2kα)2obtained in Theorem 1.1. These results are in perfect agreement with the other convergence results in Fig. 10(b): 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 2kα≈0.9717. Moreover, the numerical rate of convergence of the variance of the normalized profiles Fn/∥Fn∥L1(Rd)is much faster: it is approximately ≈0.4721 <0.5, again close to 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. 7.2. 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 Fig. 11(a) we note that the solution resembles a Gaussian distribution after a couple of iterations. After 15 iterations we obtain 54 V. Calvez, T. Lepoutre and D. Poyato Nonlinear Analysis 238 (2024) 113392 Fig. 8. Relaxation of the normalized profiles Fn/∥Fn∥L1(Rd)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α. that the normalized profile Fn/∥Fn∥L1(Rd)has approximate mean and variance respectively given by ≈ 0.0017 and ≈0.9221. In Fig. 11(b) we observe the convergence of the growth rates ∥Fn∥L1(Rd)/∥Fn−1∥L1(Rd) 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 Fig. 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 (2kα)2obtained in Theorem 1.1. Similarly, the numerical rate of convergence of the growth 55 V. Calvez, T. Lepoutre and D. Poyato Nonlinear Analysis 238 (2024) 113392 Fig. 9. Relaxation of the growth rates ∥Fn∥L1(Rd)/∥Fn−1∥L1(Rd)towards the eigenvalue λαalong time iterations 0≤n≤150 for a step function (7.1) as initial datum and weak selection parameter α= 0.015. Fig. 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:= DKL (Fn/∥Fn∥L1Fα), and (b) semi-log plot of the errors εmass n:= |∥Fn∥L1/∥Fn−1∥L1−λα|. Fig. 11. (a) Zoom near the origin of the relaxation of the normalized profiles Fn/∥Fn∥L1(Rd)towards the eigenfunction Fαalong the time iterations n= 0,1,2and 15 for the step function (7.1) as initial datum and strong selection parameter α= 0.4. The vertical dotted line represents the location of the mean of the equilibrium profile Fα. (b) Relaxation of the growth rates ∥Fn∥L1(Rd)/∥Fn−1∥L1(Rd) towards the eigenvalue λα. 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 2kα≈0.6311 in Theorem 1.1. 56 V. Calvez, T. Lepoutre and D. Poyato Nonlinear Analysis 238 (2024) 113392 8. Conclusions and perspectives In this paper, we have proven asynchronous exponential growth in a quantitative genetics model for the evolution of the distribution of traits in a population governed by sexual reproduction and multiplicative effect of selection. Our model assumes time-discrete non-overlapping generations, which rule out an eventual mixing with previous generations of ancestors. In addition, our non-linear sexual reproduction operator is set in agreement with Fisher’s infinitesimal model, and we have chosen selection to act on 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 [7], 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 [3] 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 [3], 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. Acknowledgments 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. Appendix. 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 Gand 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 57