1 Introduction

The Boussinesq equations are a standard approximate model of heat transfer in (viscous) fluids and are given by a coupled system of the Navier–Stokes equations and a dissipative transport equation for the temperature density:

$$\begin{aligned} \begin{aligned} \partial _tv+ v \cdot \nabla v + \nabla p&= (\nu _x \partial _x^2 + \nu _y \partial _y^2 )v + \theta e_2, \\ \partial _t\theta + v \cdot \nabla \theta&= (\mu _x \partial _x^2 + \mu _y \partial _y^2) \theta ,\\ \nabla \cdot v&=0. \end{aligned} \end{aligned}$$

Here \(v \in {\mathbb {R}}^2\) denotes the velocity, \(p \in {\mathbb {R}}\) is the pressure, \(\theta \in {\mathbb {R}}\) is the temperature and we consider the domain \({\mathbb {T}}\times {\mathbb {R}}\ni (x,y)\). The \(\theta e_2\) term models buoyancy which causes hotter fluid to rise and colder fluid to sink.

In Sects. 2 and 3 of this article, we consider the setting with only vertical dissipation of the velocity,

$$\begin{aligned} \nu _x=\mu _x=\mu _y=0, \nu _y=:\nu >0. \end{aligned}$$

We refer to this setting as vertical dissipation. In Sect. 4, we assume vertical dissipation in both velocity and temperature, which we refer to as full vertical dissipation.

One readily observes that at least formally any pair of functions of the form

$$\begin{aligned} v= (\beta y, 0), \theta =T(y), \end{aligned}$$

with \(\beta \in {\mathbb {R}}\) and T smooth are automatically stationary solutions of the vertical dissipation problem (choosing \(p=p(y)\) suitably). Here a particular focus in existing results has been on the case when T is affine and increasing, that is, hotter fluid is on top of colder fluid, which is known as hydrostatic balance. A main aim of this article is to study more general profiles T(y), and in particular answer how much T may oscillate if both shear and viscosity are available to counteract thermal instability.

More generally, the problem of partial dissipation has been an area of extensive research, where we in particular mention the recent works (Elgindi and Widmayer 2015; Widmayer 2018; Doering et al. 2018; Yang and Lin 2018; Wu et al. 2019; Deng et al. 2020; Wu et al. 2020; Dong et al. 2020). The question of global well-posedness has been addressed in series of works by Chae et al. (1999), Chae (2006).

In this article, we will focus on questions of asymptotic stability close to specific families of solutions and how the interaction of mixing and temperature stratification may counteract instability.

In Yang and Lin (2018), Yang and Lin studied the stability of the linearized inviscid problem around the case where \(T(y)=\alpha y\) is affine and showed that for some stability results it is necessary that \(\alpha >0\), and thus, T is increasing. We recall these results in Sect. 2 and emphasize that the threshold with respect to \(\alpha \) depends on whether one studies

  • the vorticity \(\omega \), which is always unstable,

  • the horizontal component of the velocity \(v_1\), which is stable if \(\alpha >0\) and unstable if \(\alpha <0\), or

  • the vertical component of the velocity \(v_2\), which is stable if \(\alpha > -2\) and unstable if \(\alpha <-2\).

Thus, already in this case in a specific sense one may allow \(\alpha \) to be negative if it is sufficiently small.

Recently, Masmoudi et al. (2020) showed that the associated nonlinear problem near \(T(y)=\alpha y\), \(\alpha >0\) without thermal diffusion but with viscous diffusion is asymptotically stable in Gevrey regularity. These results in particular show that this partial dissipation problem behaves similarly to the Euler equations (Bedrossian and Masmoudi 2015) instead of the Navier–Stokes equations (Bedrossian et al. 2018). If one instead considers full dissipation, in Zillinger (2020) we adapted the methods of Bedrossian et al. (2018), Liss (2020) to establish nonlinear stability in Sobolev regularity.

This article extends the results of Zillinger (2020) to the case of negative \(\alpha \) and partial dissipation. More precisely, we show for the linearized problem with vertical dissipation that the evolution is asymptotically stable provided

$$\begin{aligned} \alpha > -\frac{1}{100}\root 3 \of {\nu }. \end{aligned}$$

Similar results hold for T(y) non-affine. Thus, mixing enhanced dissipation can suppress Rayleigh–Bénard instability with an enhanced dependence on \(\nu \). We remark that beneficial interaction of shear and (in)stability in the context of reaction–diffusion and turbulence has previously been observed in Spiegel and Zaleski (1984), Doering and Horsthemke (1993), Castaing et al. (1989).

For the nonlinear problem, we further show that for affine T and full vertical dissipation the same stability results hold. As shown in Masmoudi et al. (2020) in the case of vertical dissipation only in the vorticity, a more careful analysis is required to control resonances, reminiscent of echoes in the Euler equations (Deng and Masmoudi 2018; Deng and Zillinger 2019).

Our main results concerning the linearized problem are summarized in the following theorem.

Theorem 1.1

Let \(T: {\mathbb {R}}\rightarrow {\mathbb {R}}\) be a given temperature profile. Let \(N \in {\mathbb {N}}\) and suppose that \(T'(y) \in L^\infty \). We then consider the linearized Boussinesq equations with vertical dissipation in the velocity only around \(v=(y,0)\), \(\theta =T(y)\) in coordinates \((x+ty, y)\):

$$\begin{aligned} \partial _t\omega&= \nu (\partial _y-t\partial _x)^2 \omega + \partial _x \theta , \\ \partial _t\partial _x \theta&= -T'(y)\partial _x v_2, \\ (t,x,y)&\in (0,\infty )\times {\mathbb {T}}\times {\mathbb {R}}\end{aligned}$$

Then, if the Fourier transform of \(T'\) satisfies the estimate

$$\begin{aligned} \sup _{\xi } \int |{\mathcal {F}}(T')(z-\xi )| \left( \frac{1+|z|}{1+|\xi |} + \frac{1+|\xi |}{1+|z|}\right) ^{N} (1+ \min (\nu ^{-\frac{2}{3}}, |z-\xi |^{\frac{2}{3}}))\mathrm{d}z < \frac{1}{100} \nu ^{1/3}, \end{aligned}$$
(1)

the initial value problem is stable in \(H^N \times H^N\) in the sense that there exists a constant \(C>0\) such that for any initial data \((\omega _{in}, \partial _x \theta _{in}) \in H^N \times H^N\) the solution satisfies

$$\begin{aligned}&\nu ^{1/3}\Vert \omega (t)\Vert _{H^N}^2 + \Vert \partial _x \theta \Vert _{H^N}^2 \le C (1+\nu ^{-\frac{2}{3}})^2 (\nu ^{1/3} \Vert \omega _{in}\Vert _{H^N}^2 +\Vert \partial _x \theta _{in}\Vert _{H^N}^2\\&+\nu ^{-2/3}\Vert \theta _{in}\Vert _{H^N}^2). \end{aligned}$$

In particular, if \(T'(y)=\alpha \), stability holds if \(\alpha > -\frac{1}{100} \nu ^{1/3}\).

The condition (1) is a sufficient condition to control commutators involving \(T'(y)\) and is probably not optimal in its dependence on N. In the case where \(T(y)=\alpha y\) is affine, it reduces to the condition \(|\alpha |< C_N \nu ^{1/3}\) and thus allows for \(\alpha \) to be negative. See Theorem 3.2 for further discussion. We remark that in our proof of Theorem 1.1 we consider a Fourier multiplier (20)

$$\begin{aligned} \sigma = {\mathcal {F}}^{-1} \sqrt{k^2+ \min ((\xi -kt)^2, \nu ^{-2/3})}{\mathcal {F}} \end{aligned}$$

applied to \(\theta \) in place of \(\partial _x\theta \). For simplicity of presentation, we estimate this multiplier in terms of \(\partial _x \theta \) and \(\nu ^{-1/3}\theta \) in the statement of the theorem.

For the nonlinear problem with full vertical dissipation, we obtain similar results.

Theorem 1.2

Let \(T: {\mathbb {R}}\rightarrow {\mathbb {R}}\) be a given temperature profile and consider the (forced) nonlinear problem around \(v=(y,0)\), \(\theta =T(y)\) with vertical dissipation \(\nu _y=\mu _y=:\nu >0\) in coordinates \((x+ty,y)\):

$$\begin{aligned} \partial _t\omega + v \cdot \nabla _t \omega&= \nu (\partial _y-t\partial _x)^2 \omega + \partial _x \theta , \\ \partial _t\theta + v \cdot \nabla _t \theta&= \nu (\partial _y-t\partial _x)^2 \theta - T'(y) v_2,\\ (t,x,y)&\in (0,\infty )\times {\mathbb {T}}\times {\mathbb {R}}, \end{aligned}$$

where \(\nabla _t= \begin{pmatrix} \partial _x\\ \partial _y-t\partial _x \end{pmatrix}\) denotes the gradient in these coordinates. Further suppose that \(T'\) satisfies the assumptions of Theorem 1.1. Then, this problem is stable in Sobolev regularity. More precisely, for any \(N \in {\mathbb {N}}\), \(N\ge 5\) there exists \(\epsilon _N=c_N \nu \) such that if initially

$$\begin{aligned} \Vert \omega \Vert _{H^N}^2 \le \epsilon ^2&< \epsilon _N^2,\\ \nu ^{-1} \Vert \partial _x \theta \Vert _{H^N}^2 +\nu ^{-5/3}\Vert \theta \Vert _{H^N}^2< \epsilon ^2&< \epsilon _{N}^2, \end{aligned}$$

then the solution remains bounded by \(10 \nu ^{-2/3} \epsilon ^2\) for all times.

We also obtain time integrability results for v, \((\partial _y-t\partial _x) \omega \) and \((\partial _y-t\partial _x)\theta \), which are stated in Sects. 3 and 4 and omitted here for brevity.

  • In the special case when \(T(y)=\alpha y\) is affine, the assumption reduces to \(|\alpha |\le \frac{1}{100}\nu ^{1/3}\).

  • Since we consider the small data regime with full vertical dissipation, the weights of \(\theta \) and \(\omega \) use different powers of \(\nu \) than in the viscous, vertical dissipation case.

  • We stress that \(\alpha \) here is allowed to be negative. As a related result in Lemma 2.2, we remark that the inviscid results of Yang and Lin (2018) extend to \(0\ge \alpha >-2\) when considering the vertical component of the velocity \(v_2\).

  • If there is no shear, then partial dissipation is not sufficient to restore stability of the vorticity for \(\alpha <0\) (see Lemma 2.1).

  • A combination of shear and vertical dissipation suffices to restore stability of the vorticity. Moreover, in that case we obtain an enhanced threshold in terms of \(-\nu ^{1/3}\).

  • These results further extend to the case of a non-affine, oscillating temperature profile T(y). In particular, we do not rely on cancellations or conserved quantities available in the hydrostatic balance case.

  • In addition to the linearized Boussinesq equations, we obtain results for the nonlinear small data problem, however, only with full vertical dissipation (considering T(y) non-affine as a solution of the forced problem). As recently shown in Masmoudi et al. (2020), this stronger assumption is probably necessary for stability in Sobolev regularity, since otherwise resonance chains may yield norm inflation.

The remainder of the article is structured as follows:

  • In Sect. 2 we recall some results for the inviscid problem, first obtained in Yang and Lin (2018), to introduce instability mechanisms and to discuss in which sense (partial) dissipation is necessary for stability results. With these motivations, we formulate four main questions Q1–Q4, which we address throughout the article.

  • In Sect. 3.1, we begin by studying the special case when T(y) is affine, where arguments are more transparent. In particular, we show that here the slope of the temperature profile can be allowed to be negative (colder fluid on top of hotter fluid) and that the size of the threshold depends on \(\nu \) with an enhanced rate.

  • In Sect. 3.2, we extend these linear results to the case of a general temperature profile T(y) satisfying suitable smallness conditions. In particular, T is allowed to oscillate.

  • Building on the linearized results, in Sect. 4 we study the nonlinear small data problem. Due to possible resonance chains, we here instead consider full vertical dissipation and consider T(y) as a solution of the forced problem. This extends previous nonlinear results in Zillinger (2020) for the affine, increasing case to possibly oscillating profiles.

1.1 Notation

Throughout this article, we consider solutions of the Boussinesq equations with vertical dissipation near the stationary solution

$$\begin{aligned} v&= \begin{pmatrix} y \\ 0 \end{pmatrix}, \theta = T(y),\\ (t,x,y)&\in (0,\infty ) \times {\mathbb {T}}\times {\mathbb {R}}. \end{aligned}$$

In this setting, it is natural to work in coordinates moving with the flow

$$\begin{aligned} (x+ty,y) \end{aligned}$$

and consider the equations satisfied by the perturbations in these coordinates.

If there is no possibility of confusion, these perturbations are again denoted as \(\omega \) and \(\theta \) and the linearized problem studied in Sect. 3 is given by

$$\begin{aligned} \partial _t\omega&= \nu (\partial _y-t\partial _x)^2 \omega + \partial _x \theta , \\ \partial _t\theta&= -T'(y)v_2,\\ v&= \nabla _t^{\perp }\Delta _t^{-1}\omega ,\\ (t,x,y)&\in (0,\infty ) \times {\mathbb {T}}\times {\mathbb {R}}. \end{aligned}$$

where

$$\begin{aligned} \nabla _t = \begin{pmatrix} \partial _x \\ \partial _y-t \partial _x \end{pmatrix}, \Delta _t= \partial _x^2 + (\partial _y-t\partial _x)^2 \end{aligned}$$

are the gradient and Laplacian in these coordinates.

In the nonlinear problem considered in Sect. 4, we additionally assume vertical dissipation also in the temperature and interpret T(y) as a solution of the forced problem. The system satisfied by the perturbation in coordinates moving with the shear is then given by

$$\begin{aligned} \partial _t\omega&= \nu (\partial _y-t\partial _x)^2 \omega +\partial _x \theta - v \cdot \nabla _t \omega , \\ \partial _t\theta&= \nu (\partial _y-t\partial _x)^2 \theta + T'(y)v_2 - v \cdot \nabla _t \theta ,\\ v&= \nabla _t^{\perp }\Delta _t^{-1}\omega ,\\ (t,x,y)&\in (0,\infty ) \times {\mathbb {T}}\times {\mathbb {R}}. \end{aligned}$$

We denote the Fourier transform of a function \(u(x,y) \in L^2 ({\mathbb {T}}\times {\mathbb {R}})\) by \({\tilde{u}}(k,\xi ) \in L^2 ({\mathbb {Z}}\times {\mathbb {R}})\) or \({\mathcal {F}} u\). Furthermore, we study several Fourier multipliers (see (18)), including

$$\begin{aligned} A(t, k,\xi )&= \exp (-2\int _0^t\frac{1}{1+(\frac{\xi }{k}-s)^2} \mathrm{d}s), \\ B(t,k,\xi )&= \exp (-2 \int _0^t\frac{1}{\sqrt{1+(\frac{\xi }{k}-s)^2}}1_{I}(s) \mathrm{d}t), \end{aligned}$$

with I being a prescribed time interval/Fourier region (see (9)):

$$\begin{aligned} I= \{t\ge 0 : |\frac{\xi }{k}-t|\le C\}, \end{aligned}$$

and C proportional to \(\nu ^{-1/3}\).

In Sect. 4, we study energy estimates on a given time interval (0, T). Since T is fixed throughout this section, we omit it from our notation and, for instance, write

$$\begin{aligned} \Vert u\Vert _{L^pH^N}:= \Vert \Vert u(t, \cdot )\Vert _{H^N}\Vert _{L^p((0,T))}. \end{aligned}$$

We write \(a \lesssim b\) if there exists a universal constant \(C>0\) such that \(|a|\le C |b|\).

2 Model Cases of Instability

In order to introduce ideas and mechanisms, in this section we recall results available in the literature for two special settings:

  • The linearized viscous problem without shear around \(\theta =\alpha y\), \(v=0\).

  • The linearized inviscid problem with shear around \(\theta =\alpha y\), \(v=(y,0)\).

Here, for simplicity we consider viscous dissipation in both horizontal and vertical directions, but no thermal dissipation. As a reference for the isolated mechanisms, the interested reader is referred to the textbook by Frisch and Yaglom [Yaglom (2012), Section 2.8.3]. We emphasize that the results of this section are not new, but serve to motivate our questions Q1–Q4 stated at the end of this section, which we address in this article. Furthermore, they show that under weaker assumptions instabilities may form and that the conditions in Theorem 1.1 are in this sense optimal.

In the case without shear, explicit solutions are available and it is known that the slope of \(\theta \) yields a sharp dichotomy between stability and exponential instability. The following basic lemma is reproduced from (Zillinger 2020, Proposition 2.6).

Lemma 2.1

Consider the linearized Boussinesq equations in vorticity formulation around

$$\begin{aligned} v=(0,0), \quad \theta = \alpha y, \end{aligned}$$

where \(\alpha \in {\mathbb {R}}\):

$$\begin{aligned} \begin{aligned}&\partial _t\omega = \nu \Delta \omega + \partial _x \theta , \\&\partial _t\theta + \alpha v_2 = \mu \Delta \theta . \end{aligned} \end{aligned}$$
(2)

Here \(v_2\) denotes the vertical component of the velocity field. Further suppose that at least one of \(\nu \) or \(\mu \) is zero. The evolution is stable if \(\alpha >0\) in the sense that for every \(N \in {\mathbb {N}}\) the energy

$$\begin{aligned} \alpha \Vert \omega \Vert _{H^N}^2 + \Vert \nabla \theta \Vert _{H^N}^2 \end{aligned}$$

is decreasing. In contrast, if \(\alpha <0\), there exist solutions which grow exponentially in time.

As we show in Lemma 2.2, when adding shear the instability for \(\alpha <0\) is significantly reduced and the evolution of \(v_2\) is even asymptotically stable if \(\alpha \) is not too large.

Proof

In the interest of accessibility, we reproduce the main steps of the proof from Zillinger (2020).

We observe that system (2) is a constant coefficient PDE, and hence, we obtain a decoupled system of ODEs for each Fourier mode with respect to x and y:

$$\begin{aligned} \partial _t\begin{pmatrix} {\tilde{\omega }} \\ {\tilde{\theta }} \end{pmatrix} = \begin{pmatrix} -\nu (k^2+\xi ^2) &{} \quad ik \\ \frac{ik \alpha }{k^2+\xi ^2} &{} \quad -\mu (k^2+\xi ^2) \end{pmatrix} \begin{pmatrix} {\tilde{\omega }} \\ {\tilde{\theta }} \end{pmatrix}, \end{aligned}$$
(3)

where we use \({\tilde{\omega }}\) to denote the Fourier transform of the vorticity and \(k \in {\mathbb {Z}}, \xi \in {\mathbb {R}}\) to denote the Fourier variables. In particular, we may study the problem at each frequency.

The case \(\alpha>\)0 Let \((k,\xi )\) and \(\alpha >0\) be given. Then, we may reformulate the problem as

$$\begin{aligned} \partial _t\begin{pmatrix} \sqrt{\alpha } {\tilde{\omega }} \\ \sqrt{k^2+\xi ^2} {\tilde{\theta }} \end{pmatrix} = \begin{pmatrix} -\nu (k^2+\xi ^2) &{} \frac{ik\sqrt{\alpha }}{\sqrt{k^2+\xi ^2}} \\ \frac{ik \sqrt{\alpha }}{\sqrt{k^2+\xi ^2}} &{} -\mu (k^2+\xi ^2) \end{pmatrix} \begin{pmatrix} \sqrt{\alpha } {\tilde{\omega }} \\ \sqrt{k^2+\xi ^2} {\tilde{\theta }} \end{pmatrix}. \end{aligned}$$

Note that the off-diagonal entries are equal and purely imaginary. Therefore, if we denote the matrix by M, it holds that \(M+{\overline{M}}^T\) is a real-valued, negative semi-definite diagonal matrix. Hence, it follows that

$$\begin{aligned} \frac{\mathrm{d}}{\mathrm{d}t} \left| \begin{pmatrix} \sqrt{\alpha } {\tilde{\omega }} \\ \sqrt{k^2+\xi ^2} {\tilde{\theta }} \end{pmatrix} \right| ^2 = \overline{\begin{pmatrix} \sqrt{\alpha } {\tilde{\omega }} \\ \sqrt{k^2+\xi ^2} {\tilde{\theta }} \end{pmatrix}} \cdot (M+{\overline{M}}^T) \begin{pmatrix} \sqrt{\alpha } {\tilde{\omega }} \\ \sqrt{k^2+\xi ^2} {\tilde{\theta }} \end{pmatrix} \le 0. \end{aligned}$$

Integrating this estimate with respect to \(\xi \) and k (possibly with respect to a weight \(\langle (k,\xi )\rangle ^N\)), it follows that

$$\begin{aligned} \alpha \Vert {\tilde{\omega }}\Vert _{L^2}^2 + \Vert \sqrt{k^2+\xi ^2}{\tilde{\theta }}\Vert _{L^2}^2 \end{aligned}$$

is non-increasing. The claimed result thus follows by Plancherel’s theorem.

The case \(\alpha<\) 0 Let \((k,\xi )\) with \(k\ne 0\) and \(\alpha <0\) be given. Then, the eigenvalues of the matrix

$$\begin{aligned} \begin{pmatrix} -\nu (k^2+\xi ^2) &{} ik \\ \frac{ik \alpha }{k^2+\xi ^2} &{} -\mu (k^2+\xi ^2) \end{pmatrix} \end{aligned}$$

are given by

$$\begin{aligned} \lambda _{1,2}&= -\frac{\nu +\mu }{2} (k^2+\xi ^2) \pm \sqrt{(\frac{\nu +\mu }{2} (k^2+\xi ^2))^2 - \nu \mu (k^2+\xi ^2)^2 - \alpha \frac{k^2}{k^2+\xi ^2}}\\&=-\frac{\nu +\mu }{2}(k^2+\xi ^2) \pm \sqrt{\left( \frac{\nu -\mu }{2}(k^2+\xi ^2) \right) ^2 - \alpha \frac{k^2}{k^2+\xi ^2}}, \end{aligned}$$

where we used the binomial formula \((a+b)^2-4ab=(a-b)^2\) in the last step.

We recall that by assumption (at least) one of \(\nu \), \(\mu \) vanishes. Therefore, we define \(C=\max (\nu , \mu )\) and observe that

$$\begin{aligned} (\mu +\nu )^2 = (\mu -\nu )^2 =:C^2 \end{aligned}$$

and that

$$\begin{aligned} \lambda _1 = - C (k^2+\xi ^2) + \sqrt{C^2 (k^2+\xi ^2) + (-\alpha )\frac{k^2}{k^2+\xi ^2}}, \end{aligned}$$

is strictly positive, since \((-\alpha )\frac{k^2}{k^2+\xi ^2}\) is positive. This matrix thus has a positive eigenvalue and there exist solutions of (3) which grow exponentially in time. Given these exponentially growing solutions on single Fourier modes, we next construct exponentially growing solutions in \(H^N\). We may pick a compact set in Fourier space, e.g., a ball, and construct initial data \((\omega _0, \theta _0) \in H^N\times H^{N+1}\) by prescribing the Fourier transform of the initial data to match these solutions (and vanish outside the ball). The corresponding solution then also exhibits exponential growth in time. \(\square \)

We remark that in the inviscid case, \(\nu =\mu =0\), these eigenvalues further simplify to

$$\begin{aligned} \pm \sqrt{-\alpha \frac{k^2}{k^{2}+\xi ^2}}, \end{aligned}$$

which are either purely imaginary if \(\alpha >0\) or positive and negative if \(\alpha <0\). Thus, if \(\alpha >0\) (hotter fluid is above), the evolution is not exponentially unstable. One speaks of hydrostatic equilibrium. The stability of this solution in the inviscid setting has recently been studied in Elgindi and Widmayer (2015), Widmayer (2018).

In contrast, if \(\alpha <0\) (that is, the fluid is hotter below), then one eigenvalue is positive and the solution is exponentially unstable. This phenomenon is known as Rayleigh–Bénard instability. One main question in the following will then be whether a shear flow can suppress this instability.

Having discussed the effects of dissipation without shear, we next consider the effects of an affine shear in the inviscid problem, where again explicit solutions are available. The following results have been previously obtained in Yang and Lin (2018), Zillinger (2020), Masmoudi et al. (2020) for \(\alpha >0\). By minor modifications of the proof, the results further extend to negative \(\alpha \) and higher Sobolev norms.

Lemma 2.2

Consider the linearized inviscid Boussinesq equations in vorticity formulation around

$$\begin{aligned} v=(y,0), \theta =\alpha y \end{aligned}$$

and coordinates \((x+ty,y)\) moving with the shear. Furthermore, define \(c= \frac{1}{2} \mathfrak {R}(\sqrt{1-4\alpha }) \in [0,\infty )\). Then, the velocity and temperature satisfy the following estimates

$$\begin{aligned} \Vert \theta \Vert _{H^N}&\lesssim t^{-1/2+c}(\Vert \omega _0\Vert _{H^{N+1}}+ \Vert \theta _0\Vert _{H^{N+2}}),\\ \Vert v_1 - \int _{{\mathbb {T}}} v_1 \mathrm{d}x\Vert _{H^N}&\lesssim t^{-1/2+c}( \Vert \omega _0\Vert _{H^{N+1}}+ \Vert \theta _0\Vert _{H^{N+2}}),\\ \Vert v_2\Vert _{H^N}&\lesssim t^{-3/2+c}( \Vert \omega _0\Vert _{H^{N+2}}+ \Vert \theta _0\Vert _{H^{N+3}}). \end{aligned}$$

The evolution of the temperature is hence stable if \(c<\frac{1}{2}\) (\(\alpha >0\)) and unstable if \(c>\frac{1}{2}\) (\(\alpha <0\)), while the vertical component of the velocity is stable if \(c<\frac{3}{2}\) (\(\alpha >-2\)).

The evolution of the vorticity in contrast is unstable for all \(\alpha \) in the sense that there exists non-trivial initial data such that

$$\begin{aligned} \Vert \omega (t)\Vert _{H^N}&\ge C t^{1/2+c}( \Vert \omega _0\Vert _{H^N}+ \Vert \partial _x \theta _0\Vert _{H^{N}}) \end{aligned}$$

as \(t\rightarrow \infty \).

We emphasize that for \(v_2\) we may allow \(0>\alpha >-2\) to be negative and that the evolution of the vorticity \(\omega \) is unstable for any \(\alpha \).

We remark that this combination of stability and instability is consistent with the Orr mechanism. More precisely, by an integration by parts argument it holds that

$$\begin{aligned} \Vert v_1-\int v_1 \mathrm{d}x \Vert _{L^2}\le C t^{-1} \Vert \omega (t)\Vert _{H^{1}}. \end{aligned}$$

Hence, if the velocity is asymptotically stable with a sharp decay rates of for instance \(t^{-1/2}\), this implies that the vorticity is algebraically unstable in \(H^{1}\) with a growth rate at least \(t^{-1/2+1}=t^{1/2}\).

Proof

As in the proof of Lemma 2.1 we consider the Fourier formulation, now in coordinates \((k, \xi +kt)\) moving with the shear:

$$\begin{aligned} \partial _t\begin{pmatrix} {\tilde{\omega }}\\ {\tilde{\theta }} \end{pmatrix} = \begin{pmatrix} 0 &{} ik \\ \frac{ik \alpha }{k^2+(\xi -kt)^2} &{} 0 \end{pmatrix} \begin{pmatrix} {\tilde{\omega }}\\ {\tilde{\theta }} \end{pmatrix}. \end{aligned}$$

Due to the vanishing diagonal structure, we may decouple this problem as

$$\begin{aligned} \partial _t^2 {\tilde{\omega }}&= - \frac{\alpha k^2}{k^2+(\xi -kt)^2} {\tilde{\omega }}, \\ \partial _t{\tilde{\omega }}&= ik {\tilde{\theta }}. \end{aligned}$$

After relabeling and shifting time by \(\frac{\xi }{k}\), we observe that the first equation corresponds to a Schrödinger equation with potential:

$$\begin{aligned} \left( \partial _t^2 + \frac{\alpha }{1+t^2}\right) u =0. \end{aligned}$$

As observed in Yang and Lin (2018), this problem can be solved explicitly in terms of hypergeometric functions:

$$\begin{aligned} \begin{aligned} u(t)&= c_1 {{\,\mathrm{_2F^1}\,}}\left( -\frac{1}{4}-\frac{1}{4}\sqrt{1-4\alpha }, -\frac{1}{4}+\frac{1}{4}\sqrt{1-4\alpha }, \frac{1}{2}, -t^2\right) \\&\quad + c_2 \ t\ {{\,\mathrm{_2F^1}\,}}\left( \frac{1}{4}-\frac{1}{4}\sqrt{1-4\alpha }, \frac{1}{4}+\frac{1}{4}\sqrt{1-4\alpha }, \frac{3}{2}, -t^2\right) . \end{aligned} \end{aligned}$$
(4)

As \(t\rightarrow \infty \), it holds that \({{\,\mathrm{_2F^1}\,}}(a,b,c,-t^2)\sim C t^{-2a}\) (see [DLMF, 15.8(ii)]). The same asymptotic behavior is exhibited by the approximate problem

$$\begin{aligned} (\partial _t^2 + \frac{\alpha }{t^2}) f=0, \end{aligned}$$

which we use to simplify discussion in the following. Making the ansatz \(f=t^{\beta }\), we obtain that

$$\begin{aligned} \begin{aligned} f&= c_1 t^{\beta _1} + c_2 t^{\beta _2}, \\ \beta _{1,2}&= \frac{1}{2}(1\pm \sqrt{1-4\alpha }), \end{aligned} \end{aligned}$$
(5)

which matches the asymptotic behavior of the hypergeometric functions in (4). In particular, we observe that for any \(\alpha \), \(\beta _1\) has positive real part which results in an algebraic instability of f and hence \({\tilde{\omega }}\). When considering the velocity and temperature, we recall that

$$\begin{aligned} ik \theta = \partial _t \omega \sim \partial _t f \end{aligned}$$

and that the Biot–Savart law combined with the shear by (y, 0) provides a gain of \(t^{-1}\) for \(v_1-\langle v_1 \rangle \) and by \(t^{-2}\) for \(v_2\) by the Orr mechanism. Hence, we deduce that

$$\begin{aligned}&\Vert v_1-\langle v_1 \rangle \Vert _{H^N} + \Vert \theta \Vert _{H^N} \sim t^{\beta _1 -1}, \\&\Vert v_2\Vert _{H^N}\sim t^{\beta _1-2} \end{aligned}$$

with \(\beta _1\) as in (5). In particular, we observe that

$$\begin{aligned} \mathfrak {R}(\beta _1 -1)&= c - \frac{1}{2}<0 \text { if } \alpha>0, \\ \mathfrak {R}(\beta _1 -1)&= c - \frac{1}{2}>0 \text { if } \alpha<0, \\ \mathfrak {R}(\beta _1 -2)&= c - \frac{3}{2}<0 \text { if } \alpha>-2, \\ \mathfrak {R}(\beta _1 -2)&= c - \frac{3}{2} >0 \text { if } \alpha <-2, \end{aligned}$$

where we used that \(\sqrt{1-4\alpha }=1 \Leftrightarrow \alpha =0\) and \(\sqrt{1-4\alpha }=3 \Leftrightarrow \alpha =-2\). \(\square \)

Given these (in)stability results, our main questions in this article are the following:

  1. Q1

    How much dissipation (and in which directions) needs to be added to restore linear stability?

  2. Q2

    Can we allow \(\alpha \) to be negative and how does the threshold depend on the dissipation?

  3. Q3

    When considering the problem without thermal dissipation, it is natural to consider the more general problem around \(v=(y,0)\), \(\theta =T(y)\). Under which conditions on T are such solutions linearly stable? For instance, can we allow T to oscillate?

  4. Q4

    Do these results extend to the nonlinear small data regime and if so how do stability regions depend the dissipation coefficients (that is, what perturbations can be considered “small”)?

In this paper, we focus on the case without thermal dissipation and \(v=(y,0)\), \(\theta =T(y)\). The converse problem without viscous dissipation and \(v=(U(y),0)\), \(\theta =\alpha y\) or time-dependent shear and temperature profile could be of future interest. We address questions Q1 and Q2 in Sect. 3.1 and Q3 in Sect. 3.2. The question Q4 of nonlinear stability is addressed in Sect. 4.

3 Shear can Counteract Hydrostatic Imbalance

Building on the results of Lemma 2.2 for a combination of Couette flow and an unstable affine temperature profile, in this section we consider the problem with partial dissipation.

More precisely, we consider the nonlinear Boussinesq equations with vertical dissipation of the velocity and without thermal diffusion:

$$\begin{aligned} \partial _tv + v \cdot \nabla v + \nabla p&= \nu \partial _y^2 v + \begin{pmatrix} 0 \\ \theta \end{pmatrix}, \\ \partial _t\theta + v \cdot \nabla \theta&= 0. \end{aligned}$$

As remarked in introduction, in Sect. 4 we additionally impose vertical thermal diffusion, but do not require it for the linear stability results of this section.

We observe that for any \(\beta \in {\mathbb {R}}\) and any function T(y), the collection

$$\begin{aligned} v&= \begin{pmatrix} \beta y \\ 0 \end{pmatrix},\\ \theta&= T(y),\\ p&=\int ^y T(s)\mathrm{d}s, \end{aligned}$$

is a stationary solution of these equations. As remarked in Sect. 2, it is natural to ask about the stability of such solutions.

In Sect. 2, we studied some related special cases when T(y) is affine:

  • In Lemma 2.1, we studied the problem with trivial shear, that is \(\beta =0\). In this setting the flow turned out to be linearly stable if T is increasing and linearly exponentially unstable if T is decreasing, even if the slope is very small.

  • In Lemma 2.2, we instead considered the case with shear but with trivial dissipation and saw that while the exponential instability is reduced to an algebraic one, the evolution of the vorticity is unstable.

In this section, we study the linearized problem first for the case of T affine (answering questions Q1, Q2) and then for general T in Sect. 3.2 (answering Q3). The nonlinear problem with full vertical dissipation is discussed in Sect. 4, which answers Q4. The author would like to thank Charlie Doering for raising the question of the stability of pairs \(v=(U(y),0), \theta = T(y)\) in a discussion.

3.1 Affine Temperature

In order to introduce ideas and mechanisms, we first study the case

$$\begin{aligned} T(y)=\alpha y, \end{aligned}$$

where we allow \(\alpha \in {\mathbb {R}}\) to be negative with a threshold depending on \(\nu \). More precisely, it turns out that for this special linearized problem we may allow \(\alpha \) to be arbitrarily large, but for the nonlinear setting of Sect. 4 and the non-affine problem we require a bound by \(\nu ^{1/3}\). Shear enhanced dissipation suppresses Rayleigh–Bénard instability in this case, thus answering questions Q1 and Q2 of Sect. 2.

We remark that results for \(\alpha \) positive have been previously established in Zillinger (2020). As the main novelties of this article, we show that even if \(\alpha \) is negative (but small) stability holds and that we may further allow T to be non-affine (see Sect. 3.2).

Theorem 3.1

Consider the linearized Boussinesq equations around \(v=(y,0), \theta =\alpha y\) in coordinates

$$\begin{aligned} (t, x-ty, y) \end{aligned}$$

moving with the shear flow:

$$\begin{aligned} \partial _t\omega&= \nu (\partial _y-t\partial _x)^2 \omega + \partial _x\theta , \\ \partial _t\theta&= -\alpha v_2, \end{aligned}$$

on the domain \({\mathbb {T}}\times {\mathbb {R}}\). Then, there exists \(\alpha _*=\frac{1}{100} \root 3 \of {\nu }\) such that the linearized evolution is stable at the level of the vorticity for any \(\alpha \) with \(|\alpha |< \alpha _*\). More precisely, for any \(N\in {\mathbb {N}}\) there exists a constant \(0<C=C(\nu ,\alpha )\) such that for all times \(t>0\) it holds that

$$\begin{aligned} \nu ^{1/3} \Vert \omega (t)\Vert _{H^N}^2+ \Vert \partial _x\theta (t)\Vert _{H^{N}}^2\le C \nu ^{1/3} \Vert \omega _0\Vert _{H^N}+ C \Vert \partial _x \theta _0\Vert _{H^{N}}^2, \end{aligned}$$

where \(\omega _0, \theta _0\) denote the initial data.

We stress that here we can allow \(\alpha \) to be negative and that for \(0< \nu <1\), the threshold \(\nu ^{1/3}\) is improved compared to the dissipative scale.

In particular, for this special setting we may even consider \(\alpha \in {\mathbb {R}}\) arbitrary, but in view of later results focus on the case of small negative \(\alpha \).

Proof of Theorem 3.1

Similar to the proof of Lemma 2.1, we may equivalently express the linearized Boussinesq equations around the affine temperature profile in Fourier variables as:

$$\begin{aligned} \partial _t\begin{pmatrix} {\tilde{\omega }}\\ k{\tilde{\theta }} \end{pmatrix} = \begin{pmatrix} -\nu (\xi -kt)^2 &{} i \\ \alpha \frac{i}{1+(\frac{\xi }{k}-t)^2} &{} 0 \end{pmatrix} \begin{pmatrix} {\tilde{\omega }}\\ k{\tilde{\theta }} \end{pmatrix}, \end{aligned}$$
(6)

where we consider coordinates \((k,\xi +kt)\) moving with Couette flow. Since the evolution of the x-averages of \(\omega \) and \(\theta \) decouples, in the following we without loss of generality only consider \(k\ne 0\).

We stress that the coefficients here are time dependent, and hence, this ODE system cannot anymore be explicitly solved in terms of a matrix exponential. However, a main advantage of the affine setting is that various estimates completely decouple, restrictions become trivial and operators commute, which makes this problem much simpler than the general profile case of Sect. 3.2 or the nonlinear problem of Sect. 4.

We note that the problem (6) decouples with respect to k and \(\xi \), which we thus in the following treat as arbitrary but fixed. We then claim that for any \(C>1\), \(\nu >0\) and any \(\alpha \in {\mathbb {R}}\) it holds that

$$\begin{aligned} |{\tilde{\omega }}(t)|^2 + k^2 |{\tilde{\theta }}(t)|^2&\le \left( 1+\frac{1}{|\alpha |}\right) (1+C^2) \exp \left( \frac{|\alpha |}{\nu C^2}\right) (|{\tilde{\omega }}(0)|^2 \nonumber \\&\quad + (k^2+\nu ^{-2/3}) |{\tilde{\theta }}(0)|^2) \end{aligned}$$
(7)

and thus the solution at time t is controlled in terms of the initial data. We observe some special cases for the exponential:

  • If we choose \(C=1\) we obtain a bound by

    $$\begin{aligned} \exp \left( \frac{|\alpha |}{\nu }\right) . \end{aligned}$$

    This bound holds for all \(\alpha \), but suggests a threshold \(|\alpha |<\nu \).

  • If we choose \(C=\nu ^{-1/2}\), we obtain a bound by

    $$\begin{aligned} \exp (|\alpha |), \end{aligned}$$

    where only the algebraic prefactors depends on \(\nu \).

  • If we choose \(C=\nu ^{-1/3}\), we obtain a bound by

    $$\begin{aligned} \exp (|\alpha | \nu ^{-1/3}), \end{aligned}$$

    where the exponent becomes uniformly bounded if we assume that \(|\alpha |<\nu ^{1/3}\).

The results of the theorem follow from the third case, where estimates in \(H^N\) are obtained by integrating the frequency-wise bound (7).

In order to introduce ideas and motivate the definition of C, we first discuss the case \(\alpha >0\).

Step 1 (symmetrize) Let \(\alpha >0\) be given. That is, suppose we are in the setting of hydrostatic balance. Then, one commonly exploited feature in the setting without shear is cancellation of the purely imaginary off-diagonal entries [compare Zillinger (2020), Doering et al. (2018)].

Indeed, consider the rescaled problem

$$\begin{aligned} \partial _t\begin{pmatrix} \sqrt{\alpha } {\tilde{\omega }}\\ \sqrt{k^2+(\xi -kt)^2}{\tilde{\theta }} \end{pmatrix} = \begin{pmatrix} -\nu (\xi -kt)^2 &{} \frac{i\sqrt{\alpha }}{\sqrt{1+(\frac{\xi }{k}-t)^2}} \\ \frac{i\sqrt{\alpha }}{\sqrt{1+(\frac{\xi }{k}-t)^2}} &{} \frac{t- \frac{\xi }{k}}{1+(\frac{\xi }{k}-t)^2} \end{pmatrix} \begin{pmatrix} \sqrt{\alpha }{\tilde{\omega }}\\ \sqrt{k^2+(\xi -kt)^2} {\tilde{\theta }} \end{pmatrix}. \end{aligned}$$

We observe that the off-diagonal entries are then exactly equal and imaginary and thus cancel under the matrix-valued map \(M \mapsto M + {\overline{M}}^{T}\).

Therefore, if we denote the square of the Euclidean norm of the vector as

$$\begin{aligned} E(t):= \alpha |{\tilde{\omega }}|^2 + (k^2+(\xi -kt)^2)|{\tilde{\theta }}|^2 \end{aligned}$$

it holds that

$$\begin{aligned} \partial _tE(t)&= - \nu (\xi -kt)^2\alpha |{\tilde{\omega }}|^2 + \frac{t- \frac{\xi }{k}}{1+(\frac{\xi }{k}-t)^2}(k^2+(\xi -kt)^2)|{\tilde{\theta }}|^2\\&\le \frac{\min (0,t- \frac{\xi }{k})}{1+(\frac{\xi }{k}-t)^2} E(t). \end{aligned}$$

Integrating in time and using that

$$\begin{aligned} \int _0^t \frac{\min (0,t- \frac{\xi }{k})}{1+(\frac{\xi }{k}-t)^2} \le \ln (1+t^2) \end{aligned}$$

it follows that

$$\begin{aligned} E(t)\le (1+t^2)E(0). \end{aligned}$$
(8)

Thus, irrespective of the size of \(\alpha >0\) and of \(\nu \ge 0\) we have shown that the evolution in \(H^N\) is at most algebraically unstable.

Step 2 (Using dissipation) Compared to our desired result, the estimate by (8) is not yet sufficient, since it is not uniform in time.

In the following, we hence modify the definition of E to also make use of the dissipation. More precisely, we introduce a cut-off

$$\begin{aligned} C> 1 \end{aligned}$$

to be specified later and define the resonant time interval

$$\begin{aligned} I:= \{t\ge 0 : |\frac{\xi }{k}-t|\le C\}. \end{aligned}$$
(9)

Then, it holds that

$$\begin{aligned} \begin{aligned}&\quad \partial _t\begin{pmatrix} \sqrt{\alpha } {\tilde{\omega }}\\ \sqrt{k^2+\min ((\xi -kt)^2, C^2)}{\tilde{\theta }} \end{pmatrix} \\&= \begin{pmatrix} -\nu (\xi -kt)^2 &{} \frac{ik\sqrt{\alpha }}{\sqrt{k^2+\min ((\xi -kt)^2, C^2)}} \\ \frac{ik\sqrt{\alpha }}{\sqrt{k^2+\min ((\xi -kt)^2, C^2)}}\frac{\sqrt{k^2+\min ((\xi -kt)^2, C^2)}}{\sqrt{k^2+(\xi -kt)^2}} &{} \frac{t- \frac{\xi }{k}}{1+(\frac{\xi }{k}-t)^2} 1_{I}(t) \end{pmatrix}\\&\begin{pmatrix} \sqrt{\alpha }{\tilde{\omega }}\\ \sqrt{k^2+\min ((\xi -kt)^2, C^2)} {\tilde{\theta }} \end{pmatrix}, \end{aligned} \end{aligned}$$
(10)

where \(1_I(t)\in \{0,1\}\) denotes the indicator function of I. We then define the modified energy as

$$\begin{aligned} E(t):= \alpha \Vert \omega \Vert _{H^N}^2 + \Vert \partial _x\theta \Vert _{H^N}^2 + \Vert \min (\xi -kt, C){\tilde{\theta }}\Vert _{H^N}^2. \end{aligned}$$
(11)

Step 2a (resonant region) If \(t \in I\) and \(\alpha >0\), the problem and the definition of E(t) are identical to the one considered in Step 1 and it follows that

$$\begin{aligned} \partial _tE(t) \le \frac{\min (0,t- \frac{\xi }{k})}{1+(\frac{\xi }{k}-t)^2} E(t). \end{aligned}$$
(12)

However, by definition of the interval I it holds that

$$\begin{aligned} \int _{I} \frac{\min (0, t - \frac{\xi }{k})}{1+(\frac{\xi }{k}-t)^2} \mathrm{d}t \le \ln (1+C^2), \end{aligned}$$

and thus, the growth of E during the resonant time is bounded by \((1+C^2)\).

Step 2b (non-resonant region) Next suppose that \(t \not \in I\) and thus \(|\frac{\xi }{k}-t|\) is large. In particular, \((\xi -kt)^2 \le k^2+ (\xi -kt)^2 \le 2 (\xi -kt)^2\), and thus, vertical dissipation is comparable to full dissipation.

Then, the off-diagonal entries in (10) can be estimated as

$$\begin{aligned} \left| \frac{i\sqrt{\alpha }}{\sqrt{1+\min ((\frac{\xi }{k}-t)^2, C^2)}} \right|&= \frac{\sqrt{\alpha }}{C}, \\ \left| \frac{i\sqrt{\alpha }}{\sqrt{1+\min ((\frac{\xi }{k}-t)^2, C^2)}}\frac{\sqrt{k^2+\min ((\xi -kt)^2, C^2)}}{\sqrt{k^2+(\xi -kt)^2}} \right|&\le \frac{\sqrt{\alpha }}{C}. \end{aligned}$$

Thus, using Young’s inequality with

$$\begin{aligned} \frac{\sqrt{\nu } \sqrt{k^2+(\xi -kt)^2}}{\sqrt{\nu }\sqrt{k^2+(\xi -kt)^2}}, \end{aligned}$$

we deduce that

$$\begin{aligned} \begin{aligned} \partial _tE(t)&\le -\frac{\nu }{2} (\xi -kt)^2 \alpha |{\tilde{\omega }}|^2 + \frac{\alpha }{\nu C^2} \frac{1}{k^2+(\xi -kt)^2} |\sqrt{k^2+\min ((\xi -kt)^2, C^2)} {\tilde{\theta }}|^2 \\&\le \frac{\alpha }{\nu C^2} \frac{1}{k^2+(\xi -kt)^2} E(t). \end{aligned} \end{aligned}$$
(13)

We note that the factor on the right-hand side is integrable in time.

Step 2c (Conclusion for \(\alpha>\)0) Combining the resonant estimate (12) and the non-resonant estimate (13), we deduce that

$$\begin{aligned} \begin{aligned} \partial _tE(t)&\le (1_{I}(t)\frac{\min (0,t- \frac{\xi }{k})}{1+(\frac{\xi }{k}-t)^2} + (1-1_{I}(t))\frac{\alpha }{\nu C^2} \frac{1}{k^2+(\xi -kt)^2}) E(t) \\ \Rightarrow E(t)&\le (1+C^2) \exp (\frac{\alpha }{\nu C^2}) E(0), \end{aligned} \end{aligned}$$

with E(t) defined in (11). The claimed estimate (7) for \(\alpha >0\) then follows by comparing E(t) with the squares of the \(H^N\) norms. It remains to discuss the case of negative \(\alpha \).

Step 3 (negative \(\alpha \) ) Let now \(\alpha <0\) be given and consider the problem rescaled by \(\sqrt{|\alpha |}\) instead. Then, the evolution equation (10) reads

$$\begin{aligned} \begin{aligned}&\quad \partial _t\begin{pmatrix} \sqrt{|\alpha |} {\tilde{\omega }}\\ \sqrt{k^2+\min ((\xi -kt)^2, C^2)}{\tilde{\theta }} \end{pmatrix} \\&= \begin{pmatrix} -\nu (\xi -kt)^2 &{} \frac{i\sqrt{|\alpha |}}{\sqrt{1+\min ((\frac{\xi }{k}-t)^2, C^2)}} \\ - \frac{i\sqrt{|\alpha |}}{\sqrt{1+\min ((\frac{\xi }{k}-t)^2, C^2)}}\frac{\sqrt{1+\min ((\frac{\xi }{k}-t)^2, C^2)}}{\sqrt{1+(\frac{\xi }{k}-t)^2}} &{} \frac{t- \frac{\xi }{k}}{1+(\frac{\xi }{k}-t)^2} 1_{I}(t) \end{pmatrix}\\&\begin{pmatrix} \sqrt{|\alpha |}{\tilde{\omega }}\\ \sqrt{k^2+\min ((\xi -kt)^2, C^2)} {\tilde{\theta }} \end{pmatrix}. \end{aligned} \end{aligned}$$
(14)

We thus define the energy as

$$\begin{aligned} E(t)= |\alpha | |{{\tilde{\omega }}}|^2 + (k^2+\min ((\xi -kt)^2, C^2)) |{\tilde{\theta }}|^2, \end{aligned}$$

which agrees with the previous definition if \(\alpha >0\).

Step 3a (non-resonant region) Suppose that \(t \not \in I\). We observe that in Step 2b we did not make use of the sign of \(\alpha \) but only used Young’s inequality. Furthermore, in that region \(|\xi -kt|\ge |k|\) and thus in this region vertical dissipation dominates full dissipation.

Hence, by the same argument we may deduce that also for our extended definition of E(t) it holds that

$$\begin{aligned} \partial _tE(t) \le \frac{\alpha }{\nu C^2} \frac{1}{k^2+(\xi -kt)^2} E(t). \end{aligned}$$

Step 3b (resonant region) Suppose that \(t \in I\). Then, we observe that off-diagonal terms in (14) are of the same size but have the opposite sign an hence do not cancel anymore. However, we may use Young’s inequality to still bound

$$\begin{aligned} \partial _tE(t)&\le -\nu (\xi -kt)^2|\alpha ||{\tilde{\omega }}|^2 + \frac{|\alpha |}{\sqrt{1^2+(\frac{\xi }{k}-t)^2}} E(t) + \frac{t- \frac{\xi }{k}}{1+(\frac{\xi }{k}-t)^2} E(t) \\&\le \frac{1+|\alpha |}{\sqrt{1^2+(\frac{\xi }{k}-t)^2}} E(t), \end{aligned}$$

which yields a bound on the total growth by

$$\begin{aligned} (1+C^2)^{1+|\alpha |}. \end{aligned}$$

Combining the estimates in the resonant and non-resonant region, we deduce that

$$\begin{aligned} E(t)\le (1+C^2)^{1+|\alpha |}\exp (\pi \frac{|\alpha |}{\nu C^2}) E(0). \end{aligned}$$

In particular, choosing \(C=\nu ^{-1/3}\) and supposing that \(|\alpha |<\min (\nu ^{1/3}, 1)\), this estimate reduces to

$$\begin{aligned} E(t)\le (1+\nu ^{-2/3})^2 e^{\pi } E(0), \end{aligned}$$

which implies the result.

We remark that in the proof for negative \(\alpha \) we have not relied on cancellation but only on smallness of \(\sqrt{|\alpha |}\) in combination with Young’s inequality. Hence, we may consider a modification of the energy E(t) as

$$\begin{aligned} |{\hat{\alpha }}| |{\tilde{\omega }}|^2 + (k^2+ \min ((\xi -kt)^2, C^2))|{\tilde{\theta }}|^2 \end{aligned}$$

with \({\hat{\alpha }}=\max (|\alpha |, \nu ^{1/3})\) and repeat the same proof, since \(\frac{|\alpha |}{\sqrt{{\hat{\alpha }}}}< \nu ^{1/6}\) and \(\sqrt{{\hat{\alpha }}}<\nu ^{1/6}\) satisfy the desired inequalities. \(\square \)

We remark that in the proof of this affine case we can allow \(\alpha \) to be arbitrarily large and are also free to choose C arbitrarily. As we discuss in the following, if \(T'\) is non-constant or if we study the nonlinear problem, smallness of \(\alpha \) is required in the proof. In view of resonances in the related linear inviscid damping problem (Deng and Zillinger 2019), some form of smallness condition is probably necessary.

3.2 Non-affine Temperature

Having discussed the setting of affine hydrostatic (im)balance, we next consider T(y) non-affine and address the question Q3 of Sect. 2 under which conditions on T in terms of \(\nu \) such solutions are stable. Here the problem does not decouple in frequency anymore and we thus employ a by now classical Cauchy–Kowalewskaya or ghost energy approach [compare Mouhot and Villani (2011), Bedrossian and Masmoudi (2015), Zillinger (2016)].

The linearized system around \(\theta =T(y)\) in Lagrangian coordinates is given by:

$$\begin{aligned} \begin{aligned} \partial _t\omega&= \nu (\partial _y-t\partial _x)^2 \omega + \partial _x \theta , \\ \partial _t\partial _x \theta&= -T'(y)\partial _x^2 \Delta ^{-1}_t \omega , \end{aligned} \end{aligned}$$
(15)

where we applied a derivative in x to the second equation. Since the evolution of the x-averages decouples, we assume without loss of generality that

$$\begin{aligned} \int \omega \mathrm{d}x = 0 = \int \theta \mathrm{d}x \end{aligned}$$

throughout this section.

Our main results are summarized in the following theorem.

Theorem 3.2

Let T(y) be a given temperature profile, \(N \in {\mathbb {N}}\) and consider the linearized Boussinesq equations (15) with vertical dissipation \(\nu >0\). Further suppose that the Fourier transform of \(T'\) satisfies

$$\begin{aligned} \int (1+|\xi |)^{N+5} |{\mathcal {F}}(T')(\xi )| \le 4^{-N} \nu ^{1/3}; \end{aligned}$$
(16)

then, for any initial data \(\omega _{0}, \theta _{0} \in H^N \times H^{N+1}\) it holds that

$$\begin{aligned}&\nu ^{1/3}\Vert \omega (t)\Vert _{H^N}^2 + \Vert \partial _x \theta (t)\Vert _{H^N}^2 \\&\quad \le C (1+\nu ^{-2/3})^{2} (\nu ^{1/3} \Vert \omega _0\Vert _{H^N}^2 + \Vert \partial _x\theta _0\Vert _{H^{N}}^2 +\nu ^{-2/3}\Vert \theta _0\Vert _{H^N}^2). \end{aligned}$$

We remark that (16) here is a sufficient condition to control several commutators. We expect that, in particular for large N, it is far from optimal and that it for instance would suffice to assume smallness for small N and only a finite norm for large N (compare Zillinger (2019)). Furthermore, if T happens to be strictly increasing, stability is expected also for large norms of \(T'\). The main focus of this theorem thus lies on cases where T may be oscillating. In the case \(T(y)=\alpha y\), the (tempered) Fourier transform is given by a Dirac measure and the condition (16) reduces to \(|\alpha |< \nu ^{1/3}\), as in Sect. 3.1.

Proof

In Sect. 3.1, we had seen that in the special case when \(T(y)=\alpha y\) is affine, the functions \(\theta , \omega \) satisfy the frequency-wise bound

$$\begin{aligned} \begin{aligned}&\quad \partial _t(|\alpha | |{\tilde{\omega }}|^2 + (k^2+ \min ((\xi -kt)^2, \nu ^{-2/3})|{\tilde{\theta }}|^2) \\&\le \left( \frac{1}{\sqrt{1+(\frac{\xi }{k}-t)^2}}1_{I} + \frac{1}{1+(\frac{\xi }{k}-t)^2}\right) (|\alpha | |{\tilde{\omega }}|^2 + (k^2+ \min ((\xi -kt)^2, \nu ^{-2/3})|{\tilde{\theta }}|^2). \end{aligned} \end{aligned}$$
(17)

Similar to the (linear) inviscid damping problem in the Euler equations, while this frequency-wise bounds fail in the general setting, an integrated version can be shown to hold more generally. More precisely, we define two Fourier weights

$$\begin{aligned} \begin{aligned} A(t,k, \xi )&= \exp \left( -2\int _0^{t}\frac{1}{1+(\frac{\xi }{k}-s)^2} \mathrm{d}s\right) , \\ B(t,k,\xi )&= \exp \left( -2 \int _0^{t}\frac{1}{\sqrt{1+(\frac{\xi }{k}-s)^2}}1_{I} \mathrm{d}s\right) , \end{aligned} \end{aligned}$$
(18)

where we included a factor 2 to have additional flexibility to absorb errors and the resonant time interval I is given as in (9):

$$\begin{aligned} I= \{t\ge 0 : |\frac{\xi }{k}-t|\le C\}. \end{aligned}$$

With slight abuse of notation, we identify these Fourier weights with their associated multipliers and for instance write Au in place of \({\mathcal {F}}^{-1}(A {\mathcal {F}} u)\).

Then, in this affine case the estimate (17) implies that if we define the energy

$$\begin{aligned} E(t)= \alpha \Vert AB \omega \Vert _{H^N}^2 + \Vert AB {\mathcal {F}}^{-1} \sqrt{k^2+ \min ((\xi -kt)^2, \nu ^{-2/3})} {\mathcal {F}} \theta \Vert _{H^N}^2, \end{aligned}$$
(19)

where \(\omega , \theta \) is a solution for \(T(y)=\alpha y\), then E(t) satisfies the decay estimate

$$\begin{aligned} \partial _tE(t)\le -&\sum \int \left( \frac{1}{\sqrt{1+(\frac{\xi }{k}-t)^2}}1_{I} + \frac{1}{1+(\frac{\xi }{k}-t)^2}\right) \\&\quad \quad (|\alpha | |AB {\tilde{\omega }}|^2 + (k^2+ \min ((\xi -kt)^2, \nu ^{-2/3})|AB {\tilde{\theta }}|^2). \end{aligned}$$

In particular, E(t) is non-increasing and the inequality \(E(t)\le E(0)\) implies the result of the theorem for the special case when T is affine.

Let now T(y) be given and for simplicity of notation define

$$\begin{aligned} \sigma = {\mathcal {F}}^{-1}\sqrt{k^2+\min ((\xi -kt)^2, \nu ^{-2/3})}{\mathcal {F}}. \end{aligned}$$
(20)

and introduce the constant \(\alpha \) in terms of operator norms:

$$\begin{aligned} \alpha := \Vert |\nabla _t|^{-1} B T' B^{-1} |\nabla _t|\Vert _{H^{N} \mapsto H^{N}} + \Vert (\partial _x^2 \Delta _t^{-1})^{-3/2} B T'(y) B^{-1} B (\partial _x^2 \Delta _t^{-1})^{3/2}\Vert _{H^N \rightarrow H^N}. \end{aligned}$$

As the last step of this proof, we will show that by (16) it follows that \(\alpha \le \nu ^{1/3}\). Similarly as in Theorem 3.1, we remark that in all the following estimates we may replace \(\alpha \) by \({\hat{\alpha }}=\max (\alpha , \nu ^{1/3})\) if \(\alpha < \nu ^{1/3}\).

We now claim that if E(t) is defined by the same formula as in (19) but with \(\omega ,\theta \) being solutions of the linearized problem with temperature profile T, then E(t) is non-increasing. This then implies the desired estimate by controlling B and \(\alpha \) (or \({\hat{\alpha }}\)) in terms of \(\nu \).

We hence have to estimate

$$\begin{aligned} \partial _tE(t)/2&= - \nu \alpha \Vert (\partial _y-t\partial _x) AB \omega \Vert _{H^N}^2 + \alpha \langle AB \omega , AB \partial _x \theta \rangle \\&\quad + \langle AB \sigma \theta , AB T'(y) B^{-1} B \partial _x^2 \Delta _t^{-1} \omega \rangle \\&\quad + \alpha \langle ({\dot{A}} B + A{\dot{B}})\omega , AB \omega \rangle \\&\quad + \langle ({\dot{A}} B + A{\dot{B}}) \sigma \theta , AB \sigma \theta \rangle \\&\quad + \langle AB {\dot{\sigma }} \theta , AB \sigma \theta \rangle \rangle . \end{aligned}$$

Here the dissipation term

$$\begin{aligned} - \nu \alpha \Vert (\partial _y-t\partial _x) AB \omega \Vert _{H^N}^2 \end{aligned}$$

and the terms due to the weights

$$\begin{aligned}&\alpha \langle ({\dot{A}} B + A{\dot{B}})\omega , AB \omega \rangle \\&\quad + \langle ({\dot{A}} B + A{\dot{B}}) \sigma \theta , AB \sigma \theta \rangle \end{aligned}$$

are non-positive and thus beneficial. Moreover, \({\dot{B}}\) was defined in such a way to control

$$\begin{aligned} \langle AB {\dot{\sigma }} \theta , AB \sigma \theta \rangle . \end{aligned}$$

More precisely, we note that inside the resonant interval I,

$$\begin{aligned} \partial _t\sigma ^2 = - 2k (\xi -kt) = \frac{-2 (\frac{\xi }{k}-t)}{1+(\frac{\xi }{k-t})^2} \sigma ^2 \end{aligned}$$

can be controlled by \(\dot{B} B\).

It thus remains to estimate

$$\begin{aligned} E_{\omega }:= \alpha \langle AB \omega , AB \partial _x \theta \rangle \end{aligned}$$

and

$$\begin{aligned} E_{\theta }:= \langle AB \sigma \theta , AB T'(y) B^{-1} B \partial _x^2 \Delta _t^{-1} \omega \rangle \end{aligned}$$
(21)

Estimating \(E_{\omega }\) Since the evolution equation for \(\omega \) does not involve \(T'(y)\), we may argue as in the affine case and control \(E_{\omega }\) frequency-wise. More precisely, for any given frequency \((k,\xi )\) we need to control

$$\begin{aligned} \left| \alpha (AB)^2(t,k,\xi ) \overline{{\tilde{\omega }}}(t,k,\xi ) ik {\tilde{\theta }}(t,k,\xi ) \right| \end{aligned}$$

Resonant region If \(t,k,\xi \) are such that \(|\frac{\xi }{k}-t|\le \nu ^{-1/3}\), we may bound this contribution by

$$\begin{aligned} \frac{\sqrt{\alpha }}{\sqrt{1+(\frac{\xi }{k}-t)}} (AB)^2(t,k,\xi ) \left[ \alpha |{\tilde{\omega }}|^2 + \left| \sqrt{1+(\frac{\xi }{k}-t)^2} {\tilde{\theta }}\right| ^2 \right] , \end{aligned}$$

which can be absorbed into

$$\begin{aligned} \langle A {\dot{B}} \omega , AB \omega \rangle + \langle A{\dot{B}} \sigma \theta , AB \sigma \theta \rangle \end{aligned}$$

by construction of B.

Non-resonant region If instead \(t,k,\xi \) are such that \(|\frac{\xi }{k}-t|\ge \nu ^{-1/3}\), then \(\dot{B}\) vanishes and we instead make use of the vertical dissipation. That is, we estimate

$$\begin{aligned}&\quad \left| \alpha (AB)^2(t,k,\xi ) \overline{{\tilde{\omega }}}(t,k,\xi ) ik {\tilde{\theta }}(t,k,\xi ) \right| \\&\le (AB)^2 \sqrt{\alpha } \sqrt{\nu }\sqrt{k^2+(\xi -kt)^2}|{\tilde{\omega }}| \frac{\sqrt{\alpha }}{\sqrt{\nu }\sqrt{k^2+(\xi -kt)^2}} \frac{1}{\sqrt{1+(\frac{\xi }{k}-t)^2}} \sqrt{k^2+\nu ^{-2/3}} |{\tilde{\theta }}| \end{aligned}$$

by the dissipation term

$$\begin{aligned} - \alpha (AB)^2 (\sqrt{\nu }\sqrt{(\xi -kt)^2}|{\tilde{\omega }}|)^2 \end{aligned}$$

and

$$\begin{aligned} (AB)^2\left( \frac{\sqrt{\alpha }}{\sqrt{\nu }\sqrt{k^2+(\xi -kt)^2}}\right) ^2 \frac{1}{1+(\frac{\xi }{k}-t)^2} (\sqrt{k^2+\nu ^{-2/3}} |{\tilde{\theta }}|)^2. \end{aligned}$$

Here we used that in the non-resonant region \((\xi -kt)^2\) controls the full dissipation.

The latter term can then be absorbed into

$$\begin{aligned} \langle {\dot{A}}B \sigma \theta , AB \sigma \theta \rangle \end{aligned}$$

provided

$$\begin{aligned} \frac{\sqrt{\alpha } }{\sqrt{\nu }\sqrt{k^2+(\xi -kt)^2}} \end{aligned}$$
(22)

is less than 1. Since we are in the non-resonant region, (22) can be bounded from above by

$$\begin{aligned} \sqrt{\alpha } \nu ^{-1/2+ 1/3} = (\alpha \nu ^{-1/3})^{1/2}, \end{aligned}$$

which is small since \(\alpha < \nu ^{1/3}\) by assumption.

Estimating \(E_{\theta }\) In order to estimate the contribution (21)

$$\begin{aligned} E_{\theta }= \langle AB \sigma \theta , AB T'(y) B^{-1} B \partial _x^2 \Delta _t^{-1} \omega \rangle \end{aligned}$$

we follow a similar argument as in the affine case. However, as \(T'\) is non-constant we further have to control an interaction term between the resonant and non-resonant regions.

More precisely, for any given time t we define the Fourier set

$$\begin{aligned} \Omega (t)=\{(k,\xi ): k\ne 0 , |\frac{\xi }{k}-t|\le \nu ^{-1/3}\}. \end{aligned}$$

That is, instead of time interval I associated with given frequencies, we consider frequencies for a given time t. We then split

$$\begin{aligned} \omega&= 1_{\Omega }\omega + (1-1_{\Omega }) \omega =: \omega _{\text {in}} + \omega _{\text {out}}, \\ \theta&= 1_{\Omega }\theta + (1-1_{\Omega }) \theta =: \theta _{\text {in}} + \theta _{\text {out}}. \end{aligned}$$

We then split the contributions as

$$\begin{aligned} \langle AB \sigma \theta , AB T'(y) B^{-1} B \partial _x^2 \Delta _t^{-1} \omega _{\text {out}} \rangle \\ \langle AB \sigma \theta _{\text {in}}, AB T'(y) B^{-1} B \partial _x^2 \Delta _t^{-1} \omega _{\text {in}} \rangle \\ \langle AB \sigma \theta _{\text {out}}, AB T'(y) B^{-1} B \partial _x^2 \Delta _t^{-1} \omega _{\text {in}} \rangle . \end{aligned}$$

We remark that in the affine case the third term identically vanished due to the disjoint Fourier support of \(\omega _{\text {in}}\) and \(\theta _{\text {out}}\), but that this orthogonality is lost in the general case.

\({{Step 2a }(\omega _{{out}})}\) We argue as in the affine case. Since \(\omega _{\text {out}}\) is supported in \(\Omega ^{c}\), it holds that

$$\begin{aligned} \Vert \nabla _t B \partial _x^2 \Delta _t^{-1} \omega _{\text {out}} \Vert _{H^N} \le \nu ^{2/3}\Vert (\partial _y-t\partial _x) B \omega \Vert _{H^N}. \end{aligned}$$

For \(\theta \), we do not need a further control of the support and may bound

$$\begin{aligned} \Vert |\nabla _t|^{-1} AB \sigma \theta \Vert _{H^N} \end{aligned}$$

by the time decay of A, provided

$$\begin{aligned} \Vert |\nabla _t| B T'(y) B^{-1} |\nabla _t|^{-1}\Vert _{H^N \rightarrow H^N}\le \sqrt{\alpha } \nu ^{1/6}. \end{aligned}$$

By our choice of \(\alpha \), the left-hand side is bounded by \(\alpha \) and this estimate is therefore satisfied provided \(\alpha <\nu ^{1/3}\), as assumed.

\({{Step 2b }(\theta _{{in}}, \omega _{{in}})}\) Similarly as in the proof of Theorem 3.1, we use the time decay of B to control this contribution. More precisely, we observe that

$$\begin{aligned}&\quad \langle AB \sigma \theta _{\text {in}}, ABT'(y)B^{-1}B \partial _x^2 \Delta _{t}^{-1} \omega _{\text {in}}\rangle \\&= \langle \sqrt{\partial _x^2 \Delta _{t}^{-1}}AB \sigma \theta _{\text {in}}, A (\sqrt{\partial _x^2 \Delta _{t}^{-1}})^{-1} BT'(y) B^{-1}\sqrt{\partial _x^2 \Delta _{t}^{-1}} \sqrt{\partial _x^2 \Delta _{t}^{-1}} B \omega _{\text {in}} \rangle . \end{aligned}$$

We may bound this contribution in terms of

$$\begin{aligned} \alpha \langle AB \omega _{\text {in}} ,\partial _x^2 \Delta _{t}^{-1} AB \omega _{\text {in}} \rangle \\ + \langle AB \sigma \theta _{\text {in}} ,\partial _x^2 \Delta _{t}^{-1} AB \sigma \theta _{\text {in}} \rangle , \end{aligned}$$

provided

$$\begin{aligned} \Vert (\sqrt{\partial _x^2 \Delta _{t}^{-1}})^{-1} B T'(y) B^{-1}\sqrt{\partial _x^2 \Delta _{t}^{-1}} \Vert _{H^N \rightarrow H^N} \le \sqrt{\alpha }. \end{aligned}$$

We remark that \(\sqrt{\partial _x^2 \Delta _{t}^{-1}}=|\partial _x||\nabla _t|^{-1}\) and that \(BT'(y)B^{-1}\) does not depend on x. Hence, this estimate is equivalent to the one of step 2a.

\({{Step 2c }(\theta _{{out}}, \omega _{{in}})}\) As \(T'\) is non-constant the contribution

$$\begin{aligned} \langle AB \sigma \theta _{\text {out}}, AB T'(y) B^{-1} B \partial _x^2 \Delta _t^{-1} \omega _{\text {in}} \rangle \end{aligned}$$

generally does not vanish. However, since \(\theta _{\text {out}}\) is supported away from the resonant region, we may insert an identity operator \((\partial _x^2 \Delta _t^{-1})^{3/2-3/2}\) and bound

$$\begin{aligned} \Vert AB (\partial _x^2 \Delta _t^{-1})^{3/2} \sigma \theta _{\text {out}}\Vert _{H^N} \le \nu ^{2/3} \sqrt{- \langle {\dot{A}}B \sigma \theta _{\text {out}}, AB \sigma \theta _{\text {out}}\rangle } \end{aligned}$$

and estimate

$$\begin{aligned} \Vert (\partial _x^2 \Delta _t^{-1})^{-3/2} AB T'(y) B^{-1} B \partial _x^2 \Delta _t^{-1} \omega _{\text {in}}\Vert _{H^N} \end{aligned}$$

by

$$\begin{aligned} \Vert (\partial _y-t\partial _x) AB \omega \Vert _{H^N}. \end{aligned}$$

This contribution can thus be absorbed by the same argument as in Step 2a, provided

$$\begin{aligned} \Vert (\partial _x^2 \Delta _t^{-1})^{-3/2} B T'(y) B^{-1} B (\partial _x^2 \Delta _t^{-1})^{3/2}\Vert _{H^N \rightarrow H^N}\le \nu ^{1/6}\sqrt{\alpha }, \end{aligned}$$

which by our definition of \(\alpha \) reduces to \(\alpha < \nu ^{1/3}\).

Step 3  (controlling \(\alpha )\) It remains to be shown that the estimate (16) controls \(\alpha \). Here we make use of Schur’s test, which controls the \(L^2\) operator norm of a map

$$\begin{aligned} u(x)\mapsto \int K(x,y) u(y) \mathrm{d}y \end{aligned}$$

by the square root of

$$\begin{aligned} \left( \sup _{x}\int |K(x,y)| \mathrm{d}y\right) \left( \sup _{y}\int |K(x,y)| \mathrm{d}y\right) . \end{aligned}$$

More precisely, we may express the map \(u \mapsto |\nabla _t| BT' B^{-1} |\nabla _t|^{-1} u\) as integration against a kernel on the Fourier side:

$$\begin{aligned}&{\tilde{u}}(k,\xi ) \mapsto \int \sqrt{k^2+(\xi -kt)^2} B(t,k,\xi ) \tilde{T'}(\xi -\zeta ) B^{-1}(t,k,\zeta )\\&(\sqrt{k^2+(\zeta -kt)^2})^{-1} {\tilde{u}}(k,\zeta ) d\zeta . \end{aligned}$$

Since we are further interested in a map on \(H^N\), we add an additional weight

$$\begin{aligned} \frac{1+|\xi |^{N}}{1+|\zeta |^N}. \end{aligned}$$

Then, Schur’s test asks us to control

$$\begin{aligned}&\sup _{\xi }\int \sqrt{k^2+(\xi -kt)^2} B(t,k,\xi ) \tilde{T'}(\xi -\zeta ) B^{-1}(t,k,\zeta )\\&(\sqrt{k^2+(\zeta -kt)^2})^{-1}\frac{1+|\xi |^{N}}{1+|\zeta |^N} d\zeta \le C_1 \end{aligned}$$

and

$$\begin{aligned} \sup _{\zeta }\int \sqrt{k^2+(\xi -kt)^2} B(t,k,\xi ) |\tilde{T'}(\xi -\zeta )| B^{-1}(t,k,\zeta ) (\sqrt{k^2+(\zeta -kt)^2})^{-1}\frac{1+|\xi |^{N}}{1+|\zeta |^N} d\xi \le C_2, \end{aligned}$$

which then bounds the \(L^2\) operator norm by \(\sqrt{C_1C_2}\).

We claim that this kernel can be bounded by \(|\tilde{T'}(\xi -\zeta )| (1+|\xi -\zeta |)^{N+5}\), at which point (16) implies that \(C_1=C_2=\nu ^{1/3}\), which concludes the proof.

Indeed, by construction of B, we can control

$$\begin{aligned} B(t,k,\xi )B^{-1}(t,k,\zeta ) \le \sqrt{1+|\xi -\zeta |^2}. \end{aligned}$$

Similarly, if \(|\xi |\le 3 |\zeta |\), we may simply control

$$\begin{aligned} \frac{1+|\xi |^{N}}{1+|\zeta |^N} \le 1+3^N. \end{aligned}$$

If instead \(|\xi |\ge 3 |\zeta |\), then

$$\begin{aligned} |\xi | \le |\xi -\zeta | + \frac{1}{3}|\xi | \Leftrightarrow |\xi | \le \frac{3}{2} |\xi -\zeta |. \end{aligned}$$

and thus

$$\begin{aligned} \frac{1+|\xi |^{N}}{1+|\zeta |^N} \le (\frac{3}{2})^N (1+|\xi -\zeta |^{N}). \end{aligned}$$

Finally, we need to control

$$\begin{aligned} \frac{k^2+(\xi -kt)^2}{k^2+(\zeta -kt)^2} = \frac{1+ (\frac{\xi }{k}-t)^2}{1+(\frac{\zeta }{k}-t)^2}. \end{aligned}$$

Here we may simply estimate

$$\begin{aligned} (\frac{\xi }{k}-t)^2 \le 2 (\frac{\zeta }{k}-t)^2 +2 (\xi -\zeta )^2. \end{aligned}$$

The first term cancels with the numerator, while for the second we simply bound by \(|\xi -\zeta |\).

Thus, in total it suffices to bound

$$\begin{aligned} \sup _{\zeta } \int |\tilde{T'}(\xi -\zeta )| (1+|\xi -\zeta |)^{N+5} d\xi <C_1=C_2, \end{aligned}$$

which is the assumption of our theorem.

\(\square \)

We remark that in the case when \(T'\) is increasing stronger results are possible, for instance allowing \(\alpha \) to be much larger, by using additional cancellations as in Sect. 3.1. The main advantage of this theorem hence lies in the fact that we can allow \(T'\) to be decreasing or oscillating.

4 The Nonlinear Equations with Vertical Dissipation

Given the results for the linearized problem, it is natural to ask whether they extend to the nonlinear perturbed problem:

$$\begin{aligned} \partial _t\omega + y \partial _x \omega + v \cdot \nabla \omega&= \nu \partial _y^2 \omega + \partial _x \theta , \\ \partial _t\theta + y \partial _x \theta + T'(y) v_2 + v \cdot \nabla \theta&=0,\\ (t,x,y)&\in (0,\infty )\times {\mathbb {T}}\times {\mathbb {R}}, \end{aligned}$$

and, if so, how this depends on \(\nu \). As shown recently by Masmoudi et al. (2020), this problem may exhibit an instability reminiscent of echo chains in the Vlasov-Poisson equations (Bedrossian 2020; Zillinger 2021) and Euler equations (Deng and Masmoudi 2018; Deng and Zillinger 2019). For this reason, we do not expect results in Sobolev regularity to extend (without strong modification). Therefore, in this section we instead consider the more viscous problem

$$\begin{aligned} \begin{aligned} \partial _t\omega + y \partial _x \omega + v \cdot \nabla \omega&= \nu \partial _y^2 \omega + \partial _x \theta , \\ \partial _t\theta + y \partial _x \theta + T'(y) v_2 + v \cdot \nabla \theta&= \nu \partial _y^2 \theta ,\\ (t,x,y)&\in (0,\infty )\times {\mathbb {T}}\times {\mathbb {R}}, \end{aligned} \end{aligned}$$

where we impose full vertical dissipation and view T(y) as a solution of the forced problem. Similar to results for the case of hydrostatic balance with shear studied in Zillinger (2020), our aim here is to extend the linear (asymptotic) stability results to the nonlinear equations with small data and thus answer question Q4 of Sect. 2.

Theorem 4.1

Let \(N\ge 5\) and suppose that the temperature profile T(y) satisfies the linear stability assumptions of Theorem 3.2. Then, there exists a constant \(c_N>0\) such that for any \(0< \epsilon < c_N \nu ^2\) and any initial data with

$$\begin{aligned} \Vert \omega _0\Vert _{H^N}^2&\le \epsilon ^2,\\ \Vert \partial _x \theta _0\Vert _{H^N}^2 +\nu ^{-2/3}\Vert \theta _0\Vert _{H^N}^2&\le \nu \epsilon ^2, \end{aligned}$$

the unique global solution with this initial data and in coordinates \((x-ty,y)\) satisfies

$$\begin{aligned} \Vert \omega \Vert _{L^\infty H^N}^2 + \nu \Vert (\partial _y-t\partial _x) \omega \Vert ^2_{L^2 H^N} + \Vert v_{\ne }\Vert _{L^2 H^N}^2&\le 10 \nu ^{-2/3}\epsilon ^2, \\ \Vert \partial _x \theta \Vert _{L^\infty H^N}^2&\le 10 \nu ^{1/3}\epsilon ^2, \end{aligned}$$

where \(L^p H^N := L^p((0,\infty ); H^N)\) and \(v_{\ne }=v- \int v \mathrm{d}x\) denotes the non-shear component of the velocity.

Remark 4.2

  • The nonlinear problem with vertical dissipation but without shear has been previously studied in Cao and Jiahong (2013), Li and Titi (2016), Adhikari et al. (2010).

  • The threshold \(\epsilon < c_N \nu ^2\) here is imposed to control losses of powers \(\nu ^{1/3}\) in enhanced dissipation estimates encoded in our Fourier multiplier B. The constant \(c_N\) can be obtained in terms of the constants of Sobolev embeddings and is such that for \(u,v,w \in H^N\), it holds that \(\Vert uvw\Vert _{H^N}\le \frac{1}{c_N}\Vert u\Vert _{H^N} \Vert v\Vert _{H^N}\Vert w\Vert _{H^N}\).

  • The nonlinear problem without thermal dissipation has been recently studied in Masmoudi et al. (2020). In particular, they require Gevrey regularity to control resonances, which suggests that stability in Sobolev regularity may either fail or require non-trivial modification Deng and Zillinger (2019); Deng and Masmoudi (2018).

  • In a previous work Zillinger (2020), we studied the special case where T(y) is affine (with positive slope) and with full dissipation. The present result allows for possibly oscillating profiles and only requires vertical dissipation.

  • We remark that we here estimate \(\sigma \theta \) instead of \(\nabla _t \theta \) or \(\partial _x \theta \). This is in view to the results of Sect. 3.1, for which we do not expect control of \(\nabla _t \theta \).

  • In view of the partial dissipation results of Sect. 3, we here omit questions of enhanced dissipation.

  • There has been extensive work on various partial dissipation regimes as well as on the inviscid problem. We discuss some of this literature in introductory Sect. 1.

Proof

We follow a classical bootstrap argument approach (Mouhot and Villani 2011; Bedrossian et al. 2018; Liss 2020) in the spirit of Cauchy–Kowalewskaya. As in Zillinger (2020) we here make use of multipliers constructed in Bedrossian et al. (2018), Liss (2020) for the Navier-Stokes and MHD problems, respectively, and adapt them to the problem at hand. In contrast with these works, we do not aim to derive (enhanced) dissipation estimates. However, we show that vertical dissipation is sufficient to employ these bootstrap methods (see also the discussion of echo chains Masmoudi et al. (2020) in Sect. 1). We remark that in Sect. 3.2 we have derived estimates for the associated linearized problem, which we use as a basis for our estimates in the following. A main challenge in the control of various contributions here will be that we can only control vertical dissipation and hence will have to separately consider regimes where horizontal dissipation would be large.

In our bootstrap construction, we consider \(L^p H^N\) norms on a time interval (0, T), \(T>0\), which incorporate a time-dependent Fourier multiplier M with

$$\begin{aligned} \nu ^{1/3}\le M \le 1, \end{aligned}$$

to be specified later (see equation (27)).

We then consider the maximal time \(T>0\) such that the following bootstrap estimates are satisfied:

$$\begin{aligned} \begin{aligned}&\Vert M \omega _{\ne }\Vert _{L^\infty _t H^N}^2 + \nu \Vert (\partial _y-t\partial _x) M \omega _{\ne } \Vert _{L^2 H^N}^2 + \nu ^{1/3} \Vert 1_{|\xi -kt|\le |k|} M \omega _{\ne }\Vert _{L^2 H^N}^2\\&\quad + \Vert {\mathcal {F}}^{-1}\frac{1}{\sqrt{1+(\frac{\xi }{k}-t)^2}}{\mathcal {F}} M \omega _{\ne } \Vert _{L^2 H^N}^2 \le 16 \epsilon ^2, \end{aligned} \end{aligned}$$
(23)
$$\begin{aligned} \begin{aligned}&\Vert \sigma M \theta _{\ne }\Vert _{L^\infty H^N}^2 + \nu \Vert (\partial _y-t\partial _x) \sigma M \theta _{\ne } \Vert _{L^2 H^N}^2 + \nu ^{1/3} \Vert 1_{|\xi -kt|\le |k|} \sigma M \theta _{\ne }\Vert _{L^2 H^N}^2\\&\quad + \Vert {\mathcal {F}}^{-1}\frac{1}{\sqrt{1+(\frac{\xi }{k}-t)^2}}{\mathcal {F}} \sigma M \theta _{\ne } \Vert _{L^2 H^N}^2 \le 16 \nu \epsilon ^2, \end{aligned} \end{aligned}$$
(24)
$$\begin{aligned} \begin{aligned}&\Vert \omega _{=}\Vert _{L^\infty _t H^N}^2 + \nu \Vert \partial _y \omega _{=}\Vert _{L^2 H^N}^2\le 16\epsilon ^2 , \end{aligned} \end{aligned}$$
(25)
$$\begin{aligned} \begin{aligned}&\Vert \sigma \theta _{=}\Vert _{L^\infty H^N}^2+ \nu \Vert \sigma \partial _y \omega _{=}\Vert _{L^2 H^N}^2 \le 16 \nu \epsilon ^2, \end{aligned} \end{aligned}$$
(26)

where \(\omega _{=}, \theta _{=}\) denote the x-averages and \(\omega _{\ne }, \theta _{\ne }\) their orthogonal complement and we recall that \(\sigma \) was defined in (20) as the Fourier multiplier

$$\begin{aligned} \sigma ={\mathcal {F}}^{-1}\sqrt{k^2+\min ((\xi -kt)^2, \nu ^{-2/3})}{\mathcal {F}}. \end{aligned}$$

Here the multiplier

$$\begin{aligned} {\mathcal {F}}^{-1}\frac{1}{\sqrt{1+(\frac{\xi }{k}-t)^2}}{\mathcal {F}} \end{aligned}$$

serves to control contributions by the velocity and

$$\begin{aligned} \nu 1_{|\xi -kt|\le |k|} \end{aligned}$$

is frequency localized in regions where vertical dissipation does not control full dissipation.

By local well-posedness and the assumed existence of a solution, there exists some positive time \(T>0\) such that (23)–(26) hold with \(L^2({\mathbb {R}}_{+}; \cdot )\) and \(L^\infty ({\mathbb {R}}_{+}; \cdot )\) replaced by \(L^2((0,T);\cdot )\) and \(L^{\infty }((0,T); \cdot )\). If the maximal time T with this property is infinity, this yields the results of the theorem in view of the bounds on M.

In the following, we thus assume for the sake of contradiction that \(T<\infty \) is maximal. We will then show that at the time T none of the estimates (23)–(26) attain equality. Therefore, by continuity the estimates are still satisfied for a slightly larger time, which contradicts the maximality and thus implies the result.

In order to introduce ideas, let us first consider the x-averages. We remark that in the linearized results of Sect. 3 their evolution decoupled and reduced to heat evolution. Thus, in the following we have to control the effects of the nonlinearity, where the lack of full dissipation requires us to introduce some additional splittings.

\({\textit{Estimating }\,\omega _{=}}\) We observe that \(\partial _x \theta \), \(v_{\ne }\cdot \nabla _t\omega _{=}\) and \(v_{=}\cdot \nabla _t\omega _{\ne }\) all possess a vanishing x-average and thus obtain the following evolution equation for \(\omega _{=}\):

$$\begin{aligned} \partial _t\omega _{=} +0 + (v_{\ne } \cdot \nabla _t \omega _{\ne })_{=}= \nu \partial _y^2 \omega _{=} + 0. \end{aligned}$$

Testing this equation with \(\omega _{=}\) and integrating in time, we deduce that

$$\begin{aligned} \Vert \omega _{=}(T)\Vert _{H^N}^2 + 2\nu \int _{0}^T \Vert \partial _y \omega _{=}\Vert _{H^N}^2 = \Vert \omega _{=}(0)\Vert _{H^N}^2 +2 \int _0^T \langle \omega _{=},v_{\ne } \cdot \nabla _t \omega _{\ne } \rangle . \end{aligned}$$

We recall that by assumption the initial data are of size much smaller than \(\sqrt{8}\epsilon \). Thus, if we can show that the integral on the right-hand side is bounded by \(\epsilon ^2\), this implies that equality in (25) is indeed not attained here.

As we assume only vertical dissipation, we first discuss the part involving y derivatives of \(\omega _{\ne }\):

$$\begin{aligned}&\quad \quad \quad \int _0^T \langle \omega _{=},v_{\ne }^{2} (\partial _{y}-t\partial _x) \omega _{\ne } \rangle \\&\quad \quad \le C_N \Vert \omega _{=}\Vert _{L^\infty H^N} \Vert v_{\ne }^2\Vert _{L^2 H^N} \Vert (\partial _{y}-t\partial _x) \omega _{\ne }\Vert _{L^2 H^N} \\&\underset{(25), (23), (24)}{\le } C_N \cdot 4\epsilon \cdot 4\nu ^{-1/3}\epsilon \cdot \nu ^{-1/2} 4\nu ^{-1/3} \epsilon \\&\quad \quad = \left( 4\nu ^{-7/6}C_N \epsilon \right) 16 \epsilon ^2, \end{aligned}$$

where \(v_{\ne }^2\) denotes the vertical component of the velocity field, \(C_N\) is the constant of the algebra property of \(H^N\) and the loss of factors \(\nu ^{-1/3}\) is due to the multiplier M. Since by assumption \(4\nu ^{-7/6}\epsilon \) is much smaller than 1, this term is too small to help achieve equality.

For the term involving x-derivatives, we introduce a Fourier multiplier \(\chi \) which corresponds to the projection onto the set

$$\begin{aligned} \{ (k,\xi ): |\xi -kt|\ge |k| \}. \end{aligned}$$

We note that this set is the complement of set considered in (23). Moreover, by construction it holds that

$$\begin{aligned} \Vert \chi \partial _x \omega _{\ne }\Vert _{L^2 H^N} \le \Vert \chi (\partial _{y}-t\partial _x) \omega _{\ne }\Vert _{L^2 H^N} \le \Vert (\partial _{y}-t\partial _x) \omega _{\ne }\Vert _{L^2 H^N}, \end{aligned}$$

which thus allows for an estimate of the same form as for the part involving y derivatives.

Finally, we estimate

$$\begin{aligned}&\quad \quad \quad \int _0^T \langle \omega _{=},v_{\ne }^{1} \partial _x (1-\chi ) \omega _{\ne } \rangle \\&\quad \quad = -\int _0^T \langle \omega _{=},(\partial _x v_{\ne }^{1}) (1-\chi ) \omega _{\ne } \rangle \\&\quad \quad \le C_N\Vert \omega _{=}\Vert _{L^\infty H^N} \Vert \partial _x v_{\ne }^1\Vert _{L^2 H^N} \Vert (1-\chi ) \omega _{\ne }\Vert _{L^2 H^N} \\&\underset{25, 23, 23}{\le } C_N 4\epsilon \cdot 4 \nu ^{-1/3} \epsilon \cdot 4 \nu ^{-1/6-1/3}\epsilon \\&\quad \quad = \left( C_N \nu ^{-5/6} 4\epsilon \right) 16\epsilon . \end{aligned}$$

Here the loss of powers \(\nu ^{-1/3}\) corresponds to the lower bound on M, and we observe that \(1-\chi = 1_{|\xi -kt|\le |k|}\). As this contribution is also much smaller than \(16 \epsilon ^2\), we conclude that

$$\begin{aligned} \Vert \omega _{=}(T)\Vert _{H^N}^2 + \nu \int _{0}^T \Vert \partial _y \omega _{=}\Vert _{H^N}^2 < 8 \epsilon ^2 \end{aligned}$$

and thus equality in (23) is not attained.

\({\textit{Estimating }\,\theta _{=}}\) Before discussing \(\sigma \theta _{=}\), we consider \(\theta _{=}\), where we can argue analogously to the case of \(\omega _{=}\). We may test the equation

$$\begin{aligned} \partial _t\theta _{=} + (v_{\ne }\cdot \nabla _t \theta _{\ne })_{=} = \nu \partial _y^2 \theta _{=} \end{aligned}$$

with \(\theta _{=}\) and integrate in time to again derive an integral estimate. We then estimate the contribution

$$\begin{aligned} \int _0^T \langle \theta _{=}, v_{\ne }\cdot \nabla _t \theta _{\ne } \rangle \end{aligned}$$

by a constant times

$$\begin{aligned} \Vert \theta _{=}\Vert _{L^\infty H^N} \Vert v_{\ne }^2\Vert _{L^2 H^N} \Vert (\partial _y-t\partial _x) \theta _{\ne }\Vert _{L^2 H^N}\\ + \Vert \theta _{=}\Vert _{L^\infty H^N} \Vert v_{\ne }^1\Vert _{L^2 H^N} \Vert (\partial _y-t\partial _x) \chi \theta _{\ne }\Vert _{L^2 H^N}\\ + \Vert \theta _{=}\Vert _{L^\infty H^N} \Vert \partial _x v_{\ne }^1\Vert _{L^2 H^N} \Vert (1-\chi ) \theta _{\ne }\Vert _{L^2 H^N}. \end{aligned}$$

By the bootstrap assumptions (23)–(26) this sum can be controlled in terms of \(\nu ^{-7/6} \epsilon ^3\), which is much smaller than \(\epsilon ^2\).

\({\textit{Estimating }\,\sigma \theta _{=}}\) We may extend the definition of \(\sigma \) to purely y-dependent functions as the Fourier multiplier \(\sigma = {\mathcal {F}}^{-1}\min (|\xi |, \nu ^{-1/3}) {\mathcal {F}}\). We note that the operator norm of \(\sigma \) is bounded by \(\nu ^{-1/3}\) and thus \(\sigma \theta _{=}\) could be controlled in terms of \(\theta _{=}\). However, in this way we would pass from a bound by \(\epsilon ^2\) to one by \(\nu ^{-2/3}\epsilon ^2\), which is insufficient for our bootstrap approach. Instead, we aim to show that by a similar argument as above \(\Vert \sigma \theta _{=}\Vert _{L^\infty H^N}\) can be controlled by \(\epsilon ^2\), where the loss of powers of \(\nu \) only factors into the smallness conditions on \(\epsilon \) used to control nonlinear interaction terms.

We may control

$$\begin{aligned}&\quad \Vert \sigma \theta _{=}(T)\Vert _{H^N}^2 +\nu \Vert \partial _y \sigma \theta _{=}\Vert _{L^2 H^N}^2 \\&= \Vert \sigma \theta _{=}(0)\Vert _{H^N}^2 + \int _0^T \langle \sigma \theta _{=}, \sigma v_{\ne }\cdot \nabla _t \theta _{\ne } \rangle \\&\le \Vert \sigma \theta _{=}(0)\Vert _{H^N}^2 + \Vert \sigma \theta _{=}\Vert _{L^\infty H^N} \nu ^{-1/3} C_N(\Vert v_{\ne }^2\Vert _{L^2 H^N} \Vert (\partial _y-t\partial _x) \theta _{\ne }\Vert _{L^2 H^N} \\&\quad +\Vert v_{\ne }^1\Vert _{L^2 H^N} \Vert (\partial _y-t\partial _x) \chi \theta _{\ne }\Vert _{L^2 H^N} + \Vert \partial _x v_{\ne }^1\Vert _{L^2 H^N} \Vert (1-\chi ) \theta _{\ne }\Vert _{L^2 H^N})\\&\le \Vert \sigma \theta _{=}(0)\Vert _{H^N}^2 + C_N4\nu ^{1/2} \epsilon \nu ^{-1/3} (4 \nu ^{-1/3}\epsilon \cdot 4 \nu ^{-1/2-1/3} \nu ^{1/2}\epsilon \\&\quad + C_N 4 \nu ^{-1/3}\epsilon \cdot \nu ^{-1/2/-1/3}\nu ^{1/2}\epsilon + 4\nu ^{-1/3}\epsilon \cdot 4\nu ^{-1/2-1/3} \nu ^{1/2}\epsilon )\\&\le \nu \epsilon ^2 + (16 \nu ^{-3/2}\epsilon + 16 \nu ^{-3/2}\epsilon + 16 \nu ^{-3/2}\epsilon )C_N \cdot 16\nu \epsilon ^2 \end{aligned}$$

Thus, by assumption on \(\epsilon \) and the initial data, equality in (26) is also not achieved for \(\sigma \theta _{=}\).

Estimating \(\omega _{\ne }\,{ and }\,\theta _{\ne }\) Having discussed the control of the x-averages, we now turn to control \(\omega _{\ne }, \sigma \theta _{\ne }\). Here we will first focus on contributions due to T(y) and the x-averages and finally discuss the control of the nonlinearity involving \(v_{\ne }\).

We recall that \(\omega _{\ne }\) and \(\theta _{\ne }\) satisfy the system

$$\begin{aligned} \partial _t\omega _{\ne }&= \nu (\partial _y-t\partial _x)^2 \omega _{\ne } + \partial _x \theta _{\ne } - v_{=}^1\partial _x \omega _{\ne } - v_{\ne }^2 \partial _y \omega _{=} - (v_{\ne }\cdot \nabla _t \omega _{\ne })_{\ne }, \\ \partial _t\theta _{\ne }&= \nu (\partial _y-t\partial _x)^2 \omega _{\ne } + T'(y) v_{\ne }^2 - v_{=}^1 \partial _x \theta _{\ne } - v_{\ne }^2 \partial _y \theta _{=}- (v_{\ne }\cdot \nabla _t \theta _{\ne })_{\ne }, \end{aligned}$$

where we consider \(\omega _{=}\) and \(\theta _{=}\) as given functions.

In the linearized problem of Sect. 3.2, we could without loss of generality assume that \(\omega _{=}=\theta _{=}=0\) and constructed a non-increasing energy functional. In the following, we build on these estimates and integrate them in time to show that the controls (23), (24) are stable under small nonlinear perturbations.

We recall the multipliers AB defined in (18) in Sect. 3.2:

$$\begin{aligned} A(t,k, \xi )&= \exp \left( -2\int _0^{t}\frac{1}{1+(\frac{\xi }{k}-s)^2} \mathrm{d}s\right) , \\ B(t,k,\xi )&= \exp \left( -2 \int _0^{t}\frac{1}{\sqrt{1+(\frac{\xi }{k}-s)^2}} 1_{I} \mathrm{d}s\right) , \end{aligned}$$

and for simplicity of notation write

$$\begin{aligned} M:=AB. \end{aligned}$$
(27)

Let us first study the time derivative of

$$\begin{aligned} \Vert M \omega _{\ne }\Vert _{H^N}^2. \end{aligned}$$

Then, it holds that

$$\begin{aligned}&\Vert M \omega _{\ne }(T)\Vert _{H^N}^2 - \int _0^T \langle M \omega _{\ne }, {\dot{M}} \omega _{\ne } \rangle + \nu \int _0^T \Vert M (\partial _y-t\partial _x)\omega _{\ne }\Vert _{H^N}^2 \nonumber \\&= \Vert M \omega _{\ne }(0)\Vert _{H^N}^2 + \int _0^T\langle M \omega _{\ne }, M( v_{=}^1\partial _x \omega _{\ne }) \rangle \nonumber \\&\quad + \int _0^T\langle M \omega _{\ne }, M (v_{\ne }^2 \partial _y \omega _{=}) \rangle \nonumber \\&\quad + \int _0^T \langle M \omega _{\ne }, M \partial _x \theta _{\ne }\rangle \nonumber \\&\quad + \int _0^T \langle M \omega _{\ne }, M (v_{\ne }\cdot \nabla _t \omega _{\ne }) \rangle \nonumber \\&=:\Vert M \omega _{\ne }(0)\Vert _{H^N}^2 + {\mathcal {T}}_{v_{=}^1} + {\mathcal {T}}_{\omega _{=}} + {\mathcal {T}}_{\omega _{\ne }, \theta _{\ne }} + {\mathcal {T}}_{v_{\ne }}. \end{aligned}$$
(28)

By assumption, \(\Vert M \omega _{\ne }(0)\Vert _{H^N}^2\) is much smaller than \(\epsilon ^2\), so if we can show that various terms \({\mathcal {T}}\) on the right-hand side can be controlled by the left-hand side and higher powers of \(\epsilon \), we can show that the left-hand side remains smaller than \(4\epsilon ^2\) for all times.

\({{Estimating}\,\,{\mathcal {T}}_{v^1_{=}}}\) In order to estimate \({\mathcal {T}}_{v^1_{=}}\), we make use of cancellation in an integration by parts, following a similar argument as in Zillinger (2020) with additional adjustments to account for partial dissipation. More precisely, given the multiplier M, we note that by Parseval’s identity

$$\begin{aligned}&\quad \langle M \omega _{\ne }, M (v_{=}^1 \partial _x \omega _{\ne }) \rangle \\&= \langle M \omega _{\ne }, M (v_{=}^1 \partial _x \omega _{\ne }) - v_{=}^1 \partial _x M \omega _{\ne }\rangle \\&= \sum \int \int M(t,k,\xi ) \omega _{\ne }(k,\xi )\omega _{\ne }(k,\xi +\zeta ) ik (M(t,k,\xi )- M(t,k,\xi +\zeta )) v_{=}(\zeta ). \end{aligned}$$

This cancellation is required to control \(v_{=}(\zeta )= \frac{1}{i\zeta }\omega _{=}(\zeta )\) in terms of \(\omega _{=}\). In particular, if \(|\zeta |\ge 1\) this control is trivial, while for \(|\zeta |\le 1\) we observe that M(tkz) is Lipschitz with respect to \(\frac{z}{k}\) uniformly in t:

$$\begin{aligned} |M(t,k,\xi )- M(t,k,\xi +\zeta )| \le C |\frac{\zeta }{k}|. \end{aligned}$$

Hence, we can control \({\mathcal {T}}_{v^1_{=}}\) by

$$\begin{aligned} \Vert M \omega _{\ne }\Vert _{L^2 H^N} \Vert \omega _{\ne }\Vert _{L^2 H^N} \Vert \omega _{=}\Vert _{L^\infty H^N}. \end{aligned}$$

The last factor is controlled by the preceding argument. For the first two factors, we make the observation that

$$\begin{aligned} \nu ^{1/3} \le \nu (k^2+ (\xi -kt)^2) + \frac{1}{\sqrt{1+(\frac{\xi }{k}-t)^2}} 1_{|\frac{\xi }{k}-t|\le \nu ^{-1/3}} \end{aligned}$$

and hence \(\Vert M \omega _{\ne }\Vert _{H^N}^2\) (and \(\nu ^{2/3}\Vert \omega _{\ne }\Vert _{H^N}^2\)) can be estimated in terms of the dissipation and the decay due to \({\dot{M}}\), at a loss of a factor \(\nu ^{1/3}\).

\({{Estimating}\,\,{\mathcal {T}}_{\omega _{=}}}\) We next discuss

$$\begin{aligned} {\mathcal {T}}_{\omega _{=}}= \langle M \omega _{\ne }, M (v_{\ne }^2 \partial _y \omega _{=}) \rangle . \end{aligned}$$

Here we may easily estimate by

$$\begin{aligned} \nu ^{-2/3} \Vert M \omega _{\ne }\Vert _{L^\infty H^N} \Vert M v_{\ne }^2 \Vert _{L^2 H^N} \Vert \partial _y \omega _{=}\Vert _{L^2 H^N}, \end{aligned}$$

where the factor of \(\nu ^{-2/3}\) corresponds to a rough bound of the operator norm of M. All factors are controlled in terms of the bootstrap assumption, and thus, \({\mathcal {T}}_{\omega _{=}}\) is much smaller than \(\epsilon ^2\) provided \(\epsilon ^3\) is much smaller than \(\epsilon ^2\) in terms of powers of \(\nu \).

\({{Estimating}\,\,{\mathcal {T}}_{\omega _{\ne }, \theta _{\ne }}}\) As one of the main results of Sect. 3.2, we have shown that \(M=AB\) is constructed in just such a way that

$$\begin{aligned} |\langle AB \omega _{\ne }, AB \partial _x \theta _{\ne }\rangle | \le - \langle M \omega _{\ne }, {\dot{M}} \omega _{\ne } \rangle - \alpha ^{-1} \langle M \sigma \theta _{\ne }, {\dot{M}} \sigma \theta _{\ne } \rangle \end{aligned}$$

with \(\alpha =\max (\Vert T'\Vert , \nu ^{1/3})\) (see Theorem 3.2 for the precise definition). Hence, we can absorb this contribution into the left-hand side of (28), provided we can control \({\dot{M}} \sigma \theta _{\ne }\), which will be the left-hand side of a later equation (32).

\({{Estimating}\,\,{\mathcal {T}}_{v_{\ne }, \theta _{\ne }}}\) It remains to discuss the main nonlinearity, where a key challenge is given by the lack of horizontal dissipation.

If we had full dissipation at our disposal, this estimate would reduce to controlling by

$$\begin{aligned} \Vert \omega _{\ne }\Vert _{L^\infty H^N}\Vert v_{\ne }\Vert _{L^2 H^N} \Vert \nabla _t \omega _{\ne }\Vert _{L^2 H^N}. \end{aligned}$$

However, as we only require vertical dissipation, the last factor is not easily controlled anymore. We thus have to invest additional effort to control this contribution.

As \(v_{\ne }\) is divergence-free, we observe that

$$\begin{aligned} \langle M \omega _{\ne }, M(v_{\ne } \cdot \nabla _t \omega _{\ne }) \rangle&= \langle M \omega _{\ne }, M(v_{\ne } \cdot \nabla _t \omega _{\ne }) - v_{\ne } \cdot \nabla _t M \omega _{\ne } \rangle \\&= \sum \iiint M(k,\xi ) {\tilde{\omega }}_{\ne }(k,\xi ) (M(k,\xi )- M(k-l,\xi -\zeta )) {\tilde{v}}_{\ne }(l,\zeta ) \\&\quad \cdot \begin{pmatrix} k-l \\ \xi +\zeta - (k-l)t \end{pmatrix} {\tilde{\omega }}_{\ne }(k-l,\xi -\zeta ). \end{aligned}$$

We observe that if

$$\begin{aligned} |k-l|\le \nu ^{-1} |\xi +\zeta - (k-l)t| \end{aligned}$$

the last gradient can simply be controlled by the vertical dissipation, which yields an estimate in terms of

$$\begin{aligned} \Vert \omega _{\ne }\Vert _{L^\infty H^N}\Vert v_{\ne }\Vert _{L^2 H^N} \Vert (\partial _y-t\partial _x) \omega _{\ne }\Vert _{L^2 H^N} \end{aligned}$$

and can hence be controlled. Similarly, if

$$\begin{aligned} |k-l|\le \nu ^{-1} |l| \end{aligned}$$

we can control in terms of

$$\begin{aligned} \Vert \omega _{\ne }\Vert _{L^\infty H^N}\Vert \partial _x v_{\ne }\Vert _{L^2 H^N} \Vert \omega _{\ne }\Vert _{L^2 H^N}. \end{aligned}$$

It thus only remains to discuss the region where

$$\begin{aligned} \begin{aligned} |t- \frac{\xi +\zeta }{k-l}|&\le \nu , \\ |l|&\le \nu |k|. \end{aligned} \end{aligned}$$
(29)

Here, we make use of cancellations in M. More precisely, we note that \(M(k,\xi )\) does not depend on k and \(\xi \) individually, but only on \(\frac{\xi }{k}\) and that uniformly in time

$$\begin{aligned}&|M(k,\xi )- M(k-l,\xi -\zeta )| \le C |\frac{\xi }{k}- \frac{\xi -\zeta }{k-l}|\\&= C |\frac{\xi -kt}{k}- \frac{\xi -\zeta - (k-l)t}{k-l}|\\&\le C \frac{1}{1+\nu } \frac{1}{|k-l|}(|\xi -kt| + |\xi -\zeta -(k-l)t|), \end{aligned}$$

where we used (29). We thus can control \({\mathcal {T}}_{v_{\ne }, \theta _{\ne }}\) in that region by

$$\begin{aligned} \Vert \omega _{\ne }\Vert _{L^\infty H^N}\Vert v_{\ne }\Vert _{L^2 H^N}\Vert (\partial _y-t\partial _x) \omega _{\ne }\Vert _{L^2 H^N}, \end{aligned}$$

which is controlled in terms of the bootstrap estimates.

It remains to discuss

$$\begin{aligned} \Vert {\mathcal {F}}^{-1}\frac{1}{\sqrt{1+(\frac{\xi }{k}-t)^2}}{\mathcal {F}} M \omega _{\ne } \Vert _{L^2 H^N}^2 \end{aligned}$$
(30)

and

$$\begin{aligned} \nu ^{1/3} \Vert 1_{|\xi -kt|\le |k|} M \omega _{\ne }\Vert _{L^2 H^N}^2. \end{aligned}$$
(31)

The contribution by (30) here by construction can be absorbed into the decay of A. In order to control (31), we use a combination of the decay of B and the vertical dissipation. More precisely, recalling the definition of B it suffices to show that

$$\begin{aligned} \nu ^{1/3}1_{|\xi -kt|\le |k|} \le \nu (\xi -kt)^2 + \frac{2}{\sqrt{1+(\frac{\xi }{k}-t)^2}} 1_{|\frac{\xi }{k}-t|\le \nu ^{-1/3}}. \end{aligned}$$

Indeed, if \(|\frac{\xi }{k}-t|\ge \nu ^{-1/3}\), then \(\nu (\xi -kt)^2\ge \nu \cdot \nu ^{-2/3}\ge \nu ^{1/3}\). If instead \(|\frac{\xi }{k}-t|\ge \nu ^{-1/3}\), then \(\frac{2}{\sqrt{1+(\frac{\xi }{k}-t)^2}}\ge \nu ^{1/3}\) and we can hence estimate by the latter term.

This concludes the improvement in (23).

\({\textit{Controlling }\,\sigma \theta _{\ne }}\) We next turn to controlling \(\sigma \theta _{\ne }\), where we study the time derivative of

$$\begin{aligned} \Vert M \sigma \theta _{\ne }\Vert _{H^N}^2. \end{aligned}$$

Integrating in time, we have to control

$$\begin{aligned} \begin{aligned}&\Vert M \sigma \theta _{\ne }(T)\Vert _{H^N}^2 - \int _0^T \langle M \sigma \theta _{\ne }, {\dot{M}} \sigma \theta _{\ne } \rangle + \nu \int _0^T \Vert M (\partial _y-t\partial _x)\sigma \theta _{\ne }\Vert _{H^N}^2 \\&= \Vert M \sigma \theta _{\ne }(0)\Vert _{H^N}^2 + \int _0^T \langle M \sigma \theta _{\ne }, M \sigma T'(y) v_{\ne }^2\rangle + \int _0^T \langle M \sigma \theta _{\ne }, M \sigma v_{=}^1 \partial _x \theta _{\ne } \rangle \\&\quad + \int _0^T \langle M \sigma \theta _{\ne }, M \sigma v_{\ne }^2 \partial _y \theta _{=} \rangle \\&\quad + \int _0^T \langle M \sigma \theta _{\ne }, M \sigma v_{\ne }\cdot \nabla _t \theta _{\ne }\rangle \\&=: \Vert M \sigma \theta _{\ne }(0)\Vert _{H^N}^2 + {\mathcal {T}}_{T}+ {\mathcal {T}}_{v_{=}^1} + {\mathcal {T}}_{\theta _{=}} + {\mathcal {T}}_{v_{\ne }, \sigma \theta _{\ne }}. \end{aligned} \end{aligned}$$
(32)

Here the aim again is to show that all \({\mathcal {T}}\) contributions add up to something smaller than \(\epsilon ^2\), and hence, equality is not attained.

\({{Estimating}\,\,{\mathcal {T}}_{T}}\) As one of the main results of Sect. 3.2, we have shown that \({\mathcal {T}}_{T}\) can be controlled in terms of the decay of the multipliers M and the vertical dissipation of \(\omega \) only. Thus, this contribution can estimated in terms of the left-hand side of (32) and (28).

\({{Estimating}\,\,{\mathcal {T}}_{v_{=}^1}}\) Here we may argue analogously as for \(\omega _{\ne }\), expect that M has been replaced by \(\sigma M\). We thus obtain an estimate by

$$\begin{aligned} \Vert M \sigma \theta _{\ne }\Vert _{L^2 H^N} \Vert \theta _{\ne }\Vert _{L^2 H^N} \Vert \omega _{=}\Vert _{L^\infty H^N}. \end{aligned}$$

\({{Estimating}\,\,{\mathcal {T}}_{\theta _{=}}}\) Here we may argue again analogously as for \(\omega _{\ne }\) and control by

$$\begin{aligned} \Vert M \sigma \theta _{\ne }\Vert _{L^\infty H^N}\Vert v_{\ne }^2\Vert _{L^2 H^N} \Vert \partial _y \theta _{=}\Vert _{L^2 H^N}. \end{aligned}$$

\({{Estimating}\,\,{\mathcal {T}}_{v_{\ne }, \sigma \theta _{\ne }}}\) We recall that in this theorem we assume vertical dissipation also for the temperature (in contrast with Sect. 3.2 and the problem considered in Masmoudi et al. (2020)). Therefore, in this estimate we argue largely analogously to the estimate of \({\mathcal {T}}_{v_{\ne }, \omega _{\ne }}\). However, since \(\sigma M\) also depends on k, we need some additional control in the region where the horizontal dissipation is not easily controlled.

More precisely, by the preceding arguments for \({\mathcal {T}}_{v_{\ne }, \omega _{\ne }}\) it suffices to consider

$$\begin{aligned}&\sum \iint (\sigma \theta _{\ne })(k,\xi ) \frac{1}{\sigma (k-l,\xi -\zeta )} (\sigma M (k,\xi )- \sigma M (k-l,\xi -\zeta )){\tilde{v}}_{\ne }(l,\zeta ) \\&\quad \cdot \begin{pmatrix} k-l\\ \xi -\zeta +(k-l)t \end{pmatrix} (\sigma \theta _{\ne })(k-l,\xi -\zeta ) \end{aligned}$$

in the regions where \(\xi -\zeta \) is very close to resonant and l is much smaller than k.

However, in that case we may split into differences in M and in \(\sigma \) and observe that

$$\begin{aligned} \frac{\sqrt{k^2+ (\xi -kt)^2} - \sqrt{(k-l)^2+ (\xi +\zeta - (k-l)t)^2}}{\sqrt{(k-l)^2+ (\xi +\zeta - (k-l)t)^2}}\\ \approx \frac{\sqrt{k^2}- \sqrt{(k-l)^2}}{\sqrt{(k-l)^2}} \approx \frac{l}{|k-l|}, \end{aligned}$$

where we could neglect \(\xi -kt\) and \(\xi +\zeta - (k-l)t\) since these terms could otherwise be controlled in terms of the vertical dissipation. Hence, over all we can control by

$$\begin{aligned} \Vert \sigma \theta _{\ne }\Vert _{L^\infty H^N}\Vert \partial _x v_{\ne }\Vert _{L^2 H^N} \Vert \sigma \theta _{\ne }\Vert _{L^2 H^N}, \end{aligned}$$

which concludes the proof (24). \(\square \)