Multiple scales analysis of slow--fast quasilinear systems

TL;DR

This paper applies multiple scales analysis to two quasilinear systems with fast and slow dynamics, demonstrating conservation laws and stability feedback, validated by numerical simulations and providing Python codes.

Contribution

It develops a new multiscale analysis method for quasilinear systems with coupled fast and slow modes, including instability saturation.

Findings

Wave action is conserved in the first system.

Fast mode instability is saturated through feedback on the slow variable.

Numerical simulations confirm the accuracy of the multiscale approach.

Abstract

This article illustrates the application of multiple scales analysis to two archetypal quasilinear 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 $\omega$, and the modal structure of the waves evolve, the wave action $E/\omega$ 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 multiscale 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 quasilinear 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 supplementary material.

Equations

$$\frac{\partial \mathbf{u}}{\partial t}+\left(\mathbf{u}\cdot \nabla \right)\mathbf{u}=-\nabla \left(\frac{P}{\rho }\right)+\nu △\mathbf{u},$$

$$\frac{\partial \bar{\mathbf{u}}}{\partial t}+\left(\bar{\mathbf{u}}\cdot \nabla \right)\bar{\mathbf{u}}+\overline{\left({\mathbf{u}}^{\prime }\cdot \nabla \right){\mathbf{u}}^{\prime }}=-\nabla \left(\frac{\overline{P}}{\rho }\right)+\nu △\bar{\mathbf{u}},$$

$$\frac{\partial {\mathbf{u}}^{\prime }}{\partial t}+\left(\bar{\mathbf{u}}\cdot \nabla \right){\mathbf{u}}^{\prime }+\left({\mathbf{u}}^{\prime }\cdot \nabla \right)\bar{\mathbf{u}}+\left({\mathbf{u}}^{\prime }\cdot \nabla \right){\mathbf{u}}^{\prime }-\overline{\left({\mathbf{u}}^{\prime }\cdot \nabla \right){\mathbf{u}}^{\prime }}=-\nabla \left(\frac{P^{\prime }}{\rho }\right)+\nu △{\mathbf{u}}^{\prime }.$$

$$\left|\frac{\left({\mathbf{u}}^{\prime }\cdot \nabla \right){\mathbf{u}}^{\prime }}{{\partial }_t{\mathbf{u}}^{\prime }}\right|\sim \frac{u^{\prime }k}{\omega }\sim \frac{u^{\prime }}{c},$$

$$\frac{{\partial }^{2}c}{\partial t^2}=-(c-1)+\frac{{\partial }^{3}c}{\partial x^2\partial t}-c{\left(\frac{\partial \eta }{\partial x}\right)}^2,$$

$${ϵ}^2\frac{{\partial }^2\eta }{\partial t^2}=\frac{\partial }{\partial x}\left(c(x,t{)}^2\frac{\partial \eta }{\partial x}\right),$$

$$E_{\mathrm{c}}=\frac{1}{2}{\int }_0^L\left[{\left(\frac{\partial c}{\partial t}\right)}^2+(c-1{)}^2\right]\mathrm{d}x,E_{\mathrm{w}}=\frac{1}{2}{\int }_0^L\left[{ϵ}^2{\left(\frac{\partial \eta }{\partial t}\right)}^2+c^2{\left(\frac{\partial \eta }{\partial x}\right)}^2\right]\mathrm{d}x.$$

$$\frac{\mathrm{d}\varphi }{\mathrm{d}t}=\frac{\omega (T)}{ϵ},T=t.$$

$$\frac{{\omega }^2}{{ϵ}^2}\frac{{\partial }^{2}c}{\partial {\varphi }^2}+\frac{1}{ϵ}\left[2\omega \frac{{\partial }^{2}c}{\partial \varphi \partial T}+\frac{\mathrm{d}\omega }{\mathrm{d}T}\frac{\partial \mathrm{c}}{\partial \varphi }\right]+\frac{{\partial }^{2}c}{\partial T^2}=-(c-1)-c{\left(\frac{\partial \eta }{\partial x}\right)}^2+\frac{\omega }{ϵ}\frac{{\partial }^{3}c}{\partial x^2\partial \varphi }+\frac{{\partial }^{3}c}{\partial x^2\partial T},$$

$${\omega }^2\frac{{\partial }^2\eta }{\partial {\varphi }^2}+ϵ\left[2\omega \frac{{\partial }^2\eta }{\partial \varphi \partial T}+\frac{\mathrm{d}\omega }{\mathrm{d}T}\frac{\partial \eta }{\partial \varphi }\right]+{ϵ}^2\frac{{\partial }^2\eta }{\partial T^2}=\frac{\partial }{\partial x}\left(c^2\frac{\partial \eta }{\partial x}\right).$$

$$\eta (x,\varphi ,T)={\eta }_0(x,\varphi ,T)+ϵ{\eta }_1(x,\varphi ,T)+O({ϵ}^2),$$

$$c(x,\varphi ,T)=c_0(x,\varphi ,T)+ϵc_1(x,\varphi ,T)+O({ϵ}^2),$$

$$\omega (T)={\omega }_0(T)+ϵ{\omega }_1(T)+O({ϵ}^2),$$

$${\omega }_0^2\frac{{\partial }^2{c}_2}{\partial {\varphi }^2}+\frac{{\partial }^2{c}_0}{\partial T^2}=-(c_0-1)-c_0{\left(\frac{\partial {\eta }_0}{\partial x}\right)}^2+\frac{{\partial }^3{c}_0}{\partial x^2\partial T}.$$

$$\bar{f}(x,T)={\lim_{n\to\infty}}\;\frac{1}{2n\pi}\int_{-n\pi}^{+n\pi}f(x,\phi,T)\mathrm{d}\phi,\;\;\;\;\mbox{for integer $n$}.$$

$$\frac{{\partial }^2{c}_0}{\partial T^2}=-(c_0-1)-c_0\overline{{\left(\frac{\partial {\eta }_0}{\partial x}\right)}^2}+\frac{{\partial }^3{c}_0}{\partial x^2\partial T}.$$

$$-{\omega }_0(T{)}^2\frac{{\partial }^2{\eta }_0}{\partial {\varphi }^2}+\frac{\partial }{\partial x}\left(c_0(x,T{)}^2\frac{\partial {\eta }_0}{\partial x}\right)=0.$$

$${\eta }_0(x,\varphi ,T)=A(T){\hat{\eta }}_0(x,T)\cos \varphi ,$$

$${\omega }_0^2{\hat{\eta }}_0+\frac{\partial }{\partial x}\left(c_0^2\frac{\partial {\hat{\eta }}_0}{\partial x}\right)=0,{\hat{\eta }}_0(0,T)={\hat{\eta }}_0(L,T),\frac{\partial {\hat{\eta }}_0}{\partial x}(0,T)=\frac{\partial {\hat{\eta }}_0}{\partial x}(L,T).$$

$$(f,g)={\int }_0^{L}f(x)g(x)\mathrm{d}x.$$

$${\int }_0^L{\hat{\eta }}_0^2(x,T)\mathrm{d}x=1.$$

$$(\mathcal{L}f,g)={\int }_0^{L}g(x)\left({\omega }_0^2+{\partial }_x(c_0^2{\partial }_x)\right)f(x)\mathrm{d}x={\int }_0^{L}f(x)\left({\omega }_0^2+{\partial }_x(c_0^2{\partial }_x)\right)g(x)\mathrm{d}x=(f,\mathcal{L}g).$$

$$-{\omega }_0^2\frac{{\partial }^2{\eta }_1}{\partial {\varphi }^2}+\frac{\partial }{\partial x}\left(c_0^2\frac{\partial {\eta }_1}{\partial x}\right)=2{\omega }_0\frac{{\partial }^2{\eta }_0}{\partial T\partial \varphi }+\frac{\mathrm{d}{\omega }_0}{\mathrm{d}T}\frac{\partial {\eta }_0}{\partial \varphi }+2{\omega }_0{\omega }_1\frac{{\partial }^2{\eta }_0}{\partial {\varphi }^2}-2\frac{\partial }{\partial x}\left(c_0{c}_1\frac{\partial {\eta }_0}{\partial x}\right).$$

$$\mathcal{L}{\hat{\eta }}_{1,\mathrm{c}}=F_{\mathrm{c}},$$

$$\mathcal{L}{\hat{\eta }}_{1,\mathrm{s}}=F_{\mathrm{s}},$$

$$(F_{\mathrm{c}},{\hat{\eta }}_0)=(\mathcal{L}{\hat{\eta }}_{1,\mathrm{c}},{\hat{\eta }}_0)=({\hat{\eta }}_{1,\mathrm{c}},\mathcal{L}{\hat{\eta }}_0)=({\hat{\eta }}_{1,\mathrm{c}},0)=0;$$

$${\omega }_0{\omega }_1{\int }_0^L{\hat{\eta }}_0^2\mathrm{d}x=-{\int }_0^L{\hat{\eta }}_0\frac{\partial }{\partial x}\left(c_0{c}_1\frac{\partial {\hat{\eta }}_0}{\partial x}\right)\mathrm{d}x={\int }_0^L{c}_0{c}_1{\left(\frac{\partial {\hat{\eta }}_0}{\partial x}\right)}^2\mathrm{d}x.$$

$${\omega }_1=\frac{1}{{\omega }_0}{\int }_0^L{c}_0{c}_1{\left(\frac{\partial {\hat{\eta }}_0}{\partial x}\right)}^2\mathrm{d}x.$$

$$\left(2{\omega }_0\frac{\mathrm{d}A}{\mathrm{d}T}+A\frac{\mathrm{d}{\omega }_0}{\mathrm{d}T}\right){\int }_0^L{\hat{\eta }}_0^2\mathrm{d}x=-A{\omega }_0\frac{\mathrm{d}}{\mathrm{d}T}\left({\int }_0^L{\hat{\eta }}_0^2\mathrm{d}x\right)=0;$$

$$2{\omega }_0\frac{\mathrm{d}A}{\mathrm{d}T}+A\frac{\mathrm{d}{\omega }_0}{\partial T}=0\Rightarrow \frac{\mathrm{d}}{\mathrm{d}T}\left(A^2{\omega }_0\right)=0.$$

$$\left[{\omega }_0^2+\frac{\partial }{\partial x}\left(c_0^2\frac{\partial }{\partial x}\right)\right]\frac{\delta {\hat{\eta }}_0}{\delta T}=-2\left[{\omega }_0\frac{\mathrm{d}{\omega }_0}{\mathrm{d}T}+\frac{\partial }{\partial x}\left(c_0\frac{\partial c_0}{\partial T}\frac{\partial }{\partial x}\right)\right]{\hat{\eta }}_0.$$

$$\frac{\mathrm{d}{\omega }_0}{\mathrm{d}T}=-\frac{1}{{\omega }_0}{\int }_0^L{\hat{\eta }}_0\frac{\partial }{\partial x}\left(c_0\frac{\partial c_0}{\partial T}\frac{\partial {\hat{\eta }}_0}{\partial x}\right)\mathrm{d}x=\frac{1}{{\omega }_0}{\int }_0^L{\left(\frac{\partial {\hat{\eta }}_0}{\partial x}\right)}^2{c}_0\frac{\partial c_0}{\partial T}\mathrm{d}x.$$

$$\overline{{\left(\frac{\partial {\eta }_0}{\partial x}\right)}^2}=\frac{E_0(T)}{{\omega }_0(T{)}^2}{\left(\frac{\partial {\hat{\eta }}_0}{\partial x}\right)}^2,$$

Full text

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 quasilinear (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 $\bar{\mathbf{u}}$, and of fast fluctuations, say $\mathbf{u}^{\prime}$. For instance, consider the Navier–Stokes equation for an incompressible and homogeneous fluid,

[TABLE]

where $\mathbf{u}$ is the velocity field, $P$ is the pressure, $\rho$ is the (constant) density, $\nu$ is the kinematic viscosity, and $\triangle$ 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}=\bar{\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) becomes

[TABLE]

The above set of equations is termed QL if the governing equations for the fluctuations $\mathbf{u}^{\prime}$ do not include fluctuation–fluctuation nonlinearities; i.e., if the term $\left(\mathbf{u}^{\prime}\cdot\mathbf{\nabla}\right)\mathbf{u}^{\prime}$ [and, hence, $\overline{\left(\mathbf{u}^{\prime}\cdot\mathbf{\nabla}\right)\mathbf{u}^{\prime}}$] in (3) is negligible.

Frequently, wave systems satisfy the QL approximation: denoting the wavenumber of the waves by $k$, the angular frequency by $\omega$, and the phase velocity (or “celerity”) by $c$, then

[TABLE]

which is a (generalized) Mach number. If this number is small compared to unity, then (3) reduces to QL form and the asymptotic methods introduced in this article can be applied. Since the nonlinear term $\left(\mathbf{u}^{\prime}\cdot\mathbf{\nabla}\right)\mathbf{u}^{\prime}$ is neglected, quasilinearity 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 $\bar{\mathbf{u}}$ and $\mathbf{u}^{\prime}$ is preserved: the cumulative effect of the fluctuation–fluctuation nonlinearities modifies the slow fields through the “Reynolds stress” divergence $\overline{\left(\mathbf{u}^{\prime}\cdot\mathbf{\nabla}\right)\mathbf{u}^{\prime}}$ in (2), which in turn modifies the fluctuations through terms of the form $\left(\bar{\mathbf{u}}\cdot\mathbf{\nabla}\right)\mathbf{u}^{\prime}$ and $\left(\mathbf{u}^{\prime}\cdot\mathbf{\nabla}\right)\bar{\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 $\sim 30~{}\mu\mathrm{s}$ and the streaming flow evolving on a timescale $\sim 0.1\mathrm{s}$, and subsequently shown to be described accurately by a QL dynamical system [5, 6, 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 spatiotemporally-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 $\mathbf{u}^{\prime}$ represents “eddies” and quasilinearity is realized when $\mathbf{u}^{\prime}\ll\bar{\mathbf{u}}$, hence $\left(\mathbf{u}^{\prime}\cdot\mathbf{\nabla}\right)\mathbf{u}^{\prime}\ll\left(\bar{\mathbf{u}}\cdot\mathbf{\nabla}\right)\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 quasilinear dynamical systems include strongly rotating (but weakly stratified) flows, such as open-ocean deep convection and high-latitude abyssal ocean currents[11], and rotating astrophysical disks 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[13]. 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 supplementary material.

2 Modulated Waves

2.1 Governing equations

The first system we analyze 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 $\eta$ are chosen to satisfy the following governing equations (and initial conditions):

[TABLE]

The spatial domain $x\in[0,L)$, and both $c$ and $\eta$ are $L$-periodic functions of $x$. Equation (6) is the one-dimensional (1D) wave equation with inhomogeneous celerity, where $\epsilon^{2}$ is the square of a small parameter. (Note that all variables are dimensionless.) Equation (5) 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_{\mathrm{tot}}=E_{\mathrm{c}}+E_{\mathrm{w}}$, where $E_{\mathrm{c}}$ and $E_{\mathrm{w}}$ are, respectively, the “energy” of the celerity field and of the waves, defined by

[TABLE]

For the given initial conditions, $c$ and $\eta$ are $\mathit{O}(1)$, implying the right-hand sides of (5) and (6) also are of order unity. Owing to the prefactor $\epsilon^{2}$ on the left-hand side of (6), which can be interpreted as a measure of fluctuation inertia, we expect the wave field $\eta$ to evolve on a much faster time scale than that characterizing the evolution of the celerity $c$.

2.2 Leading-order equations

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

[TABLE]

Crucially, $\phi$ and $T$ are treated as independent variables in the subsequent analysis. The introduction of the phase $\phi$ through its non-constant time derivative captures the slow evolution of the wave angular frequency $\omega$ and is referred to as the WKBJ approximation[14]. Using the chain rule, $\partial_{t}\rightarrow(\omega(T)/\epsilon)\partial_{\phi}+\partial_{T}$, the governing equations become

[TABLE]

To proceed, $c$, $\eta$, and $\omega$ are expanded as asymptotic power series in $\epsilon$,

[TABLE]

and substituted into (9)–(10). Since $\epsilon$ is small but variable, the resulting equations then can be solved order by order. At $O(1/\epsilon^{2})$, (9) yields $\partial_{\phi}^{2}c_{0}=0$, and since linear growth of $c_{0}$ with $\phi$ 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 $\partial_{\phi}c_{0}=0$, (9) at order $O(1/\epsilon)$ requires $\partial_{\phi}c_{1}=0$; i.e., $c_{1}$ also is a function only of $x$ and $T$. Finally, at $O(1)$, we obtain

[TABLE]

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

[TABLE]

This averaging operation provides a suitable definition of the overbar used to define the Reynolds decomposition $\mathbf{u}=\bar{\mathbf{u}}+\mathbf{u}^{\prime}$ in the introduction. Applying this averaging procedure to (14), and recalling $\partial_{\phi}c_{0}=0$, yields a necessary condition for the sublinear growth of $\partial_{\phi}c_{2}$:

[TABLE]

In fact, it is readily confirmed using the resulting equation for $c_{2}$ [obtained by subtracting (16) from (14)] together with the expression for $\eta_{0}$ given below in (18) that (16) ensures the boundedness of $c_{2}$ on the slow time scale.

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

[TABLE]

Equations (16)–(17) comprise a two time-scale system, for which the fluctuation equation (17) does not involve nonlinear fluctuation–fluctuation terms (e.g., $\eta_{0}^{2}$) and thus is quasilinear. (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 (17) is of the form

[TABLE]

since $c_{0}$ does not depend on $\phi$ and noting that the solution component proportional to $\sin{\phi}$ vanishes as a result of the initial conditions (6) on $\eta$. 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 (21). Crucially, at this stage of the analysis, the slow evolution of $A$ is unknown. Indeed, the ansatz (18) converts the initial-value problem (17) 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 $\mathrm{d}A/\mathrm{d}T$, which for this system is found by examining the equation for the correction field $\eta_{1}$.

2.3 Inner product and adjoint operator

Multiple scales analysis, in its most general form, requires certain notions from linear algebra that we introduce here. Substituting (18) into (17), $\hat{\eta}_{0}$ is found to satisfy

[TABLE]

This eigenvalue problem is linear, and can be expressed more compactly as $\mathcal{L}\hat{\eta}_{0}=0$, where the linear operator $\mathcal{L}=\omega_{0}^{2}+\partial_{x}(c_{0}^{2}\partial_{x})$, with periodic boundary conditions, acts on the function $\hat{\eta}_{0}$. For spatially-periodic and real-valued functions $f(x)$ and $g(x)$, we define the $L^{2}$ inner product

[TABLE]

An inner product defines a norm, and we choose to set $(\hat{\eta}_{0},\hat{\eta}_{0})=1$; that is,

[TABLE]

This relation renders the decomposition (18) unique. Moreover, after integrating by parts and making use of the periodic boundary conditions, we obtain

[TABLE]

This equality reveals that $\mathcal{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 $\mathcal{L}$ but the following procedure can be appropriately modified [14]. 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 $\eta_{1}$, which, in turn, ensures that the asymptotic expansion (11) 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[14].

2.4 Conservation of wave action

The slow temporal evolution of $A$ is obtained by collecting terms of order $O(\epsilon)$ in (10):

[TABLE]

Note that the left-hand sides of (23) and of (17) involve the same linear operator. By setting $\eta_{1}(x,\phi,T)=\hat{\eta}_{1,\mathrm{c}}(x,T)\cos{\phi}+\hat{\eta}_{1,\mathrm{s}}(x,T)\sin{\phi}$, (23) yields

[TABLE]

The notation $\delta/\delta T$ is a short-hand for $(\partial_{T}c_{0})\,\delta/\delta c_{0}$, where $\delta/\delta 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 (24) with $\hat{\eta}_{0}$ and using the fact that $\mathcal{L}$ is self-adjoint, we obtain a solvability condition:

[TABLE]

i.e., the right-hand-side of (24) must be orthogonal to the null eigenvector of the (adjoint) linear operator. Consequently,

[TABLE]

Using the normalization condition (21) then yields the following result:

[TABLE]

This expression for the angular frequency correction $\omega_{1}$, together with a relation for the evolution of the leading-order frequency $\mathrm{d}\omega_{0}/\mathrm{d}T$ (see Sec. 2.5), is required only to close the system for the correction fields $\omega_{1}$, $c_{1}$, $\eta_{1}$ arising at higher order, should they be desired. In contrast, the solvability condition $(F_{\mathrm{s}},\hat{\eta}_{0})=0$ yields a slow evolution equation for $A(T)$. Specifically, we obtain

[TABLE]

the last equality is a consequence of the normalization (21), following an interchange of the order of functional differentiation and integration. Thus, we obtain the amplitude equation

[TABLE]

The latter form of (30) 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 $E_{0}=(\omega_{0}A)^{2}/2$, as computed from (7), the adiabatic invariant derived in (30) is the so-called wave action $E_{0}/\omega_{0}$. Conservation of wave action is a generic property of all slowly-modulated, non-dissipative linear waves[15].

2.5 Slow evolution of the angular frequency

With the formalism introduced in Sec. 2.3, 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 analyzed in Sec. 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 (19) with respect to the slow time $T$:

[TABLE]

This system is of the form of $\mathcal{L}(\delta\hat{\eta}_{0}/\delta T)=G$, where $G$ is the right-hand side of (31). Thus, a solvability condition is obtained by taking the inner product of (31) with $\hat{\eta}_{0}$, yielding

[TABLE]

This result provides an explicit evolution equation for the angular frequency of the waves.

2.6 Numerical implementation

To assess the fidelity of the predicted slow dynamics to the actual system dynamics for small but finite values of $\epsilon$, direct numerical simulations (DNS) of the governing equations (5)–(6) for various values of $\epsilon$ are compared to a numerical simulation of the reduced system (16), (19), and (30) obtained from the multiple scales analysis. All simulations are performed in a domain of size $L=2\pi$ with 32 grid points, use a second-order Runge-Kutta scheme, and output the time series of the energies $E_{\mathrm{c}}$ and $E_{\mathrm{w}}$ defined in (7). Although not documented here, the convergence of each simulation with increasing spatiotemporal resolution has been confirmed.

The DNS are performed by solving the system (5)–(6) using a Fourier pseudospectral method, the timestep being set in accord with a CFL condition (e.g., $\mathrm{d}t\simeq 0.01$ for $\epsilon=0.1$). The reduced model is discretized on a Chebyshev grid, a requirement of the Dedalus eigenvalue solver; note that only (16) is explicitly time-advanced. In this equation, the term involving the fluctuations can be expressed as

[TABLE]

where $\omega_{0}(T)$ and $\hat{\eta}_{0}$ are obtained by numerically solving the eigenvalue problem (19) at each slow timestep rather than time-advancing (17) on the fast time scale, and the leading-order wave energy $E_{0}(T)=E_{0}(0)\omega_{0}(T)/\omega_{0}(0)$ owing to the conservation of wave action. For consistency with the initial conditions imposed in the DNS, we select $E_{0}(0)=\omega_{0}(0)=1$. Since the waves are not explicitly time-resolved, the timestep can be increased by an order of magnitude ($\mathrm{d}t=0.1$) relative to that required for the DNS!

The time series of the energies are reported in Fig. 1. Excellent agreement is achieved even for modest values of $\epsilon$, and, for $\epsilon=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 supplementary material. The files also can be used to generate space–time diagrams of the celerity field, should they be desired.

3 Unstable Fluctuations

3.1 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 $\Theta(x,t)$ and the fast fluctuation field $\eta(x,t)$ satisfy

[TABLE]

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 $\eta(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 $\Theta$, which evolves according to (34), the terms on the right-hand side respectively representing an external forcing $P(x,t)>0$, a linear damping ($\nu>0$ is a damping coefficient), and a quadratic feedback from the fluctuations. The small parameter $\epsilon$ controls the temporal scale separation. The coupled system (34)–(35) crudely mimics, for instance, the QL reduction of the equations governing Rayleigh–Bénard convection first proposed by Herring[16].

The coupling between the two fields $\Theta$ and $\eta$ is such that there exists an energy for this system, $E=E_{\Theta}+E_{\eta}$, where

[TABLE]

It transpires that at leading order in $\epsilon$ the nonlinear term in (35) is negligible and, thus, the total energy $E$ is conserved in the absence of forcing ($P=0$) and dissipation (requiring both $\nu=0$ in (34) and zero diffusion in (35)). In the absence of the nonlinear term in (35), the system (34)–(35) becomes quasilinear.

3.2 Leading-order equations

As for the slowly-modulated wave system considered in Sec. 2, the small parameter $\epsilon$ motivates the introduction of two scales $\phi$ and $T$, treated hereafter as independent variables and defined according to

[TABLE]

The notation $\sigma$ rather than $\omega$ is used in this example to signify that the fluctuations may grow (or decay) exponentially rather than oscillate rapidly. Like $\omega$, $\sigma$ is defined to be an eigenvalue of a certain linear operator. A crucial distinction, however, is that while $\omega$ is the slowly-varying angular frequency of any one of a countable infinity of wave modes (see Sec. 2), $\sigma$ is the slowly-evolving instantaneous growth rate of the most unstable (or least stable) fluctuation mode.

Using (37), the governing equations (34)–(35) become

[TABLE]

This set of equations is solved order by order in $\epsilon$ after substituting the following asymptotic expansions for $\sigma$, $\Theta$, and $\eta$:

[TABLE]

Equation (38) at $O(1/\epsilon)$ yields $\partial_{\phi}\Theta_{0}=0$; i.e., $\Theta_{0}$ is a function strictly of $x$ and $T$. At $O(1)$, the fast average introduced in (15) removes the term proportional to $\partial_{\phi}\Theta_{1}$, yielding

[TABLE]

The leading-order fluctuation equation [i.e. Eq. (39) at $\mathit{O}(1/\epsilon)$] is

[TABLE]

Equations (43)–(44) 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 $\partial_{T}\to(\sigma/\epsilon)\partial_{\phi}$. This system is also quasilinear: 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 quasilinearity enables a modal solution of (44) to be expressed as

[TABLE]

where $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 $\phi$ are real-valued and that a suitable normalization condition on $\hat{\eta}_{0}$ is required to render the decomposition (45) unique; see (51). Upon substituting (45) into (44), $\sigma_{0}$ is identified as the maximum eigenvalue of the operator $(\Theta_{0}+\partial_{x}^{2})$. The leading-order fluctuation equation (44) seemingly places no constraint on the evolution of $A(T)$, confirming that the system (43)–(44) is not closed on the slow time scale.

Of course, the potential exponential growth (or decay) of the leading-order fluctuation field $\eta_{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 $\Theta_{0}$:

[TABLE]

Note that, as for the fast temporal average introduced in (15), 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) $\sigma(T)=\epsilon\mathrm{d}\phi/\mathrm{d}t=0$ (implying the integrand in (46) is constant). These conditions thus require $A(T)=0$ if the maximum instantaneous growth rate $\sigma(T)$ is non-zero, a natural consequence of exponentially-fast damping for $\sigma<0$ and a requirement, were $\sigma>0$, to preclude the blow-up of the fluctuations (and, hence, of the force (46)) on the fast time scale. With $A\equiv 0$, the leading-order dynamics reduces to (43) with $\overline{\eta_{0}^{2}}=0$. Conversely, the fluctuations can have finite amplitude only if $\sigma=0$; i.e, only if the fluctuations are marginally stable. Although $\sigma$ 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) $\Theta_{0}$ is not an arbitrary function of $x$ owing to the feedback from the fluctuations, and (ii) scale-separated quasilinear 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 $\sigma_{0}$ differs from zero, then the prescription is made that the modal amplitude $A=0$. If, however, $\sigma_{0}=0$, then the coupled system (43)–(44) 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 $\sigma=0$ in the following analysis, we formally retain $\sigma_{0}$ explicitly in the derivation for clarity and generality of exposition. Since $\phi$ is constant for $\sigma=0$, the amplitude $A$ can be re-defined so that the leading-order fluctuation field

[TABLE]

and, accordingly, the fast-time average $\overline{\eta_{0}^{2}}=A^{2}\hat{\eta}_{0}^{2}$. The leading-order system therefore becomes

[TABLE]

3.3 Solvability condition

As in Sec. 2.2.3, the eigenvalue problem (48) can be cast into the form $\mathcal{L}\hat{\eta}_{0}=0$, now with $\mathcal{L}=-\sigma_{0}+\Theta_{0}+\partial^{2}_{x}$. As a result of periodicity, the linear operator $\mathcal{L}$ is self-adjoint with the inner product defined in (20); i.e.,

[TABLE]

This inner product is used to disentangle $A$ and $\hat{\eta}_{0}$ by imposing $(\hat{\eta}_{0},\hat{\eta}_{0})=1$; that is,

[TABLE]

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 Sec. 2 2.4; namely, by ensuring that the higher-order fluctuation equations are solvable. Specifically, collecting terms at $O(1)$ in (39) leads to

[TABLE]

where $\delta/\delta T=(\partial_{T}\Theta_{0})\delta/\delta\Theta_{0}$. Forming the inner product of (52) 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

[TABLE]

In the wave problem, two solvability conditions emerge from the equation for $\eta_{1}$: a closure, (28), required to evolve the dynamics of the correction fields $\omega_{1}$, $\eta_{1}$, and $c_{1}$ and an amplitude equation, (30), prescribing $\mathrm{d}A/\mathrm{d}T$ as a function of the leading-order variable $\omega_{0}$. In contrast, in the present problem, only the single constraint (53), relating the time evolution of $A$ to the higher-order field $\Theta_{1}$, is obtained. To procure a closed set of equations, it is natural to attempt to include $\Theta_{1}$ in the set of unknown fields. An equation for $\Theta_{1}$ can be derived from (38) at $O(\epsilon)$, yielding $\partial_{T}\Theta_{1}$ as a function of $\eta_{0}$ and, disappointingly, of $\eta_{1}$, another unknown variable. Thus, similarly to (28), (53) is merely a constraint required to obtain the dynamics of the correction fields $\eta_{1}$ and $\Theta_{1}$.

3.4 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 Sec. 22.5. Specifically, we differentiate the leading-order fluctuation equation (48) with respect to $T$, yielding

[TABLE]

A solvability condition is obtained by taking the inner product of (54) with $\hat{\eta}_{0}$, again utilizing the normalization condition on $\hat{\eta}_{0}$; viz.,

[TABLE]

Upon substituting (49) into (55), this solvability condition reduces to

[TABLE]

where the real coefficients $\alpha$ and $\beta$ are defined by

[TABLE]

Equation (56) prescribes the slow evolution of $\sigma_{0}$ as a function of the variables $\Theta_{0}$ and $\hat{\eta}_{0}$ and, crucially, of the fluctuation amplitude $A$. Since, for non-zero $A$, $\sigma_{0}$ must be zero, the consistency of the multiple scales expansion excludes strictly positive values of $\mathrm{d}\sigma_{0}/\mathrm{d}T$, which would immediately lead to positive growth rates; i.e., to exponentially-growing fluctuations. This scenario arises only when $\alpha>0$ and $\sigma_{0}=0$, in which case the amplitude $A(T)$ must be set to enforce $\mathrm{d}\sigma_{0}/\mathrm{d}T=0$:

[TABLE]

Equation (58) describes the effectively instantaneous (on the slow time scale) saturation of the instability via the feedback of the fluctuation field $\eta_{0}$ on the slow variable $\Theta_{0}$.

3.5 Numerical implementation

To demonstrate the efficacy of the proposed multi-scale algorithm, DNS of the governing equations (34)–(35) are compared to simulations of the multiple scales reduction given by (48), (49), and (58). All simulations are performed in a domain of size $L=2\pi$ that is discretized with 32 grid points. The damping coefficient $\nu=1$, and the external forcing

[TABLE]

A second-order Runge–Kutta time-stepping scheme is used with fixed time step $\mathrm{d}t=0.01$ for the reduced model and $\mathrm{d}t=0.005$ for the DNS. In the DNS, the governing equations (34)–(35) are solved using a Fourier pseudospectral method for $\epsilon=\{0.1,0.01\}$, while the reduced equations (valid as $\epsilon\to 0$) are discretized using Chebyshev polynomials on a non-uniformly spaced grid. For simulations of the reduced system, only the mean equation (49) is explicitly evolved in time (recall that $T=t$). At every slow timestep, however, the eigenvalue problem (48) is solved to obtain $\sigma_{0}$ and $\hat{\eta}_{0}$. Finally, the amplitude of the fluctuations $A$ is set according to

[TABLE]

where the slowly-varying coefficients $\alpha$ and $\beta$ have been defined in (57).

As evident in Fig. 2, this procedure prevents $\sigma_{0}$ from becoming positive, thereby constraining the system to evolve along a marginal stability manifold. In practice, discretization and round-off errors preclude $\sigma_{0}$ from identically equalling zero. Accordingly, for the results reported here, a tolerance of 0.01 has been enforced: $0\leq\sigma_{0}\leq 0.01$. (Smaller tolerances can be achieved by decreasing the time step $\mathrm{d}t$). The time series of the normalized fluctuation energy (half the norm of $\eta$),

[TABLE]

is also plotted in this figure. (Although not plotted here, space–time diagrams of both $\eta$ and $\Theta$ can be generated using the source files provided as supplementary material.) Excellent agreement between the DNS and the reduced model is observed for$\epsilon\simeq 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.

3.6 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.

3.6.1 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 $\sigma$ could become positive, regardless of the value and sign of the imaginary part. For the system analyzed 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:

[TABLE]

In this example, (63) describes fast oscillatory motions (waves) of frequency $O(1/\epsilon)$ that are either damped or exponentially amplified by the slow field $\Theta$. Generally, the growth rate $\sigma$ is complex, but the amplitude of the fluctuations must be set, in the multiple scales analysis, by the requirement that the real part of $\sigma$ should not increase once a state of marginal stability is reached.

3.6.2 Stabilizing nature of the fluctuation-induced feedback

Omitting the subscript ‘0’, the time-derivative of the (real) growth rate $\sigma$ obtained in (56) is given by

[TABLE]

where, again, $T$ is the slow time variable, $A$ is the amplitude of the fluctuations, and $\alpha$ and $\beta$ are two real parameters. The multiple scales approach developed in this section only applies if $\beta>0$, regardless of the value of $\alpha$. (If $\alpha>0$ and $\sigma=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 $\beta>0$ in the model analyzed here, see (57), this property is not generic; a slight modification of the set of equations (34)–(35) provides a counter-example:

[TABLE]

An analysis similar to that applied to (34)–(35) confirms that, for small values of $\epsilon$, this system is of multi-scale QL form. The corresponding expression for $\beta$, however, is modified,

[TABLE]

so that, crucially, $\beta$ is no longer sign-definite. In the scenario $\sigma=0$, $\alpha>0$, and $\beta<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 (66) 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:

[TABLE]

With these expansions, the following set of equations is obtained at leading order in $\epsilon$:

[TABLE]

where the leading-order growth rate $\sigma_{0}$ is still defined by the linear eigenvalue problem obtained from (72) 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 –$\nu\Theta_{0}$ do not contribute to the leading-order dynamics. Numerically, both (71) and (72) must be co-evolved on the fast time scale until $\sigma_{0}$ reaches zero again, and the dynamics relaxes to a slow manifold describable by the asymptotic QL reduction.

3.6.3 Degeneracy of the marginal eigenvalue

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

[TABLE]

Equation (48) then becomes

[TABLE]

subject to periodic boundary conditions. Thus, the eigenvector of largest growth rate $\sigma_{0}$ formally is equivalent to the ground state of the one-dimensional stationary Schrödinger equation, which is provably non-degenerate[17]. 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 Korteveg-de-Vries equation). Quasilinear (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 analyzed 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[18, 19], 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 [15, 20] and variational-based Hamiltonian methods[21, 22] – 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 modeling of spatially anisotropic turbulent shear flows (e.g., in wall-bounded engineering flows and in geo- and astrophysical turbulence[8, 9, 10, 11, 12]), 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 ($\epsilon\to 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 time step 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, 23, 24, 25, 26]. 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 Sec. 3 should prove valuable even for modest values of $\epsilon$.

\dataccess

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

\aucontribute

GM and GPC conceived the mathematical models, interpreted the results, and wrote the paper. GM performed the numerical simulations. Both authors gave final approval for publication.

\competing

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.

\ack

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.