{'fraction_non_alphanumeric': 0.04982753449310138, 'fraction_numerical': 0.020752099580083984, 'mean_word_length': 4.045875898600076, 'pattern_counts': {'":': 0, '<': 17, '<?xml version=': 0, '>': 18, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 18, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Rigorous Validation of Stochastic Transition Paths

Global dynamics in nonlinear stochastic systems is often difficult to analyze rigorously. Yet, many excellent numerical methods exist to approximate these systems. In this work, we propose a method to bridge the gap between computation and analysis by introducing rigorous validated computations for stochastic systems. The first step is to use analytic methods to reduce the stochastic problem to one solvable by a deterministic algorithm and to numerically compute a solution. Then one uses fixed-point arguments, including a combination of analytical and validated numerical estimates, to prove that the computed solution has a true solution in a suitable neighbourhood. We demonstrate our approach by computing minimum-energy transition paths via invariant manifolds and heteroclinic connections. We illustrate our method in the context of the classical M\"uller-Brown test potential.


Introduction
Steady states, connecting orbits and more general invariant sets are key objects to understand the long time behavior of a deterministic dynamical system. There are many theoretical existence results, but these results are often non-quantitative or difficult to apply, specially when the system is highly nonlinear. These limitations motivated the development of a field now called rigorous numerics. Its goal is to use computer-assisted verification in combination with analytic estimates to prove theorems. For example, it has turned out that the existence of many chaotic attractors could only be established with computerassisted methods [25,45]. For detailed reviews about the history and applications of rigorous numerics, we refer to [38,42,46,52] and the references therein.
Understanding the global behavior of a stochastic nonlinear dynamical system is even harder. The difficulties dealing with global dynamics induced by the deterministic nonlinear part remain, and the noise usually adds another layer of technical challenges. Nevertheless, it is frequently possible to reduce particular questions to the study of associated deterministic systems. Some examples include action functionals in large deviation theory [23,55], variance estimates for sample paths [7], pathwise descriptions of stochastic equations [37], moment closure methods [34] and partial differential equation (PDE) methods for stochastic stability [3,29]. However, if the original stochastic system was nonlinear, an associated deterministic reduced problem one will usually be also nonlinear.
In this work, we propose to combine the two strategies outlined above, i.e., to first use deterministic reduction and then employ rigorous a-posteriori validation to obtain information about the global behavior of stochastic dynamical systems. To illustrate the feasibility of our approach, we study a particular example. Here we only outline our strategy and we refer for more technical background to Section 2.
Consider a metastable stochastic gradient system [9] with potential V modelled by a stochastic ordinary differential equation (SODE) with additive noise. The minima of V correspond to deterministically stable, and for small noise stochastically metastable, states. Noise-induced transitions between minima are known to occur with high probability along minimum energy paths (MEPs) of a deterministic action functional. Many efficient numerical methods have been developed, particularly in the context of applications in chemistry, to compute MEPs [20,26,44]. For our example, MEPs can also be characterized as concatenations of heteroclinic connections [19,24] between the minima and intermediate saddle points. The main steps that are needed to rigorously validate a MEP are: (I) locate and validate saddles and minima, (II) introduce an equivalent polynomial vector field, and (III) validate heteroclinic orbits for this polynomial vector field.
Step (III) is more involved and can be decomposed in three sub-steps: (IIIa) compute and validate the local unstable manifold of the saddles, (IIIb) validate trapping regions around the minima, and (IIIc) compute and validate orbits connecting unstable manifolds and trapping regions. The output of the steps (I)-(III) is a theorem providing the existence and detailed quantitative location of the MEP. We use the standard Müller-Brown potential as a test case to present the requires steps but emphasize that our framework has many different potential applications. Our main contributions in this work are the following: • We establish the first link between stochastic dynamics and rigorous computer-assisted proofs.
• We apply our methodology to the case of MEPs in metastable stochastic system.
• We demonstrate the practical feasibility using the standard test case of the Müller-Brown potential.
• We also improve validation methods themselves by developing a priori almost optimal weights for the function spaces employed (see e.g. Proposition 4.2).
The paper is structured as follows: In Section 2 we collect different background from large deviation theory, from dynamical systems, from available algorithmic approaches to compute MEPs, and from rigorous numerics. Although these techniques are certainly well-known within certain communities, it seems useful to us to briefly introduce them here as we establish new links between several fields. In Section 3 we provide the basic setting for the rigorous validation to be carried out. In Sections 4-5 we carry out the technical estimates and required computations for step (III) regarding the local invariant manifolds and connecting orbits. The required code for the validated numerics can be found at [10].

Large Deviation Theory and Minimum Energy Paths
Consider a smooth potential V : R n → R with two local minima m 1 * and m 2 * , separated by a local maximum s * . The gradient flow ordinary differential equation (ODE) induced by the potential is Observe that m 1 * and m 2 * are locally asymptotically stable equilibria for (1) while s * is unstable. Since the Jacobian D(−∇V ) of the vector field is the (negative) Hessian of V , it follows that −D 2 V has only real eigenvalues. For the exposition here, we assume that n ≥ 2 and s * is a saddle point with a single unstable direction, i.e., −D 2 V (s * ) has a single positive eigenvalue. Denote by B 1 and B 2 the basins of attraction of m 1 * and m 2 * . If x 0 ∈ B 1 (resp. x 0 ∈ B 2 ) then x(t) converges to m 1 * (resp. to m 2 * ) as t → +∞. In particular there is no solution that connects m 1 * and m 2 * . For systems under the influence of noise, the situation is different. Let (Ω, F, F t , P) be a filtered probability space and consider a vector W (t) = W = (W 1 , W 2 , . . . , W n ) of independent identically distributed (iid) Brownian motions. Introducing additive noise in (1) gives the SODE where ε is a small parameter (0 < ε 1). Noise-induced transitions occur for (2), e.g., fixing small neighbourhoods of m 1 * and m 2 * , there exist sample paths γ = γ(t) solving (2) meeting both neighbourhoods in finite time. For the gradient system we considered here, it is well-known that the most probable path γ, also called minimum energy path (MEP), passes through the saddle point s * . Also the mean-first passage time is well-understood via the classical Arrhenius-Eyring-Kramers law [4,6,22,31], which states that transition probabilities are exponentially small in the noise level ε as ε → 0. A general framework for the small-noise case is provided by large deviation theory [23,55], which can be formulated for more general SODEs such as dX = f (X) dt + ε dW, where f is a sufficiently smooth vector field. Let φ : [0, T ] → R n be an absolutely continuous path over fixed time T > 0 and define the action (or action functional ) where the Lagrangian L is given by We can think of S T (φ) as the energy spent following a solution φ against the deterministic (ε = 0) flow for (3). Then Wentzell-Freidlin theory [23] states that the probability of a solution X = X(t) to pass through a sufficiently small tube around φ(t) for t ∈ [0, T ] is given by More generally, one can just take any set A ⊂ F of random events and obtain the large deviation principle Note carefully that this formulation has converted a stochastic problem of calculating/estimating a probability into a deterministic optimization problem. If a minimizer of the action functional exists one refers to φ * as a minimum action path (MAP). We reserve MEP for the case of a minimizer in gradient systems for a transition between states. If we want to specify the transition problem between two locally deterministically asymptotically stable points m 1 * and m 2 * , we should set where T 1 = −∞, T 2 = +∞ (or vice versa) in a theoretical framework while we select finite times for any numerical algorithm. In fact, it is usually not possible to calculate MAPs/MEPs analytically due to the need to solve a nonlinear deterministic problem. However, there are many numerical methods, mainly motivated by chemical applications. For MEPs in gradient systems, there is the classical nudged elastic band method [26] as well as the string method [56]. Both methods start with an initial piecewise-linear path in the phase space R n . The path is refined by decreasing the action. Both methods show often quite similar performance as just the re-meshing strategies differ [44]. Of course, one can also directly discretize the action and try to minimize it, which leads to another natural set of numerical methods [21]. In fact, one can also numerically compute MAPs/MEPs by noticing that in many situations they are (concatenations of) orbits for the deterministic system for ε = 0. For example, consider the gradient system ODE (1) with the three critical points m 1 * , s * , m 2 * . Then consider a path where γ s * m2 * is a heteroclinic orbit from s * to m 2 * and γ s * m1 * is a heteroclinic orbit from s * to m 1 * , i.e., Re-parametrizing γ as a path suitably, it is easy to see that it must be a minimizer of the action (4)-(5) since it is entirely composed of deterministic orbits, which have zero cost. Furthermore, it does satisfy the boundary conditions for an MEP from m 1 * to m 2 * if we reverse the time on γ s * m1 * . Hence, a MEP algorithm can also be built around calculating saddles and heteroclinic connections, which is the strategy we pursue here using rigorous numerics, i.e, we turn a numerical computation into a rigorous proof including error estimates for the location of the MEP. Nevertheless, it is very important to mention that the type of deterministic numerical algorithm used is secondary. In principle, most numerical algorithms can be combined with rigorous numerical error estimates. Therefore, our general strategy does not only apply to MEPs computed by heteroclinic connections, or more generally MAP computations, but the methodology can find applications in many other stochastic contexts. For example, there are recently developed methods, which numerically calculate the local fluctuations in stochastic systems around deterministic equilibria by using matrix equations [32,33], where rigorous numerics would be directly applicable.

Rigorous A-Posteriori Validation Methods
First, let us briefly mention that there are roughly two types of rigorous numerics techniques. The first kind are often described as geometric or topological methods. They aim at controlling the numerical solution directly in phase space, and are based on shadowing techniques [30], covering relations [61], cone conditions [60] and rigorous integration of the flow (see for instance [8,59]). The second kind are referred to as functional analytic methods and aim at validating a posteriori the numerical solution by a fixed point argument (see [16,47] for some seminal works in this direction, and [18,39,40] for some more recent references).
In this paper, we are going to use a functional analytic approach to validate numerically obtained MEPs by rigorously computing a chain of heteroclinic orbits. Let us mention that this is far from being the first time that rigorous numerics are used to study connecting orbits, and that there is a rapidly growing literature on the subject, see e.g. [2,12,17,27,43,49,57].
Given a numerical solutionX, the main idea of this approach is to construct an operator T such that: • fixed points of T correspond to genuine solutions; • T is locally contracting, in a neighborhood ofX.
If such an operator can be defined, one can then hope to apply a fixed point argument, to show that T has a fixed point in a neighborhood ofX. If this procedure is successful, we say that we have validated the numerical solutionX, since we have proven the existence of a genuine solution close toX. The following theorem gives some sufficient conditions to validate a numerical solutionX. Similar statements can be found in [1,18,40,58] and in many subsequent works.
Theorem 2.1. Let (X , · ) be a Banach space and let T : X → X be a map of class C 1 . Consider X ∈ X and assume there exists a non negative constant Y and a non negative, increasing and continuous function Z such that If there existsr > 0 such that then T has a unique fixed point in B(X,r), the closed ball of radiusr centered atX.
Proof. We show that the conditions in (13) imply that T is a contraction on B(X,r), which allows us to conclude by Banach's fixed point theorem. For X ∈ B(X,r), we estimate where we used to fact that Z is increasing to get the last inequality. Therefore (13a) yields that T maps B(X,r) into itself, and using (13b) we have that T is indeed contracting on B(X,r).
Before going further, let us make some comments about Theorem 2.1. While it is a purely theoretical statement, we use it in a context where its application requires the help of the computer at several levels. Firstly, a computer is often mandatory to obtain an approximate solution. Secondly, the definition of T will depend onX, as one can only hope to get a contraction close to the solution, and thus the bounds Y and Z that we are going to derive will depend onX (this will become evident in Sections 4 and 5). To be precise, obtaining formulas for Y and Z satisfying (12) is done by classical pen and paper estimates, but these formulas will involveX in various ways, and thus their evaluation in order to check (13) requires a computer. To make a mathematically rigorous existence statement out of a numerical computation, one must necessarily control two types of errors: • discretization or truncation errors, between the original problem and the finite dimensional approximation that the computer is working with; • round-off errors, coming from the fact that a computer can only use a finite subset of the real numbers.
The first issue is the most crucial one, and truncation errors must be estimated explicitly when deriving the bounds (12) (again see Sections 4 and 5 for explicit examples of such bounds). The issue of round-off errors can be handled quite easily thanks to the existence of many software packages and libraries which support interval arithmetic. In this work, we make use of Matlab with the Intlab package [41]. We point out that one of the advantages of functional analytic methods is that they usually allow for the numerical solution to be computed with standard floating-point arithmetic, and only require interval arithmetic to evaluate the validation bounds, i.e. to numerically but rigorously check (13). Finally, let us point out that Theorem 2.1 provides a qualitative existence statement, in the sense that the solution is proven to exists within an explicit distance, given by the validation radiusr (which in practice is very small), of the known numerical solutionX.
Remark 2.2. In Theorem 2.1, assumption (13a) is in fact enough to ensure the existence of a locally unique fixed point of T in a ball B(X,r), for somer ≤r. Indeed, first notice that (13a) implies that Z(0) < 1. We then distinguish two cases. Either Z(r) < 1 and we are done. Or Z(r) ≥ 1, in which case there exists a unique r 0 ∈ (0,r] such that Z(r 0 ) = 1. We then consider the continuous function P defined (for r ≥ 0), by which is negative at r =r by assumption (13a). Since P reaches its minimum at r = r 0 , we also have P (r 0 ) ≤ P (r) < 0. By continuity, there existsr < r 0 such that P (r) < 0, and by definition of r 0 we also have Z(r) < 1. Therefore, the assumptions (13) do indeed hold withr instead ofr, which yields the existence of a unique fixed point of T in B(X,r).

The main steps
We recall the main steps we use to validate a MEP: In this section, we outline what is required for each step. While the method introduced is rather general, we feel that the exposition is made clearer by focusing on an explicit example. Therefore, we consider for the remainder of the paper the Müller-Brown potential [28], which is a well-known test case for numerical methods in theoretical chemistry. It is given by

Finding and validating saddles points and local minima
For the Müller-Brown potential, numerically finding saddles and minima is rather straightforward, but for more complicated potentials this is where the numerical methods mentioned in Section 2.1 come into play.
Once numerical approximations of saddles and minima have been found, one can readily validate such points, for instance by using the built-in Intlab function verifynlss, that provides a tight validated enclosure of solutions of a (finite dimensional) nonlinear system. Alternatively, one could use Theorem 2.1.
We are going to detail, how this can be done for one saddle point. This example is rather trivial because the problem is already finite dimensional: we want to validate a zero of ∇V : R 2 → R 2 , and hence there are no truncation errors to take care of. However, we believe that presenting this simple example in details will help the reader to better understand the more complicated applications of Theorem 2.1 that are to come in Sections 4 and 5. Using Newton's method, we get the following numerical approximation for the first saddleX = (−0.822001558732732, 0.624312802814871).
Our aim is now to prove, using Theorem 2.1, that there is a zero of ∇V in a neighborhood ofX. To do so, we introduce the operator where I 2 is the 2 × 2 identity matrix, and is a numerically computed approximate inverse of the Hessian D 2 V (X). Therefore, T is a Newton-like operator and should be a contraction aroundX. We now consider on R 2 the supremum norm (x, y) = max(|x|, |y|) and compute the bounds Y and Z satisfying (12). We point out that, while the choice of the norm is rather unimportant for this simple example, it will be crucial in the more involved cases presented later.
To get a Y bound satisfying (12a), we simply define Notice that the evaluation of this formula must be done with interval arithmetic, because we need to be sure that T (X) −X ≤ Y . We get Next, we define a Z bound satisfying (12b). In fact, the condition (12b) only needs to hold in a neighborhood ofX, therefore we introduce an a priori upper bound r * = 10 −5 on the validation radius, and aim at obtaining a function Z such that We will only have to check that the validation radiusr that we eventually find satisfiesr ≤ r * . Based on the following splitting By the triangle inequality and the mean value theorem, we indeed have that (14) holds. Again, the evaluation of Z 1 and Z 2 must be done using interval arithmetic. Notice that (an upper bound of) sup X∈B(X,r * ) D 3 V (X) can be easily obtained with Intlab, by evaluating D 3 V (X) , whereX = X + [−r * , r * ] 2 . We obtain Z 1 = 6.606610122939668 × 10 −15 and Z 2 = 23.179491548050574.
To conclude, we must find r > 0 satisfying (13). The first condition may be re-written as and the second simply as We check, using interval arithmetic, that satisfies these two inequalities, and thus Theorem 2.1 yields the existence of a unique zero X of ∇V such that X −X ≤r. We have validated the numerical zeroX.

An equivalent formulation with a polynomial vector field
We are now almost ready to validate a heteroclinic orbit between a saddle and a minimum, for the vector field −∇V . Before doing so, we extend the dimension of the phase space, to obtain an equivalent vector field with only polynomial nonlinearities. As we will see in Sections 4 and 5, this reformulation is a sort of trade-off. On the one hand, having polynomial nonlinearities makes it easier to obtain of some estimates, whereas on the other hand, having more dimensions increases the computational cost. We point out that it is only the dimension of the phase space that is increased, and not the dimension of the object we have to validate, hence the increase in computational cost is not dramatic. To obtain this new vector field with only polynomial nonlinearities, we introduce variables for the non polynomial terms in −∇V ; that is we define consider z (i) = ψ (i) (x, y) for i = 1, . . . , 4, and compute the differential equations satisfied by these new quantities, assuming This leads us to consider the following extended variables where the z (i) are now independent variables, and the associated extended system is where the vector field is given by for j = 1, . . . , 4, and we define The formulations (15) and (17) are equivalent in the following sense. The system (17) was designed so that, if (x, y) is a solution of (15) then X defined as in (16) with z = ψ(x, y) solves (17). Conversely, by the uniqueness statement of the Picard-Lindelöf Theorem, we get that a solution of (17) with a suitable initial (or asymptotic) condition gives a solution of (15).
) be a solution of (17). Assume one of the following conditions holds: Then, (x, y) solves (15).
We can then write and given a solution X = (x, y, z (1) , z (2) , z (3) , z (4) ) of (17), (x, y) solves (15) if and only if z(t) = ψ(x(t), y(t)) for all time t for which the solution is defined. However, since X solves (17) we have and thus either (i) or (ii) ensures that z(t) = ψ(x(t), y(t)) for all time t for which the solution is defined.
We also have the following statement which shows that for any equilibrium point of (15) going to the extended system (17) only adds center directions.
Let (x 0 , y 0 ) be a zero of ∇V , and X 0 be defined as in (16) In particular, Df (X 0 ) has the same spectrum as −D 2 V (x 0 , y 0 ), up to an additional 0 of algebraic multiplicity 4.
Proof. With the same notations as in the previous proof, we have −∇V (x, y) = g(x, y, ψ(x, y)) and f (X) = I 2 Dψ(x, y) g(x, y, z).
If (x 0 , y 0 ) is such that ∇V (x 0 , y 0 ) = 0, and z 0 = ψ(x 0 , y 0 ), then g(X 0 ) = 0 and hence The result follows from the identity The polynomial reformulation is based on ideas from automatic differentiation, see [35] and the references therein.

Validation of heteroclinic orbits
We are now ready to validate heteroclinic orbits. Since these orbits features boundary conditions at infinity, we use the method of projected boundaries [5]. That is, given a saddle point S and a minimum M , we are going to look for an orbit X : Therefore, we first have to compute and validate a local unstable manifold for S and a local stable manifold for M , and then compute and validate an orbit that connects them. The fact that M is a minimum makes the determination of a local stable manifold rather straightforward, as will be explained just below. The validation a local unstable manifold for the saddle, and of an orbit connecting this unstable manifold to the minimum are more involved. We briefly present the main ideas in subsequent paragraphs, and postpone the technicalities to Sections 4 and 5 where all the needed estimates are derived.

Validation of trapping regions around the minima
In our example, each considered minimum is strict, and therefore its stable manifold contains a whole neighborhood of the minimum. We only have to explicitly find such a neighborhood. Using Intlab, we can prove the existence of a small square around each minimum, such that 1. V is strictly convex in the square, 2. −∇V is pointing strictly inward, on the boundary of the square.
The second property ensures that an orbit entering the square never leaves it again, and the first proves that there is no other minimum in the square, and hence that any orbit entering the square converges to the minimum. Such a square is thus part of the stable manifold of the minimum.

Remark 3.3.
A more precise description of the dynamics near the minimum could be obtained by applying the techniques described in Section 4 to compute a validated parameterization of the local stable manifold. This approach could also be used in the case where the minimum is not strict.

Computation and validation of a local unstable manifold for the saddles
We now explain the main ideas that we use to compute and validate a local unstable manifold for each saddle. Our technique is based on the parameterization method, introduced in [13,14,15]. We want to obtain a parameterization that conjugates the dynamics on the unstable manifold with the unstable dynamics of the linearized system. To be more precise, let us consider a saddle point S, in the extended phase space R 6 , and denote by λ the unstable eigenvalue, together with an associated eigenvector v. As in Section 3.1, λ and v can first be computed numerically and then validated rigorously, using Theorem 2.1 or via the built-in verifyeig function of Intlab. We look for a parameterization p : where φ t is the flow generated by f . However, this formulation is not the most convenient one to work with, because it involves the flow. To get rid of it, one can take a time derivative of the above equation and evaluate at t = 0, to obtain the following invariance equation One can check that, if p is such that p(0) = S and solves (19), then p satisfies (18), therefore p([−1, 1]) is a local unstable manifold of S. The invariance equation (19) is the one we are going to solve, first numerically and then rigorously by applying a variation of Theorem 2.1. Let us already mention that, since Theorem 2.1 is based on a contraction argument, it can only be used to validate locally unique solutions. However, the invariance equation (19) together with the condition p(0) = S has a one parameter family of solutions (corresponding to different rescalings) and to isolate the solution one has to add a constraint on the norm of the derivative, i.e. fix a scalar γ such that p (0) = γv.
Since f is analytic, there is an analytic solution to (19) (see [13]), therefore we can look for a power series representation of p. Starting from p(0) = S and p (0) = γv, and taking advantage of the fact that our extended vector field f is polynomial, it becomes rather straightforward to recursively compute the coefficients of the Taylor expansion. Once sufficiently many coefficients have been computed, the remaining tail can be controlled by Theorem 2.2. A detailed exposition of this procedure is the subject of Section 4. Before proceeding further, let us mention that the approach presented here, i.e., combining the parameterization method with a posteriori validation techniques, was first used in [53] and then further developed in several subsequent works (among which [11] where the role of the scaling of the eigenvectors is investigated in details, and [51], where an extension to treat resonance conditions is introduced).

Computation and validation of orbits connecting unstable manifolds and trapping regions
Our last objective is to validate an orbit starting on the unstable manifold of a saddle and ending in the attracting set of a minimum, that is we want validate a solution X : [0, τ ] → R 6 of the problem for a time τ large enough so that X (1) (τ ), X (2) (τ ) lies in the validated trapping region of a minimum of the potential V . Notice that, X(0) = p(1) ensures that the orbit is on the unstable manifold of the saddle (in the extended phase space), and by Lemma 3.1 (ii) we get that X (1) , X (2) indeed solves (15) and is on the unstable manifold of the saddle in the original two-dimensional phase space. To validate a solution X of (20), we use an approach based on Chebyshev series, introduced in [36]. The method is similar in spirit to the one described in the previous subsection to solve the invariance equation, the main difference being that we look for a solution represented as a Chebyshev series and not as a power series. Again, our approximate solution will only have a finite number of coefficients, and we control the error with the true solution thanks to a variation of Theorem 2.1. One difference with the situation described in the previous subsection is that the coefficients cannot be computed recursively, as the obtained system is fully coupled. Therefore, we have to consider a fixed point operator on the whole sequence of coefficients rather than just of the tail, but this operator has no reason to be a contraction. We remedy to this by introducing a Newton-like reformulation. In practice, the integration time τ for which we have to validate the orbit until it reaches the trapping region of the minimum can be quite long, and maybe too long for the orbit to be validated using a single Chebyshev series. To validate the orbit for longer times, we use the domain decomposition principle introduced in the context of a posteriori validation using Chebyshev series in [54]. The procedure outlined in this subsection is presented in full details in Section 5.

Validation of local manifolds 4.1 Background
The material introduced in this subsection is mostly standard in rigorous numerics. It is mainly included to fix the notation and to select/refine suitable operator norms for our computations. The reader familiar with the notions presented here might skip Section 4.1 on first reading, proceed to Section 4.2, and only refer back when needed.

Notations
For X ∈ R 6 , we denote its components by X (i) , i = 1, . . . , 6. For u ∈ R N , we denote its components by u n , n ∈ N. Given two sequences u, v ∈ R N , we denote by u v their Cauchy product For p ∈ R N 6 , using the isomorphism between R 6 N and R N 6 , we can think of p either as a sequence indexed on N and with values in R 6 , that is p = (p n ) n∈N , p n ∈ R 6 , or as 6 sequences indexed on N and with values in R, that is To a sequence p ∈ R 6 N can be associated a power series p with values in R 6 , defined by We always use this convention of denoting a (possibly multivariate) function like p with an unbold character, whereas we use a bold character like p for the associated sequence of coefficients. We introduce in the next subsection a subspace of R 6 N for which such series is guaranteed to converge, at least for θ ∈ [−1, 1]. Similarly, we denote by f : R N 6 → R N 6 the map such that, for any power series p, f (p) are the coefficients of the power series f (p), that is and so on.

Norms and Banach spaces
For η ∈ R 6 >0 , we define the following weighted 1-norm on R 6 : We denote the usual 1 -norm on R N by · 1 , that is We recall that To highlight the key space we are going to use, we state it in the next definition: Notice that we have and with a slight abuse of notation where p 1 must be understood as p (1) 1 , . . . , p (6) 1 .

Operator norms
1≤i,j≤6 is a 6 × 6 matrix, we still denote by |A| η the associated operator norm, that is We recall that For a given matrix A with positive coefficients, the weight η can be chosen in a optimal way to minimize |A| η . This is the content of the following proposition. Besides Proof. The existence ofη is given by the Perron-Frobenius Theorem. Identity (21) then yields and thus |A|η = ρ(A). Finally, since the spectral radius is bounded from above by any operator norm,η is indeed a weight that minimizes |A| η .
Remark 4.3. The above proposition allows us to choose weights in an optimal way, at least with respect to one crucial estimate needed to control the contraction rate and apply Theorem 2.1 (see Sections 4.3.2 and 5.3.3). An alternative approach would be to include the weights in the validation radiusr of Theorem 2.1 (see for instance [48]).
If A = (A m,n ) m,n∈N is an "infinite matrix" representing a linear operator on 1 , we still denote by A 1 the associated operator norm, that is We recall that Finally, a linear operator A on X η can be represented as an "infinite block-matrix", that is A (1,6) . . . . . . . . .
where each block A (i,j) is a linear operator on 1 . Each of these blocks can themselves be represented as infinite matrices A (i,j) = A (i,j) m,n m,n∈N and their 1 operator norm is then given by By a slight abuse of notation, when A is a linear operator on X η , we consider that · 1 applies block-wise to A, that is we define Similarly, |·| η applies component-wise to A 1 , that is we define

We have that
Notice that, once A 1 has been computed, Proposition 4.2 can be used to find a weight η that minimize | A 1 | η . Similarly, if B is an k-linear operator on X η , we still denote by B X η its operator norm, defined as and have that

Framework
Let (x 0 , y 0 ) ∈ R 2 be saddle point of V , denote by λ > 0 the unstable eigenvalue and by µ < 0 the stable one. We want to validate a parametrization of the local unstable manifold. We consider X 0 defined as in (16) with z 0 = ψ(x 0 , y 0 ) and work with the extended system, given by f . By Lemma 3.2, λ is also the unique unstable eigenvalue of Df (X 0 ). We denote by v an associated eigenvector. We look for a parametrization p : R → R 6 of the local unstable manifold at X 0 , that satisfies p(0) = X 0 , p (0) = γv, for some γ > 0 to be chosen later, and the invariance equation We look for a power series representation of p, that is The invariance equation (25) can be rewritten as an equation on p, yielding: We denote π n (p) = (p 0 , . . . , p n ) ∈ R 6 n+1 , which we also may identify with the element (p 0 , . . . , p n , 0, . . . , 0, . . .) ∈ R 6 N .
Notice that, for any u, v ∈ R N , the n-th coefficient of a Cauchy product (u v) n only depends on the coefficients (u m ) 0≤m≤n and (v m ) 0≤m≤n , and therefore f n (p) in fact only depends on π n (p), i.e.
Next, we want to extract from f n (p) the contribution of p n . To do so, we consider the Taylor series representing f (π n (p)) that is p k θ k p n θ n + higher order terms.
Looking at the coefficient of degree n in each of the Taylor series above, we get f n (π n (p)) = f n (π n−1 (p)) + Df (p 0 )p n , ∀n ≥ 1, and there is no contribution from the higher order terms since they have only coefficients of degree 2n or more. Introducing the 6 × 6 matrices M n = nλI 6 − Df (X 0 ), which are invertible for all n ≥ 2, the invariance equation (27) for the coefficients p n , n ≥ 2, rewrites nλp n − Df (X 0 )p n = f n (π n−1 (p)) p n = M −1 n f n (π n−1 (p)).
Therefore, starting from p 0 = X 0 and p 1 = γv, the coefficients p n can be computed recursively. However, in practice one can only compute finitely many coefficients. Assume that the coefficients p 2 , . . . , p N −1 have been computed, for a given N ∈ N. For any p ∈ X η , definep := (p 0 , . . . , p N −1 ), p := (p N , p N +1 , . . .). If γ is chosen appropriately (see Remark 4.13), and if the truncation mode N is large enough, the finite series given byp should be a good approximate parameterization of the local unstable manifold. To verify this claim, we introduce the operator T : where T n (p) = M −1 n f n (π n−1 (p)), ∀n ≥ N, with p = (p,p). For convenience we also introduce the function g defined by g n (p) = f n (π n−1 (p)), so that T n (p) = M −1 n g n (p). We note that the coefficients ofp are now simply parameters that have already been computed, rather than variables. Ifp is a fixed point of T , then p = (p,p) solves (27) and hence corresponds to an exact parameterization of the local unstable manifold. Our goal is to show that there exists a fixed point of T in a neighborhood of 0, which would prove thatp is indeed an approximate parameterization. Our main tool is going to be the following corollary of Theorem 2.1.
Corollary 4.4. Let X be a Banach space, and T : X → X a polynomial map of degree 4. LetX ∈ X and Y , Z 1 , Z 2 , Z 3 , Z 4 be non negative constants such that If there existsr > 0 such that P (r) < 0, then for all r in the non empty interval (r min , r max ), T has a unique fixed point in B(X, r), where r min is the smallest positive root of P and r max is the unique positive root of Q.
Proof. Just notice that Z defined by satisfies (12b). Condition (13a) is then equivalent to having P (r) < 0, while condition (13b) is equivalent to having Q(r) < 0, and Theorem 2.1 can then be applied to conclude.
Remark 4.5. At this stage, having that T is a polynomial is convenient but absolutely not mandatory. If T does not have a finite Taylor expansion (or if it is too cumbersome to use the full Taylor expansion), one can proceed is in Section 2.2 and consider an a priori radius r * > 0, together with non negative constants Y , Z 1 and Z 2 such that We point out that, in practice, such a Z 2 can be obtained by bounding sup X∈B(X,r * ) Then, defining Z(r) = Z 1 + Z 2 r, condition (12b) holds for all X ∈ B(X, r * ). Therefore, if assumptions (13) are satisfied for somer ≤ r * , we can still conclude that T has a unique fixed point in B(X,r).
Remark 4.6. P and Q are often called radii polynomials in the literature.
Our goal is to apply this corollary to the operator T defined just above, withX = 0, to validate the approximate parameterization defined byp.

The bounds needed for the validation
In this subsection, we obtain computable bounds Y and Z i , i = 1, . . . , 4, satisfying (29).

The bound Y
Then Proof. Since f is quartic andp only has coefficients up to order N − 1, T n (0) = 0 for all n > 4N − 4.
Notice that Y can be rigorously evaluated on a computer (or more precisely upper-bounded, using interval arithmetic).

The bound Z 1
First, we need to control the norm of M −1 n . Lemma 4.8. Let n ∈ N and η ∈ R 6 >0 such that the following holds Then we can conclude that Proof. This follows directly from a standard Neumann series argument.
Remark 4.9. In practice, the bound on M −1 n η ends up being one of the critical factors influencing the success or failure of the validation procedure. Therefore, we choose η so as to minimize |Df (X 0 )| η (see Proposition 4.2), in order to get the smallest possible estimate out of Lemma 4.8.
The downside of Lemma 4.8 is that it can only be used for n larger than |Df (X0)| η λ . If one wishes to choose the truncation level N smaller than this threshold, the following alternative bound can be used for the low order modes, i.e., for all modes below the truncation level.
Remark 4.11. In practice, such matrix Q can be obtained numerically and then validated using the function verifyeig from Intlab, combined with the fact that we know a priori by Lemma 3.2 that the eigenvalues of Df (X 0 ) can only be λ, µ and 0.
Combining the two above lemmas, we define Next, we introducing the linear operatorM : We are now ready to define the Z 1 bound.
We emphasize that the bound Z 1 defined just above is computable, since the 1 operator norms D j g (i) (p) 1 are very easy to evaluate. Indeed, the action of the linear operator D j f (i) (p) on 1 is nothing but a convolution with the sequence representing the function D j f (i) (p), and therefore these 1 operator norms are simply given by the 1 norm of the corresponding vector. For instance, since we have To go back to Dg(p), notice that, for any u, v ∈ 1 , where π 0 (u) := (0, u 1 , u 2 , . . .). Finally, using that π n−1 (p) =p for all n ≥ N , we get that Dg n (p) = Df n (p)π n−1 , and thus the 1 operator norm of D j g (i) (p) is nothing but the 1 norm of π 0 u, where u is the vector representing D j f (i) (p). For instance, we have and sincep only has finitely many non zero coefficients, this quantity can be evaluated (or more precisely, upper-bounded using interval arithmetic).
Remark 4.13. About the role of γ. The approximate parameterizationp is uniquely determined by the choice of the scaling γ in (27). In this remark, we notep =p(γ) to highlight this dependency. One crucial observation is that we havep n (γ) = γ np n (1), ∀ n ∈ N.
Therefore, D j g (i) (p(γ)) 1 goes to 0 when γ goes to 0, which means that we can have a Z 1 bound that is arbitrarily small (and in particular strictly less than 1) by taking γ small enough. In other words, up to considering a small enough patch of the local manifold, we can always get a contraction.

The bounds
Proposition 4.14. Consider T defined in (28), η ∈ R 6 >0 and define for 2 ≤ k ≤ 4 Then Proof. We have and we estimate the norm of the k-linear operator D k g(p) as in (24).
Again, this bound can be evaluated since we can compute the 1 operator norms D k (j1,...,j k ) g (i) (p) 1 .
Indeed, such multilinear operator norm is equal to the 1 norm of the vector of Taylor coefficients representing D k (j1,...,j k ) g (i) (p), and thus can be computed using interval arithmetic. For instance

Summary
We recall that we explained in Sections 3.1 and 3.3.2 how to obtain a rigorous enclosure of saddles points and of the associated eigenvalues and eigenvectors. We now show how all the estimates derived up to now can be combined with Corollary 4.4 to rigorously validate a local stable manifold around each saddle.
has been computed recursively rigorously (i.e. with interval arithmetic) so that p n solves the invariance equation (27) for all n < N . Consider the operator T defined in (28) and the bounds Y , Z 1 and Z k , 2 ≤ k ≤ 4 defined in (32), (33) and (34) respectively. Finally, assume there existsr > 0 such that P (r) < 0 and Q(r) < 0, with P and Q defined in (30). Then there exists p n ∈ R 6 , n ≥ N , such that and such that the function p defined as in (26) satisfies the invariance equation (25), that is p([−1, 1]) ⊂ R 6 is a local unstable manifold of X 0 for the vector field f . Besides is a local unstable manifold of (x 0 , y 0 ) for the vector field −∇V .
This theorem proves that the functionp defined bŷ is an approximate parameterization of the local unstable manifold of X 0 in the sense that We give in Section 6 several examples (with explicit values of the parameters) of applications of Theorem 4.15 to validate local manifolds for the Müller-Brown potential.

Background
The material introduced in this subsection is standard, and mainly included for the sake of completeness and to fix some notations. The reader familiar with the notions presented here might skip Section 5.1 on first reading, and only refer back when needed.
For k ∈ N, we then denote by T k the Chebyshev polynomial of order k, defined by T k (cos(θ)) = cos(kθ). We consider τ > 0 and a partition 0 = t (0) < t (1) < . . . < t (M ) = τ . For all k ∈ N and m = 1, . . . , M , we also introduce the rescaled Chebyshev polynomial T We associate to a sequence X ∈ R N 6 M the function X : [0, τ ] → R 6 defined as a piece-wise Chebyshev series by Notice that, for all k ∈ N, X (m) k = X (m,1) k , . . . , X (m,6) k ∈ R 6 . We also introduce in the next subsection a subspace of R 6 N M for which such series is guaranteed to converge. Given two sequences u, v ∈ R N , we denote by u * v their convolution product In this Section, f : R N 6 → R N 6 denotes the map such that, for any function X represented by a piecewise Chebyshev series X, f (X) are the coefficients of the piece-wise Chebyshev series representation of f (X), that is , and so on.

Norms and Banach spaces
Definition 5.1. For ν > 1, we define the following weighted 1 norm. We first introduce the weights For u ∈ R N , we then define We recall that we have u * v ν ≤ u ν v ν .
Next, we consider the product space X η,ν = M m=1 X η (m) ,ν endowed with the supremum norm. That is, for X ∈ R N 6M we define

Operator norms
We use the same block-representation as in Section 4. We recall that We also introduce a notation for weighted 1-operator norms with two different sets of weights: A (m,i;n,j) ν η (m,i)

Framework
We want X to satisfy X = f (X), X(0) = p(1), X(τ ) ∈ W s loc (M ). If we look for a piece-wise Chebyshev series representation of X as in (36), we obtain the following equations on the coefficients X The derivation of these equations follows readily (see for instance [36,54]) from well known properties of the Chebyshev polynomials, namely We define F by with the convention X (1), and we want to validate a zero X of F in X η,ν . The main difference with the situation considered in Section 4 is that we cannot solve first for a finite number of modes and then get a contraction for the tail, because the system is now fully coupled (since we use Chebyshev series) instead of triangular (since we used Taylor series for the manifold). To be precise, one can still solve numerically for a finite number of modes by considering a truncated system, but the obtained numerical solution and the tail must both be validated simultaneously. To do so, we introduce a Newton-like reformulation, based on a truncated system, that is suitable for our a posteriori validation procedure. We denote by π K : X η,ν → R 6M K the finite dimensional projection obtained by truncating the Chebyshev modes of order K and higher, that is π K (X) = (X 0 , . . . , X K−1 ).
We also denote by ı K the natural injection from R 6M K to X η,ν . First, we compute an approximate solutionX by solving numerically the finite dimensional problem F [K] = π K • F • ı K = 0. We use the same notation to denoteX ∈ R 6M K and its injection in X η,ν . Then, we define the linear operator A † by A † π K X = DF [K] (X)π K X, Next, we compute A [K] an approximate inverse of DF [K] (X), and define the linear operator A by The operator A enables us to recover an equivalent fixed-point formulation which should give a contraction aroundX, by considering the Newton-like operator

The bound Y
Since f is polynomial and X only has finitely many non zero coefficients, so does F (X), and hence can be evaluated on a computer using interval arithmetic.
Remark 5.4. This bound requires the evaluation of the parameterization p, since it involves F (0) 0 (X) which is defined as In practice, we only have an approximate parameterizationp, together with a validated error boundr, see Section 4.4. Thus, whenever we have to evaluate p(1) we use
We point out that, by construction of A † and A, each block B (m,i;n,j) is finite (only the first K × K coefficients are potentially non zero), and hence the 1 ν operator norms B (m,i;n,j) ν can all be evaluated on a computer.
Let u (n;l,j) be the vector of Chebyshev coefficients representing the partial derivative D (n,j) f (n,l) (X) and define and C (m,i;n,j) tail with the convention that A Then, for all 1 ≤ m ≤ M Proof. We start from and estimate each 1 ν operator norm C (m,i;n,j) ν . We do so by computing explicitly (with a computer and interval arithmetic) the norm of a finite number of columns of C (m,i;n,j) , and by estimating by hand the norm of the remaining columns, as emphasized by the following splitting: A (m,i;n,l) (DF (X) − A † ) (n,l;n,j) the second term being 0 if n = M . This yields Estimate (48) comes from a meticulous but rather straightforward analysis of the various terms appearing in (47) for the columns of indices l > 4K − 2, using that Estimates similar to (48) were obtained previously in [50,54].
Again, we point out that each row and column of A and DF (X) − A † only have a finite number of non-zero coefficients, hence the quantities C and thus we choose η (m) so as to minimize this quantity, as explained in Proposition 4.2.

The bounds
>0 andX ∈ R 6M K (again identified with its injection in X η,ν ). Consider F , A † and A as in (39), (40) and ( Then, for all 2 ≤ k ≤ 4 and all 1 ≤ m ≤ M Proof. This estimate follows readily from writing down what the norm of the multilinear operator is and simply using triangle inequalities, which is a straightforward computation albeit being quite lengthy due to the intricate norm used on X η,ν .
Again, each of the involved 1 ν operator norms can be evaluated explicitly and rigorously on a computer, using interval arithmetic.
Remark 5.9. A slightly less sharp bound, but faster to evaluate on a computer, is given by

Summary
We recall that we explained in Section 3.1 how to obtain a rigorous enclosure of minima and saddle points of V , in Section 3.3.1 how to validate trapping regions around minima, and in Sections 3.3.2 and 4 how to validate local unstable manifold for the saddles. We now show how all the estimates derived in this section can be combined with Corollary 5.3 to rigorously validate a connecting orbit between a saddle and a minimum.
Then there exists X ∈ X η,ν satisfying such that the function X defined as in (36) satisfies X = f (X) on [0, τ ] and such that X(0) belongs to the unstable manifold of the saddle. Besides, if (X (1) (τ ), X (2) (τ )) ∈ U , this yields that the existence of a connecting orbit from the saddle to the minimum for the vector field −∇V .
Notice that (50) gives a control of the error in phase space, namely We give in Section 6 several examples (with explicit values of the parameters) of applications of Theorem 5.10 to validate a MEP for the Müller-Brown potential.

Results
We use Theorem 4.15 and Theorem 5.10 to validate a MEP for the Müller-Brown potential. This MEP connects two minima through two saddles and one extra minimum (see Figure 2), and is thus composed of four connecting orbits: To validate this MEP, we follow the main steps described in Section 3. We now make precise for each step the values of the various parameters that we use and present the output of the whole procedure.
(I) Locate and validate saddles and minima. Using the method exposed in Section 3.1, we obtain the following validated enclosure for the saddles and minima. (II) Introduce an equivalent polynomial vector field. This step was already fully detailed in Section 3.2.
(IIIa) Compute and validate the local unstable manifold of the saddles. We use here the method outlined in Section 3.3.2 and detailed in Section 4. First, we obtain a validated enclosure for the unstable eigenvalue of each saddle, and for an associated eigenvector (in the extended phase space): Then, we fix a scaling for the eigenvector, namely γ 1 = 5 for the first one and γ 2 = 15, and for each saddle point solve recursively the invariance equation (27) for a finite number of coefficients. Using interval arithmetic, we compute N 1 = 20 (resp. N 2 = 30) coefficients for the parameterization of the unstable manifold of Sad 1 (resp. Sad 2 ). The values of the coefficients can be found in the file coeff para. Finally, we use Theorem 4.15 to validate these two truncated parameterization. As explained in Remark 4.9, we use Proposition 4.2 to find weights η ∈ R 6 >0 giving an optimal Z 1 bound: .
With these parameters we manage to validate both manifolds, that is we find validation radiuses r 1 = 1.024225153462953 × 10 −17 andr 2 = 1.975153579406591 × 10 −18 for which the radii polynomials P and Q are negative. See Figure 3 for a representation of the validated manifolds in (the original two-dimensional) phase space.
(IIIb) Validate trapping regions around the minima. This step is more straightforward. Using the procedure described in Section 3.3.1, we prove using Intlab that for each of the three minimum, the square of length 0.02 centered at the minimum is included in the stable manifold of the minimum. See Figure 3 for a representation of this trapping regions in phase space.
(IIIc) Compute and validate orbits connecting unstable manifolds and trapping regions. First, we numerically compute orbits starting from the end of the validated unstable manifold, to find after which integration time the orbit enters the trapping region validated at the previous step. This yields τ 11 = 0.015, τ 12 = 0.025, τ 22 = 0.012 and τ 23 = 0.018. In all this paragraph, indices ij are used to denote a quantity related to the orbit from Sad i to Min j . Then, we compute a piece-wise Chebyshev series representation of these approximate solutions. We use a uniform subdivision of size M 11 = M 12 = M 22 = M 23 = 10, and on each subdivision use K 11 = 30 (resp. K 12 = K 22 = K 23 = 20) Chebyshev coefficients. The values of the coefficients can be found in the file coeff orbit. Then, we use Theorem 5.10 to validate these four approximate orbits. As explained in Remark 5.7, we use Proposition 4.2 to find weights η ∈ R 6M >0 giving an optimal Z 1 bound. The exact value of these weights can be found in the file weights orbit. With the weights ν 11 = 1.4, ν 12 = 1.5, ν 22 = 1.7 and ν 23 = 1.5 for the 1 ν norms, we manage to validate the four orbits, that is we find validation radiusesρ 11 = 8.568033048134257 × 10 −11 , ρ 12 = 1.493891911294753×10 −11 ,ρ 22 = 5.071535766518771×10 −10 andρ 11 = 6.190811910154036×10 −10 for which the radii polynomials P (m) and Q (m) are all negative. See Figure 3 for a representation of the validated orbits in (the original two-dimensional) phase space. Finally, as suggested by Figure 3, we check using interval arithmetic and taking into account the validation radiuses ρ ij that each orbit indeed ends in the validated trapping region, which concludes the validation of each connecting orbit, and thus of the MEP. All the validation estimates can be obtained by running the script script main available at [10].