EQUATION-FREE COMPUTATION OF COARSE-GRAINED CENTER MANIFOLDS OF MICROSCOPIC SIMULATORS

. An algorithm, based on the Equation-free concept, for the approximation of coarse-grained center manifolds of microscopic simulators is addressed. It is assumed that the macroscopic equations describing the emergent dynamics are not available in a closed form. Appropriately initialized short runs of the microscopic simulators, which are treated as black box input-output maps provide a polynomial estimate of a local coarse-grained center manifold; the coeﬃcients of the polynomial are obtained by wrapping around the microscopic simulator an optimization algorithm. The proposed method is demonstrated through kinetic Monte Carlo simulations, of simple reactions taking place on catalytic surfaces, exhibiting coarse-grained turning points and Andronov-Hopf bifurcations.


1.
Introduction. The computation of center manifolds is central in model reduction and control of nonlinear systems at non-hyperbolic equilibrium points [1,2,3,10,11,24]. Restricting the dynamics of large-scale continuous and discrete systems on the center manifold allows the stability analysis of the full system based on a low-scale system whose order is determined by the dimension of the generalized eigenspace associated to the eigenvalues with zero real part. Doing so, around non-hyperbolic points one can capture the dynamics of the full system with just a few ordinary differential equations. Thus, this approach has important implications in designing controllers for stabilizing systems possessing uncontrollable modes. However, the computation of center manifolds is not trivial: it involves the extraction of a set of functional equations from the full set of the continuum models (appearing in the form of Ordinary (ODEs) and/or Partial Differential Equations (PDEs)) and consequently their solution, which most of the times is not an easy task. Furthermore, for many contemporary complex systems, coarse-grained macroscopic equations are not available in a closed form. Instead, in many situations, the level at which the physics description might be available is a more detailed one: for example, the evolution rules may be known in the form of molecular dynamics, kinetic Monte Carlo, agent-based models etc. In these circumstances, conventional continuum algorithms cannot handle the given detailed information. Bridging systematically the gap between the microscopic and macroscopic space and time scales 378 CONSTANTINOS SIETTOS of a complex system is the key for the systematic analysis and deeper understanding of the emergent responses. Over the past few years, it has been demonstrated that the so-called Equation-free approach [14,15] can establish a link between traditional continuum numerical analysis and microscopic/ stochastic simulation of complex/multiscale systems. This is a computer-assisted methodology, inspired from continuum numerical analysis, system identification and large scale iterative linear algebra, that enables microscopic-level codes to perform system-level analysis directly, without the need to pass through an intermediate, coarse-grained, macroscopic-level, description of the emergent system dynamics. The backbone of the method is the "on-demand" identification of the quantities required for continuum numerical analysis (coarse residuals, the action of coarse slow Jacobians, eigenvalues, Hessians, etc). These are obtained by repeated, appropriately initialized calls of the detailed microscopic simulators, which are treated as black boxes.
The key assumption is that macroscopic, coarse-grained models exist and close for the expected behavior of a few macroscopic system variables. These coarse-grained variables are typically a few low moments of the microscopically evolving distributions (e.g. surface coverage or the first moments of population distributions). Based on this concept, Gear and Kevrekidis [5] addressed an approach to compute steady 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. Following up, Gear et al. [6] proposed a method for an Equation-free "on the fly" approximation of slow manifolds by restricting the (m+1) derivatives of the "fast" variables to zero. The authors show that if these (m+1)derivatives are zero then the solution lies on a m-th order approximate solution of the slow manifold.
Here, a new systematic approach based on the Equation-free concept for the computation of coarse-grained center manifolds of microscopic simulators is presented. The approximation of a local coarse-grained manifold is provided in a polynomial form whose coefficients are computed by wrapping around the microscopic simulator an optimization algorithm. This is a three-step procedure including the Equation-free (a) detection of the coarse-grained non-hyperbolic equilibrium, (b) stability analysis of the critical point, and (c) approximation of a local coarsegrained manifold by identifying the quantities required for the optimization algorithm (coarse residuals, Hessians, etc). This is the first time that such a procedure is addressed. The proposed approach is demonstrated through two illustrative kinetic Monte Carlo stochastic simulations of a simple surface reaction scheme describing the dynamics of NO oxidation by H2 on Pt and Rh surfaces.

2.
Mathematical and computational background: The Equation-free framework. Let us assume that due to the complexity of the underlying dynamics evolving across temporal and spatial scales, explicit model equations for the macroscopic (emergent) level are not available in a closed form. We furthermore assume that the description of the problem is plausible in a more detailed level, meaning that we have a microscopic/ individual-based/ stochastic computational model that, given a microscopic/ detailed state at time t k = kT U , will report the values of the evolved microscopic/detailed state after a microscopic time interval T U : C T U : R N × R m → R N is the time-evolution microscopic operator, p ∈ R m is the vector of the complex system parameters. We also assume that after some time t > T U the emergent coarse-grained dynamics are governed by a few variables, say, x ∈ R n , n << N . Usually these "few" observables are the first few moments of the underlying microscopic distribution. This implies that there is a slow coarsegrained manifold that can be parametrized by x. Under this perspective and under the assumptions of the Fenichel's theorem [4] let us define the coarse-grained map where F T : R n × R m → R n having x k as initial condition and T >> T U . The above coarse-grained map which describes the system dynamics on the slow coarse-grained manifold can be obtained by simply restricting U k to x k . The assumption of the existence of a slow coarse-grained manifold asserts that the higher order moments, say, y ∈ R N −n , of the microscopic distribution U become, relatively fast over T functionals of x. Under the implicit function theorem this relation can be represented as where χ is a smooth continuous function defining the relation between the slow x and the fast y variables after a short (in the macroscopic sense) time horizon. The Equation-free approach through the concept of the coarse timestepper provides such relations "on demand": relatively short calls of the detailed simulation provide this closure (refer to [14,15,18,20] 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) Prescribe the coarse-grained initial conditions x(t 0 ) ≡ x 0 . (b) Transform them through a lifting operator µ to consistent microscopic distributions U (t 0 ) = µx(t 0 ).
(c) Evolve these distributions in time using the microscopic/detailed simulator for a short macroscopic time T to get U (t 0 + T ). The choice of T is associated with the (estimated) spectral gap of the linearization of the unavailable closed macroscopic equations.
(d) Obtain again the values of the coarse-grained variables using a restriction operator M : The above steps, constitute the black box coarse timestepper, that, given an initial coarse-grained state of the system {x k , p}, at time t k 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). Now one can "wrap" around the coarse timestepper (3) computational superstructures like Newton-Raphson (for low-order systems) to perform coarse-grained fixed point computations and stability analysis. For large-scale systems one can employ matrix-free methods like Newton-GMRES [12] to find the coarse-grained fixed points and Arnoldi iterative algorithms [21] to estimate the dominant eigenvalues of the coarse linearization, that quantify the stability of the fixed points of the unavailable macroscopic evolution equations.
For the above procedure to be accurate, one should perform the required perturbations 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 very small compared to T ; hence the coarse-grained computations will not be affected for any practical means. Nevertheless, one can enhance the computing accuracy by forcing the system to start on the slow manifold (using for example the algorithms presented in [5] and [6]).
3. Equation-free computation of center manifolds of microscopic simulators. When explicit macroscopic equations are not available in a closed form, the proposed Equation-free approximation of the coarse-grained center manifold is achieved in three parts: (a) the detection of the co-dimension one bifurcation, (b) the stability analysis of the non-hyperbolic point, and (c) the approximation of the invariant manifold. These parts can be implemented as follows: Step 1. Drive the microscopic/detailed simulator to its own open-loop bifurcation.
Within the the Equation-free framework, the detection of the co-dimension one bifurcations can be achieved in an (a) indirect or (b) direct way. In the indirect way the bifurcation is detected during the continuation procedure using methods of interpolation. The direct way requires the solution of an augmented with the socalled criticality condition (or bifurcation/test function)system: where c(x, p) is the test function. For example for a co-dimensional one turning point bifurcation [23]: The second equation simply guarantees that the k − th component of q takes a nonzero value, i.e. we seek for a nontrivial solution.
For Andronov-Hopf bifurcations several test functions have been proposed [9,17,22]. For example one might use the following where 1 ± ωi are the eigenvalues and u ± wi are the associated eigenvectors.
Step 2. Find the eigenmodes of the coarse-grained Jacobian at the criticality. The coarse-grained Jacobian is computed by appropriately perturbing the coarsegrained initial conditions fed to the coarse timestepper (3) . The i − th column of the Jacobian matrix can be evaluated numerically as where e i is the unit vector with one at the i − th component and zero in all other components.
Let us denote by v j ∈ R n the j − th eigenvector corresponding to the j − th eigenvalue λ j of the coarse-grained Jacobian ∇ x F T (x, p) computed at the criticality, say (x * , p * ). Let us denote by V the matrix whose columns are the generalized eigenvectors v j . Let us also assume that l of these eigenvectors correspond to l eigenvalues with zero real parts.
Step 3. Compute the coarse-grained center manifold. 3.1 Discretize the domain D n ⊆ R n of the state-space around the criticality where a numerical approximation of the coarse-grained center manifold is sought, in a mesh of, say, Introducing (9) in (3) we get If F T is sufficiently smooth, it can be expanded around x * in a Taylor series and eq. (10) can be written as where f r (z k , p) are higher order polynomials. Eq.(11) can be written as where At this point is it convenient to write the coarse-grained Jacobian in a block form as where J c is the k × k matrix whose eigenvalues are the l eigenvalues with zero real parts. Hence eq. (12) can be re-written as z s,k+1 = J s z s,k + g s (z s,k , z c,k , p) (15) where z c ∈ R l and z s ∈ R n−l ; g c and g s are the projections of the nonlinear part of eq.(12) in R l and R n−l respectively.
Note that z c and z s are linearly uncoupled. 3.3 If a center manifold exists, there is a smooth function: satisfying h = 0 and ∇ zc h j = 0; h j is the j − th component of h [3,19]. Then, an approximation of (16) can be constructed as polynomial function of z c . 3.4 Substitute (16) in (15) 3.4 Substitute (14) in (17) to get 3.5 Find the polynomial coefficients of h such that eq.(18) is satisfied up to an order q. This can be achieved by wrapping around the coarse timestepper an optimization method such as the unconstrained algorithm of Broyden, Fletcher, Goldfarb, Shanno (BFGS) or the Nelder-Mead method [13].
Having detected the criticality (from step 1) and computed the associated coarsegrained eigenmodes (from step 2), Step 3 above can be summarized within the Equation-free framework by the following pseudo-code: Do while {Minimize the Objective function: % For example let us assume a 3-dimensional system with two imaginary eigenvalues % λ 1,2 = ±i and the third eigenvalue with non-zero real part; O(i)) with respect to the polynomial coefficients α i , i = 1, 2, ..., q End While At this point it should be noted that as dictated from the invariant manifold theory, with the above procedure one approximates a local center manifold of a family of l-dimensional smooth local center manifolds that satisfy (18) around the equilibrium; clearly the computed center manifolds are not necessarily unique ( [3]).

4.
The illustrative examples: Computation of coarse-grained center manifolds of kMC simulations. The proposed approach is illustratred through two kMC microscopic models [18] describing the dynamics of CO oxidation on a catalyst. The species react, are 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, y) is the probability for the transition from state y to x 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 [7,8]. The reaction mechanism can be schematically described by the following elementary steps: (1) CO gas + * i ↔ CO ads,i (2) O 2,gas + * i + * j ↔ O ads,i + O ads,j (3) CO ads,i + O ads,j → CO 2,gas + * i + * j where i, j are sites on the square lattice, * denotes a site with a vacant adsorption site, while "ads" denotes adsorbed particles. For this problem a first "good" mean field approximation (which we will call model 1) is given by the following system of equations [18]: where θ i represent the coverage of species (i = A, B, corresponding to CO and O) on the catalytic surface. The parameters α, β, γ are associated with CO adsorption, O dissociative adsorption and CO desorption rates, while k r is a reaction rate constant. For α = 1.6 , γ = 0.04, k r = 1 and treating β as the bifurcation parameter, the mean field model (20) where Adding an inert site-blocking adsorbate with a reversible adsorption step the mean field approximation is given by the following equations (which we will call model 2) [18]: where θ i represent the coverages of species (i = A, B, C, resp. CO, O and inert species C) on the catalytic surface; µ denotes C adsorption and η C desorption rate. For α = 1.6, γ = 0.04, k r = 1 , η = 0.016, µ = 0.36 and treating β as the bifurcation parameter the mean field model ( It can be shown (see Appendix) that a local center manifold that is tangent to the center eigenspace is approximated by For the kMC simulations the number of the sites (system size) and the number of -consistent to the mean values of the distribution on the lattice-realizations were chosen to be N size = 1000 × 1000 and N r = 2000, respectively. The value of the time horizon was selected as T = 0.05 . The coarse timesteppers for model 1, and Applying the proposed Equation-free approach with d = 0.2 and disretization step 0.002, the coarse-grained local center manifold was computed by wrapping around the coarse timestepper the Nelder-Mead method as The residual of the mean square errors was around 0.04. Figure 1 shows the computed coarse-grained center manifold. For the kMC simulation of the kinetics of the second model the coarse timestepper converged to the "correct" coarse-grained Andronov-Hopf bifurcation (θ * A , θ * B , θ * C β * ) 1 ; the computed coarse-grained eigenvalues and corresponding eigenvectors are: Applying the proposed Equation-free approach with d = 0.2 and discretization steps of 0.005 for θ A , 0.001 for θ B and 0.01 for θ C , the coarse-grained local center manifold was computed as The residual of the mean square errors was around 1.0. Figure 2 shows the computed coarse-grained center manifold.
As it is shown for both cases the computed eigenvalues, eigenvectors and center manifolds are in a good agreement with the ones derived analytically using the mean field ODEs.

5.
Conclusions. The purpose of this paper is to demonstrate how one can compute, under certain assumptions, an approximation of coarse-grained center manifolds of microscopic simulators when good macroscopic models in a closed form are not available. In particular, it is shown how the Equation-free approach can be exploited to construct coarse-grained center manifolds in a polynomial form whose coefficients are computed by wrapping around the microscopic simulators well established optimization algorithms. The proposed Equation-free approach consists of a three tier procedure: (a) detection of the coarse-grained non-hyperbolic equilibria, (b)"on the fly" stability analysis and (c) computation of the polynomial coefficients. The underlying key assumption is that there is a separation of time scales between the higher-order and lower-order moments of the evolving microscopic distributions. This implies that around the coarse-grained fixed points there is a gap between the critical eigenvalues (the eigenvalues with zero real part for non-hyperbolic fixed points) and the rest of the spectrum. The proposed approach was illustrated through two kinetic Monte Carlo simulations of a simple reaction taking place on catalytic surfaces; for these kMC models good macroscopic models exist in the form of ODEs and are used for comparison purposes. The corresponding bifurcation diagrams exhibit turning points and Andronov-Hopf bifurcations. The coarse-grained center manifolds computed with the proposed approach approximate well the ones computed by solving a truncated set of functional equations resulting from the expansion of the corresponding macroscopic models (see Appendix). In its current form the approach requires the discretization of the phase space in N 1 × N 2 × ... × N n points (n being the dimension of the coarse-grained system) as well as the computation of the eigevectors of the coarse-grained Jacobian matrix at the non-hyperbolic points. While this is the "way to go" for low to medium scale coarse-grained systems, analogous treatment for large-scale systems is computationally demanding and thus inefficient. When this is the case the proposed approach can be modified in an alternative way: one would not need to compute all the eigenmodes and project the system to its full eigenspace but with a iterative procedure such as the Arnoldi eigensolver provide estimations of the leading eigenmodes. Then, these critical eigevectors can be used to project the system into the critical eigenspace and its compelement making use of the Fredholm Alternative Theorem (see [17]). Another point that requires further investigation is the analysis of the convergence properties of the proposed algorithm with respect to issues such as the level of stochasticity and the initial guess for the polynomial coefficients. A related work concerning the convergence properties of the algorithm presented in [6] for the Equation-free projection on slow manifolds of systems with multiple scales can be found in [25]. Finally it should be noted, that the approach requires that a "good" set of coarse-grained observables (i.e. the level n at which the macroscopic unavailable equations would close and the macroscopic variables per se) is known beforehand. However, for many complex systems such knowledge is not available a-priori. In this case, one could integrate the Equation-free framework with datamining techniques such as Principal Component Analysis, or Diffusion Maps [16] that can suggest the right coarse-grained observables.
6. Appendix: Center manifolds for the models of ODEs. Let us assume a continuous model in the form (f is considered to be sufficiently smooth) To determine a local center manifold of a non-hyperbolic fixed point (x * , p * ), the following linear transformation is introduced: where V is the matrix with columns the generalized eigenvectors v j of the Jacobian ∇ x f (x, p) computed at (x * , p * ). Expanding (28) in Taylor series around (x * , p * ) and introducing (29) we get: r(V z, p) is also expanded as a polynomial of x up to order q >= 2. The Jacobian J ≡ V −1 ∇ x f (x, p)V can be written in a block form as J = J c 0 0 J s , where J c is the l × l matrix whose eigenvalues are the l eigenvalues with zero real parts. Then a local center manifold of the form z s = h(z c ) can be calculated as a solution of the following equation [19]: For the model given by eq. (28), h(z c ) is approximated with a polynomial, which in turn is substituted in (31); then the coefficients of the terms of the same order in eq. (31) are equated. This leads to a system of algebraic equations to be solved for the unknown polynomial coefficients.