## 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 $\overline{\mathbf{u}}$, and of fast fluctuations, say **u**′. For instance, consider the Navier–Stokes equation for an incompressible and homogeneous fluid

**u**is the velocity field,

*P*is the pressure,

*ρ*is the (constant) density,

*ν*is the kinematic viscosity and $\mathrm{\u25b3}$ 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 $\mathbf{u}=\overline{\mathbf{u}}+{\mathbf{u}}^{\prime}$, where $\overline{{\mathbf{u}}^{\prime}}=0$ and the bar refers to a (running) time average over many wave periods. With this expression, (1.1) becomes

**u**′ do not include fluctuation–fluctuation nonlinearities; i.e. if the term (

**u**′ ·

**∇**)

**u**′ (and, hence, $\overline{({\mathbf{u}}^{\prime}\cdot \mathbf{\nabla}){\mathbf{u}}^{\prime}}$) 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**′ ·

**∇**)

**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 $\overline{\mathbf{u}}$ and

**u**′ is preserved: the cumulative effect of the fluctuation–fluctuation nonlinearities modifies the slow fields through the ‘Reynolds stress’ divergence $\overline{({\mathbf{u}}^{\prime}\cdot \mathbf{\nabla}){\mathbf{u}}^{\prime}}$ in (1.2), which in turn modifies the fluctuations through terms of the form $(\overline{\mathbf{u}}\cdot \mathbf{\nabla}){\mathbf{u}}^{\prime}$ and $({\mathbf{u}}^{\prime}\cdot \mathbf{\nabla})\overline{\mathbf{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 [5–7].

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 ${\mathbf{u}}^{\prime}\ll \overline{\mathbf{u}}$, hence $({\mathbf{u}}^{\prime}\cdot \mathbf{\nabla}){\mathbf{u}}^{\prime}\ll (\overline{\mathbf{u}}\cdot \mathbf{\nabla}){\mathbf{u}}^{\prime}$. 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):

*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

*E*

_{tot}=

*E*

_{c}+

*E*

_{w}, where

*E*

_{c}and

*E*

_{w}are, respectively, the ‘energy’ of the celerity field and of the waves, defined by

*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

*ϕ*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. [13–15]). Using the chain rule, ${\mathrm{\partial}}_{t}\to (\omega (T)/\u03f5){\mathrm{\partial}}_{\varphi}+{\mathrm{\partial}}_{T}$, the governing equations become

*c*,

*η*and

*ω*are expanded as asymptotic power series in

*ϵ*,

*ϵ*is small but variable, the resulting equations then can be solved order by order. At

*O*(1/

*ϵ*

^{2}), (2.5) yields ∂

^{2}

_{ϕ}

*c*

_{0}= 0, and since linear growth of

*c*

_{0}with

*ϕ*would result in unbounded growth of

*c*

_{0}, we conclude that

*c*

_{0}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 ∂

_{ϕ}

*c*

_{0}= 0, (2.5) at order

*O*(1/

*ϵ*) requires ∂

_{ϕ}

*c*

_{1}= 0; i.e.

*c*

_{1}also is a function only of

*x*and

*T*. Finally, at

*O*(1), we obtain

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

_{ϕ}

*c*

_{0}= 0, yields a necessary condition for the sublinear growth of ∂

_{ϕ}

*c*

_{2}:

*c*

_{2}(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

*c*

_{2}on the slow time scale.

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

*η*

^{2}

_{0}) 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

*c*

_{0}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 ${\hat{\eta}}_{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 d

*A*/d

*T*, 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), ${\hat{\eta}}_{0}$ is found to satisfy

*f*(

*x*) and

*g*(

*x*), we define the

*L*

^{2}inner product

*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 $\mathcal{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. [13–15]).

#### (d) Conservation of wave action

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

*δ*/

*δT*is a short-hand for (∂

_{T}

*c*

_{0})

*δ*/

*δc*

_{0}, where

*δ*/

*δc*

_{0}denotes

*functional differentiation*with respect to

*c*

_{0}(

*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

*c*

_{0}drives slow,

*O*(1) changes in the

*leading-order*wave eigenfunction ${\hat{\eta}}_{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 ${\hat{\eta}}_{0}$ and using the fact that $\mathcal{L}$ is self-adjoint, we obtain a solvability condition:

*ω*

_{1}, together with a relation for the evolution of the leading-order frequency d

*ω*

_{0}/d

*T*(see §2e), is required only to close the system for the

*correction*fields

*ω*

_{1},

*c*

_{1},

*η*

_{1}arising at higher order, should they be desired. By contrast, the solvability condition $({F}_{\mathrm{s}},{\hat{\eta}}_{0})=0$ yields a slow evolution equation for

*A*(

*T*). Specifically, we obtain

*amplitude equation*

*c*), termed an

*adiabatic invariant*. Given that the leading-order energy of the waves

*E*

_{0}= (

*ω*

_{0}

*A*)

^{2}/2, as computed from (2.3), the adiabatic invariant derived in (2.26) is the so-called

*wave action*

*E*

_{0}/

*ω*

_{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*:

*G*is the right-hand side of (2.27). Thus, a solvability condition is obtained by taking the inner product of (2.27) with ${\hat{\eta}}_{0}$, yielding

#### (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 *E*_{c} and *E*_{w} 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. d*t*≃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

*ω*

_{0}(

*T*) and ${\hat{\eta}}_{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

*E*

_{0}(

*T*) =

*E*

_{0}(0)

*ω*

_{0}(

*T*)/

*ω*

_{0}(0) owing to the conservation of wave action. For consistency with the initial conditions imposed in the DNS, we select

*E*

_{0}(0) =

*ω*

_{0}(0) = 1. Since the waves are not explicitly time-resolved, the timestep can be increased by an order of magnitude (d

*t*= 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.

### 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

*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

*ϵ*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

*σ*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

*ϵ*after substituting the following asymptotic expansions for

*σ*,

*Θ*and

*η*:

*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

*O*(1/

*ϵ*)) is

*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

*A*(

*T*) is the slowly varying modal amplitude, and ${\hat{\eta}}_{0}$ describes the spatial structure of the fluctuation mode. Note that

*A*, ${\hat{\eta}}_{0}$ and

*ϕ*are real-valued and that a suitable normalization condition on ${\hat{\eta}}_{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}+ ∂

^{2}

_{x}). 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}:

*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

*ϕ*/d

*t*= 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 $\overline{{\eta}_{0}^{2}}=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

#### (c) Solvability condition

As in §2c, the eigenvalue problem (3.15) can be cast into the form $\mathcal{L}{\hat{\eta}}_{0}=0$, now with $\mathcal{L}=-{\sigma}_{0}+{\mathit{\Theta}}_{0}+{\mathrm{\partial}}_{x}^{2}$. As a result of periodicity, the linear operator $\mathcal{L}$ is self-adjoint with the inner product defined in (2.16); i.e.

*A*and ${\hat{\eta}}_{0}$ by imposing $({\hat{\eta}}_{0},{\hat{\eta}}_{0})=1$; that is,

*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

*δ*/

*δT*= (∂

_{T}

*Θ*

_{0})

*δ*/

*δΘ*

_{0}. Forming the inner product of (3.19) with ${\hat{\eta}}_{0}$, we find that $(\mathcal{L}{\eta}_{1},{\hat{\eta}}_{0})=({\eta}_{1},\mathcal{L}{\hat{\eta}}_{0})=0$ and, hence, that $(F,{\hat{\eta}}_{0})=0$. Noting that the integral involving the functional derivative in $(F,{\hat{\eta}}_{0})$ vanishes as a result of the normalization condition, we obtain

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 *c*_{1} and an amplitude equation, (2.26), prescribing d*A*/d*T* 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

*viz.*,

*α*and

*β*are defined by

Equation (3.23) prescribes the slow evolution of *σ*_{0} as a function of the variables *Θ*_{0} and ${\hat{\eta}}_{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}/d*T*, 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}/d*T* = 0:

*η*

_{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

*t*= 0.01 for the reduced model and d

*t*= 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 ${\hat{\eta}}_{0}$. Finally, the amplitude of the fluctuations

*A*is set according to

*α*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 d*t*.) The time series of the normalized fluctuation energy (half the norm of *η*)

*η*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.

#### (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:

*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

*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=\sqrt{\alpha /\beta}$; 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:

*ϵ*, this system is of multi-scale QL form. The corresponding expression for

*β*, however, is modified

*β*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:

*ϵ*:

*σ*

_{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}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 [8–12], 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,24–27]. 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>.

### References

- 1.
Lindzen RS, Holton JR . 1968 A theory of the quasi-biennial oscillation.**J. Atmos. Sci.**, 1095–1107. (doi:10.1175/1520-0469(1968)025<1095:ATOTQB>2.0.CO;2) Crossref, Web of Science, Google Scholar**25** - 2.
Plumb RA . 1977 The interaction of two internal waves with the mean flow: implications for the theory of the quasi-biennial oscillation.**J. Atmos. Sci.**, 1847–1858. (doi:10.1175/1520-0469(1977)034<1847:TIOTIW>2.0.CO;2) Crossref, Web of Science, Google Scholar**34** - 3.
Plumb RA, McEwan AD . 1978 The instability of a forced standing wave in a viscous stratified fluid: a laboratory analogue of the quasi-biennial oscillation.**J. Atmos. Sci.**, 1827–1839. (doi:10.1175/1520-0469(1978)035<1827:TIOAFS>2.0.CO;2) Crossref, Web of Science, Google Scholar**35** - 4.
Dreeben TD, Chini GP . 2011 Two-dimensional streaming flows in high-intensity discharge lamps.**Phys. Fluids**,**23**056101 . (doi:10.1063/1.3584816) Crossref, Web of Science, Google Scholar - 5.
Chini GP, Malecha Z, Dreeben TD . 2014 Large-amplitude acoustic streaming.**J. Fluid Mech.**, 329–351. (doi:10.1017/jfm.2014.61) Crossref, Web of Science, Google Scholar**744** - 6.
Michel G, Chini GP . 2018 Strong wave–mean-flow coupling in baroclinic acoustic streaming.**J. Fluid Mech.**, 536–564. (doi:10.1017/jfm.2018.785) Crossref, Web of Science, Google Scholar**858** - 7.
Karlsen JT, Qui W, Augustsson P, Bruus H . 2018 Acoustic streaming and its suppression in inhomogeneous fluids.**Phys. Rev. Lett.**,**120**054501 . (doi:10.1103/PhysRevLett.120.054501) Crossref, PubMed, Web of Science, Google Scholar - 8.
Bouchet F, Nardini C, Tangarife T . 2013 Kinetic theory of jet dynamics in the stochastic barotropic and 2D Navier–Stokes equations.**J. Stat. Phys.**, 572–625. (doi:10.1007/s10955-013-0828-3) Crossref, Web of Science, Google Scholar**153** - 9.
Tobias S, Marston JB . 2013 Direct statistical simulation of out-of-equilibrium jets.**Phys. Rev. Lett.**,**110**104502 . (doi:10.1103/PhysRevLett.110.104502) Crossref, PubMed, Web of Science, Google Scholar - 10.
Rocha C . 2015 Coupled reduced equations for strongly stratified flows. Woods Hole Oceanographic Institution Technical Report, WHOI-2016-05. (doi:10.1575/1912/8601) Google Scholar - 11.
Nieves D, Grooms I, Julien K, Weiss JB 2016 Investigations of non-hydrostatic, stably stratified and rapidly rotating flows.**J. Fluid Mech.**, 430–458. (doi:10.1017/jfm.2016.443) Crossref, Web of Science, Google Scholar**801** - 12.
Squire J, Bhattacharjee A . 2015 Statistical simulation of the magnetorotational dynamo.**Phys. Rev. Lett.**,**114**085002 . (doi:10.1103/PhysRevLett.114.085002) Crossref, PubMed, Web of Science, Google Scholar - 13.
Hinch EJ . 1991**Perturbation methods**. Cambridge, UK: Cambridge University Press. Crossref, Google Scholar - 14.
Bender CM, Orszag SA . 1999**Advanced mathematical methods for scientists and engineers**. New York, NY: Springer. Crossref, Google Scholar - 15.
- 16.
Andrews DG, Mcintyre ME . 1978 On wave–action and its relatives.**J. Fluid Mech.**, 647–664. (doi:10.1017/S0022112078002785) Crossref, Web of Science, Google Scholar**89** - 17.
Herring JR . 1963 Investigation of problems in thermal convection.**J. Atmos. Sci.**, 325–338. (doi:10.1175/1520-0469(1963)020<0325:IOPITC>2.0.CO;2) Crossref, Web of Science, Google Scholar**20** - 18.
- 19.
Kruskal M . 1962 Asymptotic theory of Hamiltonian and other systems with all solutions nearly periodic.**J. Math Phys.**, 806–828. (doi:10.1063/1.1724285) Crossref, Web of Science, Google Scholar**3** - 20.
Kuzman GE . 1959 Asymptotic solutions of nonlinear second-order differential equations with variable coefficients.**J. Appl. Math. Mech.**, 730–744. (doi:10.1016/0021-8928(59)90164-9) Crossref, Google Scholar**23** - 21.
Bühler O . 2009**Waves and mean flows**. Cambridge, UK: Cambridge University Press. Crossref, Google Scholar - 22.
- 23.
Salmon R . 1998**Lectures on geophysical fluid dynamics**. Oxford, UK: Oxford University Press. Crossref, Google Scholar - 24.
Thomas VL, Lieu BK, Jovanović MR, Farrell BF, Ioannou PJ, Gayme DF . 2014 Self-sustaining turbulence in a restricted nonlinear model of plane Couette flow.**Phys. Fluids A**,**26**105112 . (doi:10.1063/1.4898159) Crossref, Google Scholar - 25.
Bretheim JU, Meneveau C, Gayme DF . 2015 Standard logarithmic mean velocity distribution in a band-limited restricted nonlinear model of turbulent flow in a half-channel.**Phys. Fluids**,**27**011702 . (doi:10.1063/1.4906987) Crossref, Web of Science, Google Scholar - 26.
Tobias SM, Marston JB . 2017 Three-dimensional rotating Couette flow via the generalised quasilinear approximation.**J. Fluid Mech.**, 412–428. (doi:10.1017/jfm.2016.727) Crossref, Web of Science, Google Scholar**810** - 27.
Pausch M, Yang Q, Hwang Y, Eckhardt B . 2019 Quasilinear approximation for exact coherent states in parallel shear flows.**Fluid Dyn. Res.**,**51**011402 . (doi:10.1088/1873-7005/aaadcc) Crossref, Web of Science, Google Scholar