Mixing Rates of the Geometrical Neutral Lorenz Model

The aim of this paper is to obtain polynomial decay of correlations of a Lorenz-like flow where the hyperbolic saddle at the origin is replaced by a neutral saddle. To do that, we take the construction of the geometrical Lorenz flow and proceed by changing the nature of the saddle fixed point at the origin by a neutral fixed point. This modification is accomplished by changing the linearised vector field in a neighbourhood of the origin for a neutral vector field. This change in the nature of the fixed point will produce polynomial tails for the Dulac times, and combined with methods of Araújo and Melbourne (used to prove exponential mixing for the classical Lorenz flow) this will ultimately lead to polynomial upper bounds of the decay of correlations for the modified flow.


Introduction
The study of flows on surfaces and higher-dimensional manifolds has caught the interest of many scientists because of its numerous applications such as Hamiltonian flows, geodesic and horocycle flows, billiard flows or flows from meteorological models. These flows are usually equipped with a natural invariant measure µ, for instance the SRB-measure.
The main goal is to have a better understanding of the properties of these flows, such as hyperbolicity, ergodicity, mixing (or at least weak mixing) and, in chaotic settings, rates of mixing; that is, we would like to investigate the asymptotic behaviour of the correlation coefficients where f t : M → M is a flow acting on a manifold M and µ its SRB-measure, and for observables v, w chosen from an appropriate Banach space. Knowing the rates of mixing is very helpful for proving other ergodic and statistical properties since mixing is one of the strongest statistical properties.
Obtaining good mixing rates for flows, even for hyperbolic flows, is far more difficult than for maps. Some seminal ideas were provided by Liverani [28] and Dolgopyat [19,20], with applications of these methods in e.g. [9,10,16]. To obtain sharp estimates in the polynomial setting, the operator renewal theory techniques developed by Sarig [34] and Gouëzel [23] are the only ones available.
The model we would like to study is probably one of the most emblematic ones, the Lorenz flow. In the mid seventies Afraȋmovič, Bykov and Shilnikov [1] and independently Guckenheimer and Williams [24] introduced the geometric Lorenz attractor to model the original Lorenz attractor. Our research focuses on a modified version of this geometrical model and study its rate of mixing, based on the precise estimates in [17] of Dulac times associated to a neutral saddle.
Recently, Araújo and Melbourne in [5] proved that the geometrical Lorenz flow (and hence the classical Lorenz flow), also enjoys exponential mixing. It is techniques from their papers, specifically C 1+α smoothness of the stable foliation, that leads eventually to the claimed mixing rates.

The framework
The geometrical Lorenz flow can be seen as the natural extension of a suspension semiflow built over a certain type of one-dimensional expanding map f Lor . We first consider the cross-section Σ transversal to the flow and the Poincaré map P Lor : Σ → Σ, which is decomposed in two parts. The first one is the Dulac map, denoted by P 1 , deals with the local behaviour near the origin and is obtained by considering a linear system in a neighbourhood of the origin; that is, we take the flow X t obtained from the linear systeṁ where λ u , λ s and λ ss denote the unstable, stable and strong stable eigenvalues of the original Lorenz system, respectively. Then we let points in Σ flow under X t until flow time τ ′ Lor := min{t > 0 : X t ∈ S} = −λ −1 u ln(|x|) + O(ln(|x|)) as x → 0. Thus we have that X τ ′ Lor = P 1 : Σ → S, where S ± is the image of Σ ± under P 1 and has a cusp-like shape, see Figure 1.
The second part, denoted by P 2 , consists of the return of S to Σ and simulates the random turns of a regular orbit around the origin and describes a butterfly-like shape. This is done by a composition of a rotation, expansion and translation with hitting time τ 2 (x) ∈ C ϵ . Thus, the full return time of the Poincaré map P Lor = P 2 • P 1 is given by As we will see later, the lines in the y-direction (i.e.,, parallel to the y axis) in Σ form the stable foliation, invariant under P Lor ; that is, for any leaf γ of this foliation, its image P Lor (γ) is contained in a leaf of the same foliation. By quotienting out the stable direction we can rewrite the Poincaré map as a skew-product; that is, P Lor (x, y) = (f Lor (x), g Lor (x, y)).
The geometric Lorenz flow is the couple (W, X t W ), where W = {X t (x) |x ∈ Σ, t ∈ R + }. Consider U = x∈Σ X [0,rLor(x)] W (x), then the geometric Lorenz attractor (of the corresponding vector field) is given by Λ Lor = t>0 X t W (U ). In [5], exponential mixing for the geometrical Lorenz flow was proven under two conditions: the stable foliation has to be C 1+α and a uniform non-integrability (UNI) condition needs to be satisfied.
The modified version is obtained by changing the local behaviour near the origin. We achieve this by replacing the linear system for the following system; where a 0 , a 1 , a 2 , b 0 , b 1 , b 2 , c 0 , c 2 and ℓ > 0, a 2 b 0 < 9a 0 b 2 , ∆ := a 2 b 0 − a 0 b 2 ̸ = 0 and O(4) refers to terms of order four or higher, under the condition that they are of the form x 2 O(2) near the yz-plane and z 2 O(2) near the xy-plane. This system has a polynomial Dulac time (see (8) and Figure 3) given by; as x → 0 and β 2 = a2+b2 2b2 . To obtain the flow time τ ′ N eu , we make use of the estimates of the Dulac map and the tails of the return map obtained by Bruin and Terhesiu in [17]. This change of flow time, from logarithmic to polynomial, will ultimately allow us to deduce the polynomial decay of correlations.
We denote by N t the flow obtained from the system given by (4). This change in the local behaviour near the origin leads to a change on the map P 1 ; that is, we have now N τ ′ Neu = D 1 : Σ → S. For the second part, the return of S to Σ, we consider the same diffeomorphism P 2 with same hitting time. In this way, we obtained the modified Poincaré map P Neu = P 2 • d 1 with return time given by, Similarly, we define the geometric neutral Lorenz flow as the couple (W, N t W ), where W = {N t (x) |x ∈ Σ, t ∈ R + }. We consider again U = x∈Σ N [0,rNeu(x)] (x), the geometric neutral Lorenz attractor (of the corresponding vector field) is given by Λ Neu = t>0 N t (U ).
As we will see in more detail in Section 2, the geometrical neutral Lorenz flow will be split into three models. Model 1 is obtained when we take the parameters c 0 = c 2 = 0 in (4). Model 2 when we consider a 1 = b 1 = 0. Finally, Model 3, the most general, will be given by taking all parameters strictly positive.

Main results
Let C η be the space of functions that are η-Hölder in the space direction, and C m,η be the space of functions that are m + η-Hölder (i.e., m time differentiable with an η-Hölder m-th derivative) in the flow direction, see Section 4 for the precise definitions. The main result in this paper is the following theorem: 1.1 Theorem. Let N t : Λ Neu → Λ Neu be the geometrical neutral Lorenz flow for Model 1 and Model 2 obtained from the neutral form given by (4), with corresponding parameters. Λ Neu its attractor and its SRB measure µ. Then N t has polynomial decay of correlations (with exponent β 2 = a2+b2 2b2 ); that is, there exist m ≥ 1 and a constant C > 0 such that for observables v ∈ C η (M ) ∩ C 0,η (M ), w ∈ C m,η (M ), and time t > 1 we have A first question that presents itself is of course if these bounds are sharp, and if current operator renewal theory methods [23,34] cannot prove that. We say more on this at the end of Section 4.
For the proof of Theorem 1.1, we obtain an explicit form of the Poincaré map, since we can solve the differential equation in the y component. Thus we are able to prove polynomial decay of correlations by using the results on non-uniformly hyperbolic flows established by Bálint et al. in [11].
For the third model the situation is more subtle since, to our knowledge, finding the solution of the differential equation in the y component is next to impossible. To overcome this problem we will analyse and compare, with numerical methods, the limit behaviour of the Dulac maps obtained in [14] and [17] and adapted to our framework. More precisely, we will analyse the limit behaviour of the maps D 1 : Σ → S obtained for each Neutral model. This is sufficient since the Poincaré maps considered in this work are given by P Neu = P 2 • D 1 , where P 2 is a diffeomorphism and the map D 1 is the Dulac map from the cross-section Σ to the cusps S, which depends on the differential equation being considered. Therefore, the behavioural changes exhibited by the map P Neu are represented by the changes of the map D 1 . Dulac in [21] made a significant contribution to solving Hilbert's 16th problem by incorporating his map as an element to establish that polynomial vector fields in the plane possess a limited number of limit cycles, demonstrating that they cannot have an infinite number of such cycles.
The numerical analysis on the behaviour of the Dulac maps will give us the plausibility of the following conjecture.
1.2 Conjecture. Let N t : Λ Neu → Λ Neu be the geometrical neutral Lorenz flow for Model 3 obtained from the neutral form given by (4), with the corresponding parameters. Λ Neu its attractor and its SRB measure µ. Then N t has polynomial decay of correlations (with exponent β 2 = a2+b2 2b2 ). The organization of this paper is as follows: In Section 2 we will give the construction of the Poincaré maps of the Neutral Model 1 and 2. In Section 3 we will be devoted to the proof that the stable foliation for the geometrical neutral models is C 1+α and the UNI condition is satisfied by adapting the existing proofs for the geometrical Lorenz model. Section 4 contains the framework of non-uniformly hyperbolic flows and the proof of Theorem 1.1. Finally, in Section 5 we will present the numerical analysis and results we obtained for the Dulac map and the tails of the return map.
For the remaining of this paper we will adopt the following notation.
Data availability statement: The methods of how the numerical graphics were computed are given in the last section of this paper. For further details, please contact the authors.

Neutral model 1
To create the modified models we will apply local surgery in a neighbourhood of the hyperbolic saddle equilibrium of the geometrical Lorenz model, namely the origin, and transform it into a neutral equilibrium. We do this because we aim to slow down the orbit and thus increase the time that orbits take to flow from the cross-section Σ to the cusps S, see Figure 1, and see the changes this new motion produces in the decay of correlations. The flow obtained from this modification will be an almost Anosov flow [26,27]. Existence of a finite or infinite SRB measure for two-dimensional almost Anosov diffeomorphisms was already proven in [26,27]. Bruin and Terhesiu in [17] proved mixing rates in the infinite SRB measure setting for almost Anosov diffeomorphism and established the required spectral properties for the transfer operator (acting on an appropriate anisotropic Banach space of distributions) of an induced map so as to obtain optimal rates of mixing. Furthermore, they gave more precise tail estimates for the inducing scheme. We will take advantage of these methods and estimates and use them to deduce the rates of mixing of our almost Anosov flow.
We consider now Figure 2. The section Σ is transversal to the flow and every trajectory eventually crosses Σ in the direction of the negative axis z. Then for each (x, y, 1) ∈ Σ, the time τ ′ Neu such that N τ ′ Neu (x, y, 1) ∈ S is determined by the estimates of the Dulac map provided in [17], as we will explain now. Let us start with the Neutral model 1; that is, we consider a neighbourhood U of the fixed point0 = (0, 0, 0) where the vector field has the form (4) (in local Euclidean coordinates) with c 0 = c 2 = 0 and the other parameters satisfying the constraints given before. Note that U is taken much smaller than the scale of Figure 2. This vector field, denoted by Z, is cubic at0 in the direction transversal to the stable manifold of0, but this is the only source of nonhyperbolicity. The y-axis is invariant and all solutions tend to0. The divergence is given by Since we want the flow to shrink volume exponentially fast, as does the Lorenz flow, we need Div(Z) ≤ −c < 0. Therefore, we let ℓ be large enough such that The solution for the y-component is given by y(t) = y 0 e −ℓt . Thus we obtain a non-autonomous system of differential equations, since the contribution of the y-component to the x and z-component is exponentially small as time increases, these terms are of smaller order than the higher order terms. Thus we are left with the two-dimensional system studied in [17]: Now let W s and W u be two mutually tranversal foliations of the interior of U ∩ {positive quadrant} that is invariant under the flow of (7) and such that • the leaves of W s accumulate in C 1 topology on the stable manifold of (0, 0) and are transversal to the unstable manifold of (0, 0), and • the leaves of W u accumulate in C 1 topology on the unstable manifold of (0, 0) and are transversal to the stable manifold of (0, 0).
One would like to use the stable and unstable foliation of the flow for W s and W u , but as long as we defined the flow only locally, the above properties suffice.
We make the following observations: has the shape of a cusp at (±1, 0, 0) and (with some abuse of notation) we will denote these images as S ± and S = S + ∪ S − .
where c is a constant, the line segments in Σ parallel to the y-axis and by ℓ h (c) = {(±1, y, z) ∈ S | z = c}, the line segments in S parallel to the y-axis. Then D 1 (ℓ v (c 0 )) = ℓ h (c 1 ); that is, the map D 1 takes the y-direction lines in Σ to the horizontal line segments in S as illustrated in Figure 2.
The return of the cusps S to the cross-section Σ is described by the map P 2 = T • E a • R θ , where R θ is a rotation by an angle of θ = 3π 2 and the rotation axis are the boundaries of the cross-section Σ parallel to the y-axis, E a is an expansion by a factor of a > 1 in the x-direction and a translation T such that the unstable direction which starts from the origin is sent to the boundary of Σ; that is, we want to send the cusp points C ± to A ± , see Figure 1. Thus, the full Poincaré map P Neu = P 2 • D 1 : Σ → Σ, shown in Figure 4, is given by, In the positive quadrant the matrix DP Neu has eigenvalues λ 1 = a β0 β2 |x| β 0 β 2 −1 and λ 2 = e −ℓA2(x,β2) . By restricting 1 2 < β 0 < 2 we have that 0 < β 0 β 2 < 1 since β 2 > 2. Then λ 1 → ∞ as x → 0. For the other eigenvalue we have that λ 2 < 1 and λ 2 → 0 as x → 0 since A 2 (x, β 2 ) → ∞ as x → 0. Thus the modified Poincaré map P Neu is hyperbolic when x approaches the origin. The foliation given by the lines ℓ v (c) is invariant under the map P Neu ; that is, given any leaf γ of this foliation its image P Neu (γ) is contained in a leaf of the same foliation (See Figure 4). Therefore, we can express P Neu as P Neu (x, y) = (f Neu (x), g Neu (x, y)), where f Neu : I\{0} → I is a Lorenz-like map with exponent β = β0 β2 with β 0 = a0+b0 2a0 and β 2 = a2+b2 2b2 ; that is, f Neu is given by, g1.-The map g Neu is piecewise C 2 and for fixed x 0 , the map g Neu (x 0 , y) is a contraction in the y-direction, i.e., where d is the Euclidean distance in I and 0 < c < 1.
g2.-DP Neu has the following bound on its partial derivatives: a) For all (x, y) ∈ Σ we have ∂ y g Neu (x, y) = e −ℓA2(x,β2) . Since β 2 > 2 and |x| ≤ 1, there is 0 < η < 1 such that and |y|, |x| ≤ 1 we get that |∂ x g Neu (x, y)| is bounded. In fact, it tends to zero exponentially fast as x approaches the origin.
g3.-From g2.a) above follows the uniform contraction of the foliation given by the lines ℓ v (c); in other words, there is C > 0 such that, for any given leaf γ of the foliation and for y 1 , y 2 ∈ γ, we have d(P n Neu (y 1 ), P n Neu (y 2 )) ≤ Cη n d(y 1 , y 2 ), when n → ∞.

Neutral model 2
Now we will consider the Neutral model 2; that is, we consider the same neighbourhood U of the origin where the flow has the local form given by (4) with a 1 = b 1 = 0 and the remaining parameters satisfying the same constraints stated in this framework. Again0 = (0, 0, 0) is the only neutral periodic orbit and the vector field is cubic in the direction transversal the stable manifold of0, but this is the only source of non-hyperbolicity. If x = 0 and z = 0, then we see thaṫ x = 0,ẏ = −ℓy andż = 0. Hence, the y-axis is invariant and all solutions tend to the origin as in the previous model. Moreover, sinceẋ andż are decoupled from y, we have (7). Thus the asymptotics for the Dulac map and the flow time given in [17] follow. Also Div Since we want a flow that shrinks volume exponentially fast as before, we take ℓ large enough so that (3a0−b0) We consider the same cross-section Σ as before and proceed to construct the Poincaré map in the same way. Denote by N t the flow obtained from (4) with the pertinent constraints in the parameters; that is, N t (x, y, z) = (x(t), y(t), z(t)). By (8) we obtain the following form for the flow, N t(x,z) (x, y, z) = (x 0 , y(t(x, z)), ω(z, t(x, z)).
Note thatẏ = y(−ℓ(1 + c 0 x 2 + c 2 z 2 )), applying Grönwall's Lemma we get, By the estimates of the Dulac map and since z ∈ Σ, we obtain that the time t(x, z) becomes a function of the variable x and the integral t 0 (c 0 x 2 + c 2 z 2 )ds can be expressed as a function q of the variable x. Observe that q(x) > 0 for every x. Therefore we get that, Hence y(t) decreases exponentially fast as before but with a faster rate. All together, we get that the map D 1 : Σ → S is given by where β = β0 β2 , compare this with (10). The statements form Observation 2.2 for this new version of the map D 1 will also hold. To finish the construction of the Poincaré map we compose now with the map P 2 . Therefore, the full return map P Neu : Σ → Σ of Σ is given by where β = β0 β2 ∈ (0, 1). The matrix DP Neu has eigenvalues λ 1 = aβ|x| β−1 and λ 2 = e −ℓ(A2(x,β2)+q(x)) . Since β ∈ (0, 1) we have that λ 1 → ∞ as x → 0. For the other eigenvalue we have that λ 2 < 1 and λ 2 → 0 as x → 0 since (A 2 (x, β 2 ) + q(x)) → ∞ as x → 0. Thus the modified Poincaré map P Neu is hyperbolic when x approaches the origin.
The properties stated before remain true for this new modified return map P Neu like the invariance of the stable foliation given by the vertical lines ℓ v (c) under the map P Neu . Hence, we can express again P Neu as P Neu (x, y) = (f Neu (x), g Neu (x, y)), where f Neu : I \ {0} → I is again a Lorenz-like map with exponent β (see (12)) and g Neu : (I \ {0}) × I → I satisfy the same properties stated in the previous section. The existence of a unique a.c.i.p and statistical properties such as exponential decay of correlations for observables with bounded variation for the Lorenz-like map f Neu are well established, see for example [35].

Existence and regularity of the strong stable foliation
In this subsection we will study the properties of the strong stable foliation F ss for the neutral geometrical Lorenz model we built in Section 2.
For the neutral geometrical Lorenz attractor, denoted by Λ Neu , we consider the Lorenz attractor Λ Lor in an open neighbourhood U of the origin. Instead of considering the linearised vector field we consider the vector field given by (4) with the parameters corresponding for model 1 and 2 described in Section 2.1 and 2.2, respectively. More precisely, we take an open neighbourhood U in which the cross-section Σ is contained. Then the Dulac map from Σ to S has the form given by (10) and (16) for the models 1 and 2, respectively. This modification yields a different flow time from the cross-section Σ to S. In the original Lorenz construction we have a logarithmic Poincaré return time but for these modifications we have a polynomial Poincaré return time given by (5). The rest of the construction, however, remains unchanged; that is, the flow constructed from S to Σ is made by a composition of an expansion, a rotation and a translation. Therefore we have the same hitting time τ 2 (x) and thus the full return time for the modified Poincaré map P Neu is given in (6). The modified Poincaré map P Neu : Σ → Σ is given by (11) and (17) for the model 1 and 2, respectively. We saw that the lines in the y-direction, denoted by ℓ v (c), in the cross-section Σ form the stable foliation which is preserved by the return map P Neu . Thus by quotienting out the stable direction we can rewrite the Poincaré map as a skew-product; that is, In finding the eigenvalues, we get λ ss = −ℓ for Model 1 (i.e., c 0 = c 2 = 0) and λ ss = −ℓ(1+c 0 x 2 +c 2 z 2 ) for Model 2 (i.e., a 1 = b 1 = 0), In both cases, the other two eigenvalues are 1 2 To ensure that these eigenvalues are real and of opposite sign, it suffices to check that The worst case is when y = 0, so we consider the terms not including y (and divide by 3 for simplicity): Divide by z 4 and introduce the new coordinate u = x 2 /z 2 : The left hand side is indeed negative for u = 0, and it is negative for all u if the equation −P u 2 + Qu − R = 0 has no real solution, i.e., if the discriminant is negative: Divide by a 2 b 0 and use the new coordinate γ = a0b2 a2b0 . Then we get the inequality which is equivalent to |γ − 5 9 | < 4 9 . That is, it fails if γ ≤ 1 9 or γ ≥ 1. Now consider equality in (18) and we divide by a 2 b 0 , which brings it to the form If γ ≥ 1, then these solutions are negative, and since u = x 2 /z 2 , this means that there are no solutions (x, y, z) ∈ U . The remaining case γ ≤ 1 9 is exactly the excluded case in the lemma. This concludes the proof.
Using our assumption a 2 b 0 < 9a 0 b 2 and Lemma 3.1, we obtain that the origin is the only point where we have λ ss = −ℓ and λ s = λ u = 0. Before continuing, we recall the definitions of a partially hyperbolic set and strongly dissipativity.
3.2 Definition. Let Λ be a compact invariant set of X ∈ X r (M ), c > 0 and 0 < λ < 1. We say that Λ has a (c, λ)−dominated splitting if the tangent bundle T Λ M has a DX t -invariant splitting of sub-bundles such that for all t > 0 and x ∈ Λ, we have We say that Λ is partially hyperbolic if it has a (c, λ)−dominated splitting such that E 1 is uniformly contracting; that is, for some c > 0 and all t > 0 and every x ∈ Λ it holds In this case we will denote E 1 by E s and call it the contracting direction. Also E 2 will be denoted by E cu and called the center-unstable direction.
3.3 Definition. Let G be a C ∞ vector field on R 3 with a Lorenz-like equilibrium; that is, the eigenvalues of DG p are real and satisfy 0 < −λ s < λ u < −λ ss . We say that G is strongly dissipative if the divergence of the vector field G is strictly negative; that is, there is c > 0 such that Div(G)(x) ≤ −c for every x and the eigenvalues of the singularity at p satisfy the constraint Figure 2 shows how the flows given by the Neutral model 1 and 2 send the lines in the y-direction in Σ to lines in the y-direction in S. Thus, its derivative D N t preserves the lines in the y-direction. Furthermore, by the way the flow from S to Σ was constructed ( Figure 1) we notice that horizontal lines in S; that is, parallel to the y-axis, are taken to parallel lines to the same axis in Σ. In other words, the flow from S to Σ preserves parallel lines to the y-axis. Since this flow is a composition of a rotation, an expansion and a translation, the derivative of the flow also preserves planes orthogonal to the y-axis. From this we can deduce that the splitting where N t is the flow obtained from Equation (4) with the corresponding parameters for model 1 and 2. Since we have uniform contraction along E ( D N t | Ex ≤ e λsst with λ ss = −ℓ < 0 for every x ∈ U ) and domination of the x ∈ U ) we can conclude that U is partially hyperbolic. It is worth noticing that the origin is the only point that spoils the singular hyperbolicity condition (a set A is singular hyperbolic if all its singularities are hyperbolic and it has volume expanding central direction) since J cu t (0) = | det D N t | F0 | = e (λu+λs)t = 1 (recall λ s = 0 = λ u ); that is, there is no area expansion along the subbundle F . Hence, Λ Neu is a partially hyperbolic attractor. Theorem 6 in [4] provides us with local strong-stable and center-unstable laminations W ss ϵ (x) and W cu ϵ (x), respectively, through the points x ∈ Λ Neu \ {0}. We note that both, W ss ϵ (x) and W cu ϵ (x), are embedded discs and hence submanifolds of M . Also W ss ϵ (x) is uniquely determined since E s is uniformly contracting. Corollary 6 in [4] shows us that the local strong-stable lamination can be extended to an invariant foliation F ss (x) of a open neighbourhood of Λ Neu with C 2 leaves and whose foliated charts are C 1 . Moreover, the leaves are uniformly contracted by the action of the flow.
We note that Σ is a C 2 embedded compact disk transversal to the flow N t . Furthermore, Σ is contained in the open neighbourhood V of Λ Neu . By Theorem 6 and Corollary 6 from [4] we know that local strong-stable lamination W ss ϵ (x) extends to an invariant foliation F ss (x). In this way, for x ∈ Σ we define W ss (x, Σ) to be the connected component of is the center-stable leaf. Since the flow (N t ) t∈R is C 2 , W ss (x, Σ) is a C 2 onedimensional embedded curve for every x ∈ Σ and their leaves form a C 1 foliation F ss Σ of Σ. Given a pair of embedded disks D 1 and D 2 in Σ intersecting transversally a set {W ss (x, Σ)} x∈Σ of stable leafs, the holonomy map H : Figure 5. From the developments in the partially hyperbolicity theory by Brin-Pesin [13] and Pugh-Shub [32] we have that the projection along leaves, also known as holonomies, between pair of transversal surfaces to F ss have a Hölder continuous Jacobian with respect to Lebesgue measure. This implies a similar conclusion for the holonomies transversal to F cs . It follows that the holonomy between pairs of transversal curves to F ss Σ along the lines of F ss Σ can be thought of as interval maps having a Hölder Jacobian. Hence these holonomies are C 1+α for some 0 < α < 1. In this setting, the leaves W s s(x, Σ), with x ∈ Σ, determine a foliation F ss Σ of Σ with transversal smoothness C 1+α . Therefore, we can assume that Σ is the image of the unit square I × I under the action of a C 1+α diffeomorphism h, for some 0 < α < 1. Furthermore, h sends vertical lines inside the leaves of F ss Σ . The next step is to prove that the strong stable foliation F ss (x) is not only C 1 but also C 1+α , for some 0 < α < 1. This was done by Araújo, Melbourne and Varandas in [6], stated as Lemma 2.2. This result is a consequence of domination of the splitting (Equation (19)), uniform contraction along the stable direction (Equation (20)) and strong dissipativity (Definition 3.3). Therefore, we can conclude that the neutral geometrical flow has a strong-stable foliation F ss (x) which is C 1+α . Furthermore, the modified return map also has a strong-stable foliation F ss Σ , whose transversal smoothness is C 1+α , with 0 < α < 1.
A final remark concerning this subsection, Theorem 6 in [4] is stated for a singular hyperbolic attractor. However, the conclusions and arguments still hold true if we consider a compact partially hyperbolic invariant set instead of a singular hyperbolic set. We will also like to mention that the situation for the origin is slightly different, since the splitting of the tangent bundle is given by a one-dimensional strong-stable direction E s and a two-dimensional center direction E c . However, we are only concerned with the existence of the strongstable manifold W ss (0), which follows from the theory of partial hyperbolicity.

The UNI condition
The main goal of this subsection is to show that the stable and unstable manifolds of the modified geometrical model are jointly nonintegrable. The joint nonintegrability of stable and unstable foliations can be interpreted as follows: The stable and unstable foliation of an Anosov flow are always transversal, therefore, if they are jointly integrable, this provides us with a codimension one invariant foliation which is transversal to the flow direction. In contrast, if there exists a codimension one submanifold transversal to the flow direction, then this foliation must be subfoliated by both the stable and unstable foliations. Thus they must be jointly integrable. In this situation it is known [22,Proposition 3.3] that the flow is semiconjugate to a suspension with a locally constant rooffunction over a subshift of finite type. Such a flow need not mix! From the work of Araújo, Butterley and Varandas [3], we know that the joint nonintegrability of the stable and unstable manifolds is equivalent to the uniformly nonintegrability (UNI) condition. As we will see later in this section, the UNI condition, stated formally in Definition 3.7, will ensure that the roof function of the suspension flow, which we will use to model the Neutral geometrical Lorenz flow, is not cohomologous to a constant function; that is, the time it takes for a point in the base to reach the roof of the suspension flow is not the same for every point. This property will guarantee the mixing properties.
In this chapter we will aim to prove the UNI condition. In [8] and [6] Araújo et al. prove the UNI condition by exploiting the properties obtained by using hyperbolic times. From now on, when we talk about the neutral geometrical model, we will refer to both models we constructed in Section 2.
Let I = [−1, 1] and f Neu : I → I be the one-dimensional Lorenz-like map obtained from the neutral geometrical model. We notice that {0} is a nondegenerate 1 exceptional set for f Neu . From Theorem 4.3 in [8] we know that there areX a neighbourhood of the singularity 0, a countable partitionQ of X Lebesgue modulo zero into subintervals, a function τ :X → Z + constant on partition elements and the induced mapF = f τ Neu :Q →X which is a C 2 uniformly expanding diffeomorphism with bounded distortion for anyQ ∈Q.
The motivation behind taking the inducing scheme is that we aim to extend the Gibbs-Markov mapF to a two-dimensional Gibbs-Markov map F and build the suspension flow F t over the map F with roof function R (see Equation (22)). This will allow us to use the results given in [11] to deduce the decay of correlations for the suspension flow F t and ultimately for the geometrical neutral Lorenz flow N t . Furthermore, we want to use the properties provided by the hyperbolic times and the bounded distortion of the mapF , to obtain the bound stated in Proposition 3.9 for the roof function R, which will help us prove the UNI condition for the geometrical neutral flow. We make now some observations regarding the induced mapF .
3.4 Observation. The mapF is obtained by inducing f Neu on the interval X, and the inducing time is given by the sum of a hyperbolic time with a nonnegative integer bounded by N , where N is such that Furthermore,F is a full branch Markov map ontoX since 0 has dense preimages under f Neu . For more details we refer the reader to [8].
In addition to the bounded distortion and uniform expansion ofF , we have the following inequalities for f Neu as a consequence of hyperbolic times. Given σ ∈ (0, 1) and c > 0, there is a constant b > 0 such that: 1.-(Backward contraction) LetQ(x) denotes the element of the partitionQ containing x, for y ∈Q(x) 2.-(Slow recurrence to the singular point) Following [6], our next step is to construct a piecewise uniformly hyperbolic map F with infinitely many branches, which coversF . First, let W ss PNeu (x) denote the stable leaf under the Poincaré map P Neu containing the point x, π : Σ → I the projection map. We define X = {W ss PNeu (x) | x ∈X} as the union of stable leaves alongX. We also extend the induced time τ to a function on X denoted also by τ and given by τ (x) = τ (π(x)). We make the following important observation on the tails of τ .
3.5 Observation. The tails of the return time τ and its extension also denoted by τ are exponential (see [33]); i.e., there exists a constant c > 0 such that µX (τ > n) = O(e −cn ) and µ X (τ > n) = O(e −cn ), where µX and µ X are the SRB-measures of the Gibbs-Markov maps F andF , see below.
Now, we construct F : X → X the Poincaré map by setting F (x) = P τ (π(x)) Neu (x) for x ∈ X. Furthermore, let Q be the measurable partition of X by taking {W ss PNeu (x) | x ∈Q} as its elements, withQ ∈Q. We will make use of X and F when making the model of the neutral geometrical Lorenz flow by a suspension flow.
It is standard [33] that the mapF has a unique a.c.i.p measure µX onX. Furthermore, r Neu ∈ L 1 (µX ) and there exists a unique invariant measure µ X for F , µ Σ for P Neu and µ I for f Neu satisfying π * (µ X ) = µX , π * (µ Σ ) = µ I , and also Moreover, µ X ≪ µ Σ and µ X (X) = 1, thus µ Σ (X) > 0. We take the induced roof function R : X → R + given by Notice that R is constant along stable leaves because r neu is constants along stable leaves. We also call R the quotient induced roof function R :X → R + .
With this in mind we can state the definition of the UNI condition. First, we give the definition of suspension flow.
3.6 Definition. Let (Σ, ν) be a probability space and P : Σ → Σ an ergodic measure-preserving transformation. Let r : Σ → R + be a measurable (Hölder continuous) roof function. We define the suspension space as The suspension flow f t : Σ r → Σ r is given by f t (x, u) = (x, u + t) computed modulo identifications and the measure µ = ν ×λ, where λ is the Lebesgue measure, is ergodic and f t -invariant.
In the finite measure case, we normalise byr = Σ rdµ so that µ = ν×λ r is a probability measure.
3.7 Definition. Let R : X → R + be a roof function as above, F t : X R → X R the suspension flow built over F : for h 1 , h 2 ∈ H n ; that is, inverse branches of F n . Then the UNI condition holds if there exist D > 0 and h 1 , h 2 ∈ H n0 for some sufficiently large n 0 ≥ 1, such that |ψ ′ h1,h2 | ≥ D. We saw in the previous subsection that the Poincaré map P Neu : Σ → Σ has a strong stable foliation. We observe that the leaves of this foliation cross Σ, hence the induced map F has a strong stable manifold W ss F (x) = W ss PNeu (x) that crosses X. Araújo et al. in [6,Proposition 2.4] provide us with local unstable manifolds of uniform size for F and defined almost everywhere, and by [6, Proposition 2.4] we obtain a local product structure, see Figure 6. To use the arguments given in [6] for the temporal distortion function we need to adapt the proofs of the uniform bound of the derivative of the induced roof function R and a bound on the flow time r Neu . An adjustment is required because we are changing the flow time from Σ to S. Recall that the original flow time for the geometrical Lorenz model is given by (3), whereas the modified flow time is given by (6); that is, we change from logarithmic to polynomial. More precisely, we have the following propositions. Proof. LetQ ∈Q and h ∈ H; that is, h :X →Q be an inverse branch ofF with inducing time τ = τ (Q) ≥ 1 and fix x ∈Q. We first observe that From the construction of the inducing partition using hyperbolic times (for more details on hyperbolic times see [2] and [8]), we have that backward contraction and slow recurrence to the singular point, see Observation 3.4; that is, there are constants σ ∈ (0, 1), b 0 > 0 and c ∈ (0, 1 2 ] such that 1.-(Backward contraction) For y ∈Q(x)

Notice that
where the inequality follows from the backward contraction.
Moreover, from the slow recurrence to the singularity we also get the following inequality; for some constant b 1 > 0. Altogether, Equations (23) and (24) imply that where s = 1 − c(1 + 1 β2 ). Since 0 < c ≤ 1 2 and β 2 > 2, we have s < 1. Therefore the sum converges, and we have that For x, y ∈X we define the separation time s(x, y) as the least integer n ≥ 0 such thatF n (x) andF n (y) are in different partition elements of Q 0 . For given 0 < η < 1, the symbolic metric is defined onX as d η (x, y) = η s(x,y) . Finally, we set r (k) Neu (x)). 3.9 Proposition. There exists B > 0 such that for all x, y ∈X with s(x, y) ≥ 1 and 0 ≤ k ≤ τ (x) = τ (y) we have |r denotes the Lipschitz constant of the quotient induced roof function R :X → R + with respect to d η . Moreover, |F (x) −F (y)| ≤ Bd η (x, y).
Proof. For convenience, in this proof f will denote f Neu . Let x, y ∈X such that s(x, y) = n ≥ 1 and 0 ≤ k ≤ τ (x) = τ (y). Thus y ∈Q n (x), whereQ n (x) = is the nth refinement ofQ(x), and so τ (F i (x)) = τ (F i (y)) for i = 0, . . . , n − 1. Hence, the choice of the cross-section assures that r Neu is constant along stable leaves and that r Neu (x) = |x| − 1 β 2 h(x) + τ 2 (x), where h(x) is bounded and bounded away from zero. In fact, h(x) is of the form h 0 + O(|x| γ ) where h 0 is a positive constant and γ > 0 depending on whether the higher order terms are consider in (7) or not. Also h(x) is differentiable for x > 0, because the Dulac map is differentiable. Then we can write We first notice that, The second inequality in (26) follows from Observation 3.4-1. Now, denote ) by A, then we obtain that, We notice that is bounded. Indeed we have the following: where C 0 > 0 and 0 < c ≤ 2 3 . The first and second inequalities follow from the Bernoulli inequality and from Observation 3.4-1 and 2, respectively. Since h(f i (x)) is bounded and by (28), we can bound the first term in the sum of (27) To finish the proof it remains to find a bound for Notice that |h( . Hence, by Observations 3.4-1 and 2, we have that (29) is bounded by Assuming that h ′ (ξ) is bounded and combining all the previous bounds we have that, . As in the previous proof, the sum converges since β 2 > 2, 0 < c ≤ 2 3 and hence 0 < s, u < 1. This establishes what we were aiming to prove.
One caveat we have to make here is that this argument is only valid if we assume boundedness for h ′ . This is true if the higher order terms in (7) are not present. If higher order terms are present, boundedness of h ′ is plausible, but since the required perturbation argument in [17] is less constructive so as to immediately derive this boundedness, we will try to convince the reader with the numeric analysis performed at the end of this work that this boundedness is indeed true. Thus, rather than a rigorous proof we give a combination of mathematical arguments and numerical verification. Now let Q ∈ Q be a partition elements for F . The temporal distortion function T : Q × Q → R is defined almost everywhere by, where [x, y] is the local product of x and y (see Figure 7). The second inequality follows from the property of r Neu of being constant along stable leaves. Now, for every x, y ∈ X in the same unstable manifold for F : X → X we define The continuity and other properties of T 0 are stated in [6, Lemma 3.1]. Furthermore, we can rewrite the temporal distortion function T (x, y) in terms of T 0 ; that is, The main result concerning the temporal distortion function establishes the joint nonintegrability of the stable and unstable foliations for the flow by proving that the temporal distortion function T is not identically zero, that is, there is Q ∈ Q and x, y ∈ Q such that T (x, y) ̸ = 0 and is stated and proven in [ For fixed x ∈ X, we define the map h : W u F (x) → R given by the map h is C 1 . Furthermore, there exists a nonempty open set U ⊂ W u F (x) such that h| U is a C 1 diffeomorphism. For the proofs of the properties of the map h see Proposition 3.6 and Corollary 4.7 in [6].
The next result will be of great help in proving the decay of correlations for the geometrical neutral Lorenz flow. The original statement involves the geometrical Lorenz flow, but the same arguments can be used to prove the same result for our setting. Before stating the result we give a definition.
3.11 Definition. Let X and F be as in the beginning of Section 6.3. A subset Z 0 ⊂ X is called a finite subsystem of X if Z 0 = n≥0 F −n Z, where Z is the union of finitely many partition elements of X.
Let Q 1 and Q 2 ∈ Q be two partition elements and consider Q = Q 1 ∪ Q 2 . We define the finite subsystem Q 0 = n≥0 F −n Q, then we have the following: 3.12 Proposition. [6, Proposition 3.8] For the finite subsystem Q 0 , the set T (Q 0 × Q 0 ) has positive lower box dimension.
We will like to end this Section with the following remark. To establish their results on decay of correlations, Bálint et al. in [11] and Melbourne in [30] assumed a very important, yet technical property namely, absence of approximate eigenfunctions. They also provide some criteria that guarantees the absence of approximate eigenfunctions. The first one, involves the temporal distortion function providing a nonintegrability condition. This criteria is given by Proposition 3.12; that is, when the temporal distortion function is not identically zero. In other words, when the UNI condition is satisfied. The second one, is a Diophantine condition on the periods of three periodic solutions [20], which is satisfied with probability one. It is important to remark that from these criteria, the UNI condition is robust while the Diophantine condition is not.

Decay of correlations
In this section we prove Theorem 1.1 for the first two neutral models. The third and more general model will be analysed in Section 5. We will use the results of Bálint et al. [11] to prove our theorem. In [11] polynomial decay of correlations for non-uniformly hyperbolic flows is proven under absence of approximate eigenfunctions. Let us start by giving the description of a nonuniformly hyperbolic flow described in [11].
First, we observe that the geometrical neutral Lorenz flow N t : Λ Neu → Λ Neu , where Λ Neu is the geometrical neutral Lorenz attractor, can be modelled as the suspension flow S t : Σ rNeu → Σ rNeu over the Poincaré map P Neu with base the cross-section Σ and roof function r Neu from (6). However, in order to use the results of [11], we take the alternative model F t : X R → X R , where X ⊂ Σ is a cross-section to the flow with nice hyperbolic structure (local product structure) and with induced roof function R : X → R + given by Neu (x)), see Section 3.2. Then the suspension flow F t built over the uniformly hyperbolic map F : X → X is identical to the suspension flow S t , thus F t is an extension of the underlying flow, namely the neutral geometrical Lorenz flow. Within this framework, N t is called in [11] a non-uniformly hyperbolic flow.
Under suitable conditions it can be shown that the suspension flow F t described above is a Gibbs-Markov flow [11,Section 6]. Therefore, the mixing rates for non-uniformly hyperbolic flows can be deduced from the corresponding results for Gibbs-Markov flows, see [11,Corollary 8.1].
For observables v and w, let ρ t (v, w) denote the decay of correlations of the geometrical neutral Lorenz flow; that is, where µ is the SRB measure of N t . Before giving the proof of Theorem 1.1 for the neutral models 1 and 2, we will give the definitions of the space of observables.
Let (M, d) be a metric space with diam(M ) ≤ 1 and define a flow T t : M → M on M . We fix η ∈ (0, 1] and for a given observable v : M → R we define and the norm ∥v∥ C η = |v| ∞ + |v| C η . We define the Banach space of Hölder functions on M by C η (M ); i.e., Furthermore, let and define ∥v∥ C 0,η = |v| ∞ + |v| C 0,η . We denote the space of Hölder observables in the flow direction by We will say that an observable w : We will denote the space of observables that are m-times differentiable in the flow direction by For a Borel set X ⊂ M we define C η (X) as above by using the restriction of the metric d to X.
Proof Theorem 1.1. We first note that the suspension flow F t projects to a quotient suspension semiflowF t :X R →X R , whereX is the quotient space obtained from X after quotienting out the stable leaves andF :X →X is a Gibbs-Markov map. Proposition 3.9 ensures that the following inequality holds.
where {Q i } i≥0 is a countable Lebesgue modulo zero partition into subintervals, 0 < γ < 1 and s(x, y) is the separation time.
Therefore, we have thatF t is a Gibbs-Markov semiflow and consequently that F t is a Gibbs-Markov flow. Then the conclusion follows from [11,Corollary 8.1].
There are still four details concerning the hypothesis in [11,Corollary 8.1] that we have not mentioned yet. The first one is regarding condition (H); for us this condition holds automatically since R is constant along stable leaves.
The second concerns the absence of approximate eigenfunctions for F t . Melbourne gave in [30,Chapter 5] sufficient conditions for the absence of approximate eigenfunctions, namely the existence of a finite subsystem with positive lower box dimension. Hence, it follows from Proposition 3.12 and Lemma 8.9 in [11] that F t has absence of approximate eigenfunctions.
The fourth and last detail concerns how to improve the estimates for µ X (R > t) and remove the logarithmic term. For this, we make use of [12,Lemma 4.1]. There the settings is made for infinite horizon planar periodic Lorentz gases, for that setting the tails of the flow time (in this work denoted by r Neu ) is of order O(t −2 ). By replacing in [12,Lemma 4.1] the order of the tails from O(t −2 ) to O(c 1 t −β2 ) we can use the same proof to remove the logarithmic term. Hence, we have that µ X (R > t) = O(t −β2 ). With this we conclude our proof.
The natural questions is about the lower bounds, i.e., if the bounds given in this theorem are sharp.
Although we definitely think they are, the currently available literature is insufficient to conclude this, although the margin is fairly narrow.
In [23], the renewal operator methods are developed to get such lower bounds, but his paper is for maps, not flows. Melbourne and Terhesiu come the closest in [31], where they consider suspension semiflows with polynomial roof functions over Gibbs-Markov base maps, and indeed, their results imply polynomial mixing for the flow F t : X R → X R , for Hölder observables that are constant on the stable fibers. The step from this suspension flow to the actual flow N t : Λ Neu → Λ Neu , however, is not trivial at all. (Here the effect of the second inducing step needs to be undone, in a way).
In [15] this step was taken for discrete suspensions (i.e., Young towers), specifically for billiard maps, but not for flows. However, the lower bounds that we obtain for F t : X R → X R as a corollary of the results in [31] are sufficient to prove stable laws (with exponent 1/β 2 ∈ (1, 2]) for the neutral Lorenz flow N t : Λ Neu → Λ Neu , cf. [14,18].

Numerical analysis and results
In this section we provide the results of the numerical approximation we obtained for the exponent β of the Dulac map (see Equation (8)) and the exponent β 2 of the tails of the return map (see Theorem 2.1) for the two-dimensional setting (the setting of [17]) and for the 3-dimensional setting concerning this work.

Numerics of the 2-dimensional case
To verify the existing theoretical asymptotics from [17], we will start the numerical analysis in the 2-dimensional case; that is, we consider the framework of [14] and [17]. There, the following neutral form was considered: where a 0 , a 2 , b 0 and b 2 > 0 and ∆ := a 2 b 0 − a 0 b 2 ̸ = 0. For simplicity we let κ = 2. For the analysis of the Dulac map close to the neutral equilibrium of Equation (34), we take an unstable leaf W u (0, y 0 ) and a stable leaf W s (x 0 , 0), then the Dulac map D : W u (0, y 0 ) → W s (x 0 , 0), shown in Figure 3, assigns the firs intersection of the integral curve through (x, y 0 ) with the stable leaf W s (x 0 , 0), where x ∈ W u (0, y 0 ) and T is the flow time; that is, the exit time.
For the setting considered in this work, we will perform the numerical experiments with x 0 = 1. In order to corroborate the estimates of the Dulac map given by Equation (8), we expect the numerical experiments to show us that We actually will see that .
From [14] and [17] we know that the constants c(y 0 ) are given by a specific formula which is not easy to compute. For this reason, we decided to use the least-squares method to calculate these constants.
For the numerical experiments we will take different values of β and 250 points x ∈ [1.0 × 10 −5 , 1.0 × 10 −4 ] at the unstable leave W u (0, y 0 ) with y 0 = 1.0. The integration method we will use for the numerical experiments concerning this work is the so-called Radau quadrature method, to deal with the numerical complications of integrating near a neutral stationary point, see [25].  We proceed in the same way to perform the numerical analysis for the exponent β 2 . From the estimates obtained in [17], we can see that, where t is the flow time that it takes a point from the unstable leaf to hit the stable leaf and x ∈ [1.0 × 10 −5 , 1.0 × 10 −4 ]. We will actually show that .
(38) Figure 9 shows the approximation of the exponent of the tail of the return map with values β 2 = 0.1333 and β 2 = 2.0 which correspond to the values β = 0.400 and β = 0.266, respectively. For this case the constants for the adjusted approximations are c(y 0 ) = ln(0.06) and c(y 0 ) = ln(0.045) for Figure 9 a) and b), respectively. From the approximations we see, as before, that the error decreases as we approach to x = 0.

Numerics of the 3-dimensional case
In this subsection we will perform the numeric experiments for the 3-dimensional models. We will start with Neutral model 1. Recall that the Neutral model 1 was given by Equation (39) where a 0 , a 1 , a 2 , b 0 , b 1 , b 2 and ℓ > 0 and ∆ := a 2 b 0 − a 0 b 2 ̸ = 0.Precise asymptotics are not available, but since y(t) decreases exponentially fast, the same asymptotics as in (34) are expected, and our numerics indeed confirm this. Note that the strong stable direction of (39) is still purely y-directed.
For the analysis of the Dulac map close to the neutral equilibrium of Equation (39), we will perform the numerical analysis on N 1 : Σ → S. For the purpose of this work, we want to show with the numeric experiment that the x and z components behave like the 2-dimensional model from the previous subsection regardless of the y value. To perform the numerical analysis we will take different unstable leaves W u (x, y 0 , z 0 ) and a stable leaf W s (1, 0, 0), where y 0 = 1.0 and z 0 = 1.0. Hence, like in the 2-dimensional analysis we expect the numerical experiments to show us that, where z is the last value of the integral curve with initial condition (x, y 0 , z 0 ) for x ∈ [1.0 × 10 −5 , 1.0 × 10 −4 ]. Again, we will actually show that As before, we will use the Radau quadrature method and take 250 points for the values of  Next, we consider the Neutral model 2 given by Equation (42), with the parameters satisfying the usual constraints, and present the numerical results obtained by performing the same experiments we did for the Neutral model 1. We consider this form since it is no longer a skew product like the previous model and poses a new challenge to deduce its asymptotics and u decay of correlations.  Until now we have performed the numerical experiments for the Neutral model 1 and 2 corresponding to the normal forms given by Equations (39) and (42), respectively. We saw there that for both models, the asymptotic behaviour of the x and z components are the same as the x and y component of the 2-dimensional model. Recall that for this two models we had explicit formulas for the map N : Σ → S and hence for the modified Poincaré map P Neu : Σ → Σ. Our goal is to see whether the asymptotic behaviour of the Neutral model 3 given by Equation (4) is similar to the other two models. Therefore, the numerical results obtained for the Neutral models 1 and 2 will be our reference and we will compare them to the numerical results obtained for the third model.  In the previous subsection we saw the numerical analysis on N 1 : Σ → S and showed, with the numeric experimentation, that the x and z components behaves like the 2-dimensional model. For the next numerical analysis, we will take an unstable leaf W u (x, y 0 , z 0 ) and a stable leaf W s (1, 0, 0), where y 0 = 1.0 and z 0 = 1.0. From the estimates obtained in [17] we can see that where t is the flow time that it takes a point from the unstable leaf to hit the stable leaf and x ∈ [1.0 × 10 −5 , 1.0 × 10 −4 ]. We will actually show that As before, we will use the Radau quadrature method and take 50 points for the values of  (44), and the theoretical value of β 2 , corresponding to the blue graph in all figures, is obtained from the parameters a 2 and b 2 ; that is, β 2 = a2+b2 2b2 . We start considering the neutral model 1. Figure 13 a) shows the approximation for β 2 = 1.333 which corresponds to the case β = 0.40, for the adjusted approximation the constant is ln(c(z 0 )) = 0.06 and b) displays the approximation for β 2 = 2.0 corresponding to the case β = 0.266 with adjustment constant ln(c(z 0 )) = 0.05. Next, we consider the neutral model 2. Figure 14 a) and b) show the approximation for β 2 = 1.333 and β 2 = 2.0, respectively. Their adjusted approximation the constant are ln(c(z 0 )) = 0.06 and ln(c(z 0 )) = 0.04, respectively. Finally, we consider the neutral model 3. Figure 15 a)   This numerical experiment has shown us a good approximation of the exponent of the decay of correlations for the 3 neutral models. From this we can deduce the same results concerning the decay of correlations, and obtaining Theorem 1.1 in its full generality.