STABLE NEARLY SELF-SIMILAR BLOWUP OF THE 2D BOUSSINESQ
AND 3D EULER EQUATIONS WITH SMOOTH DATA
JIAJIE CHEN AND THOMAS Y. HOU
Abstract.
Inspired by the numerical evidence of a potential 3D Euler singularity [54, 55],
we prove finite time blowup of the 2D Boussinesq and 3D axisymmetric Euler equations with
smooth initial data of finite energy and boundary. There are several essential difficulties in
proving finite time blowup of 3D Euler with smooth initial data. One of the essential diffi-
culties is to control a number of nonlocal terms that do not seem to offer any damping effect.
Another essential difficulty is that the strong advection normal to the boundary introduces a
large growth factor for the perturbation if we use weighted
L
2
estimates. We overcome this
difficulty by using a combination of a weighted
L
∞
norm and a weighted
C
1
/
2
norm, and
develop sharp functional inequalities using the symmetry properties of the kernels and some
techniques from optimal transport. Moreover we decompose the linearized operator into a
leading order operator plus a finite rank operator. The leading order operator is designed in
such a way that we can obtain sharp stability estimates. The contribution from the finite rank
operator can be captured by an auxiliary variable and its contribution to linear stability can
be estimated by constructing approximate solution in space-time. This enables us to establish
nonlinear stability of the approximate self-similar profile and prove stable nearly self-similar
blowup of the 2D Boussinesq and 3D Euler equations with smooth initial data and boundary.
1.
Introduction
The question whether the 3D incompressible Euler equations can develop a finite time singu-
larity from smooth initial data of finite energy is one of the most outstanding open questions in
the theory of nonlinear partial differential equations and fluid dynamics. The main difficulty is
due to the presence of the vortex stretching term in the vorticity equation:
(1.1)
ω
t
+
u
·∇
ω
=
ω
·∇
u
,
where
ω
=
∇ ×
u
is the
vorticity vector
of the fluid, and
u
is related to
ω
via the
Biot-
Savart law
. The velocity gradient
∇
u
formally has the same scaling as vorticity
ω
. Thus the
vortex stretching term has a nonlocal quadratic nonlinearity in terms of vorticity. However, the
nonlocal nature of the vortex stretching term can lead to dynamic depletion of the nonlinear
vortex stretching, which could prevent a finite time blowup, see e.g. [18, 26, 46]. The interested
readers may consult the excellent surveys [17, 35, 44, 49, 56] and the references therein.
In this paper, we study the blowup of the axisymmetric Euer equations with smooth initial
data and boundary. Denote by
ω
θ
,
u
θ
and
φ
θ
the angular vorticity, angular velocity, and angular
stream function, respectively. The 3D axisymmetric Euler equations are given below:
(1.2)
∂
t
(
ru
θ
) +
u
r
(
ru
θ
)
r
+
u
z
(
ru
θ
)
z
= 0
, ∂
t
(
ω
θ
r
) +
u
r
(
ω
θ
r
)
r
+
u
z
(
ω
θ
r
)
z
=
1
r
4
∂
z
((
ru
θ
)
2
)
,
where the radial velocity
u
r
and the axial velocity
u
θ
are given by the Biot-Savart law:
(1.3)
−
(
∂
rr
+
1
r
∂
r
+
∂
zz
)
φ
θ
+
1
r
2
φ
θ
=
ω
θ
, u
r
=
−
φ
θ
z
, u
z
=
φ
θ
r
+
1
r
φ
θ
,
with the no-flow boundary condition
φ
θ
(1
,z
) = 0 on the solid boundary
r
= 1 and a periodic
boundary condition in
z
. For 3D Euler blowup that occurs at the boundary
r
= 1, we know that
the axisymmetric Euler equations have the same scaling properties as those of the 2D Boussinesq
equations [56]. Thus, we also study the 2D Boussinesq equations on the upper half space:
ω
t
+
u
·∇
ω
=
θ
x
,
(1.4)
θ
t
+
u
·∇
θ
= 0
,
(1.5)
Date
: October 20, 2022.
1
arXiv:2210.07191v2 [math.AP] 19 Oct 2022
2
JIAJIE CHEN AND THOMAS Y. HOU
where the velocity field
u
= (
u,v
)
T
:
R
2
+
×
[0
,T
)
→
R
2
+
is determined via the Biot-Savart law
(1.6)
−
∆
φ
=
ω, u
=
−
φ
y
, v
=
φ
x
,
where
φ
is the stream function with the no-flow boundary condition
φ
(
x,
0) = 0 at
y
= 0. By
making the change of variables
̃
θ
,
(
ru
θ
)
2
,
̃
ω
=
ω
θ
/r
, we can see that
̃
θ
and ̃
ω
satisfy the 2D
Boussinesq equations up to the leading order for
r
≥
r
0
>
0.
The main result of this paper is a rigorous proof of the nearly self-similar blowup of the 2D
Boussinesq and the 3D axisymmetric Euler equations with smooth initial data and boundary.
The blowup mechanism of the 2D Boussinesq equations with boundary is essentially the same as
that of the 3D axisymmetric incompressible Euler equations with boundary. Our work is inspired
by the computation of Luo-Hou [54, 55] in which they presented some convincing numerical
evidence that the 3D axisymmetric Euler equations with smooth initial data and boundary
develop a potential finite time singularity. Inspired by the recent breakthrough of Elgindi [27]
(see also [28]) on the blowup of the axisymmetric Euler equations without swirl for
C
1
,α
initial
velocity, we have proved asymptotically self-similar blowup of the 2D Boussinesq equations and
the nearly self-similar blowup of the 3D axisymmetric Euler equations with
C
1
,α
velocity and
boundary in [12]. The blowup analysis presented in [12] takes advantage of the
C
1
,α
velocity in an
essential way and does not generalize to prove the Hou-Luo blowup scenario with smooth initial
data. The results presented in this paper provide the first rigorous proof of stable nearly self-
similar blowup of the 2D Boussinesq and 3D Euler equations with smooth data and boundary.
1.1.
Some main ingredients in our analysis.
We follow a general strategy that we have
established in our previous works [12–14]. We first reformulate the problem of studying the
finite time blowup of the 2D Boussinesq and 3D Euler equations as the problem of proving
nonlinear stability of an approximate steady state of the dynamic rescaling formulation. A
very important first step is to construct an approximate steady state of the dynamic rescaling
formulation with sufficiently small residual errors. We achieve this by decomposing the solution
into a semi-analytic part that captures the far field behavior of the solution and a numerically
computed part that has compact support. The approximate steady state gives an approximate
self-similar profile. See more discussions in Section 5. We remark that there has been some
recent exciting development of using a physics-informed neural network (PINN) to construct an
approximate steady state of the 2D Boussinesq equations, see [75].
Establishing linear stability of the approximate steady state is the most crucial step in our
blowup analysis. One essential difficulty is that the advection normal to the boundary for
smooth initial data introduces a large growth factor if we use weighted
L
2
energy estimates
similar to [12–14, 27], see more discussions in Section 2. Such difficulty is absent in the analysis
of 1D models [13, 14] and in the case of
C
1
,α
initial velocity [12] since we gain a small factor
α
for the advection term. To overcome the destabilizing effect due to advection normal to the
boundary, we choose a weighted
L
∞
norm, which allows us to extract the maximal amount of
damping from the local terms without suffering from the destabilizing effect due to advection
normal to the boundary. In order to close the energy estimates, we use a combination of the
weighted
L
∞
norm and the weighted
C
1
/
2
norm.
To estimate the nonlocal terms effectively, we derive sharp
C
1
/
2
estimates for
∇
u
using the
symmetry properties of the kernels and some techniques from optimal transport. We note that
novel functional inequalities on similar Biot-Savart laws have played a crucial role in [27, 48].
The authors identify the main term in the Biot-Savart law, which has a much simpler structure,
and use it to control the velocity effectively. The sharp H ̈older estimates play a similar role
in our work. We decompose the Biot-Savart law into two parts. The main terms capture the
most singular part of the Biot-Savart law, and the remaining terms are more regular. We can
generalize these estimates to the weighted sharp
C
1
/
2
estimate since the commutators between
the singular weights and the velocity operators are more regular. Compared with [27, 48], our
remaining terms are not small and cannot be controlled by
a-priori
estimates.
We decompose the linearized equations into three parts. The first part is the main part
and differs from the linearized equations by a finite rank operator
K
. We choose
K
to capture
the main contributions from the remaining terms in the Biot-Savart law and estimate
K
in
STABLE BLOWUP OF 2D BOUSSINESQ EQUATIONS
3
the second part. The stability analysis for the first part can be established using the new
functional inequalities. The second part accounts for the contributions from
K
to the linearized
equations. Thanks to the finite rank structure of
K
, we perform space-time estimates with
computer assistance to obtain sharp stability estimates for the second part. The third part is
small and can be treated as a small perturbation. Note that the first part is not small since the
linearized operator is not compact. By integrating these estimates, we establish linear stability.
The main results of this paper are stated by the two informal theorems below. The more
precise and stronger statement of Theorem 1 will be given by Theorem 3 in Section 2 and the
precise statement of Theorem 2 will be given Theorem 4 in Section 6.
Theorem 1.
There is a family of smooth initial data
(
θ
0
,ω
0
)
with
θ
0
being even and
ω
0
being
odd, such that the solution of the Boussinesq equations
(1.4)
-
(1.6)
develops a singularity in
finite time
T <
+
∞
. The velocity field
u
0
has finite energy. The blowup solution
(
θ
(
t
)
,ω
(
t
))
is nearly self-similar in the sense that
(
θ
(
t
)
,ω
(
t
))
with suitable dynamic rescaling is close to
an approximate blowup profile
(
̄
θ,
̄
ω
)
up to the blowup time. Moreover, the blowup is stable for
initial data
(
θ
0
,ω
0
)
close to
(
̄
θ,
̄
ω
)
in some weighted
L
∞
and
C
1
/
2
norm.
Theorem 2.
Consider the 3D axisymmetric Euler equations in the cylinder
r,z
∈
[0
,
1]
×
T
.
The solution of the 3D Euler equations
(1.2)
-
(1.3)
develops a nearly self-similar blowup (in the
sense described in Theorem 1) in finite time for some smooth initial data
ω
θ
0
,
u
θ
0
supported away
from the axis
r
= 0
. The initial velocity field has finite energy,
u
θ
0
and
ω
θ
0
are odd and periodic in
z
. The blowup is stable for initial data
(
u
θ
0
,ω
θ
0
)
that are close to the approximate blowup profile
( ̄
u
θ
,
̄
ω
θ
)
after proper rescaling subject to some constraint on the initial support size.
1.2.
A novel framework of analysis with computer assistance.
One of our main contri-
butions is to introduce a novel framework of analysis that enables us to obtain sharp stability
estimates by combining sharp functional inequalities, energy estimates, and approximate space-
time solutions constructed numerically with rigorous error control. Such errors are treated as
small perturbations in the energy estimates. Here we give a high level description of the linear
stability analysis using this new framework of analysis. We use the 2D Boussinesq equations
as an example. The same analysis also applies to 3D Euler with minor modifications. More
discussions and motivation will be provided in Section 2. Let ̄
ω
,
̄
θ
be an approximate steady
state. We denote
W
= (
ω,θ
x
,θ
y
) and decompose
W
=
W
+
̃
W
with
W
= ( ̄
ω,
̄
θ
x
,
̄
θ
y
). We further
denote by
L
the linearized operator around
W
that governs the perturbation
̃
W
in the dynamic
rescaling formulation (see Section 2.5),
(1.7)
̃
W
t
=
L
(
̃
W
)
,
where we have neglected the contributions from the third part, the nonlinear terms and the
residual error in the above linearized equation. We note that the coefficients of
L
depend on the
approximate steady state
W
. We decompose the linearized operator
L
into the main part
L
0
plus a finite rank perturbation
K
, i.e
(1.8)
L
=
L
0
+
K
.
The leading order operator
L
0
is constructed in such way that we can obtain sharp stability
estimates and extract damping effect by using sharp functional inequalities.
We first perform energy estimates to extract the damping effect from the leading order oper-
ator
L
0
in the weighted
L
∞
norm. By choosing appropriate singular weights, we can establish
stability associated with
L
0
using the weighted
L
∞
norm with a weak coupling to the weighted
C
1
/
2
norm. Obtaining sharp energy estimates in the weighted
C
1
/
2
norm is more challenging.
To obtain a sharp bound in the weighted
C
1
/
2
norm, we exploit the symmetry properties of
the kernels and apply
optimal transport
along the horizontal and the vertical directions sep-
arately to obtain relatively sharp
C
1
/
2
estimates. Since the contribution from the horizontal
direction dominates that from the vertical direction, our strategy gives a relatively tight bound
in the weighted
C
1
/
2
norm.
To estimate the perturbed operator
K
, we perform space-time estimates with computer as-
sistance. We use the following toy model to illustrate the main ideas by considering
K
as a
4
JIAJIE CHEN AND THOMAS Y. HOU
rank-one operator
K
(
̃
W
) =
a
(
x
)
P
(
̃
W
) for some nonlocal bounded linear operator
P
satisfying
the properties (i)
P
(
̃
W
) is constant in space and is only a function of time; (ii)
‖
P
(
̃
W
)
‖≤
c
‖
̃
W
‖
.
Given initial data
̃
W
0
, we decompose (1.7) as follows
(1.9)
∂
t
̃
W
1
(
t
) =
L
0
̃
W
1
,
̃
W
1
(0) =
̃
W
0
,
∂
t
̃
W
2
(
t
) =
L
̃
W
2
+
a
(
x
)
P
(
̃
W
1
(
t
))
,
̃
W
2
(0) = 0
.
It is easy to see that
̃
W
=
̃
W
1
+
̃
W
2
solves (1.7) with initial data
̃
W
0
since
L
=
L
0
+
a
(
x
)
P
by
our assumption. By construction, the leading operator
L
0
has the desired structure that enables
us to obtain sharp stability estimates. The second part
̃
W
2
is driven by the rank-one forcing
term
a
(
x
)
P
(
̃
W
1
(
t
)). Using Duhamel’s principle, the fact that
P
(
̃
W
1
(
t
)) is constant in space, we
yield
(1.10)
̃
W
2
(
t
) =
∫
t
0
P
(
̃
W
1
(
s
))
e
L
(
t
−
s
)
a
(
x
)
ds.
Since
̃
W
1
(
t
) =
e
L
0
(
t
)
̃
W
0
decays in
L
∞
(
φ
) (
φ
is a singular weight), we can control
P
(
̃
W
1
(
s
)).
If
e
L
(
t
)
a
(
x
) decays in
L
∞
(
φ
) for large
t
, we can show that
̃
W
2
(
t
) also decays in
L
∞
(
φ
) and
establish stability estimate of
̃
W
2
.
A crucial idea in the estimate of
K
is that we bridge the energy estimates and numerical
PDEs via an approximate solution in space and time. To see this, we note that
e
L
(
t
)
a
(
x
) is
equivalent to solving the linear evolution equation
v
t
=
L
(
v
) with initial data
v
0
=
a
(
x
). Due
to the rapid decay of the linearized equation, we only need to solve this initial value problem
using a numerical scheme up to a modest time
T
1
and pad the solution beyond
T
1
to be zero.
The residual error of the numerical PDE in space and time is treated as a small perturbation in
the energy estimates. We then construct an approximate solution by interpolating the solution
in time by a cubic polynomial. The stability property of
̃
W
1
allows us to control the numerical
error in computing
e
L
(
t
)
a
(
x
) and obtain sharp stability estimates for
̃
W
2
.
We remark that we only need to modify the linearized operator by a finite rank operator in
a small sector near the boundary where we have the smallest amount of damping. We also add
a few more terms to approximate the regular part of the solution in the bulk of the domain.
The amount of damping improves significantly as we move away from the boundary and in the
far field. Moreover, due to the decay of the far field solution
W
, we only need to construct a
finite rank operator
K
in a finite size domain within this narrow sector. The rank of
K
that
we use is less than 50. To estimate the contributions from the finite rank operator, we need to
solve the linear PDE in space-time with a number of initial data, which can be implemented in
full parallel. We have not explored the full potential of our method. If we try to perturb the
linearized operator by a rank-
N
operator, we just need to solve the linearized equation with
N
initial data, which can be implemented very efficiently using a powerful parallel cluster.
Our stability analysis uses some quantities involving the approximate steady state. The grid
values of these quantities are available with rigorous error bounds using interval arithmetic. Since
we use piecewise polynomials to approximate the steady state, we can obtain rigorous bounds
for the high order derivatives of the steady state. Such bounds in turn provide rigorous bounds
for lower order derivatives, the pointwise values and various integrals involving the approximate
steady state by using standard numerical analysis. The verification involves several case studies.
We will present all the essential tools for verification in the Appendix and also in later sections,
but will defer the actual verification steps of various quantities to the supplementary material.
To pass from the 2D Boussinesq equations to the 3D axisymmetric Euler equations, we follow
the same ideas presented in our previous work [12] by controlling the support of the vorticity
to be in a small region that is close to the boundary and does not intersect the symmetry axis.
The asymptotic scaling properties of the Biot-Savart kernels are exactly the same as those of
the Biot-Savart kernels for the 2D Boussinesq equations up to some asymptotically small terms
after making appropriate changes of variables. We will provide some additional estimates to
control these asymptotically small terms and prove the blowup of the 3D Euler equations.
STABLE BLOWUP OF 2D BOUSSINESQ EQUATIONS
5
The stability analysis presented in this paper can be substantially simplified if we have access
to a very powerful parallel cluster. If the rank of
K
is sufficiently large and the residual errors can
be made arbitrarily small, we can treat the third part of the linearized operator, the nonlinear
terms and the residual errors as small perturbations. Thus, we just need to focus on establishing
linear stability of
L
0
. Once we establish linear stability, nonlinear stability is relatively easy to
establish since the residual errors are very small. The
a priori
estimates of the perturbations in
the weighted
L
∞
and
C
1
/
2
norms make it easier to establish energy estimates in a higher order
norm, e.g. weighted
C
3
/
2
norm. As in [12–14, 27], by differentiating the perturbed equation for
̃
W
in time and performing
C
3
/
2
estimates, one can further prove the existence of the asymptot-
ically self-similar blowup profile. This provides a very powerful tool to give a constructive proof
of stable finite time blowup for a large class of nonlinear PDEs with smooth initial data.
1.3.
Comparison between our method of analysis and the topological argument.
Our
method of analysis shares some similarity with the recently developed blowup analysis using a
topological argument, see e.g. [58, 62–64]. In the topological argument, one also constructs a
compact perturbation operator
K
to the linearized operator
L
. After subtracting the compact
perturbation operator from the linearized operator, one can establish linear stability of the lead-
ing operator
L
0
in some Hilbert space. The compact perturbation operator can be approximated
by a finite rank operator. Typically, this compact perturbation operator contains a number of
unstable directions. A topological argument is used to show that there exists a stable trajectory
leading to a finite time blowup, which avoids a finite number of unstable directions. This method
has been used to prove unstable blowup of several nonlinear PDEs with great success.
The main difference between our method of analysis and the topological argument is in the
way we estimate the finite rank operator
K
. First of all, in our framework, we do not require
the energy space to be a Hilbert space. The main innovation of our approach is that we develop
a constructive method of analysis to establish stability of the finite rank operator by solving a
finite number of decoupled linear PDEs in space-time with rigorous error control. This allows
to establish linear stability of the original linearized operator
L
. In comparison, a typical
topological argument may only allow one to establish stability of the leading order operator
L
0
at the expenses of creating finitely many unstable directions induced by the finite rank operator
K
. Such method is ideal if we expect to have unstable blowup. However, if we expect to have
stable blowup as in the Hou-Luo blowup scenario, proving blowup with finitely many unstable
directions for the Hou-Luo blowup scenario using a topological argument is not satisfactory.
1.4.
Review of literature.
There has been a lot of effort in studying 3D Euler singularities
using various simplified models. Several 1D models, including the Constantin-Lax-Majda (CLM)
model [19], the De Gregorio (DG) model [24,25], the gCLM model [68] and the Hou-Li model [45],
have been introduced to study the effect of advection and vortex stretching in the 3D Euler
equations. Singularity formation from smooth initial data has been established for the CLM
model in [19], for the DG model in [14], and for the gCLM model with various parameters
in [5, 8, 10, 14, 29, 31]. In the viscous case, singularity formation of the gCLM model with some
parameters has been established in [8,70]. In [15], the authors proved the blowup of the Hou-Luo
model proposed in [55]. In [13], Chen-Hou-Huang proved the asymptotically self-similar blowup
of the Hou-Luo model by extending the method of analysis established for the finite time blowup
of the De Gregorio model by the same authors in [14].
In [16,38,39,50], the authors proposed several simplified models to study the Hou-Luo blowup
scenario [54, 55] and established finite time blowup of these models. In these works, the velocity
is determined by a simplified Biot-Savart law in a form similar to the key lemma in the seminal
work of Kiselev-Sverak [48]. In [30, 32], Elgindi and Jeong proved finite time blowup for the 2D
Boussinesq and 3D axisymmetric Euler equations in a domain with a corner using
̊
C
0
,α
data.
There has been some interesting recent results on the potential instability of the Euler blowup
solutions, see [51, 72]. In a recent paper [11], we showed that the blowup solutions of the 2D
Boussinesq and 3D Euler equations with
C
1
α
velocity considered in [12, 27] are also unstable
using the notion of stability introduced in [51, 72]. These two seemingly contradictory results
reflect the difference of the two approaches in studying the stability of 3D Euler blowup solutions.
6
JIAJIE CHEN AND THOMAS Y. HOU
The linear stability analysis in [51,72] is performed by directly linearizing the 3D Euler equations
around a blowup solution in the original variables. It does not take into account the changes
in the blowup time and the blowup exponent due to the change of the initial condition. Such
information has been used in establishing the nonlinear stability of the blowup profile in [12, 27].
In [42,43], Hou and Huang reported a potential two-scale traveling singularity for the 3D Euler
equations in the interior domain. Inspired by the work reported in [42,43], Hou discovered a new
class of potential Euler singularity at the origin whose scaling properties are compatible with
those of the Navier-Stokes equations [41]. Moreover, the Navier-Stokes equations with the same
smooth initial data develop potentially singular solutions at the origin with maximum vorticity
increasing by a factor of 10
7
[40]. Various blowup criteria have been applied to provide strong
support for this potentially singular behavior of the Navier-Stokes equations.
The rest of the paper is organized as follows. Sections 2 – 5 will be devoted to the blowup
analysis for the 2D Boussinesq equations and Section 6 will be devoted to the blowup analysis
for 3D Euler equations. In Section 2, we provide detailed discussions and some key ingredients in
establishing linear stability of an approximate profile using various simplified models. In Section
3, we show how to obtain sharp H ̈older estimates using optimal transport. Section 4 is devoted
to energy estimates and Section 5 is devoted to the construction of an approximate self-similar
profile using the dynamic rescaling formulation. In Section 7, we discuss the construction of the
approximate space-time solution to the linearized operator
L
. Finally we show how to estimate
the
L
∞
and H ̈older norms of the velocity in the regular case in Section 8. Some technical
estimates and derivations are deferred to the Appendix.
2.
Linear stability analysis and the main ideas
In this section, we will outline the main ingredients in our stability analysis. We will mainly
focus on the 2D Boussinesq equations since the analysis for the 3D Euler is very similar with
minor modifications, see Section 6. As in [12–14], we will use the dynamic rescaling formulation
for the 2D Boussinesq equations in an essential way. The most essential part of our analysis lies
in the linear stability. We need to use a number of techniques to extract the damping effect from
the linearized operator around the approximate steady state of the dynamic rescaling equations
and obtain sharp estimates of various nonlocal terms. Since the damping coefficients we obtain
are relatively small, we need to construct an approximate steady state with a very small residual
error. This is extremely challenging since the solution is supported on the upper half plane with
a slowly decaying tail in the far field. We use analytic estimates and numerical analysis with
rigorous error control to verify that the residual error is small in the energy norm. See more
detailed discussions in Section 5 and the supplementary material.
Passing from linear stability to nonlinear stability is relatively easier since the perturbation is
quite small due to the small residual error. Yet we need to verify various inequalities involving
the approximate steady state using the interval arithmetic [36,67,69] and numerical analysis with
computer assistance. The most essential part of the linear stability analysis can be established
based on the grid point values of the approximate steady state. The reader who is not interested
in the rigorous verification can skip the verification process in the supplementary material.
2.1.
Notations and operators.
The upper bar notation is reserved for the approximate steady
state, e.g. ̄
ω,
̄
θ
.
We introduce the notations for the nonlinear terms
(2.1)
N
1
=
−
u
·∇
ω
+
c
ω
ω,
N
2
=
−
u
·∇
η
−
u
x
η
−
v
x
ξ
+2
c
ω
η,
N
3
=
−
u
·∇
ξ
−
u
y
η
−
v
y
ξ
+2
c
ω
ξ.
Without specification,
N
i
depends on (
ω,η,ξ
). Given the approximate steady state ̄
ω,
̄
θ,
̄
c
l
,
̄
c
ω
,
we denote by
F
i
and
̄
F
ω
,
̄
F
θ
the residual error
(2.2)
̄
F
ω
=
−
( ̄
c
l
x
+
̄
u
)
·∇
̄
ω
+
̄
θ
x
+ ̄
c
ω
̄
ω,
̄
F
θ
=
−
( ̄
c
l
x
+
̄
u
)
·∇
̄
θ
+ ̄
c
θ
̄
θ,
F
1
,
̄
F
ω
,
F
2
,
∂
x
̄
F
θ
,
F
3
,
∂
y
̄
F
θ
.
STABLE BLOWUP OF 2D BOUSSINESQ EQUATIONS
7
Notations.
Denote by
C
α
x
,C
α
y
the partial H ̈older seminorms
(2.3)
[
ω
]
C
α
x
(
D
)
,
sup
x,z
∈
D,x
2
=
z
2
|
ω
(
x
)
−
ω
(
z
)
|
|
x
1
−
z
1
|
α
,
[
ω
]
C
α
y
(
D
)
,
sup
x,z
∈
D,x
1
=
z
1
|
ω
(
x
)
−
ω
(
z
)
|
|
x
2
−
z
2
|
α
.
Given a weight
g
(
h
) :
R
2
→
R
+
that is
−
α
-homogeneous , i.e.
g
(
λh
1
,λh
2
) =
λ
−
α
g
(
h
), e.g.,
g
(
h
) =
|
h
|
−
α
, we define the weighted H ̈older seminorm
(2.4)
||
ω
||
C
α
g
(
D
)
= sup
x,z
∈
D
|
(
ω
(
x
)
−
ω
(
z
))
g
(
x
−
z
)
|
.
We will mostly use
D
=
R
2
++
. In this case, we drop
D
to simplify the notations.
We define
〈
f,g
〉
the inner product in
R
2
++
(2.5)
〈
f,g
〉
=
∫
R
2
++
f
(
x
)
g
(
x
)
dx.
2.2.
Dynamic rescaling formulation.
Following [12–14], we consider the dynamic rescaling
formulation of the 2D Boussinesq equations. Let
ω
(
x,t
)
,θ
(
x,t
)
,
u
(
x,t
) be the solutions of (1.4)-
(1.6). Then it is easy to show that
(2.6)
̃
ω
(
x,τ
) =
C
ω
(
τ
)
ω
(
C
l
(
τ
)
x,t
(
τ
))
,
̃
θ
(
x,τ
) =
C
θ
(
τ
)
θ
(
C
l
(
τ
)
x,t
(
τ
))
,
̃
u
(
x,τ
) =
C
ω
(
τ
)
C
l
(
τ
)
−
1
u
(
C
l
(
τ
)
x,t
(
τ
))
,
are the solutions to the dynamic rescaling equations
(2.7)
̃
ω
τ
(
x,τ
) + (
c
l
(
τ
)
x
+
̃
u
)
·∇
̃
ω
=
c
ω
(
τ
) ̃
ω
+
̃
θ
x
,
̃
θ
τ
(
x,τ
) + (
c
l
(
τ
)
x
+
̃
u
)
·∇
̃
θ
=
c
θ
̃
θ,
where
̃
u
= ( ̃
u,
̃
v
)
T
=
∇
⊥
(
−
∆)
−
1
̃
ω
,
x
= (
x,y
)
T
,
(2.8)
C
ω
(
τ
) = exp
(
∫
τ
0
c
ω
(
s
)
dτ
)
, C
l
(
τ
) = exp
(
∫
τ
0
−
c
l
(
s
)
ds
)
, C
θ
= exp
(
∫
τ
0
c
θ
(
s
)
dτ
)
,
t
(
τ
) =
∫
τ
0
C
ω
(
τ
)
dτ
and the rescaling parameters
c
l
(
τ
)
,c
θ
(
τ
)
,c
ω
(
τ
) satisfy [12]
(2.9)
c
θ
(
τ
) =
c
l
(
τ
) + 2
c
ω
(
τ
)
.
We have the freedom to choose the time-dependent scaling parameters
c
l
(
τ
) and
c
ω
(
τ
) accord-
ing to some normalization conditions. These two free scaling parameters are related to the fact
that Boussinesq equations have scaling-invariant property with two parameters. The 3D Euler
equations enjoy the same property. See [12]. After we determine the normalization conditions
for
c
l
(
τ
) and
c
ω
(
τ
), the dynamic rescaling equation is completely determined and the solution
of the dynamic rescaling equation is equivalent to that of the original equation using the scaling
relationship described in (2.6)-(2.8), as long as
c
l
(
τ
) and
c
ω
(
τ
) remain finite.
We remark that the dynamic rescaling formulation was introduced in [52, 60] to study the
self-similar blowup of the nonlinear Schr ̈odinger equations. This formulation is also called the
modulation technique in the literature and has been developed by Merle, Raphael, Martel, Zaag
and others. It has been a very effective tool to analyze the formation of singularities for many
problems like the nonlinear Schr ̈odinger equation [47,61], compressible Euler equations [2,3], the
nonlinear wave equation [66], the nonlinear heat equation [65], the generalized KdV equation [57],
and other dispersive problems. Recently, this method has been applied to study singularity
formation in incompressible fluids [12, 27] and related models [8–10, 14].
To simplify our presentation, we still use
t
to denote the rescaled time in (2.7) and simplify
̃
ω,
̃
θ
as
ω,θ
(2.10)
ω
t
+ (
c
l
x
+
u
)
·∇
ω
=
θ
x
+
c
ω
ω,
θ
t
+ (
c
l
x
+
u
)
·∇
θ
=
c
θ
θ.
Following [13], we impose the following normalization conditions on
c
ω
,c
l
(2.11)
c
l
= 2
θ
xx
(0)
ω
x
(0)
, c
ω
=
1
2
c
l
+
u
x
(0)
, c
θ
=
c
l
+ 2
c
ω
.
8
JIAJIE CHEN AND THOMAS Y. HOU
For smooth data, these two normalization conditions play the role of enforcing
(2.12)
θ
xx
(
t,
0) =
θ
xx
(0
,
0)
, ω
x
(
t,
0) =
ω
x
(0
,
0)
for all time. In fact, we can derive the ODEs of
θ
xx
(
t,
0) and
ω
x
(
t,
0)
d
dt
ω
x
(
t,
0) = (
c
ω
−
c
l
−
u
x
(0))
ω
x
(
t,
0) +
θ
xx
(
t,
0)
,
d
dt
θ
xx
(
t,
0) = (
c
θ
−
2(
c
l
+
u
x
(0)))
θ
xx
(
t,
0)
,
where we use
v
|
y
=0
= 0
,v
x
(
t,
0) = 0. Under the conditions (2.11), the right hand sides vanish.
2.3.
Main Result.
In this section, we state our main result for the 2D Boussinesq equations.
We first introduce some notations and define our energy. Let
ψ
i
,φ
i
,ψ
i,g
,g
i
be the singular
weights defined in (C.1), (C.2), (C.3), and
μ
1
,μ
2
,τ
1
,τ
2
,μ
4
be the parameters chosen in (C.4).
We define the energy
E
on three variables
f
1
,f
2
,f
3
as follows
(2.13)
P
1
= max
1
≤
i
≤
3
||
f
i
φ
i
||
∞
, P
2
=
τ
−
1
1
max(
||
f
1
ψ
1
||
C
1
/
2
g
1
,μ
1
||
f
2
ψ
2
||
C
1
/
2
g
2
,μ
2
||
f
3
ψ
3
||
C
1
/
2
g
3
)
P
3
=
τ
2
max(
μ
2
||
f
1
φ
g
1
||
∞
,
||
f
2
φ
g,
2
||
∞
,
||
f
3
φ
g,
3
||
∞
)
,
P
4
= max(
1
65
|
c
ω
(
ω
)
|
,
1
10
|
f
2
,xy
(0)
|
,
1
5
|
f
1
,xy
(0)
|
)
, E
= max(
P
1
,P
2
,P
3
,P
4
)
.
where
u
x
(
f
)(0) =
−
4
π
∫
R
++
2
y
1
y
2
|
y
|
4
f
(
y
)
dy
.
Theorem 3.
Let
(
̄
θ,
̄
ω,
̄
c
l
,
̄
c
ω
)
be the approximate self-similar profile constructed in Section 5
and
E
∗
= 5
·
10
−
6
. For even initial data
θ
0
and odd
ω
0
of
(2.10)
with a small perturbation to
(
̄
θ,
̄
ω
)
with
E
(
ω
0
−
̄
ω,θ
0
,x
−
̄
θ
x
,
̄
θ
0
,y
−
̄
θ
y
)
< E
∗
, we have
(2.14)
||
ω
−
̄
ω
||
L
∞
,
||
θ
x
−
̄
θ
x
||
L
∞
,
||
θ
y
−
̄
θ
y
||
∞
<
10
3
E
∗
,
|
u
x
(
t,
0)
−
̄
u
x
(0)
|
,
|
̄
c
ω
−
c
ω
|
<
100
E
∗
for all time. In particular, we can choose smooth initial data
ω
0
,θ
0
∈
C
∞
c
in this class with
finite energy
||
u
0
||
L
2
<
+
∞
such that the solution to the physical equations
(1.4)
-
(1.6)
with these
initial data blows up in finite time
T
.
Remark
2.1
.
In our analysis, we decompose the nonlinear operator that governs the evolution
of the perturbation into two parts. Correspondingly, we decompose the perturbation
̃
W
=
(
ω
−
̄
ω,θ
x
−
̄
θ
x
,θ
y
−
̄
θ
y
) into two components,
̃
W
=
̃
W
1
+
̃
W
2
, see Section 2.10.3 for the precise
decomposition. At the linear level, the evolution of
̃
W
1
is governed by a leading order operator
L
0
that enjoys sharp energy estimates after we subtract a finite rank operator. The perturbation
̃
W
2
mainly captures the contributions from the finite rank perturbation operator. One advantage
of using such decomposition is that
̃
W
1
is only weakly coupled to
̃
W
2
through the nonlinear terms
and the residual errors, which are small. We can estimate the contributions of
̃
W
2
to the linear
stability by solving a finite number of linear PDEs in space and time and estimate the space-time
residual errors via energy estimates. The stability of the leading order operator
L
0
enables us to
control the residual errors arising from solving the linear PDEs in space and time rigorously. A
crucial step in proving our main theorem is to establish the estimate
E
4
(
̃
W
1
)
< E
∗
with energy
E
4
defined in (4.59). Note that
̃
W
1
=
̃
W
at
t
= 0 and the energy
E
4
agrees with
E
at
t
= 0.
See Sections 2.10, 4, and (4.8.5) for more discussions.
We will follow the framework in [12–14] to establish finite time blowup. It consists of the
following steps: (a) construct approximate steady state to (2.10) with small residual error in
suitable functional spaces; (b) perform linear and nonlinear stability analysis around the ap-
proximate steady state; (c) choose small perturbation in the energy space to obtain smooth
initial data with finite energy and obtain finite time blowup using a rescaling argument and the
stability results.
In the remaining of this section, we will outline some main ingredients in our blowup analysis
by using a number of simplified models to illustrate and motivate the main ideas behind our
method of analysis.
STABLE BLOWUP OF 2D BOUSSINESQ EQUATIONS
9
2.4.
Basic properties of the approximate steady state.
Following the ideas in [13, 14], we
construct the approximate steady state ( ̄
ω,
̄
θ,
̄
c
ω
,
̄
c
l
) of the dynamic rescaling equations (2.10),
(2.11) by solving them numerically for a long enough time. In Figure 1, we plot the approximate
steady state ̄
ω,
̄
θ
x
. We plot the variable
̄
θ
x
rather than
̄
θ
since
̄
θ
grows in the far-field. Given
the approximate steady state, we construct the numerical stream function
̄
φ
N
by solving the
Poisson equations. Then we can derive the residual (2.2) up to the error in solving the Poisson
equations. In Figure 2, we plot the piecewise rigorous bound of the weighted
L
∞
(
φ
1
) norm
of
̄
F
1
and the
L
∞
(
φ
2
) norm of
̄
F
2
. We remark that
φ
1
,φ
2
are very singular near
x
= 0 with
leading order
|
x
|
−
2
.
9
,
0
.
385
|
x
|
−
3
. This is why the weighted
L
∞
(
φ
1
) norm of
̄
F
1
and the
L
∞
(
φ
2
)
norm of
̄
F
2
are relatively large near the origin.
We observe that unweighted errors of
̄
F
1
,
̄
F
2
are very small near the origin, less than 2
·
10
−
12
since we use a uniform fine grid near the origin. In Figure 3, we plot the piecewise rigorous
bound of the unweighted
L
∞
norm of
̄
F
1
and the
L
∞
norm of
̄
F
2
in a local domain near the
origin. Note that the unweighted
L
∞
norm of
̄
F
1
and the
L
∞
norm of
̄
F
2
increase as
|
x
|
increases. On the other hand, since the singular weights have a very mild growth rate in the
far field of order
O
(
|
x
|
1
/
16
), the contributions from these weighted errors to the energy
E
3
scale
like 0
.
165
∗|
x
|
1
/
16
‖
̄
F
i
‖
L
∞
(
i
= 1
,
2). Thus the weighted
L
∞
norm of
̄
F
1
and the
L
∞
norm of
̄
F
2
only amplify the unweighted norms very mildly in the far field. Since the damping effect
is stronger in the far field, we are able to close the energy estimates with the residual errors
that we obtain. We defer the details of numerical computation to Section 5. Here, we list some
important properties of the approximate steady state.
Figure 1.
Approximate steady state in the near-field. Left figure: profile ̄
ω
;
right figure:
̄
θ
x
.
Exponents.
The exponents and the velocity near the origin satisfy
(2.15)
̄
c
l
≈
3
.
00649898
,
̄
c
ω
≈−
1
.
02942516
,
̄
u
x
(0)
≈−
2
.
532674
,
̄
v
x
(0)
≈
0
.
We remark that the ratio ̄
c
l
/
̄
c
ω
≈ −
2
.
9205600 is very close to the one reported by Hou-Luo
[54, 55].
Regularity and representation.
The variables ̄
ω,
̄
ψ
are odd in
x
and
̄
θ
is even in
x
. Denote
by
φ
the stream function. One should not confuse the stream function
φ
with singular weights
φ
1
,
φ
2
, etc. The approximate steady state ( ̄
ω,
̄
θ,
̄
φ
) is represented by piecewise fifth order
polynomials ̄
ω
2
,
̄
θ
2
,
̄
φ
2
supported in [0
,D
1
]
2
with
D
1
≈
10
15
, and semi-analytic parts ̄
ω
1
,
̄
θ
1
,
̄
φ
1
that capture the far-field behavior of the solutions
̄
ω
= ̄
ω
1
+ ̄
ω
2
,
̄
θ
=
̄
θ
1
+
̄
θ
2
,
̄
φ
=
̄
φ
1
+
̄
φ
2
.
See (5.2). In particular, we have ̄
ω,
̄
θ,
̄
φ
∈
C
4
,
1
. The solution enjoys the decay rate
̄
ω
∼
r
α
,
̄
θ
∼
r
1+2
α
, α
≈
̄
c
ω
/
̄
c
l
.
10
JIAJIE CHEN AND THOMAS Y. HOU
Figure 2.
Weighted residual errors of the approximate steady state. Left
figure: piecewise rigorous
L
∞
(
φ
1
) bound of
F
1
in the
ω
equation. Right figure:
piecewise rigorous
L
∞
(
φ
2
) bound of
F
2
in the
θ
x
equation
.
Figure 3.
Unweighted residual errors of the approximate steady state. Left
figure: piecewise rigorous
L
∞
bound of
F
1
in the
ω
equation. Right figure:
piecewise rigorous
L
∞
bound of
F
2
in the
θ
x
equation.
Anisotropic.
The solutions
̄
θ
and ̄
ω
are anisotropic in the sense that the
y
-derivative of the
profile is much smaller than the
x
-derivative, especially in the near field:
(2.16)
|
̄
θ
y
|
< c
3
|
̄
θ
x
|
, c
3
≈
0
.
16
,
|
̄
ω
y
|
< c
4
|
̄
ω
x
|
, c
4
≈
0
.
23
,
for (
x,y
)
∈
[0
,
1]
2
. Note that the anisotropic property of the solutions has been observed for the
C
1
,α
singular solution [12].
The advection.
The advection in (2.10) satisfies the following important inequalities
̄
c
l
x
+ ̄
u
(
x,y
)
≥
c
1
x, c
1
≈
0
.
47
,
̄
c
l
y
+ ̄
v
(
x,y
)
≥
c
2
y, c
2
≈
3
,
for all
x,y
∈
R
2
++
. For
x,y
∈
R
2
++
near the origin, we have
̄
c
l
x
+ ̄
u
(
x,y
)
≈
0
.
47
x, c
l
y
+ ̄
v
(
x,y
)
≈
5
.
54
y.