A numerical method for the approximation of stable and unstable manifolds of microscopic simulators

We address a numerical methodology for the approximation of coarse-grained stable and unstable manifolds of saddle equilibria/stationary states of multiscale/stochastic systems for which a macroscopic description does not exist analytically in a closed form. Thus, the underlying hypothesis is that we have a detailed microscopic simulator (Monte Carlo, molecular dynamics, agent-based model etc.) that describes the dynamics of the subunits of a complex system (or a black-box large-scale simulator) but we do not have explicitly available a dynamical model in a closed form that describes the emergent coarse-grained/macroscopic dynamics. Our numerical scheme is based on the equation-free multiscale framework, and it is a three-tier procedure including (a) the convergence on the coarse-grained saddle equilibrium, (b) its coarse-grained stability analysis, and (c) the approximation of the local invariant stable and unstable manifolds; the later task is achieved by the numerical solution of a set of homological/functional equations for the coefficients of a polynomial approximation of the manifolds.


Introduction
The computation of invariant manifolds of dynamical systems is very important for a series of system-level tasks, particularly for the bifurcation analysis and control.For example, the detection of stable manifolds of saddle points allows the identification of the boundary between different basins of attraction, while the intersection of stable and unstable manifolds most often leads to complex dynamical behavior such as chaotic dynamics [10,49].Their computation is also central to the control of nonlinear systems and especially in the control of chaos [11,33,39,49].However, their computation is not trivial: even for relatively simple low-dimensional ODEs, their analytical derivation is most of the times an overwhelming difficult task.Thus, one has to resort to their numerical approximation.However, this task is not easy; at the beginning of the 1990s only one-dimensional global invariant manifolds of vector fields could be computed.Guckenheimer and Worfolk [18] proposed an algorithm for converging on the stable manifold of saddles based on geodesics emanating from the saddle by iteratively rescaling the radial part of the vector field on the submanifold.Johnson et al. [20] introduced a numerical scheme to reconstruct twodimensional stable and unstable manifolds of saddles.The proposed method starts with the creation of a ring of points on the local-linear eigenspace and successively creates circles of points that are then connected by a triangular mesh.The appropriate points are selected through time integration so that the velocity of the vector field is similar in an arc-length sense for all trajectories.Krauskopf and Osinga [26] developed a numerical method based on geodesics; the manifold is evolved iteratively by hyperplanes perpendicular to a previous detected geodesic circle.Krauskopf et al. [27] addressed a numerical method for the approximation of two-dimensional stable and unstable manifolds which incorporates the solution of a boundary value problem; the method performs a continuation of a family of trajectories possessing the same arc-length.For a survey of methods for the numerical computation of stable and unstable manifolds see also Krauskopf et al. [27].In the above methods, the stable manifold is computed as the unstable manifold of the inverse map, i.e., by following the flow of the vector field backward in time [8].Thus, an explicit knowledge of the vector field and its inverse is required which however is not always available.England et al. [8] presented an algorithm for computing one-dimensional stable manifolds for planar maps when an explicit expression for the inverse map is not available and/or even the map is not invertible.Triandaf et al. [47] proposed a procedure for approximating stable and unstable manifolds given only experimental data based on time-delay embeddings of a properly selected data set of initial conditions.Another approach to compute invariant manifolds, the so-called parametrization method has been introduced by Cabre et al. [5][6][7].This is a numerical-assisted approach based on functional analysis tools for deriving analytical expressions of the local invariant manifolds.This involves the expansion of the invariant manifold as series and the construction of a system of homological equations for the coefficients of the series.Based on this approach, Haro et al. [19] addressed a numerical approach for the computation of the coefficients of high-order power series expansions of parametrizations of two-dimensional invariant manifolds.Breden et al. [4] employed the parametrization method to compute stable and unstable manifolds of vectors fields.For the implementation of the method it is assumed that the vector field is explicitly available in a closed form.Finally, focusing on singularly perturbed systems, Zagaris et al. [51] and Kristiansen et al. [28] extended the Computational Singular Perturbation (CSP) algorithm of Lam and Goussis [17,29] to approximate stable and unstable fiber directions on the slow manifold.
However, for many complex systems of contemporary interest, the equations that can describe adequately the dynamics at the macroscopic-continuum scale are not explicitly available in the form of ODEs or PDEs in a closed form.Take for example the case where the laws that govern the dynamics of the interactions between the subunits that constitute the complex system may be known in the form of, e.g., molecular dynamics, Brownian dynamics, agent-based modeling, and Monte Carlo, but a macroscopic description (ODEs or PDEs) is not available in a closed form.The lack of such a macroscopic description hinders the systematic numerical analysis, optimization and control.In general, two paths are traditionally followed for the numerical analysis of the emergent dynamics.The first one is the simple bruteforce simulation in time.An ensemble of many initial condition configurations would be set up; then a large enough number of microscopic runs would be performed; some of the parameters would be modified and finally the statistics of the detailed simulations would be computed.However, this practice is not appropriate for the systematic numerical analysis (for example one cannot find saddles with temporal simulations).The second path follows the statistical-mechanics/assisted approach where one aims at extracting closures between the moments of the microscopic distributions.For example, for Monte Carlo Markovian models a Master equation is usually derived for a few moments of the underlying probability distribution.However, these equations usually involve higher order moments whose evolution dynamics are functions of higher order moments.This leads to an infinite hierarchy of evolution equations and at some point these higher order moments have to be expressed as functions of the lower order moments in order to "close" the system of equations.However, the assumptions that underlie such "closures" may introduce biases to the analysis of the "actual" dynamics (see for example the discussion in [36]).
The equation-free approach addressed by Kevrekidis et al. [24,25,31,45], a multiscale numerical-assisted framework allows the bridging between traditional continuum numerical analysis methods, and microscopic/stochastic simulation of complex/multiscale systems bypassing the explicit derivation of moment closures.The equation-free approach identifies "on-demand" the necessary quantities for the numerical analysis of the emergent dynamics and has been used for the bifurcation analysis, control, optimization, rare-events analysis of a wide range of microscopic models and problems.Regarding the computation of coarse-grained invariant manifolds, Gear and Kevrekidis [14] introduced a method for the convergence on the coarse-grained slow manifolds of legacy simulators by requiring that the change in the "fast" variables (i.e., the variables that are quickly "slaved" to the variables that parametrize the slow manifold) is zero.Gear et al. [13] computed coarse-grained slow manifolds by restricting the derivatives of the "fast" variables to zero.Zagaris et al. [50] performed a systematic analysis of the accuracy and convergence of equationfree projection to the slow manifold.Frewen et al. [12] traced within the equation-free framework two-dimensional slow manifolds to get out of potential wells.Finally, Quinn et al. [34] have exploited the concept of equation-free approach to compute the one-dimensional stable manifold of a one-dimensional delay differential equation.
Here, building up on a previous work for the computation of coarse-grained center manifolds of microscopic simulators [43], we present a new numerical method based on the equation-free approach for the computation of coarse-grained stable and unstable manifolds of saddles of microscopic dynamical simulators (and in general large-scale discrete-time black-box simulators).The approximation of the coarsegrained stable and unstable manifolds is achieved using a polynomial approximation; the coefficients of the polynomials are computed iteratively by a Newton-Raphson scheme applied on a coarse-grained map of the microscopic simulator.Thus, the proposed numerical method involves a three-step procedure for the equation-free: (a) detection of the coarse-grained saddle, (b) computation of the coarse-grained Jacobian on the saddle and the computation of the corresponding eigenmodes, (c) identification of the coefficients of the polynomials that provide an approximation of the coarse-grained stable and unstable manifolds; this step involves: (i) the numerical construction of a back-box coarse-grained map for the coefficients of the polynomial approximation, and (ii) the iterative estimation of the polynomial coefficients by applying Newton's method around the constructed coarse-grained map.The method is illustrated through two examples whose stable and unstable manifolds are also approximated analytically.The first example is a simple toy discrete-time map and the second one is a Gillespie-Monte Carlo realization of a simple catalytic reaction scheme describing the dynamics of CO oxidation on catalytic surfaces.This is the first work that addresses the equation-free computation of both coarse-grained stable and unstable manifolds of saddles of microscopic simulators.

Computation of stable and unstable manifolds of saddles for discrete-time models
We will first present the basis for approximating the local stable and unstable manifolds of a saddle point of discrete-time systems when the equations are given explicitly.Then, in Section 3, we will show how this can be exploited within the equation-free framework for the approximation of coarse-grained stable and unstable manifolds, when equations are not given explicitly, thus when one has a large-scale black-box simulator.The later case includes both large-scale black-box simulators of let's say systems of ODEs and particularly microscopic/stochastic multiscale models, which is the focus of this work.
Let us begin by considering a discrete-time model given by: where F : R n × R p → R n is a smooth multivariable, vector-valued time-evolution operator that takes as initial condition at time t k = (kT ), x k ∈ R n and after some time horizon (sampling time) T reports the evolved state x k+1 , ∈ R n at time t k+1 = (k + 1)T ; p ∈ R p denotes the vector of parameters.
Regarding the computation of the local stable and unstable manifolds of a saddle point of the above system, the following theorem can be easily proven: Theorem 1 Let us denote by (x * , p * ) a saddle point of the discrete-time model (1) which satisfies x * = F (x * , p * ).Let us also assume that the Jacobian ∇F (x * , p * ) is diagonalizable/block-diagonalizable.Let V s be the n × l matrix whose l columns are the l eigevectors of ∇F (x * , p * ) that correspond to the l eigenvalues lying inside the unit circle, and V u be the n × n − l matrix whose columns are the eigevectors of ∇F (x * , p * ) that correspond to the n − l eigenvalues lying outside the unit circle.Let us also define z s ∈ R l and z u ∈ R n−l by the transformation where x k = x k − x * .Then the fixed point x * = 0 has: (A1) a C r l-dimensional local stable manifold W s (0) tangent to the subspace spanned by the columns of V s at the origin defined by: where tangent to the subspace spanned by the columns of V u at the origin defined by: where h u : R n−l → R l is a C r function which satisfies h u (0) = 0 and (B1) On the stable manifold the following system of functional equations hold: (B2) On the unstable manifold the following system of functional equations hold: In the above, (5) and (6) provide the equations of the invariant manifolds in diagonalized coordinates.Λ s is the l × l (block) diagonal matrix containing the l eigenvalues with |λ i | < 1 and Λ u is the (n − l) × (n − l) (block) diagonal matrix containing the (n − l) eigenvalues with |λ i | > 1; g s and g u are l and (n − l) vector-valued functions, respectively, obtained by where g(x , x * , p * ) corresponds to the n-vector-valued nonlinear function: containing all, but the linearization around the saddle, nonlinear terms of F (x k , p), satisfying: Proof Equations ( 5) and ( 6) are derived by re-writing (1) using the transformation given by (2) as: and Then, by taking (3), Eq. ( 10) reads: Finally, by (11) and ( 9) we obtain: with h s (0) = 0.In a similar manner, it can be shown that ( 6) in (B2) holds true on the unstable manifold.

Local approximation of the stable and unstable manifolds with truncated polynomial sequence
As by Theorem 1, the local stable and unstable manifolds are smooth nonlinear functions of z s and z u , respectively, then according to the Stone-Weierstrass theorem [37] they can be approximated by any accuracy around (x * , p * ) by a sequence of polynomial functions of z s and z u , respectively.For example, for the stable manifold (and similarly for the unstable manifold), ∀z uj = h sj (z s ), j = 1, 2, . . .n − l we can write: where P k i , i = 1, 2, ..l are polynomials (e.g., Chebyshev polynomials) of degree k i .
Truncating at a degree, say, M, we get the following truncated polynomial approximation: A simple choice would be to take as polynomials the power series expansion of z s .For example, if l = 2, M = 2, ( 14) results to an expression with (M + 1) l = 9 terms (including the constant term which under the above formulation should be zero): The existence of a local analytic solution for the form of the nonlinear functional ( 14) is guaranteed by the following theorem [46] (see also [21]): Theorem 2 Consider the following system of nonlinear functional equations: where φ : R n → R m is an unknown function.Then if: 1. f : R n → R n , w : R n × R m → R m are analytic functions such that f (0) = 0 and w(0, 0) = 0. 2. The function φ admits a formal power series solution.3. The fixed point that satisfies f (0) = 0 is a hyperbolic point, i.e., none of the eigenvalues of the Jacobian ∇ z f (z = 0) is on the unit circle.
where ẑs = {ẑ s1 , . . .ẑsl } are nonlinear functions of z s = {z s1 , . . .z sl }: Note that in general, both the left-hand side and the right-hand side of (17) contain higher order terms than M due to (5) and the nonlinearities in g uj .By equating both sides of ( 17) the terms up to an order r <= M with respect to {z s1 , . . .z sl }, we get the following coupled system of nonlinear homological equations with respect to the (n − l) The above system constitutes a nonlinear (in general) system of (n − l) • (r + 1) l unknowns with (n − l) • (r + 1) l equations that can be solved iteratively, e.g., using the Newton-Raphson algorithm.
What is described above, forms the basis for the numerical algorithm presented in Section 3 for the equation-free computation of coarse-grained stable and unstable manifolds of saddles of microscopic simulators.

An illustrative example
Let us consider the following discrete dynamical system: The above system can be written as: It can be shown that the local stable manifold of the saddle (0, 0, 0) is given by Let us choose a power series expansion up to order two (i.e., M = 2) of the stable manifold around the fixed point Hence, an approximation of the stable manifold is given by: Here, V = ⎡ ⎣ 1 0 0 0 1 0 0 0 1 ⎤ ⎦ .Hence, from (5) we get: or or Thus, from (22) we have: or By equating the coefficients of the corresponding power series up to order two, we get the following system of equations: From the above system we get: Thus, an approximation of the stable manifold around the saddle point is given by: It can be shown, that the unstable manifold is the trivial solution x 1 = 0, x 2 = 0 as follows.Let us again choose a power series expansion up to order two (i.e., M = 2) of the unstable manifold around the saddle.Hence, an approximation of the stable manifold is given by: Hence, from (6) we get: u (x 3 ) h (2)  u (x 3 ) For the above system, it can be easily verified that the unstable manifold is the one with x 1 = 0, x 2 = 0.

Numerical approximation of the stable and unstable manifolds of microscopic-stochastic multiscale and black-box simulators
Let us now assume that explicit model equations (such as the ones given by ( 1)) for the macroscopic (emergent) dynamics are not available in a closed form.Under this hypothesis, we cannot follow the procedure for the analytical approximation of the invariant manifolds as one needs to explicitly know the operator F (i.e., g s and g u in ( 7)).Thus, when explicit macroscopic equations are not available in a closed form, but a microscopic dynamical simulator is available, the approximation of the invariant stable and unstable manifolds at the macroscopic (the coarse-grained) level requires: (a) the bridging of the micro and macro scale, and (b) the numerical approximation of the coarse-grained manifolds.In what follows, we address a new multiscale numerical method for the approximation of the local stable and unstable manifolds based on the equation-free framework.Thus, let us assume that we have a microscopic (such as a Brownian dynamics, Monte Carlo, molecular dynamics, agent-based) computational model that, given a microscopic/detailed distribution of states at time t k = kT U , will report the values of the evolved microscopic/detailed distribution after a time horizon T U : is the time-evolution microscopic operator, p ∈ R p is the vector of the parameters.
A key hypothesis for the implementation of the equation-free numerical framework is that after some time t T U the emergent macroscopic dynamics can be described by a few observables, say, x ∈ R n , n N. Usually these "few" observables are the first moments of the underlying microscopic distribution.This implies that there is a slow coarse-grained manifold that can be parametrized by x.The hypothesis of the existence of a slow coarse-grained manifold dictates that the higher order moments, say, y ∈ R N−n , of the microscopic distribution U become, relatively fast over the macroscopic time, functions of the n lower order moments.At the moments-space, this dependence can be written as a singularly perturbed system of the form: where > 0 is a sufficiently small number.Under the above description and assumptions, Fenichels' theorem for continuous systems [9] can be extended to the following theorem that guarantees the existence of an invariant low-dimensional "slow" manifold on which evolve the coarse grained dynamics of the discrete system (35) (see also [3]).

Theorem 3 Let us assume that the functions
Then the dynamics of the system given by (35) can be reduced to: on a smooth manifold defined by: The manifold M is diffeomorphic and O( ) close to the M 0 manifold defined for = 0.Moreover, the manifold M is locally invariant under the dynamics given by (35).M defines the "slow" manifold on which the dynamics of the system evolve after a short (in the macroscopic scale) time horizon.Under this perspective and under the above theorem assumptions, let us define the coarse-grained map: where F T : R n ×R p → R n is a smooth multivariable, vector-valued function having x k as initial condition at time t k = kT , where T is a macroscopic reporting time horizon with T >> T U .The above coarse-grained map, which describes the system dynamics on the slow coarse-grained manifold M can be obtained by finding χ that relates the higher order moments of the microscopic distribution U k to the lower order moments x that describe the macroscopic observations/dynamics.The equation-free approach, through the concept of the coarse timestepper, bypasses the need to extract such a relation analytically which in most of the cases is an "overwhelming" difficult task and can introduce modeling biases (see the critical discussion in [36]).The equationfree approach provides such relations in a numerical way: relatively short calls of the detailed simulator provide this closure (refer to [24,25,45] for more detailed discussions).Briefly, the coarse timestepper consists of the following basic steps: Given the set of the macroscopic variables at time t 0 : (a) Set the coarse-grained initial conditions x k=0 ≡ x 0 .(b) Transform the coarse-grained initial conditions to consistent microscopic distributions U 0 = μx 0 , where μ is a lifting operator.(c) Run the microscopic/detailed simulator for a short macroscopic interval T to get the resulting microscopic distributions U k+1 .The choice of T is associated with the (estimated) gap of the eigenspectrum of the Jacobian of the unavailable closed macroscopic equations around the stationary state.(d) Compute the values of the coarse-grained variables using a restriction operator M: x k+1 = MU k+1 .
The above steps constitute the black box coarse timestepper, which, given an initial coarse-grained state of the system {x k , p}, at time t k = kT will report the result of the integration of the microscopic rules after a given time-horizon T (at time t k+1 ), i.e., x k+1 = F T (x k , p).
At this point, one can use iterative linear algebra numerical methods such as the Newton-Raphson method (for low-order systems) to converge to the coarse-grained fixed points and to perform bifurcation analysis.For large-scale systems, one can also resort to matrix-free methods such as the Newton-GMRES [23] to find the coarse-grained fixed points and the Arnoldi method [38] to analyze the stability of the coarse-grained fixed points of the unavailable macroscopic evolution equations.
The coarse-grained Jacobian ∇F T (x * , p * ) can be computed by appropriately perturbing the coarse-grained initial conditions fed to the coarse timestepper (38).For low to medium dimensions, the i-th column of the Jacobian matrix can be evaluated numerically, e.g., using central finite differences as: where e i is the unit vector with one at the i-th component and zero in all other components.Then, one can solve the eigenvalue problem with direct solvers.The tracing of solution branches of saddles around turning points can be achieved by standard continuation techniques such as the pseudo-arc-length continuation [22].
For the above procedure to be accurate, one should perform the required computations when the system lies on the slow manifold.If the gap between the fast and slow time scales is very big then the time required for trajectories starting off the slow manifold to reach the slow manifold will be small compared to T ; hence, we expect that the coarse-grained computations will not be significantly affected for any practical means.As also discussed in Kevrekidis et al. [25], under the "strong assumption" of a big separation of time scales, the fast off-the slow manifold dynamics will lead to quick "healing" of the lifting error.Nevertheless, one can relax the above "strong" (and "vague" [25]) assumption and enhance the computing accuracy by producing lifting operators that bring the system on the slow manifold (using for example the algorithms presented in [13,14,32,44,50]).
Returning back to the problem of numerical approximation of the stable and unstable manifolds, as now there are no analytical expressions, due to the nonlinear dependence of g s and g u on z s and z u , the coefficients a of the polynomial approximation of the stable z u = h(z s ) (and correspondingly of the unstable) manifold can be determined numerically by solving in an iterative way the system of nonlinear homological equations given by (19).
Here, we solve the above task through the concept of the coarse-timestepper as described in the following steps (in what follows, we present the algorithm for the computation of the stable manifold; the procedure for the computation of the unstable manifold is analogous): 1. Construct the coarse-timestepper given by the map (38) using appropriate lifting μ and restricting M operators of the microscopic evolved distributions.2. "Wrap" around the coarse-timestepper a continuation technique (e.g., the pseudo-arc-length continuation) to converge to a saddle fixed point (x * , p * ). 3. Compute the coarse-grained Jacobian ∇F T (x * , p * ) and solve the eigenvalue problem ∇F T (x * , p * )V = ΛV .Find the l stable and n−l unstable eigenmodes.Rearrange V as V = V s V u with V s being the n × l matrix whose columns are the eigevectors of the Jacobian that correspond to the l eigenvalues lying inside the unit circle, V u is a n × n − l matrix whose columns are the eigevectors of the Jacobian that correspond to the n − l eigenvalues lying outside the unit circle.4. Choose a certain set of polynomials as well as their maximum order M for the numerical approximation of the j-th element, say h js of the stable manifold in the form of: where, the variables z s , z u are defined by the transformation (2). 5. Let us denote with q (j) , j = 1, 2, . . .n−l, the column vector of dimension ((M + 1) l − 1) × 1 containing the unknown polynomial coefficients a (j) k 1 ,k 2 ,...,k l of the j-th element (h js ) of the stable manifold as approximated by (41).Set an initial guess at time t = 0, for each one of the unknown coefficients of the (expansion) contained in q (j) , say q (j) ] , i.e., the column vector of dimension ((M + 1) l − 1) • (n − l) containing all the unknown coefficients.6. Create a grid of n p initial conditions z s,0 within a certain distance B around (z s = 0, p * ) where an approximation of the stable manifold is sought, and at a certain distance from it, i.e., d < z s,0 < B. 7. Set convergence tolerance, say, tol, for the approximation of the polynomial coefficients.Set the number of time steps k max for calling the timestepper (38) and

End For
Based on all n p simulations construct the matrix A that contains all the values of each one of the polynomials P k 1 (z 1s ) • • • P k l (z ls ), k 1 , k 2 , . . .k l = 0, 1, . . .M. Thus, the above procedure will result to a matrix A of dimension (n p • (k max + 1)) × ((M + 1) l − 1) (the constant terms are set equal to zero).
Find q (j) k max , (j = 1, 2, . . ., n − l) , by solving the linear least squares problem: where, b (j) = z (1)  uj,0 z (1)  uj,1 . . .z (1)  uj,k max . . .z The optimal solution of the above linear least squares problem is given by the solution of A Aq (j) If the matrix A A is of full rank, then the above system has a unique solution given by: q (j) Note that if the initial points x (i) 0 are chosen close enough to the fixed point x * and/or the number of time-steps k max are relatively large then as z s → 0, the matrix A A will not be of full rank.In that case one could use the Moore-Penrose pseudo-inverse of A A to solve (44) and the solution reads: where, the pseudo-inverse matrix A + that is obtain by the Singular Value Decomposition (SVD) of the matrix A: where, Σ + is the inverse of sub-block diagonal matrix containing the nonzero singular values of the SVD decomposition of A.
k max , q (2) Thus the above procedure creates the map: • Update the values of the polynomials coefficients by using, e.g., the Newton-Raphson scheme around the map given by ( 42): -Compute the Jacobian ∇Ξ (Q 0 ) by perturbing appropriately Q 0 (thus, repeating the pipeline described in step 8 above, ((M + 1) l − 1) • (n − l) times, thus starting simulations using Q 0 ± e i (where e i , i = 1, 2, . . ., ((M+1) l −1)•(n−l) is the unit vector with one at the i-th component and zero in all other components of Q 0 ) for an approximation of the Jacobian Ξ (Q 0 ) with central differences).-Solve the linear system to get dQ.-Update the estimates for the polynomial coefficients: End Do }.
Note, that the above scheme can be seen as a Gauss-Newton-like scheme for the estimation of the unknown coefficients of the series expansion by the aid of the coarsetimestepper.While the existence of a unique solution is on the one hand guaranteed by Theorems 1 and 2, on the other hand, the convergence of the scheme depends on the choice of the distance B around the saddle point, where the approximation of the invariant manifolds in sought, the initial guesses for the polynomial coefficients, the level of noise due to the microscopic (stochastic) simulations, as well as the "quality" of the lifting operator (see the discussion for equation-free computations, e.g., in [13,14,25,31,32,42,44,50,51]).An initial guess for the polynomial coefficients could be the values of the coefficients of the linear terms of the polynomial expansion that reconstruct the subspace spanned by the columns of V s or V u .The convergence properties of the above scheme with respect to the above factors will be studied more thoroughly in a future work.

The illustrative examples: numerical results
The proposed approach is illustrated through two examples: (a) the toy model ( 20) and (b) a Monte Carlo simulation of a catalytic reaction on a lattice presented in [31], for which we have also derived analytically an approximation of the stable and unstable manifolds based on the mean field model.

The toy model
We have shown that the stable manifold of the discrete model given by ( 20) is given by Here, we will derive a numerical approximation of the stable manifold by assuming that the equations of the model are not explicitly known.Our assumption is that we have a black-box model that given initial conditions (x 1,0 , x 2,0 , x 3,0 ) it outputs 3 ) = (0, 0, 0).The Jacobian on the saddle is approximated by central finite differences with = 0.01 as perturbation on the initial conditions and running the simulator for one step.By doing so, the numerical approximation of the Jacobian actually coincides for any practical means with the analytical one.The eigenvalues are λ 1 = −0.5, λ 2 = −0.5, λ 3 = 2 and the eigenvectors are given by e i , i.e., the unit vectors with one at the i-th component and zero in all other components.From the above, it is clear that z 1s = x 1 , z 2s = x 2 , z 1u = x 3 .Thus, we chose a power series expansion of the manifold around the saddle as For the construction of the map (see ( 42)) we have used the following parameters: k max = 3, n p = 4, z s1,0 = (−0.1,0.1), z s2,0 = (−0.1,0.1), and central finite differences with = 0.01, for the numerical approximation of the Jacobian ∇Q that is required for the Newton-Raphson iterations; the tolerance was set to tol = 1E-04, and the initial guess of the power expansion coefficients was set as q 0 ≡ (a 0,1 , a 0,2 , a 1,0 , a 1,1 , a 1,2 , a 2,0 , a 2,1 , a 2,2 )=(0,0,0,0,0,0,0,0).
The algorithm converged in two Newton iterations to the following expression for the stable manifold: The numerical scheme outputs also a non-zero coefficient for a 2,2 which is not present in the analytical approximation.This is due to the truncation of the power expansion to second-order terms: when equating the terms on both sides of ( 27), higher order powers than three are set to zero.One can confirm the contribution of this extra term found by the numerical scheme by simple simulations.For example, by setting as initial conditions x 1,0 = 0.2, x 2,0 = 0.2 and x 3,0 = − 4 7 x 2 2,0 + 32 119 x 2 1,0 x 2,0 , we get the results shown in Table 1.Note that x 3,k goes to zero and then after k = 3 it diverges due to the (truncated) approximation of the manifold.This is what should be expected: the lower order model results in a trajectory that is a little bit off the stable manifold; thus at the beginning, the system approaches the saddle (see the first two time steps) and then diverges away from it as expected.In fact, if we add the extra term found with the numerical scheme, and start with the same initial conditions for x 1,0 and x 2,0 , but with  2, that is after k = 3 the system converges to the saddle.

Kinetic monte carlo simulation of CO oxidation on a catalyst
The proposed numerical approach is illustrated through a kMC microscopic model describing the dynamics of CO oxidation on a catalyst [31].The species react adsorbed or desorbed on a finite lattice with periodic boundary conditions.At each time instant, the sites of the lattice are considered to be either vacant or occupied by the reaction species.The system dynamics are described by the following chemical master equation: where, P (x, t) is the probability that the system will be in state x at time t and Q(y, x) is the probability for the transition from state y to x (and vice versa) per unit time.The summation runs over all possible transitions (reactions).Here, the numerical simulation of the above stochastic equation was realized using the Gillespie kMC algorithm [15,16,31].The reaction mechanism can be schematically described by the following elementary steps (for more details see in [30]): ( where i, j are neighbor sites on the square lattice, * denotes a site with a vacant adsorption site, while "ads" denotes adsorbed particles.By adding an inert siteblocking adsorbate with a reversible adsorption step the mean field approximation can be derived by the master (54) and is given by the following system of ordinary differential equations [31]: For the time integration of the mean field model (55), we used the Matlab ode suite function "ode15s" [40], a variable-step, variable-order solver with absolute and relative tolerances set to 1E-03.For the kMC simulations the number of the sites (system size) and the number of runs-consistent to the mean values of the distribution on the lattice realizations-were chosen to be N size = 2000×2000 and N r = 2000, respectively.The value of the time horizon was selected as T = 0.05 (for a discussion on the appropriate choice of the reporting horizon T and as well as the system size and time realization please refer to [31]).The coarse-grained bifurcation diagram was obtained by applying the equation-free approach upon convergence of the Newton-Raphson to a residual of O(10 −3 ) for ≈ 10 −2 .We have chosen this model as for big enough lattice realizations and runs, the coarse-grained bifurcation diagram and stability practically coincides with the one obtained from the mean filed model [31]; thus one can perform a direct comparison of the numerical approximation of the stable manifold obtained with the kMC simulator and the one derived analytically from the mean field model.
Here, we have chosen to find the stable and unstable manifolds at β = 20.7.For this value of the bifurcation parameter and the choice of the values of T , N r and N, the coarse-grained fixed point is found to be

Numerical approximation of the stable manifold
For our illustrations, we have chosen to provide a third-order approximation of the stable manifold for the mean-field model (55) and its corresponding kMC simulation.Thus we sought for an approximation of the stable manifold as: Following the approach described in Section 2, one obtains a nonlinear system of six algebraic equations (see Appendix) which was solved for the unknown coefficients with Newton-Raphson.Here, the convergence tolerance was of the order of 10 −6 , while the perturbation for computing the Jacobian matrix with finite differences was of the order of 10 −3 .The initial guess of the coefficients was set as q 0 = (0, 0, 0, 0, 0, 0) .In this case, the expression for the approximation of the stable manifold is given by: Next, we applied the proposed numerical method, treating the mean-field timestepper-resulting from the integration of the system of ODEs (55) with ode15s-as a black-box with a reporting time horizon T = 0.05.We have set n p = 6 and a grid of initial conditions z s,0 : {−0.02, −0.01, −0.005, 0.005, 0.01, 0.02} around the saddle and k max = 6.Again, the convergence tolerance was set to 10 −6 , while the perturbation for computing the Jacobian matrix with finite differences was of the order of 10 −2 .Starting the algorithm with zeros as initial guesses for all the unknown coefficients (i.e., by setting q 0 = {0, 0, 0, 0, 0, 0}), the algorithm results to the following expression: The scheme converged within 3 iterations (the sequence of the convergence criterion dq (j) ∼ was (after the first iteration) ∼ O(0.1), ∼ O(10 −7 )).By selecting a grid closer to the equilibrium as z s,0 : {−0.005, −0.003, −0.001, 0.001, 0.003, 0.005}, the algorithm results to the following expression: By comparing (57) and (59), we see that the numerical approximation of the coarse-grained manifold as obtained by "wrapping" the proposed algorithm around the black-box ODEs timestepper is almost identical with the one computed using Newton-Raphson for the solution of the analytically derived nonlinear homological equations (given in Appendix A.2).A relatively small deviation in the coefficients of the 3rd-order terms in (58) which was computed with the "wider" grid is due to the larger truncation error.Finally, we applied the proposed numerical approach to the kMC timestepper with N size = 2000 × 2000 , N r = 2000 and a reporting time horizon T = 0.05.In this case, the convergence tolerance for the kMC simulations was set 10 −3 , while the perturbation for computing the Jacobian matrix with finite differences was of the order of 10 −2 .We have set n p = 6 and a grid of initial conditions z s,0 : {−0.02, −0.01, −0.005, 0.005, 0.01, 0.02} around the coarse-grained saddle, and k max = 6.Starting again the algorithm with zeros as initial guesses for all the unknown coefficients (i.e., by setting q (0) = {0, 0, 0, 0, 0, 0}, the algorithm results to the following expression for the coarse-grained stable manifold: In this case, the scheme converged within 4 iterations (the sequence of the convergence criterion dq (j) was (after the first iteration) ∼ O(0.1), ∼ O(0.01), ∼ O(0.001)).The computation time for one Newton iteration was of the order of 50 min on a INTEL Xeon CPU E5-2630, 2.2GHz with 64GB RAM.
Taking N size = 4000 × 4000 , N r = 2000, and leaving everything else the same, the algorithm results to: In this case, the scheme converged again within 4 iterations (the sequence of the convergence criterion dq (j) was (after the first iteration) ∼ O(0.1), ∼ O(0.01), ∼ O(0.001)).At this point we should note that by taking larger lattice sizes N size and more runs N r , the convergence to the mean-field results will increase.On the contrary, for small sizes and number of runs, the algorithm fails to converge.For example if one takes, N size = 200 × 200, N r = 2000, the algorithm does not converge.
For the particular problem of kMC simulations, a more thorough discussion on the interplay between the level of noise, the selection of the amplitude of the perturbation for the estimation of the coarse-grained Jacobian, the "healing" assumption of the equation-free approach due to the lifting operation and the choice of the reporting horizon T can be found in [30,31] where the equation-free approach has been addressed to construct the coarse-grained bifurcation diagram of the same model.For example, if the time reporting horizon T is too short, the errors that are introduced due to the lifting do not have the time to "heal" down to the slow manifold.If on the other hand, one selects a reporting horizon T too long, then the quantification of the Jacobian is inaccurate, thus introducing biases in the computations model.

Numerical approximation of the unstable manifold
Here, we seek for the following approximation of the unstable manifold Following the approach described in the Appendix, we obtained analytically seven algebraic equations, which were solved for the unknown coefficients with Newton-Raphson; the convergence tolerance was set 10 −6 while the perturbation for computing the Jacobian matrices was of the order of 10 −2 .In this case, the approximation of the unstable manifold is given by: (63) Again, we first applied the proposed numerical method to the mean-field timestepperresulting from the integration of the system of ODEs (55) with ode15s-as a black-box-with a reporting time horizon T = 0.05.We have set n p = 4 and a grid of initial conditions z u1,0 , z u2,0 : {−0.02, −0.01, 0.01, 0.02} around the saddle setting k max = 6.Again, the convergence tolerance was set 10 −6 , while the perturbation for computing the Jacobian matrix with finite differences was of the order of 10 −3 .Starting the algorithm with zeros as initial guesses for all the unknown coefficients (i.e., by setting q 0 = {0, 0, 0, 0, 0, 0, 0}), the algorithm results to the following expression for the unstable manifold: The scheme converged within 3 iterations (the sequence of the convergence criterion dq (j) was ∼ O(10 −1 ), ∼ O(10 −2 ), ∼ O(10 −6 )).
By selecting a grid closer to the equilibrium as z u1,0 , z u2,0 : {−0.005, −0.003, 0.003, 0.005}, the algorithm results to the following expression for the unstable manifold: (65) By comparing (64) and ( 65) with (63), we see that the numerical approximation of the coarse-grained manifold as obtained by "wrapping" the proposed algorithm around the black-box ODEs timestepper is almost identical with the one computed using Newton-Raphson for the solution of the analytically derived nonlinear homological equations (given in Appendix A.2).
Finally, we applied the proposed numerical approach to the kMC timestepper with N size = 2000 × 2000 , N r = 200 and a reporting time horizon T = 0.05.In this case, the convergence tolerance for the kMC simulations was of the order of 10 −3 , while the perturbation for computing the Jacobian matrix with finite differences was of the order of 10 −2 .We have set n p = 4 and a grid of initial conditions z u1,0 , z u2,0 : {−0.02, −0.01, 0.01, 0.02} around the coarse-grained saddle setting k max = 6.Starting again the algorithm with zeros as initial guesses for all the unknown coefficients (i.e., by setting q (0) = {0, 0, 0, 0, 0, 0, 0}, the algorithm results to the following expression for the coarse-grained unstable manifold: In this case, the scheme converged within 3 iterations (the sequence of the convergence criterion dq (j) was ∼ O(0.1), ∼ O(0.01), ∼ O(0.001)).The computation time for one Newton iteration was of the order of 3.5 h on a INTEL Xeon CPU E5-2630, 2.2GHz with 64GB RAM.
Taking N size = 4000 × 4000 , N r = 2000, and leaving everything else the same, the algorithm results to: In this case, the scheme converged again within 3 iterations (the sequence of the convergence criterion dq (j) was ∼ O(0.1), ∼ O(0.01), ∼ O(0.001)).By comparing the expressions (66) and ( 67) with (63), we see that the numerical approximation of the coarse-grained manifold of the kMC simulator is in a fair agreement with the one obtained by the mean field model.

Conclusions
We propose a numerical method for the approximation of the local coarse-grained stable and unstable manifolds of saddle points of microscopic simulators when macroscopic models in a closed form in the form of ODEs are not explicitly available.The methodology is based on the equation-free multiscale framework.The proposed numerical algorithm consists of three steps: (a) detection of the coarse-grained saddle by constructing the coarse-timestepper of the microscopic dynamics, (b) estimation of the coarse-grained Jacobian and evaluation of its eigenvalues and eigenvectors, and (c) estimation of the coefficients of the polynomial approximation of the invariant manifolds.The later step involves the construction of a map for the coefficients of the polynomial expansion of the manifolds.The key hypothesis is that a macroscopic model in the form of ODEs can in principle describe the emerging macroscopic dynamics but it is not available in a closed form.The proposed numerical approach was illustrated through two examples, a toy model treated as a black-box time-stepper and a kinetic Monte Carlo simulator of a simple catalytic reaction.For the kMC simulator a mean field model in the form of ODEs was also given.For both models, we have also derived analytically the parametrization of the invariant manifolds for comparison purposes.As we show, the proposed numerical method approximates fairly well the analytical approximation when considering the vector fields as known.
The proposed numerical method provides a polynomial approximation of the local stable and unstable manifolds in a neighborhood of the coarse-grained saddle, based on appropriately initialized temporal simulations of the microscopic simulator.It should be noted that the task of finding even linear stable manifolds from brute-force simulations (and thus physical experiments) is itself a non-trivial task [1,47].In physical experiments the additional difficulty is that most often one cannot set the initial conditions at will; towards this aim, one could resort to control-based continuation methods [1,41,45].In a future work, we aim at extending the proposed numerical method to perform a piece-wise approximation of the global manifold.This could be done for example by coupling the proposed algorithm with parametercontinuation of the polynomial coefficients as we move far from the equilibrium.Of course the computation of global manifolds, especially of two-dimensional stable and unstable manifolds constitutes a much more difficult task.Towards this direction, Quinn et al. [34] have exploited the equation-free approach to compute the one-dimensional stable manifold of an one-dimensional delay differential equation.Other points that require further investigation are the analysis of the convergence properties of the algorithm, the numerical convergence properties of the scheme with respect to the amplitude of the noise to stochasticity, the sensitivity of the approximation with respect to the discretization of the domain around the saddle as well as the issue of finding confidence intervals for the coefficients of the polynomial expansion.Another extremely interesting problem is that of equation-free uncertainty quantification (UQ) [2,35,48,52], an approach that can be exploited to deal with the inherent stochasticity of the microscopic simulations and to provide "on demand" the appropriate coefficients of generalized polynomial chaos expansions.

Appendix. Extraction of the stable and unstable manifolds for the mean field model of ODEs
Let us assume a continuous model in the following form of ODEs: where f is considered to be sufficiently smooth.
To determine the local stable and unstable manifold of a saddle fixed point (x * , p * ), the following linear transformation is introduced: where V is the matrix with columns the eigenvectors v j of the Jacobian ∇ x f (x, p) computed at (x * , p * ).As in Section 2 expanding the right-hand side of (A.1) around (x * , p * ) and introducing (A.2) we get: g(V z, p) contains the higher order terms with respect to x.By rearranging appropriately the columns of V , the Jacobian J ≡