Nonlinear emergent macroscale PDEs, with error bound, for nonlinear microscale systems

Many physical systems are formulated on domains which are relatively large in some directions but relatively thin in other directions. We expect such systems to have emergent structures that vary slowly over the large dimensions. Common mathematical approximations for determining the emergent dynamics often rely on self-consistency arguments or limits as the aspect ratio of the `large' and `thin' dimensions becomes nonphysically infinite. Here we extend to nonlinear dynamics a new approach [IMA J. Appl. Maths, DOI: 10.1093/imamat/hxx021] which analyses the dynamics at each cross-section of the domain via a rigorous multivariate Taylor series. Then centre manifold theory supports the global modelling of the system's emergent dynamics in the large but finite domain. Interactions between the cross-section coupling and both fast and slow dynamics determines quantitative error bounds for the nonlinear modelling. We illustrate the methodology by deriving the large-scale dynamics of a thin liquid film, where the film is subject to a Coriolis force induced by a rotating substrate. The approach developed here quantifies the accuracy of known approximations, extends such approximations to mixed order modelling, and may open previously intractable modelling issues to new tools and insights.


Introduction
Many systems of interest in science and engineering occur in a domain with disparate length scales (Davis 2017, e.g.): often a fine structure is modulated on a much larger scale (Mielke 1992, e.g.). Such disparate scales often are a major challenge in computational simulations (Rüde et al. 2016, p.14, e.g.). Two classic examples are Taylor-Couette flow (Iooss and Adelmeyer 1992, e.g.) and Benard convection (Segel 1969, e.g.). Often the fine-scale detail is crucial to the accurate modelling of a multiscale system, but with multiple length scales a simulation resolving the physical fine scale is not only prohibitively inefficient and severely constrained by memory limitations, it is also an arduous task to analyse simulation data generated at a scale much smaller than the scale of interest. This article develops a unified mathematical theory to use multiple length scales to reduce the full set of nonlinear governing equations to a simplified evolution equation, with quantified error, enabling more efficient simulations and analysis. This theoretical methodology should be able to better justify and illuminate many extant longwave and homogenisation theories (Bakhvalov and Panasenko 1989;Cross and Hohenberg 1993;Oron, Davis, and Bankoff 1997;Roberts and Li 2006;Pavliotis and Stuart 2008).
We focus on a class of multiscale systems whose physical domain is 'large' in multiple directions, but have a relatively 'thin' cross-section in the other directions. As an specific example, Section 1.1 takes as prescribed a variant of the integrated boundary layer pdes for a thin liquid film of Newtonian fluid spreading over a planar rotating surface (Chang 1994), and rigorously derives a simpler lubrication model of the nonlinear flow of the film (Wilson, Hunt, and Duffy 2000;Oron, Davis, and Bankoff 1997, §II.K, e.g.). Appendix A lists computer algebra code for deriving this lubrication model, with the code written to be readily adaptable to a wide range of systems. Thin liquid flows are important in biology, physics, and engineering, as well as in the environment. They may be of common liquids such as water or oil, or of more rheologically complex fluids, and display interesting nonlinear wave patterns (Oron, Davis, and Bankoff 1997;Wilson, Hunt, and Duffy 2000;Murisic et al. 2011;Lam et al. 2015). Other examples of systems amenable to our methodology include flood and tsunami modelling (Noakes, King, and Riley 2006;Bedient and Huber 1988;LeVeque, George, and Berger 2011), pattern formation (Newell and Whitehead 1969;Cross and Hohenberg 1993;Westra, Binks, and Water 2003), wave interactions (Nayfeh and Hassan 1971;Griffiths, Grimshaw, and Khusnutdinova 2006), elastic shells (Naghdi 1973;Mielke 1986;Lall, Krysl, and Marsden 2003), and microstructured materials (Romanazzi, Bruna, and Howey 2016).
Section 2 defines the generic nonlinear pde system to which the methodology applies, and defines the 'large' but 'thin' multiscale domain on which the system evolves. In pattern evolution problems a 'thin' domain variable is the phase of the underlying small scale pattern (Roberts 2015, §3.3, e.g.). A first step in the reduction of such pdes was taken by Mielke (1992), but the analysis required solutions to exist for all R-time and for all R m -space which precludes initial/ boundary value problems. Here we analyse the dynamics of a general pde in a thin cross-section of the large domain by constructing a multivariate Taylor expansion for the local spatial structures and analysing the evolution of the coefficients. Being the union of local-space-time modelling means the approach is valid everywhere outside of boundary layers (Roberts 1992, e.g.) and initial transients (Roberts 1989, e.g.). Section 3 details how to capture the emergent behaviour of the system at every chosen cross-section via constructing a set of generalised eigenvectors which span the centre subspace on which the slow system evolves. Based upon these eigenvectors the system's centre manifold is parametrised and an emergent nonlinear pde derived for the slow evolution. The order of the multivariate Taylor expansion determines the order of accuracy of the derived slow evolution, and Lagrange's Remainder Theorem provides a novel exact expression for the error of this slow evolution.
This article extends the methodology initially developed by Roberts (2015) to derive the long, slow space evolution of nonlinear pde systems on a long onedimensional physical domain with a relatively thin cross-section. Roberts and Bunder (2017) then developed the methodology to linear systems that have a domain that is large in multiple dimensions. Here we further develop the methodology to nonlinear pdes with multiple large dimensions. Our methodology is different to many other methods which derive the large-scale evolution in that here no asymptotic limit is required for the scale separation between the large domain and the thin cross-section (Kondic 2003;Náraigh and Thiffeault 2010, e.g.). Specifically, we have no requirement that a scale separation parameter (often named ) must be asymptotically small; we only require that such a small-large scale separation exists so that we can establish a centre-stable subspace (Assumption 3), and then our approach is valid at finite scale separation.

Example of a rotating shallow fluid flow
As an example application of some of the results, consider the flow of a shallow layer of fluid on a solid flat rotating substrate, such as in spin coating (Wilson, Hunt, and Duffy 2000;Oron, Davis, and Bankoff 1997, §II.K, e.g.) or large-scale shallow water waves (Dellar and Salmon 2005;Hereman 2009, e.g.). Let x = (x 1 , x 2 ) parametrise location on the rotating substrate, and let the fluid layer have thickness h( x, t) and move with depth-averaged horizontal velocity v( x, t) = (v 1 , v 2 ). We take as given (with its simplified physics) that the governing set of pdes is the nonlinear system ∂h ∂t where b represents viscous bed drag, f is the Coriolis coefficient, g is the acceleration due to gravity, and we neglect surface tension. The pdes (1) are similar to those used by Dellar and Salmon (2005, eq. (21)), but with only one component of the Coriolis force and the addition of viscous drag, and also similar to that used by Hereman (2009, eqs. (22)-(24)), but with a flat substrate. For such a shallow fluid flow, the horizontal gradient ∇ of quantities are mostly relatively small (Davis 2017, e.g.). Then the flow driven by variations of film thickness, ∇h, is approximately balanced in (1b) by the rotation and the substrate drag, leading to the velocity field Substituting this balance in the conservation of mass equation (1a) derives the single component, 'lubrication', model (2) Having just one component, we can use the pde model (2) as a simpler description of the shallow fluid dynamics. But pde (2) is an approximation to the 'original' pde (1), and so three outstanding questions are: can we find a rigorous error? can the analysis be extended to higher order? and can such an approach apply generally? Our answer to all three questions is yes.
Returning to the original system of pdes (1) and defining the system field u( x, t) = (h, v 1 , v 2 ), we rewrite system (1) as one equation, while also combining terms of like order: where matrices The system's first line (3a) is the linear part, whereas the second line (3b) contains the nonlinear terms which encode inertial acceleration. When the nonlinear terms (3b) are negligible, such as for very viscous fluids with low Reynolds numbers, the method described by Roberts and Bunder (2017) derives a slow linear pde approximation of (3), but the analysis and results of that article do not account for nonlinear dynamics. Herein we develop the methodology and theory to account for nonlinear effects. For the shallow fluid flow described by (3), the 'large' domain is some 'physical' subset of the x 1 x 2 -plane, and the 'thin' cross-section is the three-components of u = (h, v 1 , v 2 ) . The aim is to capture the long, slow behaviour of the original u( x, t) field in a one-component slow field U( x, t) (instead of the three components which describe the thin cross-section), and to construct a pde for U( x, t) which is correct to some order N in spatial derivatives, with known error. In general, higher orders N are potentially able to capture more extreme spatial fluctuations but may not be structurally stable, and so we address up to some low-moderate order N. This restricts U( x, t) to describing long, relatively gradual, spatial variations of the original field u( x, t). Section 3.2 defines the slow field U( x, t) in the general case and proves that it describes the behaviour of the original microscale field on the slowly evolving centre manifold (or slow manifold, as is the case in this shallow fluid example).
To determine the nature of the slow field U we first seek to understand the 'slow' and 'fast' evolution of the lowest-order linearisation of the system field u (Assumption 3 elaborates the general case). The linear dynamics of u are dominantly characterised by the lowest order linear term in (3), the term L 0 u where here The eigenvalues of L 0 indicate that u evolves on a one-dimensional slow subspace (one zero eigenvalue) and a two-dimensional stable subspace (two eigenvalues, −b± if, with negative real part). The stable part of u decays relatively quickly, roughly like e −bt , while the slow component of u, namely h, evolves on the one-dimensional slow subspace. In the notation introduced by Assumption 3 and Section 3, there exists a slow subspace of dimension m = 1 with right and left eigenvectors V 0 = (1, 0, 0) and Z 0 = (1, 0, 0), and eigenvalue A 0 = 0 .
Once the lowest-order linear dynamics of the system field u are known from matrix L 0 , we construct the generalised eigenvectorsṼ n , for | n| N , which span the spatially-local slow subspace of the linear system (3a) to the specified order of accuracy N. This order of accuracy N is that of a local multivariate Taylor expansion of the field u( x, t). The advantages of such a Taylor expansion are that not only does it provide a straightforward way to increase the order of accuracy N, it also provides a rigorous error term from Lagrange's Remainder Theorem. Section 3.1 discusses the general construction of the eigenvectorsṼ n , which in essence detail local out-of-equilibrium structures, and then Section 3.2 models the full nonlinear system to derive the slow manifold pde.
Appendix A lists computer algebra Reduce 1 code which applies the general theory of Sections 2 and 3 to determine the slow pde for any order N of the Taylor expansion, thus constructing pdes which describe the slow U = h field evolution of the shallow fluid dynamics to various orders of accuracy. See Sections 2 and 3 for details that justify (2) for order N = 2, and that for order N = 3 the slow pde is the 2D advection-diffusion pde with multi-indices 1,2 ∈ N 2 0 , symmetry a 1 2 = a 2 1 when | 1 | = | 2 | , and nonzero constant coefficients 2 a (01)(01) = a (10)(10) = a (02)(00) = a (20)(00 The four c coefficients equal to bg/(b 2 + f 2 ) correspond to the N = 2 approximation (2), but the rigorous derivation from Appendix A on the slow subspace 1 Reduce [http://reduce-algebra.com/] is a free, fast, general purpose, computer algebra system.
2 Equation (3b) describes two nonlinear terms which, following the derivation in Section 3, result in nonlinear terms with constant coefficients a j 1 2 for j = 1, 2 . But in this example, since the nonlinear terms are both second order (P j = 2 for j = 1, 2), to obtain (4a) we set a 1 1 2 empowers a far richer description of the slow dynamics in the x 1 x 2 -plane. Further, Section 3.2 provides an exact expression (45) for the error of slow pdes such as (4a); Appendix A.3 details how components of this error are constructed for this shallow fluid flow. Executing the computer algebra code in Appendix A to obtain slow pdes of higher orders is straightforward. Although the computational time increases rapidly with N, in principle the code is applicable to any order N. For example, the slow pde for order N = 4 is with the same nonzero constant coefficients a 1 2 given in (4b), as well as For small enough damping, b < f, the fourth order hyperdiffusion in (5a) makes the model structurally unstable. As is often necessary, for b < f one would then regularise the model as in the Benjamin, Bona, and Mahony (1972) regularised long wave equation. Notwithstanding such practical regularisation, our derived error expression (45) applies and is useful for as long as the spatial gradients in the solutions to (5a) remain small enough.
The following sections develop theoretical support for the derivation of nonlinear slow pdes such as (2), (4) and (5), and derive the novel exact algebraic expression (45) for the error.

Local expansion of general nonlinear dynamics
Consider some multiscale spatial domain X × Y where X is some open domain of large macroscale extent, and Y is a 'relatively small' microscale domain (in some Hilbert space). We analyse the dynamics of some field u within the multiscale spatial domain X × Y and determine the emergent behaviour of this field on the macroscale; that is, we aim to derive a description, over some time interval T, of the long, slow u field dynamics on the macroscale domain X while accounting for the fine details in the microscale domain Y in an 'averaged', 'homogenised' or 'slaved' sense (Roberts 2015;Roberts and Bunder 2017). As the domain Y is a small cross-section of the full domain of the system, a description of the large-scale behaviour should not involve fluctuations across Y as dynamic variables.
We consider the field u( x, y, t) in a given Hilbert space U (finite or infinite dimensional), where u : X × Y × T → U is a function of M-dimensional position x ∈ X ⊆ R M , cross-sectional position y ∈ Y, and time t ∈ T ⊆ R. We suppose the field u( x, y, t) satisfies some specified nonlinear pde in the form where f[u] : U → U is a 'strictly' nonlinear function of field u and its derivatives, the L k are linear operators (in y), the mulrivariate (mixed) derivative is of order | k| = k 1 + k 2 + · · · + k M , and where the infinite sum in the pde (6) is notionally written as being over all possible multi-indices k ∈ N M 0 (as usual, the set of natural numbers N 0 := {0, 1, 2, . . .}), although in practice there will be only a finite number of terms in this sum.
In application to fluid or heat convection the nonlinear term is the quadratic f[u] = u· ∇ u, whereas for the Ginzburg-Landau equation f[u] = u 3 . Consequently we consider nonlinearities that are sums of products of such factors.
Assumption 1. The nonlinear function f[u] may be written as, or usefully approximated as, a sum of products of u and its derivatives: for P j the order of each nonlinear term, and for some M-dimensional index p j i . (Sometimes we detail the case when f[u] has only one term in its sum.) Roberts (2015) considered systems such as (6) on the multiscale domain X × Y for one dimensional X ⊂ R , and by analysing the dynamics of the system at some cross-sectional station X ∈ X , constructed a reduced pde for the slowly varying dynamics. The construction relied on a Taylor expansion of the field u to order N, with the expansion made exact by including the Lagrange's Remainder term in the derivation. Analysis of the Taylor coefficients then reveals the slow behaviour of the system near X ∈ X within a centre manifold. Then a projection of the u field pde onto this centre manifold, generalised to the union over all stations X ∈ X , defined the slow pde within domain X, and a projection of Lagrange's Remainder determines the error of the pde. Roberts and Bunder (2017) analogously considered linear systems on the multiscale domain X × Y , but generalised to M-dimensional X ⊂ R M . Here we further generalise these earlier developments to the class of nonlinear pdes (6) that have multiple macroscale dimensions.

Large-scales modulates local Taylor coefficients
In this section we develop the procedure of Roberts and Bunder (2017, §3.1) to the additional complication of nonlinear effects f [u]. Both the field u and the nonlinearity f[u] are written as Taylor expansions with Lagrange Remainder terms. These remainder terms ensure the analysis of the dynamics of the system is exact.
The Taylor expansions and subsequent analysis require some assumptions about the smoothness of u. For k max denoting the largest magnitude derivative in the linear term of pde (6), for p j i representing all magnitude derivatives in the nonlinear term (7), and for Taylor expansion to order N, the field u must be in differentiability class C N+max(p j i ,k max ) . At every cross-section station X ∈ X ⊂ R M we expand the field u as an Nth order Taylor multinomial about x = X : with the multi-index factorial n! := n 1 !n 2 ! · · · n M ! , the multi-index magnitude | n| = n 1 + n 2 + · · · + n M , the multi-index power x n := x n 1 1 x n 2 2 · · · x n M M , and where • in the first sum, for | n| < N , the coefficients and u ( n) : X × Y × T → U ; • and in the second sum, for | n| = N , by Lagrange's Remainder Theorem for multivariate Taylor series (Roberts and Bunder 2017, eq. 18, e.g.), the coefficients and The Taylor expansion of the nonlinear term f[u] in pde (6) is expressed in the same way as for the field u; that is, we do not expand f[u] in a series in u, but instead expand f[u( x, y, t)] in a series in x to order N. As Assumption 1 specifies that f[u] is a sum of a product of linear functions, (7), the smoothness requirements for constructing the Taylor expansion of f[u] are satisfied by u being • and, for | n| = N , The Taylor coefficients f ( n) : U → U of the nonlinearity are, in principle, functions of the u field Taylor coefficients u ( k) with | k| N . They may be obtained by substituting the Taylor expansion (8a) of the field u into equations (9b) and (9c). For example, in the case of nonlinearity f[u] being only one term, a direct substitution of (8a) into the nonlinear form (7) gives where, for sums over multi-indices such as m= k we require that k i m i i for each component i = 1, 2, . . . , M .
As in the linear case (Roberts and Bunder 2017, Eq. (19)), the multivariate Taylor multinomial (8a) of a field u gives, after some rearrangement, that the th spatial derivative where appearing for the limits of some sums, ( k) ⊕ denotes the multi-index vector with ith component max(k i , 0), thus ensuring all multi-index components are nonnegative in the sums. Now, for every , | | N, take the th spatial derivative of the nonlinear pde (6), and substitute (11) for the spatial derivatives of field u (replacing with + k where appropriate), . (12) As the multivariate Taylor multinomial (8a) is exact, for all stations X ∈ X and x ∈ X , equation (12) is exact for every x ∈ χ( X), where χ( X) is an open subset of X such that for all points x ∈ χ( X) the convex combination X + s( x − X) ∈ χ( X) for every 0 s 1 ; this condition ensures that when we take the limit x → X , x will always remain inside χ( X) ⊂ X and ( x − X) → 0 . Now set x = X in equation (12) so that all terms containing factors of ( x − X) vanish. Unless otherwise specified, hereafter u ( n) denotes u ( n) ( X, y, t) when | n| < N and denotes u ( n) ( X, X, y, t) when | n| = N . Similarly, f ( n) denotes f ( n) ( X, y, t) when | n| < N and denotes f ( n) ( X, X, y, t) when | n| = N . In addition, swap the n and multi-indices in (12). Then, where the remainder where throughout means for each component, but excluding exact equality of the two multi-indices. The second term on the right-hand side of (12) (when | k| 1) determines the remainder (13b). Since multi-index n ∈ N M 0 and | n| N, the total number of unique multi-indices, and thus the number of coupled odes (13a), is For all indices | n| = 0, . . . , N, the u ( n) terms in equation (13a) are evaluated at station X, but the spatial derivatives of u ( n) ( X, x, y, t) with | n| = N that appear in the remainder term r n (13b) couple the dynamics at station X to dynamics of the system along the line joining fixed station X to variable position x, that is, the dynamics at X are coupled to the dynamics at points in the neighbourhood χ( X). This dependence of derivatives of u ( n) ( X, x, y, t) on the dynamics at points in χ( X) is directly seen from an application of the integral mean value theorem on equation (8c). By this theorem, there exists someŝ ∈ (0, 1) such that and X +ŝ( x − X) ∈ χ( X). Spatial derivatives of u ( n) ( X, x, y, t) retain dependence onŝ, and thus on the dynamics about X, even when evaluated at x = X. In contrast, u ( n) ( X, X, y, t) = ∂ n x u X is independent ofŝ. Whereas anŝ ∈ (0, 1) must exist for each u ( n) ( X, x, y, t), theseŝ are generally not determined, and so we view gradients of u ( n) ( X, x, y, t) as 'uncertain'. We therefore classify the remainder r n as uncertain forcing which couple the local dynamics at X to the dynamics in its neighbourhood, and thereby to the global dynamics over X.
The nonlinear f ( n) may also contain 'uncertain' gradients of u ( n) ( X, x, y, t), depending on the particular nonlinearity. For example, for the case of a singleterm nonlinearity we obtain the last line of equation (10) which contains spatial derivatives up to order p i of u ( n+ p i ) ( X, x, y, t) where | n + p i | = N . So, if at least one p i > 0 , the nonlinear term contains uncertain gradients which couple the dynamics at X to the dynamics in χ( X). Section 2.2 explicitly identifies these uncertain gradients in the nonlinear f ( n) .

Generating multinomial and PDE
We now pack all the multivariate Taylor coefficients u ( n) together into a generating function (multinomial). For every station X ∈ X and time t ∈ T consider the field u in terms of a local Taylor multinomial (8a) about the cross-section x = X . In terms of the indeterminate ξ ∈ R M , define the generating multinomial where this generating multinomialũ, through its range denoted by U N , is implicitly a function of the indeterminate ξ and the cross-sectional variable y. This generating multinomialũ : X × T → U N for the vector space and where ⊗ t represents the vector space tensor product. The generating operator acts to convert the original field u( x, y, t) into the generating multinomialũ( X, t) = Gu( x, y, t) . The generating operator (16) similarly converts the nonlinear term of pde (6) into a multinomial in ξ, withf [ũ] : U N → U N appearing as the nonlinear term in (18) of the next Proposition 2. We introduce the multinomialũ because it is more conveniently compact to deal with one multinomial and one pde than the N Taylor coefficients u ( n) and the N differential equations (13) derived in the previous section. Roberts and Bunder (2017, §3.2) constructed a similar multinomialũ and pde via the generating operator G; however, here we make new special provisions for the nonlinear term f [u]. Although the compact form of multinomialũ is useful, a more important property is the equivalence of the dynamics ofũ and the original field u( x, y, t), as described by Proposition 2.
Proposition 2. Let u( x, y, t) be governed by the specified nonlinear pde (6). Then the dynamics at every locale X ∈ X ⊂ R M is equivalently governed by the nonlinear pde for the generating function multinomialũ( X, y, t) defined in (15), the 'uncertain' forcingr [u] given by (22), nonlinearf[ũ] defined by (17), and operator To establish Proposition 2, we first show that multinomialũ( X, y, t) (15) satisfies pde (18). To construct a pde forũ, take the time derivative of (15) and replace ∂u ( n) /∂t using (13a): where in the first term the n and k sums are exchanged, and we then simplify this term using the useful identity obtained from derivatives of the generating multinomial (15) with respect to ξ. The above pde (20) is precisely pde (18) of Proposition 2 with forcing 'remainder' termr upon using expression (13b) for r n . The second task for establishing Proposition 2 is to show that the generating pde (18) and the original pde (6) describe the same dynamics at every locale X ∈ X ⊂ R M . We do this by providing a more physical interpretation of the generating operator G and the generating multinomialũ( X, t), beyond just a convenient way to pack the Taylor coefficients of u( x, y, t).
Consider the Taylor expansion of some general function g( x) ∈ C N+1 at x = X + ξ about x = X : where R N (g) is the order N Lagrange remainder term of g( X+ ξ) (Roberts and Bunder 2017). So, Gg( x) evaluates g( X + ξ) correct to O | ξ| N+1 . Similarly,ũ( X, t) = Gu( x, y, t) evaluates u( X + ξ, y, t) correct to O | ξ| N+1 . We interpretũ( X, t) as the projection of u( x, y, t) at x = X + ξ onto the space U N = U ⊗ t G N , with O | ξ| N+1 interpreted not as an error but as the difference between u( X + ξ, y, t) and its projection onto U N (Roberts and Bunder 2017). As G commutes with the temporal derivative G∂u( x, y, t)/∂t = ∂ũ( X, t)/∂t and ∂ũ( X, t)/∂t is equivalent to the Taylor expansion of ∂u( x, y, t)/∂t at x = X+ ξ correct to O | ξ| N+1 . Therefore the generating pde (18) for multinomialũ( X, t) is equivalent to the pde (6) for u( x, y, t) evaluated at x = X + ξ correct to O | ξ| N+1 . Thus the dynamics of pde (18) are identical to the dynamics of pde (6) at every x = X ∈ X . This completes the proof of Proposition 2. The dynamics of the original nonlinear pde (6) for field u are equivalent to the dynamics of the nonlinear pde (18) for the N dimensional multinomialũ (15); furthermore, the two pdes are symbolically the same with u ↔ũ and x ↔ ξ , plus a forcing term. But the advantage of the multinomial form is that the derivatives ∂ ξ operate only on G N , that is, multinomials of at most degree N in ξ ∈ R M , and are thus bounded in G N . In contrast, the derivatives ∂ x in the original pde are potentially unbounded (e.g., for u rapidly oscillating or containing irrational functions). The slowly varying modelling of Section 3.2 takes advantage of the near symbolic equivalence between pde (6) and pde (18) with u ↔ũ and x ↔ ξ.
We now expand the nonlinear term (17) of pde (18) explicitly in terms of generating multinomialũ and nonlinear 'uncertain' terms involving gradients of u ( n) with | n| > N . Section 3.2 makes use of this expansion to simplify the remainder term ρ of the slow pde, and the appendix applies the expansion in the construct the slow pde for the fluid flow example discussed in Section 1.1.
The nonlinear term (17) in the generating pde (18), expanded according to (7) in Assumption 1 is The components with | m i + p j i | > N in the last term on the right hand side are 'uncertain', similar to the uncertain forcingr[u] (22), although in the special case where all p j i = 0 , no such uncertain nonlinear terms exist. For the uncertain gradients, consider expansion (11) with | | > N evaluated at x = X , Using this expansion for the uncertain terms, as well as (8b) and (21) and Assumption 1, we rewrite the nonlinear term (24) as where in the second term f i [u,ũ] is either a function of the generating multinomialũ or of the uncertain gradients of original field u, Thus in the Taylor expansion (17), where the second term contains all uncertain gradients.

A slow nonlinear model emerges
Section 3.1 constructs the eigenspace which describes the emergent slow dynamics of the generating pde (18) by analysing a linearisation of the pde (18). We then show how this eigenspace and associated eigenvalues describe the slow dynamics of original nonlinear pde (6). It is possible to determine the slow dynamics of pde (6) from the eigenspace of the linearised Taylor coefficient pdes (13a), without introducing the generating multinomial and generating function, but then one must explicitly deal with N Taylor coefficients and their coupled N pdes, as seen in the linear example of Roberts and Bunder (2017) [ §2.2]. Employing the linear eigenspace as a foundation to describe the dynamics of a nonlinear system is justified by centre manifold theory (Carr 1981; Aulbach and Wanner 2000; Haragus and Iooss 2011, e.g.) which assures us that generically the stability properties of a linear system with centre-stable dynamics persist under nonlinear perturbations and time-dependent forcing. Section 3.2 extends the analysis to the nonlinear pde (18) to describe the slow dynamics of the nonlinear system, including coupling between the centre and stable subspaces via uncertain terms in both the forcing and the nonlinear term. The slow dynamics are characterised by a set of generalised eigenvectorsṼ n which span the centre subspace on which the slow dynamics ofũ evolve. As these generalised eigenvectors are determined from the linearised pde, they are the same as those we determined (Roberts and Bunder 2017, §3.3) for linear pdes. However, there we constructed the generalised eigenvectorsṼ n via a two-step process, firstly constructing generalised eigenvectors for the linearised version of the original pde (6) (i.e., for f[u]) = 0), and then mapping these eigenvectors into U N (Roberts and Bunder 2017, eq. (37) and Lemma 6). Here we show how to construct the generalised eigenvectors in U N directly from the generating pde (18).
To analyse the eigenspace and determine the slow dynamics of the linearised generating pde (18), we apply Assumption 3 which describes the eigenspace of L 0 , the lowest order operator inL. However, Assumption 3 does not provide necessary assumptions for the extraction of a slow model; for example, here we derive the slow model after assuming the Hilbert space U is a centre-stable subspace, but an analogous derivation is possible when U is a slow-stable subspace (the shallow fluid example of Section 1.1 is on a slow-stable subspace and the code in Appendices A.2 and A.3 permit either slow-stable or centre-stable dynamics).
Assumption 3. We assume the following for the primary case of purely centrestable dynamics.
1. The Hilbert space U is the direct sum of two closed L 0 -invariant subspaces, E 0 c and E 0 s , and the corresponding restrictions of L 0 generate strongly continuous semigroups (Gallay 1993;Aulbach and Wanner 1996).
For convenience, Definition 4 packs the m eigenvectors which span the centre subspace E 0 c of L 0 into one matrix V 0 , and similarly packs the eigenvalues into the matrix A 0 .
Definition 4. Assumption 3 identifies a subset of m eigenvectors of L 0 which span the centre subspace E 0 c ⊂ U . • With these eigenvectors define • Since the centre subspace is an invariant space of L 0 , define complex matrix A 0 ∈ C m×m to be such that L 0 V 0 = V 0 A 0 (often A 0 will be in Jordan form, but it is not necessarily so).
• Use ·, · to also denote the inner product on the Hilbert space U, ·, · : U × U → C , the field of complex numbers.
Interpret this inner product when acting on two matrices/vectors with elements in U as the matrix/vector of the corresponding elementwise inner products. For example, for Z 0 , V 0 ∈ U 1×m , Z 0 , V 0 ∈ C m×m .
• Define Z 0 ∈ U 1×m to have m linearly independent columns which are the m left eigenvectors of L 0 , ordered such that the jth columns of V 0 and Z 0 have the same eigenvalue and normalised such that Z 0 , V 0 = I m .
The next Section 3.1 uses the centre subspace eigenvectors V 0 of L 0 to generate a set of generalised eigenvectors ofL which describe the slow dynamics of linear pde ∂ tũ =Lũ confined to the centre subspace.

Generalised eigenvectors span the centre subspace
We invoke Assumption 3 to construct a set of eigenvectors (possibly generalised) which span the centre subspace E N c ⊂ U N of the linear operatorL (19). These eigenvectors capture the slow behaviour of the linear pde which is the linearisation of generating pde (18), with neglected forcing terms. For 0 < | n| N , we construct the generalised eigenvectorṼ n ∈ U 1×m ⊗ t G N = U 1×m N from the following recurrence relations, beginning withṼ 0 = V 0 , The m rows of allṼ n with | n| N form a subset of U N with mN elements. Below we show that these mN elements are generalised eigenvectors ofL which span the centre subspace E N c . To do this we show that the mN elements are linearly independent and that the generalised eigenvector equationLṼ n −Ṽ n A 0 only produces linear combinations ofṼ k with 0 k < n .
The inner product (28c) ensures thatṼ n = ξ n n! V 0 +Ṽ n for someṼ n ∈ U 1×m N such that Z 0 ,Ṽ n = 0 m for all | n| > 0 . Further, since the ξ n n! V 0 part ofṼ n gives zero in the left hand side of (28b), the objective of (28b) is to determine theṼ n part ofṼ n . As the right hand side of (28b) is a function ofṼ k with 0 k < n , we conclude thatṼ n has order of ξ no larger than the order of theseṼ k . Since we know that V 0 is independent of ξ, for | n| = 1 equation (28b) ensures thatṼ n is independent of ξ and the highest order ofṼ n when | n| = 1 must be ξ n . By induction we conclude that for any n,Ṽ n is of order k in ξ, where 0 k < n , andṼ n is of order n in ξ. ThusṼ n is an nth order multinomial in U 1×m N and for all | n| N we have N linearly independentṼ n . Now consider the rows of eachṼ n . Sincẽ with linearly independent eigenvectors v 0 j for j = 1, . . . , m , and since Z 0 ,Ṽ n = 0 m , each of the m elements ofṼ n are linearly independent. Therefore, the m elements of allṼ n with | n| N, form a set of mN linearly independent elements of U N .
To show that the rows ofṼ n are generalised eigenvectors ofL in the centre subspace E N c , consider The left-hand side only producesṼ k with 0 k < n, and thus the rows ofṼ n are generalised eigenvectors of rank n with eigenvalues in matrix A 0 . Since the rows of allṼ n with | n| N provide mN linearly independent generalised eigenvectors ofL with eigenvalues contained in A 0 , these mN eigenvectors must span the centre subspace E N c .
Lemma 5. For generalised eigenvectorṼ n constructed from recurrence relations (28), derivatives ofṼ n with respect to ξ satisfy ∂ m ξṼ n =Ṽ n− m for 0 < m n .
SinceṼ n = ξ n n! V 0 +Ṽ n withṼ n of order less than n in ξ, then ∂ n ξṼ n =Ṽ 0 , in agreement with Lemma 5 when m = n . However, to prove Lemma 5 we need only prove the | m| = 1 case for general n, as the additive property of derivative powers 3 then ensure the Lemma is true for all 0 < m n .
We prove Lemma 5 by induction. For | n| = 1 , sinceṼ n = ξ n n! V 0 +Ṽ n withṼ n independent of ξ, we know that ∂ n ξṼ n = V 0 =Ṽ 0 , and thus Lemma 5 is true for all | n| = 1 . Now assume that ∂ m ξṼ k =Ṽ k− m with | m| = 1 is true for all 0 < k n . Then, for | m| = 1 , replace n with n + m in (28b) and take the mth derivative with respect to ξ, where in the third line we recall that the highest ξ order ofṼ n+ m− k is n + m − k , so to take the mth derivative we must have k n ; in the fourth line we apply the assumption that ∂ m ξṼ k =Ṽ k− m for all 0 < k n; and the fifth line comes from equation (28b) . On comparing the first and last lines we see that ∂ m ξṼ V k− m with | m| = 1 is true for all 0 < k n and 1 | n| < N , then ∂ m ξṼ n+ m =Ṽ n is also true. Since ∂ n ξṼ n =Ṽ 0 is true when | n| = 1 , ∂ m ξṼ n+ m =Ṽ n must be true for all | n| > 0 when | m| = 1 . Finally, because derivative orders are additive, Lemma 5 must be true for all 0 < m n .

Slow field and PDE
In this section we complete our primary aim, which is to model the slow dynamics of the original field u( x, y, t). To do this, we project u( x, y, t) onto the centre subspace E 0 c and define this projection as the slow field U( x, t) = Z 0 , u( x, y, t) ∈ C m . The aim of this section is to construct a pde for U( x, t) with an exact error term. For the shallow fluid flow example of Section 1.1, pdes of different order are (4a) and (5a), but the error term is new.
The slow field U( x, t) evaluated at station x = X is equivalent to Z 0 , u( X + ξ, y, t) evaluated at ξ = 0 , and sinceũ( X, t) = u( X + ξ, y, t) + O | ξ| N+1 (Section 2.2 and equation (23)) the slow field is also equivalent to Z 0 ,ũ( X, t) evaluated at ξ = 0 . We expandũ( X, t) in terms of the centre modesṼ n and the analogous stable modes, and then project this parameterisation onto the the centre subspace E 0 c to obtain the slow field U( X, t). Since we projectũ( X, t) onto E 0 c , the stable modes may at first seem superfluous in the expansion ofũ( X, t). However, while the stable modes decay exponentially rapidly, resulting in the emergence of the evolution of u( x, y, t) on the centre subspace, through the nonlinearity these stable modes are not generally negligible and their influence must be accounted for in U( X, t).
We defineW n ∈ U 1×m ⊗ t G N = U 1×m N as the generalised eigenvectors which span the stable subspace E N s ⊂ U N ofL. The full set of generalised eigenvectors, V n andW n , span U N ofL, fully parameterising the fieldũ. Many of the properties of theW n are analogous to those of theṼ n and can be established by proofs similar to those presented in Section 3.1. Therefore, here we only briefly comment on those properties ofW n which are required for the analysis of generating multinomialũ, and ultimately the slow field U( x, t).
For convenience, define matrixṼ = [Ṽ n ] where the columns ofṼ are the centre subspace eigenvectorsṼ n . The ordering of theṼ n inṼ is according to the magnitude | n|, so that the first column isṼ 0 . From (29),LṼ n = 0 k nṼ n− k A k and so we define block upper triangular matrix A such thatLṼ =ṼA . The upper block triangular matrix A consists of N × N blocks with 0 m below the diagonal, A 0 ∈ C m×m along the main diagonal, and the ( k, n) block above the diagonal (that is, for n > k) is A n− k ∈ C m×m . The centre eigenvalues ofL are the eigenvalues of the blocks along the diagonal of A, namely, the m eigenvalues of A 0 , λ 1 , . . . , λ m , repeated N times. Similarly to matrix A which satisfiesLṼ =ṼA, we define matrix B such that LW =WB. Analogous to A, B is upper block triangular with N × N blocks of size C m ×m , with block B 0 along the main diagonal, and B n− k the ( k, n) block above the diagonal. Thus the stable eigenvalues ofL are the eigenvalues of B which must be λ m+1 , λ m+2 , . . . repeated N times. Recall that these stable eigenvalues all have real part −β < −Nα , whereas the magnitude of the real part of the centre eigenvalues are α (Assumption 3).
We capture the full centre-stable dynamics of the linear generating pde (27) on U N withũ for parameters U = [U n ] and S = [S n ] with U n ∈ C m and S n ∈ C m . As r[u],f[ũ] ∈ U N and the generalised eigenvectors,Ṽ n andW n , span U N the forcing and nonlinear terms are uniquely parameterised in terms of these generalised eigenvectors,r where r c = [r n c ] and f c = [f n c ] , and similarly for r s and f s , with r n c , f n c ∈ C m and r n s , f n s ∈ C m . We substitute the expansion ofũ (30) into the nonlinear generating pde (18) and separate the forcing and nonlinear terms into centre and slow components (31). From Section 3.1 we know that the centre subspace eigenvectorsṼ n are linearly independent, and similarly the stable subspace eigenvectorsW n are linearly independent. Therefore, we separate the centre and slow components of the pde to obtain A general solution of the pde for the stable parameter S (32b) is with convolution h(t) g(t) := t 0 h(t − τ)g(τ) dτ . As all eigenvalues of B have real part −β < −Nα, for some decay rate γ ∈ (α, β) , This solution for S shows that, after a sufficiently long time, the forcing and nonlinear terms dominate S through the convolution, thus showing how the forcing and nonlinear terms couple the centre and stable solutions through f s/c and r s/c and why the influence of the stable modes are not negligible. We now construct a pde for the slow field U( X, t) by considering both its the temporal and spatial derivative in terms of the centre-stable dynamics. Since U( X, t) = Z 0 ,ũ( X, t) ξ= 0 at station x = X , From (28c), the inner product Z 0 ,Ṽ ξ=0 = I m 0 m · · · 0 m . Also, because of the upper block triangular structure of A, the k element of AU is Therefore, in the first term on the right hand side of (35) we only retain the k = 0 element of AU, and in the second and third terms we only retain f 0 c and r 0 c , respectively. So now, Now consider the order n spatial derivative of the slow field, From Lemma 5, Z 0 , ∂ n ξṼ k ξ= 0 = Z 0 ,Ṽ k− n ξ= 0 which, from (28c), equals the identity I m if k = n , but zero otherwise. Thus, in the first inner product of (37), only the U n element of U remains. Then, in the second inner product of (37), substitute the solution of the stable parameter S (34). The spatial derivative is now Combining (36) and (38), Whereas this equation symbolically resembles a pde, it is strictly a differentialintegral equation which couples the dynamics at each station X via the 'uncertain' gradient terms and the stable parameter S(t), which is dependent on the history convolution integrals (33). To obtain a slow pde without this coupling to different stations, such as pde (4a) used in the shallow fluid flow example, we retain all terms which do not couple to different stations (i.e., no dependence on derivatives ∂ x u ( n) with | n| = N and no dependence on S(t)) and regulate all other terms to a remainder. The last line of (39) contains a convolution, so is part of the remainder, and the two forcing terms r 0 c (t) and Z 0 ,W ξ= 0 r s (t) are dependent on uncertain gradients, so are also in the remainder. In contrast, the nonlinear terms f 0 c (t) and f s (t) contain parts which we want to retain in the slow pde, as well as terms which should be in the remainder. For specific cases, removing the remainder components from the nonlinear terms in (39) is achieved using (26), as shown in Appendix A.2. Here, for the general case, we show that the nonlinear terms which are retained in the slow pde must take a particular form. First, separate the nonlinear terms in the second line of (39) into two parts, where f 0 c ( X, t) ∈ C m and f 0 s ( X, t) ∈ C m contain no uncertain terms and no dependence on S (so are retained in the slow pde), and where f 0 c,r (t) and f s,r (t) contain all uncertain terms and S dependent terms. The f 0 c,r (t) and f s,r (t), as well as the convolutions and forcing terms in (39), are not in retained in the slow pde. As the nonlinear function f [u] in the original pde (6) is a sum of nonlinear terms f j [u] (Assumption 1), the nonlinear f 0 c ( X, t) and f s ( X, t) are also a sum of nonlinear terms indexed by integer j, where f j 0 c ( X, t) ∈ C m , and f j 0 s ( X, t) ∈ C m . As f j [u] is of order P j in u and its derivatives (Assumption 1), f j 0 c and f j 0 s must be of order P j in U n ∈ C m for all | n| N 4 . So, in general, each k = 1, . . . , m element of f j 0 c + f j s must have the form for some constant vector a j k 1 2 ... P j ∈ C m P j and where ⊗ represents the usual Kronecker product, for which U p ⊗ U q ∈ C m 2 and U 1 ⊗ U 2 ⊗ · · · ⊗ U P j ∈ C m P j . 5 On replacing all U with spatial derivatives of ∂ x U, as shown in (38), the kth coordinate of the m-dimensional nonlinear term retained in the slow pde must have the form Now, on replacing arbitrary station X with x ∈ X , the slow pde determined from the differential-integral equation (39) is with a jT 1 2 ... P j (∂ 1 x U) ⊗ · · · ⊗ (∂ P j x U) the m-dimensional vector with elements k = 1, 2, . . . , m defined by (43), and with remainder The nonlinear f 0 c (t) and Z 0 ,W ξ= 0 f s (t) are sums of nonlinear terms of order P j in both U n and S n for all | n| N , but as all S n dependence is contained in f 0 c,r (t) and f s,r (t) there is only U n dependence in f 0 c ( X, t) and f s ( X, t). 5 For m = 1, as in the fluid flow example of Section 1.1, the constant vector in (42) reduces to a scalar and the Kronecker products reduce to a multiplication of scalars. For P j = 2 and any value of m, each term in (42) Analogous slow pdes were derived by Roberts (2015) (equation (22)) and Roberts and Bunder (2017) (equation (51)), but without the nonlinear terms.
Simplifications of the remainder ρ (45) are possible when the order N is chosen to be higher than the order of the spatial derivatives in the original pde (6). The original pde contains linear operators L k ∂ k x for k satisfying 0 | k| < ∞, but in practice there will be an upper limit on | k|, say k max (often k max = 1, 2-the example of Section 1.1 has k max = 1). Assume that N > k max and consider the uncertain linear terms in ρ: Since Z 0 is independent of ξ, we need only consider ξ = 0 inr[u] (22). When ξ = 0 the right hand side of (22) requires | n| = 0 , | | = N and k , but since N > k max | k| we can never satisfy | | = N and k. So, when N > k max we havẽ r[u] ξ=0 = 0 and r 0 c (t) + Z 0 ,W ξ= 0 r s (t) = 0 . Similarly, consider the projection of the nonlinear term and then expandf[ũ] using (25). If N is chosen to be larger than any spatial derivative in the nonlinear term, that is N > p j i for all i = 1, 2, . . . , P j and for all j the number of nonlinear terms, theñ contains no uncertain terms. So, when separating the nonlinear terms according to (31) and (40) f 0 c,r (t) and f s,r (t) contain all S dependence and any convolution terms, but no uncertain terms.
We have shown that N > max(p j i , k max ) removes the uncertain terms from f 0 c,r (t) and f s,r (t) and sets r 0 c (t) + Z 0 ,W ξ= 0 r s (t) = 0 , but this does not remove all uncertain terms from the remainder ρ (45). Uncertain terms are still present in the remainder because of f s (t) and r s (t) which appear in the convolution in the second line of (45). Section 1.1 presents the example of a shallow fluid flow on a rotating substrate and, with computer algebra code provided in Appendix A, constructs slow pde of the form given in (44), as shown in equations (4a) and (5a), for N = 3 and N = 4 , respectively. Appendix A.3 calculates parts of the remainder ρ, such as Z 0 ,W B and A n Z 0 , ∂ n ξW , and shows that, since the order is sufficiently large (N > max(p j i , k max ) = 1) we have r 0 c (t) + Z 0 ,W ξ= 0 r s (t) = 0 . Whereas the appendix is written to support the example presented in Section 1.1, only Appendix A.1 is specific to this example, with the code in Appendices A.2 and A.3 written in a general format so as to be readily adaptable to a large number of systems.

Conclusion
This article further develops a general theory to support practical approximations of slow variations in space. This methodology was initially developed by Roberts (2015) for one dimensional space, and then extended by Roberts and Bunder (2017) to linear systems in multi-dimensional space. We here provide theoretical support and a practical example for the case of a nonlinear system of pdes in a spatial domain that is large in multiple dimensions. The significant advantages of the theoretical methodology are: • the approach is readily applicable to a wide range of systems, as illustrated by the general theory provided in Sections 2 and 3; • higher order pde are obtained in a straightforward manner by increasing the order N of the Taylor expansion; and • the resulting slow pde has a well-defined error, with a derived an algebraic form which can be bounded in applications.
In the general theory, we make some assumptions about the structure of the nonlinear microscale system and its dynamics. Assumption 1 requires that the nonlinearity in the original microscale pde should be a sum of products of the unknown field and its derivatives, and Assumption 3 requires centre-stable dynamics.
The key requirement for the presented methodology is the persistence of the centre manifold of the linear system (described in Section 3.1) when perturbed by nonlinearities and time-dependent forcing, thus justifying the importance of the eigenspace of the linearised system to the full nonlinear microscale system (Section 3.2). As other invariant manifolds are often similarly persistent, we expect that the methodology is not restricted to the centre-stable dynamics required by Assumption 3. Indeed, we show that other invariant manifolds are possible with the fluid flow example in Section 1.1, which has slow-stable dynamics. Furthermore, this fluid flow example has a more complex nonlinear structure than that required by Assumption 1. The Reduce Algebra code presented in Appendix A is designed for this fluid flow example, but is written so as to be adaptable to other systems, including those with different nonlinear structures.
Future research will aim to further generalise the methodology. Of particular interest is stochastic dynamics (Arnold and Imkeller 1998;Roberts 2008), deriving boundary conditions for the slowly varying model from microscale boundary conditions (Segel 1969;Roberts 1992;Mielke 1992), and non-local operators (i.e., beyond the local operators L k ∂ k x in (6)) (Calcagni, Montobbio, and Nardelli 2008).
A Computer algebra determines the emergent macroscale model The Reduce Algebra code presented here takes the microscale fluid flow equations (1) for field u( x, y, t) and derives the emergent slow macroscale dynamics (4a) or (5a) for the slow field U( x, t) on x ∈ X , where domain X ⊂ R M with M = 2 is specified as 'large'. The microscale system (1) decomposes into a stable subspace and an m = 1 dimensional slow subspace, with the slow field evolving on this slow subspace U ∈ R m . The Reduce code is applicable to other systems with original microscale field u ∈ U ⊂ R d on an M = 2 dimensional large domain, and described by pdes with an order P = 2 nonlinear term, provided the dynamics of these systems decompose into a d−1 stable subspace and an m = 1 dimensional centre or slow subspace, thus producing a slow field U( x, t) ∈ R . But this code is readily adaptable to systems with different specifications. Other parameters in the code are easily changed; for example, the dimension d of the field u and the order of the Taylor expansion N are variables in the computer algebra.
We firstly set some printing options and ensure that Reduce will provide us with complex number solutions. We then specify the dimensions of the u field and the desired order of the Taylor expansion N.
1 on div; off allfac; on revpri; on complex; 2 d:=3; % dimension of u field 3 nn:=3; % order of Taylor expansion Appendix A.1 defines the microscale pde of field u( x, y, t), introduced in Section 1.1, which describe fluid flow on a rotating substrate. Appendix A.2 is generic code for any system with an m = 1 dimensional slow subspace, a 'large' domain of dimension M = 2 , and nonlinearity of order P = 2 . This generic code constructs both slow eigenvectors and eigenvalues,Ṽ n and A n , and stable eigenvectors and eigenvalues,W n and B n , of the matrix operatorL (19) for | n| N . Then the code produces the slow macroscale pde for the slow field U( x, t) ∈ R m with x defined on the large domain X ⊂ R M .

A.2 Generalised eigenvectors and slow PDE
The code below derives the centre (or slow) and stable modes for some pde of the form (6), and constructs the slow pde projected onto the centre (or slow) subspace (as describe in Section 3.2). The code is applicable to any d dimensional microscale field on an M = 2 dimensional 'large' domain with a pde of nonlinear order P = 2 , where the d dimensional field separates into a m = d − 1 stable subspace and a m = 1 dimensional slow or centre subspace (although v, defined in line 63, is list of at least d dummy variables, so the number of dummy variables in this list should be increased if field u has dimension d > 7 , and similarly for w defined in line 97, is list of at least d(d − 1) dummy variables).
This code is used to derive the slow pde of the fluid flow problem described in calculate the left and right eigenvectors associated with the eigenvalues in B 0 of matrix L 0 , and then apply a recurrence relation analogous to (28) to construct the stable eigenvectorsW n and associated matrices of eigenvalues B n ofL for all 0 < | n| N .

A.3 Remainder of the slow PDE
In this section we calculate some terms which appear in the remainder ρ (45) to illustrate how this remainder is constructed for a specific model. We assume that the remainder is evaluated after some relatively long time so that terms O e −γt are negligible; as γ ∈ (α, β), this ensures that the linear effects of the stable modes are negligible (Assumption 3). We calculate all inner products Z 0 ,W n ξ= 0 and Z 0 , ∂ n ξW k ξ= 0 and use these to calculate ipzw(j,k1,k2)*sub(solr,rs(j,k1,k2))); 296 zero:=rc0+rsall; 297 298 end; Here we do not calculate the nonlinear terms which appear in the remainder ρ (45), but these terms can be calculated by editing the nonlinear calculations in Appendix A.2 to include the second line of (24). As before, equation (31) separates the centre and stable components off [ũ], and (40) separates the nonlinear components into terms which are included in the slow pde (44) and those which are in the remainder ρ (45). As the nonlinear calculations in Appendix A.2 are already fairly memory intensive, we do not extend these calculations to deriving the nonlinear remainder terms in ρ.