Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences
You have accessResearch articles

Multiple scales analysis of slow–fast quasi-linear systems

G. Michel

G. Michel

Laboratoire de Physique Statistique, École Normale Supérieure, CNRS, Université P. et M. Curie, Université Paris Diderot, Paris 75005, France

[email protected]

Google Scholar

Find this author on PubMed

and
G. P. Chini

G. P. Chini

Department of Mechanical Engineering and Program in Integrated Applied Mathematics, University of New Hampshire, Durham, NH 03824, USA

Google Scholar

Find this author on PubMed

    Abstract

    This article illustrates the application of multiple scales analysis to two archetypal quasi-linear systems; i.e. to systems involving fast dynamical modes, called fluctuations, that are not directly influenced by fluctuation–fluctuation nonlinearities but nevertheless are strongly coupled to a slow variable whose evolution may be fully nonlinear. In the first case, fast waves drive a slow, spatially inhomogeneous evolution of their celerity field. Multiple scales analysis confirms that, although the energy E, the angular frequency ω and the modal structure of the waves evolve, the wave action E/ω is conserved in the absence of forcing and dissipation. In the second system, the fast modes undergo an instability that is saturated through a feedback on the slow variable. A new multi-scale analysis is developed to treat this case. The key technical point, confirmed by the analysis, is that the fluctuation energy and mode structure evolve slowly to ensure that the slow field remains in a state of near marginal stability. These two model systems appear to be generic, being representative of many if not all quasi-linear systems. In each case, numerical simulations of both the full and reduced dynamical systems are performed to highlight the accuracy and efficiency of the multiple scales approach. Python codes are provided as electronic supplementary material.

    1. Introduction

    Many physical, chemical, biological and ecological systems evolve on disparate time scales. By explicitly accounting for this scale separation, multiple scales analysis enables reduced equations governing the slow evolution of the comparably fast processes to be derived. This flexible mathematical formalism has proven fruitful in a variety of contexts, not only enabling significant computational savings but also leading to the introduction of adiabatic invariants (e.g. wave action) and new governing equations (e.g. the nonlinear Schrödinger equation). In this article, we illustrate the derivation of reduced equations for quasi-linear (QL) dynamical systems evolving on two distinct temporal scales. Slow–fast QL systems generically arise in the mathematical description of two broad classes of phenomena: (i) slowly modulated waves and (ii) instabilities that saturate through a feedback on the slow variable. A primary theme of the present work is that there are fundamental differences in the mathematical analysis of these multi-scale QL systems, with the latter requiring the introduction of a new asymptotic formalism.

    QL systems have a characteristic mathematical structure comprising two sets of equations that govern the coupled evolution of slow variables, say u¯, and of fast fluctuations, say u′. For instance, consider the Navier–Stokes equation for an incompressible and homogeneous fluid

    ut+(u)u=(Pρ)+νu,1.1
    where u is the velocity field, P is the pressure, ρ is the (constant) density, ν is the kinematic viscosity and is the Laplacian operator, and assume that the dynamics is characterized by fast waves coupled with a mean velocity field that evolves comparably slowly in time. To account for this scale separation, we introduce a Reynolds-like decomposition u=u¯+u, where u¯=0 and the bar refers to a (running) time average over many wave periods. With this expression, (1.1) becomes
    u¯t+(u¯)u¯+(u)u¯=(P¯ρ)+νu¯1.2
    and
    ut+(u¯)u+(u)u¯+(u)u(u)u¯=(Pρ)+νu.1.3
    The above set of equations is termed QL if the governing equations for the fluctuations u′ do not include fluctuation–fluctuation nonlinearities; i.e. if the term (u′ · )u′ (and, hence, (u)u¯) in (1.3) is negligible.

    Frequently, wave systems satisfy the QL approximation: denoting the wavenumber of the waves by k, the angular frequency by ω and the phase velocity (or ‘celerity’) by c, then

    |(u)utu|ukωuc,1.4
    which is a (generalized) Mach number. If this number is small compared with unity, then (1.3) reduces to QL form and the asymptotic methods introduced in this article can be applied. Since the nonlinear term (u′ · )u′ is neglected, quasi-linearity also guarantees that the fluctuations cannot interact directly to generate new harmonics, thereby limiting the range of temporal scales that must be resolved in simulations. Nevertheless, the dynamics remains rich and the two-way nonlinear coupling between u¯ and u′ is preserved: the cumulative effect of the fluctuation–fluctuation nonlinearities modifies the slow fields through the ‘Reynolds stress’ divergence (u)u¯ in (1.2), which in turn modifies the fluctuations through terms of the form (u¯)u and (u)u¯. These attractive attributes account for the increasing prevalence of QL models in the atmospheric, oceanic and astrophysical sciences (although in these models the averaging operation frequently is generalized to incorporate a spatial mean) and, more generally, in various branches of fluid mechanics, some examples of which are described below.

    One well-documented QL phenomenon is the quasi-biennial oscillation of the winds in the lower equatorial stratosphere, which undergo reversals approximately every 14 months. This slowly evolving flow results from a two-way coupling with internal gravity waves, which have a period of only a few tens of minutes [1]: the waves drive a shear flow through their Reynolds stress divergence, and the shear flow in return refracts the waves. Given the evident separation in temporal scales distinguishing the waves and the shear flow, along with the small Mach numbers characterizing the waves, these slow reversals can be captured in an idealized QL model [2,3]. As a second illustration, acoustic waves in a fluid exhibiting strong stable density inhomogeneities can nonlinearly interact to generate a mean Eulerian (streaming) flow that advects the density inhomogeneities, thereby modifying the wave frequency, amplitude and modal structure. This baroclinic acoustic streaming was first observed in high-intensity discharge lamps [4], with the wave period being 30 μs and the streaming flow evolving on a time scale of order 0.1 s, and subsequently shown to be described accurately by a QL dynamical system [57].

    The application of multi-scale QL models is not restricted to dynamically stable wave systems. Indeed, the QL reduction also has proved surprisingly effective for the investigation of spatio-temporally chaotic and even turbulent dynamical systems dominated by spectrally non-local energy transfers. In this latter scenario, e.g. for turbulent fluid flows, the fluctuation term u′ represents ‘eddies’ and quasi-linearity is realized when uu¯, hence (u)u(u¯)u. In the atmospheres of the Earth and gas giants, for example, turbulent small-scale eddies are known to drive slowly evolving zonal jets that in turn undergo shear instabilities and spawn small-scale eddies. Both the QL character and the time-scale separation manifest in this coupled eddy/jet system can be derived in the limit of small forcing and dissipation [8], thus providing a consistent framework to explain the spontaneous generation of these jets [9]. Similarly, in strongly (stably) stratified turbulence, anisotropic layers of horizontally moving fluid (oriented orthogonally to the direction of the imposed density gradient) spontaneously emerge that, owing to their relative motion, are susceptible to small-scale instabilities. A QL system has been derived in the asymptotic limit of strong stratification and shown to be capable of describing the dynamics of these anisotropic layers [10]. Other examples of turbulent yet approximately QL dynamical systems include strongly rotating (but weakly stratified) flows, such as open-ocean deep convection and high-latitude abyssal ocean currents (K. Julien 2018, private communication and see also [11]), and rotating astrophysical discs in which magnetic fields are generated [12].

    In this work, we describe the multiple scales analysis of two systems that are intended to be representative of the two broad classes of QL dynamics, involving waves and instabilities, respectively. These examples are sufficiently simple that the derivations are reasonably concise, allowing us to highlight the differences: specifically, for QL systems exhibiting fast exponentially growing instabilities rather than stable high-frequency wave motions, we show that the amplitudes of the ‘most dangerous’ fluctuation modes are slaved to the slowly evolving mean fields. The main novelty of our study is the new procedure we introduce to capture this slaving, which is not associated with the usual dissipative contraction to a slow manifold but rather with a marginal stability constraint. Furthermore, to illustrate the utility of the multiple scales approach in each case, direct numerical simulations of the master equations are compared to numerical simulations of the reduced models we derive. The simulations are coded in Python, using the open-source framework Dedalus, which allows for an easy and efficient implementation of initial-value, boundary-value and eigenvalue problems.1 Another attractive feature is that the user need only explicitly enter the equations, and the numerical scheme is then automatically generated. For completeness, all of the associated source files are commented and provided as electronic supplementary material.

    2. Modulated waves

    (a) Governing equations

    The first system we analyse is a one-dimensional model of fast, linear waves that are coupled strongly to their celerity field. For simplicity, only a bounded spatial domain is considered, although the analysis can be generalized to wave propagation in spatially extended domains. This system can be viewed as a generic model of wave/mean-flow interaction, crudely mimicking the phenomenology of, e.g. baroclinic acoustic streaming noted in the introduction. More generally, the analysis detailed in this section can be applied to any QL system in which the fluctuations are slowly modulated waves that are not susceptible to a rapidly amplifying instability. Here, the celerity field c and the wave field η are chosen to satisfy the following governing equations (and initial conditions):

    2ct2=(c1)+3cx2tc(ηx)2,c(x,0)=1,ct(x,0)=02.1
    and
    ϵ22ηt2=x(c(x,t)2ηx),η(x,0)=2Lcos(2πxL),ηt(x,0)=0.2.2
    The spatial domain x∈[0, L), and both c and η are L-periodic functions of x. Equation (2.2) is the one-dimensional wave equation with inhomogeneous celerity, where ϵ2 is the square of a small parameter. (Note that all variables are dimensionless.) Equation (2.1) describes the evolution of an initially uniform celerity distorted by the wave-induced Reynolds stress and subject to both restoring and dissipative processes. This set of equations conserves, in the absence of dissipation, the total energy Etot = Ec + Ew, where Ec and Ew are, respectively, the ‘energy’ of the celerity field and of the waves, defined by
    Ec=120L[(ct)2+(c1)2]dxandEw=120L[ϵ2(ηt)2+c2(ηx)2]dx.2.3
    For the given initial conditions, c and η are O(1), implying the right-hand sides of (2.1) and (2.2) also are of order unity. Owing to the prefactor ϵ2 on the left-hand side of (2.2), which can be interpreted as a measure of fluctuation inertia, we expect the wave field η to evolve on a much faster time scale than that characterizing the evolution of the celerity c.

    (b) Leading-order equations

    The dynamics thus evolves on both a fast, O(ϵ), and a slow, O(1), time scale. To exploit this temporal scale separation, we introduce a fast phase ϕ and a slow time T, where

    dϕdt=ω(T)ϵ,T=t.2.4
    Crucially, ϕ and T are treated as independent variables in the subsequent analysis. The introduction of the phase ϕ through its non-constant time derivative captures the slow evolution of the wave angular frequency ω and is referred to as the WKBJ approximation (e.g. [1315]). Using the chain rule, t(ω(T)/ϵ)ϕ+T, the governing equations become
    ω2ϵ22cϕ2+1ϵ[2ω2cϕT+dωdTcϕ]+2cT2=(c1)c(ηx)2+ωϵ3cx2ϕ+3cx2T2.5
    and
    ω22ηϕ2+ϵ[2ω2ηϕT+dωdTηϕ]+ϵ22ηT2=x(c2ηx).2.6
    To proceed, c, η and ω are expanded as asymptotic power series in ϵ,
    η(x,ϕ,T)=η0(x,ϕ,T)+ϵη1(x,ϕ,T)+O(ϵ2),2.7
    c(x,ϕ,T)=c0(x,ϕ,T)+ϵc1(x,ϕ,T)+O(ϵ2)2.8
    andω(T)=ω0(T)+ϵω1(T)+O(ϵ2),2.9
    and substituted into (2.5) and (2.6). Since ϵ is small but variable, the resulting equations then can be solved order by order. At O(1/ϵ2), (2.5) yields ∂2ϕc0 = 0, and since linear growth of c0 with ϕ would result in unbounded growth of c0, we conclude that c0 is a function of x and T only. This result confirms that at leading order c is a field that evolves strictly on the slow time scale. With ∂ϕc0 = 0, (2.5) at order O(1/ϵ) requires ∂ϕc1 = 0; i.e. c1 also is a function only of x and T. Finally, at O(1), we obtain
    ω022c2ϕ2+2c0T2=(c01)c0(η0x)2+3c0x2T.2.10

    We next introduce the fast average over ϕ (denoted with an overbar) of any function f(x, ϕ, T) bounded in ϕ:

    f¯(x,T)=limn12nπnπ+nπf(x,ϕ,T)dϕ,for integer n.2.11
    This averaging operation provides a suitable definition of the overbar used to define the Reynolds decomposition u=u¯+u in the Introduction. Applying this averaging procedure to (2.10), and recalling ∂ϕc0 = 0, yields a necessary condition for the sublinear growth of ∂ϕc2:
    2c0T2=(c01)c0(η0x)2¯+3c0x2T.2.12
    In fact, it is readily confirmed using the resulting equation for c2 (obtained by subtracting (2.12) from (2.10)) together with the expression for η0 given below in (2.14) that (2.12) ensures the boundedness of c2 on the slow time scale.

    We now turn to the leading-order reduction of the fluctuation equation (2.6):

    ω0(T)22η0ϕ2+x(c0(x,T)2η0x)=0.2.13
    Equations (2.12) and (2.13) comprise a two time-scale system, for which the fluctuation equation (2.13) does not involve nonlinear fluctuation–fluctuation terms (e.g. η20) and thus is QL. (Of course, for this system, no such nonlinearity was included in the master set of equations.) The key question that arises is how to accurately and efficiently integrate this two time-scale system? We proceed by observing that the solution of (2.13) is of the form
    η0(x,ϕ,T)=A(T)η^0(x,T)cosϕ,2.14
    since c0 does not depend on ϕ and noting that the solution component proportional to sinϕ vanishes as a result of the initial conditions (2.2) on η. Here, A(T) is a slowly evolving amplitude, and η^0 is a function that characterizes the spatial structure of the wave mode. A normalization condition must be prescribed to make this decomposition unique; see (2.17). Crucially, at this stage of the analysis, the slow evolution of A is unknown. Indeed, the ansatz (2.14) converts the initial-value problem (2.13) into a linear eigenvalue problem that places no constraint on the modal amplitude. Consequently, the reduced system is not closed on the slow time scale. The objective of the remainder of the analysis is to derive an equation for dA/dT, which for this system is found by examining the equation for the correction field η1.

    (c) Inner product and adjoint operator

    Multiple scales analysis, in its most general form, requires certain notions from linear algebra that we introduce here. Substituting (2.14) into (2.13), η^0 is found to satisfy

    ω02η^0+x(c02η^0x)=0,η^0(0,T)=η^0(L,T),η^0x(0,T)=η^0x(L,T).2.15
    This eigenvalue problem is linear, and can be expressed more compactly as Lη^0=0, where the linear operator L=ω02+x(c02x), with periodic boundary conditions, acts on the function η^0. For spatially periodic and real-valued functions f(x) and g(x), we define the L2 inner product
    (f,g)=0Lf(x)g(x)dx.2.16
    An inner product defines a norm, and we choose to set (η^0,η^0)=1; that is,
    0Lη^02(x,T)dx=1.2.17
    This relation renders the decomposition (2.14) unique. Moreover, after integrating by parts and making use of the periodic boundary conditions, we obtain
    (Lf,g)=0Lg(x)(ω02+x(c02x))f(x)dx=0Lf(x)(ω02+x(c02x))g(x)dx=(f,Lg).2.18
    This equality reveals that L is self-adjoint, a feature that, while not essential, simplifies the following analysis. For other boundary conditions (and/or other differential operators), the adjoint operator may not be equal to L but the following procedure can be appropriately modified. As shown in the next subsection, computing the adjoint linear operator enables a solvability condition to be imposed on the equation for the correction field η1, which, in turn, ensures that the asymptotic expansion (2.7) remains uniformly valid on the slow time scale T. This procedure generalizes the operation of ‘removing resonant or secular forcing terms’ usually introduced in early lectures on multiple scales analysis (e.g. [1315]).

    (d) Conservation of wave action

    The slow temporal evolution of A is obtained by collecting terms of order O(ϵ) in (2.6):

    ω022η1ϕ2+x(c02η1x)=2ω02η0Tϕ+dω0dTη0ϕ+2ω0ω12η0ϕ22x(c0c1η0x).2.19
    Note that the left-hand side of (2.19) and of (2.13) involve the same linear operator. By setting η1(x,ϕ,T)=η^1,c(x,T)cosϕ+η^1,s(x,T)sinϕ, (2.19) yields
    Lη^1,c=Fc,Fc=2ω0ω1Aη^02Ax(c0c1η^0x)2.20
    and
    Lη^1,s=Fs,Fs=2ω0(dAdTη^0+Aδη^0δT)Adω0dTη^0.2.21
    The notation δ/δT is a short-hand for (∂Tc0) δ/δc0, where δ/δc0 denotes functional differentiation with respect to c0(x, T), as arises here because of the tight, two-way coupling between the waves and the slowly evolving celerity field; specifically, the evolution of c0 drives slow, O(1) changes in the leading-order wave eigenfunction η^0, rendering the analysis non-standard. Nevertheless, we demonstrate below that the required functional derivatives can be evaluated explicitly, obviating the need for costly sensitivity computations.

    Taking the inner product of (2.20) with η^0 and using the fact that L is self-adjoint, we obtain a solvability condition:

    (Fc,η^0)=(Lη^1,c,η^0)=(η^1,c,Lη^0)=(η^1,c,0)=0;2.22
    i.e. the right-hand side of (2.20) must be orthogonal to the null eigenvector of the (adjoint) linear operator. Consequently,
    ω0ω10Lη^02dx=0Lη^0x(c0c1η^0x)dx=0Lc0c1(η^0x)2dx.2.23
    Using the normalization condition (2.17) then yields the following result:
    ω1=1ω00Lc0c1(η^0x)2dx.2.24
    This expression for the angular frequency correction ω1, together with a relation for the evolution of the leading-order frequency dω0/dT (see §2e), is required only to close the system for the correction fields ω1, c1, η1 arising at higher order, should they be desired. By contrast, the solvability condition (Fs,η^0)=0 yields a slow evolution equation for A(T). Specifically, we obtain
    (2ω0dAdT+Adω0dT)0Lη^02dx=Aω0ddT(0Lη^02dx)=0;2.25
    the last equality is a consequence of the normalization (2.17), following an interchange of the order of functional differentiation and integration. Thus, we obtain the amplitude equation
    2ω0dAdT+Adω0T=0ddT(A2ω0)=0.2.26
    The latter form of (2.26) reveals the existence of a conserved quantity in the limit of slow external processes (here, slow variation of the celerity c), termed an adiabatic invariant. Given that the leading-order energy of the waves E0 = (ω0A)2/2, as computed from (2.3), the adiabatic invariant derived in (2.26) is the so-called wave action E0/ω0. Conservation of wave action is a generic property of all slowly modulated, non-dissipative linear waves [16].

    (e) Slow evolution of the angular frequency

    With the formalism introduced in §2c, it is possible to derive an additional equation that, although not necessary for this problem, is in fact crucial for predicting the slow dynamics of the system analysed in §3. Accordingly, this additional condition is introduced here to facilitate comparison between the wave and instability problems and is derived by differentiating the leading-order eigenvalue equation (2.15) with respect to the slow time T:

    [ω02+x(c02x)]δη^0δT=2[ω0dω0dT+x(c0c0Tx)]η^0.2.27
    This system is of the form of L(δη^0/δT)=G, where G is the right-hand side of (2.27). Thus, a solvability condition is obtained by taking the inner product of (2.27) with η^0, yielding
    dω0dT=1ω00Lη^0x(c0c0Tη^0x)dx=1ω00L(η^0x)2c0c0Tdx.2.28
    This result provides an explicit evolution equation for the angular frequency of the waves.

    (f) Numerical implementation

    To assess the fidelity of the predicted slow dynamics to the actual system dynamics for small but finite values of ϵ, direct numerical simulations (DNS) of the governing equations (2.1) and (2.2) for various values of ϵ are compared to a numerical simulation of the reduced system (2.12), (2.15) and (2.26) obtained from the multiple scales analysis. All simulations are performed in a domain of size L = 2π with 32 grid points, use a second-order Runge–Kutta scheme, and output the time series of the energies Ec and Ew defined in (2.3). Although not documented here, the convergence of each simulation with increasing spatio-temporal resolution has been confirmed.

    The DNS are performed by solving system (2.1) and (2.2) using a Fourier pseudospectral method, the timestep being set in accordance with a CFL condition (e.g. dt≃0.01 for ϵ = 0.1). The reduced model is discretized on a Chebyshev grid, a requirement of the Dedalus eigenvalue solver; note that only (2.12) is explicitly time-advanced. In this equation, the term involving the fluctuations can be expressed as

    (η0x)2¯=E0(T)ω0(T)2(η^0x)2,2.29
    where ω0(T) and η^0 are obtained by numerically solving the eigenvalue problem (2.15) at each slow timestep rather than time-advancing (2.13) on the fast time scale, and the leading-order wave energy E0(T) = E0(0)ω0(T)/ω0(0) owing to the conservation of wave action. For consistency with the initial conditions imposed in the DNS, we select E0(0) = ω0(0) = 1. Since the waves are not explicitly time-resolved, the timestep can be increased by an order of magnitude (dt = 0.1) relative to that required for the DNS!

    The time series of the energies are reported in figure 1. Excellent agreement is achieved even for modest values of ϵ, and, for ϵ = 0.1, the simulations of the reduced equations provide a quantitatively accurate representation of the energy exchanges between the waves and the external medium. Of course, for this unforced and dissipative system, the total energy decays toward zero and a rest state is eventually reached. As noted in the Introduction, the associated source files for these simulations are commented and provided as electronic supplementary material. The files also can be used to generate space–time diagrams of the celerity field, should they be desired.

    Figure 1.

    Figure 1. Energy of the waves Ew (Ew(0) = 1) and of the celerity field Ec (Ec(0) = 0) obtained from direct numerical simulations (DNS) for ϵ = {0.1, 0.5} and from the multiple scales model (derived in the limit ϵ → 0) in which the waves are computed via the solution of an eigenvalue problem (EVP) instead of being explicitly time-resolved. (Online version in colour.)

    3. Unstable fluctuations

    (a) Governing equations

    The second example provides an illustration of a QL system in which the fluctuations would grow (or decay) exponentially fast in the absence of feedback on the slow variable. Nevertheless, we confirm that scale separation can be preserved as a result of the two-way coupling. In this example, the slowly evolving field Θ(x, t) and the fast fluctuation field η(x, t) satisfy

    Θt=PνΘη2,Θ(x,0)=13.1
    and
    ϵηt=Θη+2ηx2ϵη3,η(x,0)=cos(2πxL),3.2
    respectively, subject to periodic boundary conditions in a domain of spatial period L. It may be noted that the fluctuations satisfy a Ginzburg–Landau (GL) equation, e.g. describing the evolution of a non-uniform pattern amplitude η(x, t), although neither this specification nor interpretation is a requirement of the formalism. Moreover, unlike standard GL models, the ‘distance’ to the instability threshold is controlled by the (slow) field variable Θ, which evolves according to (3.1), the terms on the right-hand side respectively representing an external forcing P(x, t) > 0, a linear damping (ν > 0 is a damping coefficient), and a quadratic feedback from the fluctuations. The small parameter ϵ controls the temporal scale separation. The coupled system (3.1) and (3.2) crudely mimics, for instance, the QL reduction of the equations governing Rayleigh–Bénard convection first proposed by Herring [17].

    The coupling between the two fields Θ and η is such that there exists an energy for this system, E = EΘ + Eη, where

    EΘ=120LΘ2dxandEη=ϵ20Lη2dx.3.3
    It transpires that at leading order in ϵ the nonlinear term in (3.2) is negligible and, thus, the total energy E is conserved in the absence of forcing (P = 0) and dissipation (requiring both ν = 0 in (3.1) and zero diffusion in (3.2)). In the absence of the nonlinear term in (3.2), system (3.1) and (3.2) becomes QL.

    (b) Leading-order equations

    As for the slowly modulated wave system considered in §2, the small parameter ϵ motivates the introduction of two scales ϕ and T, treated hereafter as independent variables and defined according to

    dϕdt=σ(T)ϵ,T=t.3.4
    The notation σ rather than ω is used in this example to signify that the fluctuations may grow (or decay) exponentially rather than oscillate rapidly. Like ω, σ is defined to be an eigenvalue of a certain linear operator. A crucial distinction, however, is that while ω is the slowly varying angular frequency of any one of a countable infinity of wave modes (see §2), σ is the slowly evolving instantaneous growth rate of the most unstable (or least stable) fluctuation mode.

    Using (3.4), the governing equations (3.1) and (3.2) become

    σϵΘϕ+ΘT=PνΘη23.5
    and
    σϵηϕ+ηT=Θηϵ+1ϵ2ηx2η3.3.6
    This set of equations is solved order by order in ϵ after substituting the following asymptotic expansions for σ, Θ and η:
    σ(T)=σ0(T)+ϵσ1(T)+O(ϵ2),3.7
    Θ(x,ϕ,T)=Θ0(x,ϕ,T)+ϵΘ1(x,ϕ,T)+O(ϵ2)3.8
    andη(x,ϕ,T)=η0(x,ϕ,T)+ϵη1(x,ϕ,T)+O(ϵ2).3.9
    Equation (3.5) at O(1/ϵ) yields ∂ϕΘ0 = 0; i.e. Θ0 is a function strictly of x and T. At O(1), the fast average introduced in (2.11) removes the term proportional to ∂ϕΘ1, yielding
    Θ0T=PνΘ0η02¯.3.10
    The leading-order fluctuation equation (i.e. equation (3.6) at O(1/ϵ)) is
    σ0η0ϕ=Θ0η0+2η0x2.3.11
    Equations (3.10) and (3.11) comprise a two time-scale system that (as for the modulated waves example) can be closed on the fast time scale upon making the non-asymptotic replacement ∂T → (σ/ϵ)∂ϕ. This system is also QL: fluctuation–fluctuation nonlinearities are of higher order and have been consistently neglected in the fluctuation equation (3.11). Typically, the nonlinearity in (3.2), the equation from which (3.11) is derived, is crucial for the saturation of amplifying solutions that obey Ginzburg–Landau dynamics; its systematic omission in (3.11) implies that a distinct saturation mechanism is involved here.

    Once again, the question that arises is how to accurately and efficiently integrate this slow–fast QL system when two formally independent time scales are retained? The combination of temporal scale separation and quasi-linearity enables a modal solution of (3.11) to be expressed as

    η0(x,ϕ,T)=A(T)η^0(x,T)eϕ,3.12
    where A(T) is the slowly varying modal amplitude, and η^0 describes the spatial structure of the fluctuation mode. Note that A, η^0 and ϕ are real-valued and that a suitable normalization condition on η^0 is required to render the decomposition (3.12) unique; see (3.18). Upon substituting (3.12) into (3.11), σ0 is identified as the maximum eigenvalue of the operator (Θ0 + ∂2x). The leading-order fluctuation equation (3.11) seemingly places no constraint on the evolution of A(T), confirming that system (3.10) and (3.11) is not closed on the slow time scale.

    Of course, the potential exponential growth (or decay) of the leading-order fluctuation field η0 on the fast time scale would, if not properly addressed, break the posited scalings, as evidenced by the force acting on the slow field Θ0:

    η02¯=limnA2η^022nπnπ+nπe2ϕdϕ.3.13
    Note that, as for the fast temporal average introduced in (2.11), the fast time variable ranges from negative to positive infinity (rather than from, e.g. zero to positive infinity) while the slow time T is fixed. To ensure the convergence of this fast-time average for all x, one of the following two conditions must be satisfied: (i) A(T) = 0 or (ii) σ(T) = ϵdϕ/dt = 0 (implying the integrand in (3.13) is constant). These conditions thus require A(T) = 0 if the maximum instantaneous growth rate σ(T) is non-zero, a natural consequence of exponentially fast damping for σ < 0 and a requirement, were σ > 0, to preclude the blow-up of the fluctuations (and, hence, of the force (3.13)) on the fast time scale. With A≡0, the leading-order dynamics reduces to (3.10) with η02¯=0. Conversely, the fluctuations can have finite amplitude only if σ = 0; i.e. only if the fluctuations are marginally stable. Although σ initially need not equal zero, it continuously evolves under the forcing of the slow field and may eventually pass through zero. Heuristically, the persistence of one or more zero eigenvalues may be understood by recognizing that (i) Θ0 is not an arbitrary function of x owing to the feedback from the fluctuations and (ii) scale-separated QL systems with unstable fast fluctuations must self-tune to a state of near marginal stability. Indeed, as we demonstrate through this example, it is this regime that characterizes the long-term dynamics of the forced system.

    Thus, if the largest eigenvalue σ0 differs from zero, then the prescription is made that the modal amplitude A = 0. If, however, σ0 = 0, then the coupled system (3.10) and (3.11) is not closed on the slow time scale: a slow evolution equation for A(T) is required. In the next two subsections, we show how to derive such an equation. Although we enforce σ = 0 in the following analysis, we formally retain σ0 explicitly in the derivation for clarity and generality of exposition. Since ϕ is constant for σ = 0, the amplitude A can be re-defined so that the leading-order fluctuation field

    η0(x,ϕ,T)=A(T)η^0(x,T),3.14
    and, accordingly, the fast-time average η02¯=A2η^02. The leading-order system therefore becomes
    (σ0+Θ0)η^0+2η^0x2=03.15
    and
    Θ0T=PνΘ0A(T)2η^02.3.16

    (c) Solvability condition

    As in §2c, the eigenvalue problem (3.15) can be cast into the form Lη^0=0, now with L=σ0+Θ0+x2. As a result of periodicity, the linear operator L is self-adjoint with the inner product defined in (2.16); i.e.

    (Lf,g)=0Lg(x)(σ0+Θ0+x2)f(x)dx=0Lf(x)(σ0+Θ0+x2)g(x)dx=(f,Lg).3.17
    This inner product is used to disentangle A and η^0 by imposing (η^0,η^0)=1; that is,
    0Lη^02dx=1.3.18
    Guided by the analysis of the wave system, we first attempt to derive a slow-evolution equation for the amplitude A through the ‘usual procedure’ introduced in §2d; namely, by ensuring that the higher-order fluctuation equations are solvable. Specifically, collecting terms at O(1) in (3.6) leads to
    Lη1=F,F=η0TΘ1η0+η03=(dAdTAΘ1+A3η^02)η^0+Aδη^0δT,3.19
    where δ/δT = (∂TΘ0)δ/δΘ0. Forming the inner product of (3.19) with η^0, we find that (Lη1,η^0)=(η1,Lη^0)=0 and, hence, that (F,η^0)=0. Noting that the integral involving the functional derivative in (F,η^0) vanishes as a result of the normalization condition, we obtain
    1AdAdT+A20Lη^04dx=0LΘ1η^02dx.3.20

    In the wave problem, two solvability conditions emerge from the equation for η1: a closure, (2.24), required to evolve the dynamics of the correction fields ω1, η1, and c1 and an amplitude equation, (2.26), prescribing dA/dT as a function of the leading-order variable ω0. By contrast, in the present problem, only the single constraint (3.20), relating the time evolution of A to the higher-order field Θ1, is obtained. To procure a closed set of equations, it is natural to attempt to include Θ1 in the set of unknown fields. An equation for Θ1 can be derived from (3.5) at O(ϵ), yielding ∂TΘ1 as a function of η0 and, disappointingly, of η1, another unknown variable. Thus, similar to (2.24), (3.20) is merely a constraint required to obtain the dynamics of the correction fields η1 and Θ1.

    (d) Slow evolution of the growth rate

    In this example, a suitable constraint on the amplitude A can be derived by employing an operation analogous to that performed (solely for illustrative purposes) in §2e. Specifically, we differentiate the leading-order fluctuation equation (3.15) with respect to T, yielding

    Lδη^0δT=η^0(dσ0dTΘ0T).3.21
    A solvability condition is obtained by taking the inner product of (3.21) with η^0, again using the normalization condition on η^0; viz.,
    dσ0dT=0Lη^02(Θ0T)dx.3.22
    Upon substituting (3.16) into (3.22), this solvability condition reduces to
    dσ0dT=0L(PνΘ0A(T)2η^02)η^02dxαβA(T)2,3.23
    where the real coefficients α and β are defined by
    α=0L(PνΘ0)η^02dxandβ=0Lη^04dx.3.24

    Equation (3.23) prescribes the slow evolution of σ0 as a function of the variables Θ0 and η^0 and, crucially, of the fluctuation amplitude A. Since, for non-zero A, σ0 must be zero, the consistency of the multiple scales expansion excludes strictly positive values of dσ0/dT, which would immediately lead to positive growth rates; i.e. to exponentially growing fluctuations. This scenario arises only when α > 0 and σ0 = 0, in which case the amplitude A(T) must be set to enforce dσ0/dT = 0:

    A(T)=αβifσ0=0andα>0.3.25
    The above equation describes the effectively instantaneous (on the slow time scale) saturation of the instability via the feedback of the fluctuation field η0 on the slow variable Θ0.

    (e) Numerical implementation

    To demonstrate the efficacy of the proposed multi-scale algorithm, DNS of the governing equations (3.1) and (3.2) are compared to simulations of the multiple scales reduction given by (3.15), (3.16) and (3.25). All simulations are performed in a domain of size L = 2π that is discretized with 32 grid points. The damping coefficient ν = 1, and the external forcing

    P=1+12(cos(t)cos(x)+sin(0.6t)cos(2x)).3.26
    A second-order Runge–Kutta time-stepping scheme is used with fixed timestep dt = 0.01 for the reduced model and dt = 0.005 for the DNS. In the DNS, the governing equations (3.1) and (3.2) are solved using a Fourier pseudospectral method for ϵ = {0.1, 0.01}, while the reduced equations (valid as ϵ → 0) are discretized using Chebyshev polynomials on a non-uniformly spaced grid. For simulations of the reduced system, only the mean equation (3.16) is explicitly evolved in time (recall that T = t). At every slow timestep, however, the eigenvalue problem (3.15) is solved to obtain σ0 and η^0. Finally, the amplitude of the fluctuations A is set according to
    A(T)={0ifσ0<0orα<0,α/βotherwise,}3.27
    where the slowly varying coefficients α and β have been defined in (3.24).

    As evident in figure 2, this procedure prevents σ0 from becoming positive, thereby constraining the system to evolve along a marginal stability manifold. In practice, discretization and round-off errors preclude σ0 from identically equalling zero. Accordingly, for the results reported here, a tolerance of 0.01 has been enforced: 0≤σ0≤0.01. (Smaller tolerances can be achieved by decreasing the timestep dt.) The time series of the normalized fluctuation energy (half the norm of η)

    Eηϵ=0Lη22dx3.28
    is also plotted in this figure. (Although not plotted here, space–time diagrams of both η and Θ can be generated using the source files provided as electronic supplementary material.) Excellent agreement between the DNS and the reduced model is observed for ϵ≃0.01 after a transient that persists for a few slow time units. This ‘bursting’ regime, evident in the DNS, is very sensitive to the initialization of the fluctuations, because the amplitude of the fluctuations is extremely small when positive growth rates are first attained. (Recall that, before this instant, the fluctuations experience approximately exponential decay.) This transient is absent from the reduced model, since the energy of the fluctuations is instantaneously adjusted to a finite value once a state with zero growth rate is reached.
    Figure 2.

    Figure 2. Normalized energy of the fluctuations Eη/ϵ computed from direct numerical simulations (DNS) for ϵ = {0.1, 0.01} and from a novel multiple scales algorithm in which the fluctuations are not time-evolved but computed from the solution to an eigenvalue problem (EVP). The maximum instantaneous growth rate, σ0, obtained from the solution to the EVP is also reported. (Online version in colour.)

    (f) Additional comments

    The multiple scales analysis described in this section applies to QL systems in which fluctuations have the potential to be amplified via a linear instability mechanism. These systems may have other attributes that warrant further discussion, as detailed below.

    (i) Oscillatory instabilities

    The key requirement of the algorithm—tuning the amplitude of the fluctuations to prevent their exponential growth—results from the possibility that the real part of the growth rate σ could become positive, regardless of the value and sign of the imaginary part. For the system analysed in this section, the imaginary part of the growth rate is strictly zero; i.e. the instability is stationary. Nevertheless, a similar procedure can be applied when a bifurcation gives rise to an oscillatory instability, as, for instance, can occur in the following modified system:

    Θt=PνΘ+ν2Θx2ϵ2(ηt)2,Θ(x,0)=13.29
    and
    ϵ22ηt2=ϵΘηt+2ηx2,η(x,0)=cos(2πxL),ηt(x,0)=0.3.30
    In this example, (3.30) describes fast oscillatory motions (waves) of frequency O(1/ϵ) that are either damped or exponentially amplified by the slow field Θ. Generally, the growth rate σ is complex, but the amplitude of the fluctuations must be set, in the multiple scales analysis, by the requirement that the real part of σ should not increase once a state of marginal stability is reached.

    (ii) Stabilizing nature of the fluctuation-induced feedback

    Omitting the subscript ‘0’, the time-derivative of the (real) growth rate σ obtained in (3.23) is given by

    dσdT=αβA(T)2,3.31
    where, again, T is the slow time variable, A is the amplitude of the fluctuations, and α and β are two real parameters. The multiple scales approach developed in this section only applies if β > 0, regardless of the value of α. (If α > 0 and σ = 0, then A=α/β; otherwise A = 0.) That is, the fluctuations must provide a restoring force, not drive the slow field toward an even more unstable state. Although β > 0 in the model analysed here, see (3.24), this property is not generic; a slight modification of the set of equations (3.1) and (3.2) provides a counter-example:
    Θt=PνΘη2eΘ23.32
    and
    ϵηt=ΘηeΘ2+2ηx2ϵη3.3.33
    An analysis similar to that applied to (3.1) and (3.2) confirms that, for small values of ϵ, this system is of multi-scale QL form. The corresponding expression for β, however, is modified
    β=0Lη^04dxβ=0Lη^04(12Θ02)e2Θ02dx,3.34
    so that, crucially, β is no longer sign-definite. In the scenario σ = 0, α > 0 and β < 0, the fluctuations cannot prevent positive values of the growth rate from being realized, thereby invalidating the posited asymptotic scalings. In fact, the resulting exponential growth would be saturated by the nonlinear term in (3.33) by a physical process not captured by the QL system. Nevertheless, asymptotic methods still can be applied to this regime via a suitable rescaling of the unknowns:
    σ(T)=σ0(T)+ϵσ1(T)+O(ϵ2),3.35
    Θ(x,ϕ,T)=Θ0(x,ϕ,T)+ϵΘ1(x,ϕ,T)+O(ϵ2)3.36
    andη(x,ϕ,T)=η0(x,ϕ,T)ϵ+ϵη1(x,ϕ,T)+O(ϵϵ).3.37
    With these expansions, the following set of equations is obtained at leading order in ϵ:
    σ0Θ0ϕ=η02eΘ023.38
    and
    σ0η0ϕ=Θ0η0eΘ02+2η0x2η03,3.39
    where the leading-order growth rate σ0 is still defined by the linear eigenvalue problem obtained from (3.39) without the nonlinear term. Thus, in this regime, the system evolves strictly on the fast time scale, with the dynamics taking the form of punctuated bursting. The asymptotic analysis also reveals that, in this regime, the slow external forcing P and damping –νΘ0 do not contribute to the leading-order dynamics. Numerically, both (3.38) and (3.39) must be co-evolved on the fast time scale until σ0 reaches zero again, and the dynamics relaxes to a slow manifold describable by the asymptotic QL reduction.

    (iii) Degeneracy of the marginal eigenvalue

    Our last point concerns the assumption, implicitly made throughout this section, that the eigenvalue problem (3.15) has non-degenerate eigenvalues, at least for the eigenvalue having the largest growth rate (i.e. σ = 0 is a simple eigenvalue), in accord with the ansatz (3.12). In fact, the non-degeneracy of the leading eigenvalue can be demonstrated for this particular eigenvalue problem. Consider (3.15), and rewrite the variables as follows:

    η^0Ψ,σ0E,Θ0V.3.40
    Equation (3.15) then becomes
    EΨ=HΨ,H=2x2+V,3.41
    subject to periodic boundary conditions. Thus, the eigenvector of largest growth rate σ0 formally is equivalent to the ground state of the one-dimensional stationary Schrödinger equation, which is provably non-degenerate (e.g. item 21 of [18]). This property, however, clearly is very specific. More generally, it may be necessary to consider two (or more) independent marginal modes, each having its own amplitude. In this scenario, the corresponding number of equations describing the slow-time evolution of each associated growth rate can be derived, and the algorithm described in this section adapted accordingly.

    4. Conclusion

    Multiple scales analysis is usually introduced to capture the dynamics of a single variable or field evolving over disparate (spatio-)temporal scales. Canonical examples arise in the study of ordinary differential equations governing nonlinear oscillators (e.g. the van der Pol oscillator) and of partial differential equations governing nonlinear waves (e.g. the Korteweg–de-Vries equation). QL partial differential equations constitute an important class of dynamical systems that requires a generalization of this approach to treat two (or more) tightly coupled fields evolving on different time scales.

    Most (if not all) multi-scale QL systems can be categorized into one of the two model systems analysed in this article. For systems involving wave/mean-flow interactions, and more generally for any QL system in which the fluctuations are slowly modulated waves, an amplitude equation can be derived by imposing an appropriate solvability condition on the dynamics of the correction fields. As is well known [19,20], the amplitude equation reduces to the conservation of wave action if neither dissipation nor external forcing acts directly on the waves, in which case there exist complementary and, arguably, more elegant mathematical approaches—including the generalized Lagrangian mean formalism [16,21] and variational-based Hamiltonian methods [22,23]—to derive the appropriate conservation law(s) for the slow dynamics. Nevertheless, these alternative methods cannot readily incorporate forcing and dissipation, which are naturally and straightforwardly included in the methodology presented here. Moreover, the tight coupling between the wave and mean (e.g. celerity) fields renders our non-variational approach non-standard, as evidenced by the occurrence of functional derivatives in the coefficients arising in the amplitude equation that account for changes in the leading-order wave eigenfunction resulting from the slow evolution of the leading-order mean field; nevertheless, these functional derivatives can be explicitly evaluated, so that costly sensitivity analyses are not required.

    For many QL systems, particularly those arising in the modelling of spatially anisotropic turbulent shear flows (e.g. in wall-bounded engineering flows and in geo- and astrophysical turbulence [812], K. Julien, private communication), the fluctuations can exhibit various forms of dynamical instability. Generically, the QL reduction for such systems is only valid in the limit in which there is temporal scale separation and, in this limit (ϵ → 0 in our notation), the fluctuations can grow (or decay) exponentially on the fast time scale. Importantly, for these systems, any attempt to derive an amplitude equation for the slow evolution of the fluctuations by employing the ‘usual’ multiple scales approach (i.e. seeking a solvability condition by examining the equations for the correction fields) fails, instead merely providing closures needed for evaluation of the corrections fields, should they be desired. A primary contribution of the present article is the introduction of a new multiple scales formalism that is applicable to QL systems that exhibit dynamically unstable fluctuation fields. The essential point is that the slowly evolving amplitude of the fluctuations must be instantaneously prescribed—i.e. slaved to the slow mean field(s)—to prevent positive growth rates from being realized once a state of marginal stability is attained. In ongoing work, we are applying this new formalism to strongly stratified turbulent shear flows.

    We conclude by emphasizing that multiple scales analysis can yield a sizable improvement in computational efficiency, as the timestep employed by the numerical time-integrator can be increased from a fraction of a fast-time unit to a fraction of a slow-time unit. Clearly, this reduction in computational expense is crucial for the study of several realistic physical systems exhibiting strong scale separation. (In the quasi-biennial oscillation, for instance, the slow and fast time scales are separated by five orders of magnitude!) Even for systems in which the scale separation is less extreme, QL models have proven useful in many applications [9,12,2427]. We emphasize that when dynamical instabilities are possible these QL systems must self-tune toward a state of marginal stability, at least in a statistical (i.e. time-averaged) sense. Thus, the analysis and algorithm developed in §3 should prove valuable even for modest values of ϵ.

    Data accessibility

    All data included in this article can be generated from the Python codes provided as electronic supplementary material.

    Authors' contributions

    G.M. and G.P.C. conceived the mathematical models, interpreted the results, and wrote the paper. G.M. performed the numerical simulations. Both authors gave final approval for publication.

    Competing interests

    We declare we have no competing interests.

    Funding

    This research was supported in part by the National Science Foundation under grant no. NSF PHY-1748958.

    Acknowledgements

    The authors are grateful to S. Tobias, C. Doering and K. Julien for fruitful discussions, and acknowledge the hospitality of the Kavli Institute for Theoretical Physics at the University of California, Santa Barbara, where much of this research was completed.

    Footnotes

    1 See the entry in the Astrophysical Source Code Library, <http://ascl.net/1603.015>, and the Dedalus project homepage <http://dedalus-project.org>.

    Electronic supplementary material is available online at https://dx.doi.org/10.6084/m9.figshare.c.4413425.

    Published by the Royal Society. All rights reserved.

    References