Large time decay and growth for solutions of a viscous Boussinesq system

In this paper we analyze the decay and the growth for large time of weak and strong solutions to the three-dimensional viscous Boussinesq system. We show that generic solutions blow up as $t\to\infty$ in the sense that the energy and the $L^p$-norms of the velocity field grow to infinity for large time, for $1\le p<3$. In the case of strong solutions we provide sharp estimates both from above and from below and explicit asymptotic profiles. We also show that solutions arising from $(u_0,\theta_0)$ with zero-mean for the initial temperature $\theta_0$ have a special behavior as $|x|$ or $t$ tends to infinity: contrarily to the generic case, their energy dissipates to zero for large time.


Introduction
In this paper we address the problem of the heat transfer inside viscous incompressible flows in the whole space R 3 . Accordingly with the Boussinesq approximation, we neglect the variations of the density in the continuity equation and the local heat source due to the viscous dissipation. We rather take into account the variations of the temperature by putting an additional vertical buoyancy force term in the equation of the fluid motion.
x ∈ R 3 , t ∈ R + Here u : R 3 × R + → R 3 is the velocity field. The scalar fields p : R 3 × R + → R and θ : R 3 × R + → R denote respectively the pressure and the temperature of the fluid. Moreover, e 3 = (0, 0, 1), and β ∈ R is a physical constant. For the decay questions that we address in this paper, it will be important to have strictly positive viscosities in both equations: ν, κ > 0. By rescaling the unknowns, we can and do assume, without loss of generality, that ν = 1 and β = 1. To simplify the notation, from now on we take the thermal diffusion coefficient κ > 0 such that κ = 1.
Like for the Navier-Stokes equations, obtained as a particular case from (1.1) putting θ ≡ 0, weak solutions to (1.1) do exist, but their uniqueness is not known. The global existence of weak solutions, or strong solutions in the case of small data has been studied by several authors. See, e.g. [1], [11], [17], [18], [34]. Conditional regularity results for weak solutions (of Serrin type) can be found in [9]. The smoothness of solutions arising from large axisymmetric data is addressed in [2] and [24]. Further regularity issues on the solutions have been discussed also [20,15].
The goal of this paper is to study in which way the variations of the temperature affect the asymptotic behavior of the velocity field. We point out that several different models are known in the literature under the name of "viscous (or dissipative) Boussinesq system". The asymptotic behaviour of viscous Boussinesq systems of different nature have been recently addressed, e.g., in [3,12]. But the results therein cannot be compared with ours.
Only few works are devoted to the study of the large time behavior of solutions to (1.1). See [21,26]. These two papers deal with self-similarity issues and stability results for solutions in critical spaces (with respect to the scaling). On the other hand, we will be mainly concerned with instability results for the energy norm, or for other subcritical spaces, such as L p , with p < 3.
The above estimate for the temperature looks optimal, since the decay agrees with that of the heat kernel. On the other hand the optimality of the estimate for the velocity field is not so clear. For example, in the particular case θ 0 = 0, the system boils down to the Navier-Stokes equations and in this simpler case one can improve the bound for the velocity into u(t) 2 ≤ u 0 2 . In fact, u(t) 2 → 0 for large time by a result of Masuda [28]. Moreover, in the case of Navier-Stokes the decay of u(t) 2 agrees with the L 2 -decay of the solution of the heat equation. See [35,25,39] for a more precise statement.
The goal of the paper will be to show that the estimate of weak solutions u(t) 2 ≤ C(1 + t) 1/4 can be improved if and only if the initial temperature has zero mean. To achieve this, we will establish the validity of the corresponding lower bounds for a class of strong solutions.
In particular, this means that very nice data (say, data that are smooth, fast decaying and "small" in some strong norm) give rise to solutions that become large as t → ∞: our results imply the growth of the energy for strong solutions : (1.2) c(1 + t) 1/4 ≤ u(t) 2 ≤ C(1 + t) 1/4 , t > > 1.
The validity of the lower bound in (1.2) (namely, the condition c > 0), will be ensured whenever the initial temperature is sufficiently decaying but We feel that is important to point out here an erratum to the paper [23]. Unfortunately, the lower bound in (1.2) contredicts a result in [23,Theorem 2.3], where the authors claimed that u(t) 2 → 0 under too general assumptions, weaker than those leading to our growth estimate. The proof of their theorem (in particular, of inequalities (5.6) and (5.8) in [23]) can be fixed by putting different conditions on the data, including θ 0 = 0. This is essentially what we will do in part (b) of our Theorem 2.2 below. Similarity, the statement of Theorem 2.4 in [23] contredicts our lower bound (1.6) below (inequalities (5.18)- (5.20) in their proof do not look correct). This will be also corrected by our Theorem 2.2. We would like to give credit to the paper [23] (despite the above mentioned errata), because we got from there inspiration for our results of Sections 3 and 4.
Our main tool for establishing the the lower bound will be the derivation of exact pointwise asymptotic profiles of solutions in the parabolic region |x| > > √ t. This will require a careful choice of several function spaces in order to obtain as much information as possible on the pointwise behavior of the velocity and the temperature. A similar method has been applied before by the first author in [6] in the case of the Navier-Stokes equations, although the relevant estimates were performed there in a different functional setting.
Even though several other methods developped for Navier-Stokes could be effective for obtaining estimates from below (see, e.g., [13,22,30]), our analysis has the advantage of putting in evidence some features that are specific of the Boussinesq system: in particular, the different behavior of the flow when |x 3 | → ∞ or when x 2 1 + x 2 2 → ∞, due to the verticality of the bouyancy forcing term θe 3 (see Theorem 2.6 below). Moreover, the analysis of solutions in the region |x| > > √ t and our use of weighted spaces completely explains the phenomenon of the energy growth: the variations of the temperature push the fluid particles in the far field; even though in any bounded region the fluid particles slow down as t → ∞ (this effect is measured e.g. by the decay of the L ∞ -norm established in Proposition 2.5), large portions of fluid globally carry an increasing energy during the evolution. Our result thus illustrates the physical limitations of the Boussinesq approximation, at least for the study of heat convection inside fluids filling domains where Poincaré's inequality is not available, such as the whole space.
In fact, our method applies also to weighted L p -spaces, so let us introduce the weighted norm Then we will show that strong solutions starting from suitably small and well decaying data satisfy, for t > 0 large enough, As before in (1.2), these lower bounds hold true with a constant c > 0 as soon as the initial temperature has non-zero mean. Notice that in this case the L p -norms asymptotically blow up for large time if and only if p < 3. The above restriction r + 3 p < 3 on the parameters is optimal. Indeed, under the same conditions yielding to (1.3) we will also prove that when r ≥ 0, 1 ≤ p < ∞ and r + 3 p ≥ 3, then one has u(t) L p r = ∞, for all t > 0. The fact that the L p r -norm becomes infinite instanteneously for this range of the parameters is related to that the velocity field immediately spatially spreads out and cannot decay faster than |x| −3 as |x| → ∞ for t > 0, and this even if u 0 ∈ C ∞ 0 (R 3 ).
As we observed before, the lower bound in (1.3) brakes down when θ 0 = 0. In such case, our decay estimates can be improved. We will establish the upper bound for weak solutions by the Fourier splitting method. This method was first introduced in [35]. We will need no smallness assumption on u 0 to prove (1.4). We have to put however a smallness condition of the form We do not know if it is possible to get rid of (1.5) to establish (1.4). Such smallness condition however looks natural as it respects the natural scaling invariance of the system (1.1).
As before, the decay estimate (1.4) is optimal for generic solutions satisfying θ 0 = 0. Indeed, when we start from localized and small velocity, we can establish the upper-lower bounds for strong solutions Similarily as before, the validity of the lower bound (the condition c ′ > 0) now requires a non vanishing condition in the first moments of θ.
1.1. Notations. We denote by C ∞ 0 the space of smooth functions with compact support. The L p -will be denoted by · p and in the case p = 2 we will simply write · . Moreover, σ denotes the completion of V under the norm · p , and V the closure of V in H 1 0 We will adopt the following convention for the Fourier transform of integrable functions: Ff (ξ) = f (ξ) = f (x)e −2πix·ξ dx. Here and throughout the paper all integral without integration limits are over the whole R 3 .
The notations L p,q and L p r have a different meaning. For 1 < p < ∞ and 1 ≤ q ≤ ∞, L p,q denotes the classical Lorentz space. For the definition and the basic inequalities concerning Lorentz spaces (namely the generalisation of the straightforward L p -L q Hölder and Young convolution inequalities) the reader can refer to [27,Chapter 2]. Here we just recall that L p,p agrees with the usual Lebesgue space L p , and L p,∞ agrees with the weak Lebesgue space L p w (or Marcinkiewicz space) 1 p is equivalent to the natural norm on L p,∞ , for 1 < p < ∞. Here as it is usual we defined where λ denotes the Lebesgue measure. On the other hand, for 1 ≤ p ≤ ∞ and 0 ≤ r ≤ ∞, L p r is the weighted Lebesgue space, consisting of the functions f such that (1+ |x|) q f ∈ L p . Notice that the bold subscript σ introduced above is not a real parameter: in the notation L p σ , this subscript simply stands for "solenoidal".
To denote general constants we use C which may change from line to line. In certain cases we will write C(α) to emphasize the constants dependence on α.
We denote by e t∆ the heat semigroup. Thus, We denote by E(x) the fundamental solution of −∆ in R 3 . The partial derivatives of E are denoted by E x j , E x j ,x k , etc.
1.2. Organization of the paper. All the main results are stated without proof in Section 2. The statement of our theorems are splitted into two parts: part (a) is devoted to the properties of solutions in the general case (where the integral θ 0 is not necessarily zero). In Part (b) of our theorems we are concerned with the special case θ 0 = 0.
The rest of the paper is organized as follows. Sections 3-4 are devoted to the proof of our results on weak solutions and Sections 5-6 to strong solutions. In Section 7 we collect a few technical remarks.

Statement of the main results
2.1. Results on weak solutions. We will consider weak solutions to the viscous Boussinesq equations and establish existence and natural decay estimates in L p spaces, 1 ≤ p < ∞, for the temperature, together with bounds for the growth of the velocities. Such estimates rely on the fact that in both equations of the system (1.1) we have a diffusion term. They complete those of obtained in [16] where there was no temperature diffusion.
We start recalling the basic existence result of weak solutions to the system (1.1). See, e.g., [9,16].
σ . There exists a weak solution (θ, u) of the Boussinesq system (1.1), continuous from R + to L 2 with the weak topology, with data (u 0 , θ 0 ) such that, for any T > 0, Such solution satisfies, for all t ∈ [0, T ], the energy inequalities for all t ≥ 0, and some absolute constant C > 0.
One can improve the growth estimate on the velocity as soon as θ 0 belongs to some L p space, with p < 2. For simplicity we will consider only the case θ 0 ∈ L 1 ∩ L 2 .
Moreover, it is natural to ask under which supplementary conditions on the initial data one can insure that the energy of the fluid u(t) 2 remains uniformly bounded. Theorem 2.2 provides an answer. Beside a smallness assumption on θ 0 1 , we need to assume that θ 0 = 0.
In addition, it is possible to prove that u(t) not only remains uniformly bounded, but actually decays at infinity (without any rate). Explicit decay rates for u(t) can also be prescribed provided the linear part e ∆ u 0 decays at the appropriate rate.
Under the additional condition θ 0 ∈ L 1 the estimates on the weak solution constructed in Proposition 2.1 can be improved into

3)
Moreover, if θ 0 ∈ L 1 ∩ L p , for some 1 ≤ p < ∞, then (b) (The θ 0 = 0 case) In this part we additionally assume θ 0 ∈ L 1 1 and θ 0 = 0. Then there exists an absolute constant ε 0 > 0 such that if then the weak solution of the Boussinesq system (1.1) constructed in Proposition 2.1 satisfies, for some constant C > 0 and all t ∈ R + , and Moreover under the additional condition u 0 ∈ L 3/2 ∩ L 2 σ , we have The proof of part (a) of Theorem 2.2 is straightforward. Part (b) is more subtle: its proof relies on an inductive argument: assuming that the approximate velocity u n−1 grows in L 2 at most like t 1/8 (which is actually better than provided by Proposition 2.1) we can prove that the same growth estimate remains valid for u n . A smallness condition on θ 0 is needed to insure that in the inequality u n (t) ≤ C(1 + t) 1/8 the constant can by taken independently on n. With this improved control on the growth of the velocity, applying the Fourier splitting method (introduced in [35]) and a few boot-strapping we can improve our estimates up to the rates given in (2.5) and (2.7). Remark 2.3. Solutions of the Boussinesq system (1.1) with energy decaying faster than t −1/2 might exist, but they are likely to be highly non-generic. Indeed, it seems difficult to construct such solutions assuming only that the data belong to suitable function spaces (possibly with small norms). The main obstruction is that one would need stringent cancellations properties on the data, which however turn out to be non-invariant under the Boussinesq flow. A possible way to obtain such fast decaying solutions would be to start with data satisfying some special rotational symmetries, as those described in [5].

2.2.
Results on strong solutions. The best way to prove the optimality of the estimates contained in Theorem 2.2 is to establish the corresponding lower bound estimates at least for a subclass of solutions. For the study of the estimates from below we will limit our considerations to a class of strong solutions. This is not a real restriction as lower bound estimates established for solutions emanating from well localized, smooth and small data are expected to remain valid in the larger class of weak solutions. Studying strong solutions has also the advantage of better putting in evidence some interesting properties specific of the Boussinesq system, such as the influence of the vertical buoyancy force on the pointwise behavior of the fluid in the far-field.
The existence of strong solutions to the system (1.1) will be insured by a fixed point theorem in function spaces invariant under the natural scaling of the equation. Thus, if u ∈ X , where X is a Banach space to be determined, we want to have, for all λ > 0, is the rescaled velocity. A suitable choice for norm of the space X , inspired by [10], is This choice for X is quite natural. Indeed, whenever |u 0 (x)| ≤ C|x| −1 , the linear evolution e t∆ u 0 belongs to X and this paves the way for the application of the fixed point theorem in such space. More precisely, we define X as the Banach space of all locally integrable divergence-free vector fields u such that u X < ∞, and continuous with respect to t in the following usual sense: u(t) → u(0) in the distributional sense as t → 0 and ess sup x∈R 3 In the same way, if θ belongs to a Banach space Y, we want to have is the rescaled temperature. We then define a Banach space Y of scalar functions through the norm and the natural continuity condition on the time variable as before.
The starting point of our analysis will be the following proposition providing a simple construction of mild solutions (u, θ) ∈ X ×Y. We will refer to them also as strong solutions. Indeed, one could prove that such solutions turn out to be smooth, as one could check by adapting to the system (1.1) classical regularity criteria for the Navier-Stokes equations like that of Serrin [37]. See the paper [9]. The smoothness of these solutions, however, plays no special role in our arguments.
Moreover, these conditions define u and θ uniquely.
Next Proposition shows it is possible to obtain better space-time decay estimates provided one starts with suitably decaying data.

Proposition 2.5.
(a) Let u 0 and θ 0 as in Proposition 2.4, and satisfying the additional decay estimates, for some 1 ≤ a < 3, b ≥ 3, and a constant C > 0, Then the solution constructed in Proposition 2.4 satisfies, for another constant C > 0 independent on x and t, (b) (the θ 0 = 0 case) Assume now 2 ≤ a < 4, a = 3, and b ≥ 4 and let u 0 and θ 0 satisfying the previous assumptions. If, in addition, then the decay of u and θ is improved as follows: Thus, in the following asymptotic expansions, ∇E x 3 and ∇E x h ,x 3 are vectors whose components are homogeneous functions of degree −3 and −4 respectively. In particular |∇E We are now in the position of stating our main results on strong solutions. The first theorem describes the asymptotic profiles of solutions in the parabolic region |x| > > √ t. Roughly, it states that all sufficiently decaying solutions (u, θ) of (1.1) behave in such region like a potential flow.
(a) Let a > 3 2 and b > 3. Let (u, θ) be a (mild) solution of (1.1) satisfying the decay estimates (2.13)-(2.14). Then the following profile for u holds: Let (u, θ) be a solution satisfying the decay condition (2.16). Then the following profiles for u j (j = 1, 2, 3) hold: where R is a lower order term for |x| > > √ t > > 1, namely The following remark should give a better understanding of the theorem.
(a) Let a > 3 2 , b > 3 and let (u, θ) be a solution as in Part (a) of Theorem 2.6. Then for all r, p such that there exists t 0 > 0 such that the solution satisfies the upper and lower estimates in the weighted-L p -norm we have for another suitable continuous function φ : R + → R such that φ(0) = 0 and φ(σ) > 0 for σ > 0.
Remark 2.9. When θ 0 = 0, we thus get by (2.23) the sharp large time behavior u(t) L p r ≃ t . The condition m = 0 is satisfied for generic solutions. It prevents θ to have oscillations at large times.

The mollified Boussinesq system and existence of weak solutions
The existence of weak solutions to the Boussinesq system is well known, see [9]. Their uniqueness, however, is an open problem. Moreover, we do not know if any weak solutions satisfy the energy inequality and the decay estimates stated in Proposition 2.1. For this reason, we now briefly outline another construction of weak solutions, which is well suited for obtaining all our estimates.
We begin by introducing a mollified Boussinesq system. As the construction below is a straightworward adaptation of that of Caffarelli, Kohn and Nirenberg, [8], we will be rather sketchy. For completeness we recall the definition of the "retarded mollifier" as given in [8]. Let ψ(x, t) ∈ C ∞ such that Consider, for n = 1, 2, . . . and δ = T /n, the mollified Cauchy problem The iteration scheme starts with u 0 = 0. Note that since div u = 0, we also have div (Ψ δ (u n )) = 0, for t ∈ R + . At each step n, one solves recursively n + 1 linear equations: first one solves the transport-diffusion equation (with smooth convective velocity) for the temperature; after θ n is computed, solving the second of (3.1) amounts to solving a linear equation on each strip R 3 × (mδ, (m + 1)δ), for m = 0, 1, . . . , n − 1.
For solutions to (3.1) we have the following existence and uniqueness result, For each n ∈ {1, 2, . . . } there exists a unique weak solution (θ n , u n , p n ) of the approximating equations with data (3.2) such that, for any T > 0, Moreover, for all t > 0, θ n and u n satisfy the energy inequalities as in (2.1) and (2.2).
In particular, the sequences θ n , u n and p n , n = 1, 2, . . . are bounded in their respective spaces.
Proof. This can be proved using the Faedo-Galerkin method. As the the argument is standard (see, e.g. [38, Theorem 1.1, Chapter III], or [8,Appendix]), we skip the details. We only prove the condition on the pressure since this is the only change that we have to make in [8].
Taking the divergence of the second equation in (3.1) we get p = p n 1 + p n 2 , where ∆p n and Thus, p n 1 ∈ L 5/3 (0, T ; L 5/3 ) uniformly with respect to n, by the energy inequality for u n , interpolation, and the Calderon-Zygmund theorem, as proved in [8]. On the other hand, belongs to the Lorentz space L 3/2,∞ (R 3 ). But θ n ∈ L ∞ (0, T ; L 2 ) uniformly with respect to n, hence Young convolution inequality in Lorentz spaces (see [27,Chapter 2]) yields p n 2 ∈ L ∞ (0, T ; L 6 ) uniformly with respect to n.
It follows from Proposition 3.1 that, extracting suitable subsequences, p n = p n 1 + p n 2 , where (p n 1 ) converges weakly in L 5/3 (0, T ; L 5/3 ) and p n 2 converges in L ∞ (0, T ; L 6 ) in the weak-* topology. Moreover, (u n ) is convergent with respect to the topologies listed in [8, p. 828-829]. On the other hand (θ n ) will be convergent with respect to the same topologies, because all estimates available for u n hold also for θ n .
No additional difficulty in the passage to the limit in the nonliner terms arises, in the equation of the temperature other those already existing for the Navier-Stokes equations. Hence, the distributional limit (θ, u, p) of a convergent subsequence of (θ n , u n , p n ) is a weak solutions of the Boussinesq system. This establishes Proposition 2.1.
We finish this section by establishing the natural L p -estimates for the approximating temperatures: σ and let (θ n , u n , p n ) be the solution of the mollified Boussinesq system (3.1) for some n ∈ {1, 2, . . .}. Let also 1 ≤ p < ∞. If θ 0 ∈ L 1 ∩ L p , then where A = A(p, θ 0 1 , θ 0 p ) and c > 0 is an absolute constant.
Proof. First notice that, for each n, θ n is the solution of a linear transport-diffusion equation with smooth and divergence-free velocity Ψ(u n−1 ). The L p decay estimates for these equations are well known. We reproduce the same proof as in [14,19]) expliciting better the constants, as we will need the expressions of such constants later on.
See [14, Corollary 2.6] for a nice proof of (3.4) that remains valid in the much more general case of trasport equations with (or without) fractional diffusion. We start with the case 2 ≤ p < ∞. Multiplying the equation for θ n by p|θ n | p−2 θ n and integrating we get By the Sobolev embedding theorem,Ḣ 1 ⊂ L 6 , hence The interpolation inequality yields Combining these two inequalities with the basic estimate θ n (t) 1 ≤ θ 0 1 , we obtain the differential inequality Integrating this we get
The case 1 ≤ p < 2 is deduced by interpolation.
Proof. We only have to estimate the L 2 norm of the velocity. We make use of the identity that can be justified exactly as for the mollified Navier-Stokes equations, see [8] and [35].
Multiplying the velocity equation in the Boussinesq system (3.1) by u n and integrating we get Dividing by u n (t) , d dt u n (t) ≤ θ n (t) .
Now we use the decay of θ n (t) obtained in Lemma 3.2. Integrating we obtain the second of (3.5).
4. Improved bounds for weak solution in the case θ 0 = 0 The estimates obtained in Lemma 3.3 can be considerably improved provided we additionally assume θ 0 = 0 and the moment condition θ 0 ∈ L 1 1 . First of all, from an elementary heat kernel estimate one easily checks that in this case where A 1 > 0 depends only on the data, through its L 2 -norm and the |x| |θ 0 (x)| dx integral.
We will now see that in this case the approximated temperature θ n (t) also decays at the faster rate (t + 1) −5/4 in the L 2 -norm. Using this new decay rate for θ n it is possible to show that the velocities are uniformly bounded in L 2 . Once we have a solution such that u n (t) remains bounded as t → ∞ one can go further and prove that u n actually decays at infinity in L 2 at some algebraic decay rate, depending on the decay of the linear evolution e t∆ u 0 . In view of the passage to the limit from the mollified system (3.1) to the Boussinesq system (1.1), all the estimates must be independent on n.  for all n ∈ N and t ∈ R + . Here A > 0 is some constant depending on the data u 0 and θ 0 , and independent on n, and t.
Proof. We denote by C a positive absolute constant, which may change from line to line. We also denote by A 1 , A 2 , . . . positive constants that depend only on the data. More precisely, A j = A j θ 0 , θ 0 L 1 1 , u 0 .
We make use of the Fourier splitting method introduced in [35]. The first step consists in multiplying the temperature equation by θ n and to integrate by parts. Using the Plancherel theorem in the energy inequality for θ n , we get Now split the integral on the right hand side into S ∪ S c , where (4.5) S = ξ : |ξ| ≤ k 2(t + 1)

1/2
, and k is a constant to be determined below. Noting that for ξ ∈ S c one has −|ξ| 2 ≤ − k 2(t+1) , it follows that Multiplying by (1 + t) k we obtain where S is as in (4.5).
Taking the Fourier transform in the equation for θ n in (3.1) we get Replacing this in (4.6) and applying the Plancherel theorem we get From where it follows, letting k = 7/2,  Recalling the estimate (4.1) for the linear evolution we get, for all n ∈ N, We now use the following inequality, deduced from estimate (3.5), Putting this inside (4.9) we obtain a new bound for θ n 2 , namely (4.11) Step 2: The inductive argument.
We will now prove by induction that, for all positive integer n we have u n−1 (t) ≤ u 0 + M t 1/8 (4.12) where M > 0 is some constant independent on n (but possibly dependent on the data θ 0 , u 0 ) to be determined. Notice that estimate (4.12) is actually better than what we have so far (compare with the second of (3.5)).
For n = 1 the inductive condition (4.12) is immediate since u 0 = 0. Let us now prove that u n ≤ u 0 + M t 1/8 assuming that (4.12) holds true.
We get from (4.11) and the induction assumption (4.12) This implies This last inequality will be used to estimate u n as follows. First recall that (4.14) d dt u n (t) ≤ θ n (t) .
After an integration in time we get, using (4.13), For the induction argument we need to prove u n ≤ u 0 + M t 1/8 . Hence we need that Choosing M large enough, for example, our condition then boils down to the inequality The validity this last inequality is insured by assumption (4.2). This concludes the induction argument and establishes the validity of the estimate (4.12) for all n.
Step 3: Uniform bound for the L 2 -norm of the velocities u n .
The result of Step 2 implies the existence of a constant A 2 > 0 such that, for all n ≥ 1, In the proof of the previous Step we also deduced that, for some A 3 > 0, Combining such two estimates with inequality (4.9) we easily get Now using this improved estimate for θ n 2 with (4.15) in (4.9) arrive at θ n (t) 2 ≤ A 5 (1 + t) −9/4 .
Going back to the differential inequality (4.14) we finally get, for some constant A 6 > 0 independent on n, and t ∈ R + , u n (t) 2 ≤ A 6 .
Replacing in (4.9) we can further improve the decay of θ n up to We now would like to improve the result of the previous Proposition by establishing decay properties for u n (t) 2 . Specifically, if we assume in addition that the linear part of the velocity satisfies e t∆ u 0 2 ≤ C(1 + t) −1/2 (this happens e.g. when u 0 ∈ L 3/2 ∩ L 2 σ ), then the same decay holds for the approximate velocities u n .
Putting this inside (4.18) and integrating on an interval of the form [t ǫ , t], with ǫ > 0 arbitrary and t ǫ chosen in a such way e t∆ u 0 2 < ǫ for t ≥ t ǫ , we obtain On the other hand, in the case

Now going back to (4.18) and integrating in time we finally get
We are now in the position of deducing our result on weak solutions to the Boussinesq system (1.1).
Proof of Theorem 2.2. Now this is immediate: passing to a subsequence, the approximate solutions θ n and u n converge in L 2 loc (R + , R 3 ) to a weak solution (θ, u) of the Boussinesq system (1.1). Moreover, the previous Lemmata imply that θ n and u n satisfy estimates of the form v n (t) ≤ f (t), for all t > 0, where f (t) is a continuous function independent on n. Then the same estimate must hold for the limit θ and u, except possibly points in a set of measure zero. But since weak solutions are necessarily continuous from [0, ∞) to L 2 under the weak topology, θ(t) and u(t) are lower semi-continuous and hence they satisfy the above estimate for all t > 0. This observation on the weak semi-continuity is borrowed from [25].

Strong solutions: preliminary lemmata
The integral formulation for the Boussinesq system, formally equivalent to (1.1) reads ∇ · u 0 = 0 The above system will be solved applying the following abstract lemma, which slightly generalizes that of G. Karch  Lemma 5.1. Let X and Y be two Banach spaces, let B : X × X → X and B : Y × X → Y be two bilinear maps and L : Y → X a linear map satisfying the estimates B(u, v) X ≤ α 1 u X v X , B(θ, v) Y ≤ α 2 θ Y u X and L(θ) X ≤ α 3 θ Y , for some positive constants α 1 , α 2 and α 3 .
Let 0 < η < 1 be arbitrary. For every (U, Θ) ∈ X × Y such that the system has a solution (u, θ) ∈ X × Y. This is the unique solution satisfying the condition Proof. In the case 0 < α 3 < 1, one can take η = α 3 . In such particular case, this lemma is already known, see [26,Lemma 2.1]. Therefore, we only have to prove that we can get rid of the restriction 0 < α 3 < 1. This is straightforward. We introduce on the space Y an equivalent norm, defined by θ Y ′ = α 3 η θ Y . By Lemma 2.1 of Karch and Prioux, applied in the space (X , · X ) and (Y, · Y ′ ), we have that if then the system (5.2) has a unique solution such that The conclusion of Lemma 5.1 is now immediate.
Remark 5.2. Using our improved version of Lemma 2.1 of [26], it is be possible to get rid of the smallness assumption |β| < 1 in the main results of Karch and Prioux [26].
Remark 5.3. The proof of Lemma 2.1 in [26] relies on the contraction mapping theorem. In particular, the solution can be obtained passing to the limit with respect to the X ×Y-norm in the iteration scheme (k = 1, 2, . . .) See [26] for more details. See also [33,Lemma 4.3] for similar abstract lemmata.
Let a ≥ 1. We define X a as the Banach space of divergence-free vector fields u = u(x, t), defined and measurable on R 3 × R + , such that, for some C > 0, In the same way, for b ≥ 3 we define the space Y b of functons θ ∈ L ∞ t (L 1 ) satisfying the estimates Such spaces are equipped with their natural norms. They are obviously decreasing with respect to inclusion as a and b grow. Recalling the definition of X and Y in Section 2.1, we see that X 1 = X ∩ L ∞ x,t with equivalence of the norms and that Y 3 = Y ∩ L ∞ x,t . The estimates (2.13)-(2.14) in Proposition 2.5 are thus equivalent to the conditions u ∈ X a and θ ∈ Y b .
We start with some elementary embeddings.
Lemma 5.4. Let L p,q be the Lorentz space, with 1 < p < ∞ and 1 ≤ q ≤ ∞. Then the following four inequalities hold: for some constants C depending only on p and q. In particular, choosing p = q one gets the correponding estimates for the classical L p -spaces.
Proof. The above estimates for the weak-Lebesgue spaces L p,∞ are simple. Indeed, if u ∈ X , then |u(x, t)| ≤ C|x| −3/p t 1 2 ( 3 p −1) and one has only to recall that any function bounded by |x| −3/p belongs to L p,∞ . The other L p,∞ -estimates for u and θ contained in (5.6) In the case 1 ≤ q < ∞, we use that L p,q is a real interpolation space between L p−ε,∞ and L p+ε,∞ , for all 1 < p − ε < p < p + ε < ∞. Therefore estimates (5.6) for all 1 ≤ q ≤ ∞ follow from the corresponding estimates in the particular case q = ∞ via the interpolation inequality.
The first useful estimate in view of the application of Lemma 5.1 is the following.
Lemma 5.5. Let 1 ≤ a < 3 and θ ∈ Y 3 . We have, for some constant C > 0 depending only on a, Moreover, Proof. We prove only (5.7) since the proof of (5.8) is essentially the same. By a renormalization, we can and do assume that θ Y 3 = 1. Let K(x, t) be the kernel of the operator e t∆ P. Then we can write K(x − y, t − s)θ(y, s)e 3 dy ds.
We have the well known estimates for K (see, e.g., [7, Prop. 1]) where C > 0 is come constant independent on x, t and 0 ≤ a ≤ 3. We also recall the scaling relation and the fact that K(·, t) ∈ C ∞ (R 3 ) for t > 0. The usual L p estimates for K are Using the L 2 -L 2 convolution inequality, we get Owing to this estimate, the conclusion L(θ) ∈ X a will follow provided we prove the pointwise inequality, This leads us to decompose where which is even better in the region {(x, t) : |x| ≥ 2 √ t} than what we need (recall that 1 ≤ a < 3). Using now |θ(x, t)| ≤ C|x| −3 and the scaling properties of K we obtain by a change of variables for |x| ≥ 2 √ t, and 1 ≤ a < 3. Next, using again |θ(x, t)| ≤ C|x| −3 and |K(x, t)| ≤ C|x| −3 , Therefore, Lemma 5.5 in now established.
We collect in the following Lemma all the estimates on B(u, v) that we shall need. (We will apply estimate (5.15) in the proof of Proposition 2.4, estimate (5.14) for Proposition (2.4) and estimate (5.14) in Theorem 2.6. Lemma 5.6. Let 1 ≤ a < 3. For some constant C > 0, depending only on a we have where (2a) * = min{2a, 4}. Moreover, Proof. We begin with the proof of the first estimate. As before, we can assume u X = v Xa = 1. We start writing where F (x, t) is the kernel of the operator e t∆ P∇.
The conclusion follows combining this with estimate (5.20). The proof of estimate (5.14) is similar. Notice the limitation (2a) * ≤ 4, which is due to the restriction on η in inequality (5.17). The proof of (5.15) also follows along the same lines and is left to the reader.
Lemma 5.7. Let a ≥ 1, b ≥ 3. For some constant C > 0 depending only on a, b, we have Moreover, Proof. As before, we give details only for the first estimate. Denoting F (x, t) the kernel of e t∆ ∇, we can write Notice that F rescales exactly as F . Moreover, These are the same estimates as for F , but there is now no limitation to the spatial decay rate (i.e., the restriction η ≤ 4 appearing in (5.17) can be removed). Therefore, the space-time pointwise decay estimates for B(θ, u)(x, t) can be proved essentially in the same way as in the previous Lemma.
The L 1 -estimate (useful for estimating the Y-norm is straightforward: This allows us to conclude. Proof of Proposition 2.4. We need two elementary estimates on the linear heat equation. Namely, Both estimates immediately follow from direct computations on the heat kernel g t (x) = (4π t) −3/2 e −|x| 2 /(4t) . (See, e.g. [4,29]) Here one only needs to use |g t (x)| ≤ C|x| −3 and the usual L 1 -L ∞ estimates for g t . Letting U = e t∆ u 0 and Θ = e t∆ θ 0 , by assumption (2.10) we get, for some C > 0, U X + Θ Y ≤ Cǫ. The system (5.1) can be written in the abstract form (5.2). By inequalities ( Proof of Part (a) of Proposition 2.5. By construction, the solution (u, θ) of Proposition 2.4 is obtained as the limit in X × Y of the sequence (u k , θ k ) defined in (5.3). By the first of assumptions (2.12), and applying straightforward estimates on the heat kernel (see also [4,29]), |e t∆ u 0 (x)| ≤ C(1 + |x|) −a and |e t∆ u 0 | ≤ C(1 + t) −a/2 . (Here we need 0 ≤ a < 3). These two conditions imply in particular that e t∆ u 0 ∈ X a . Similarily one deduces from the second inequality in (2.12) that e t∆ θ 0 ∈ Y b (when b = 3 one here needs also θ 0 ∈ L 1 ).
By estimate (5.21) and assumption (2.10), for all k = 1, 2 . . . we get If ǫ > 0 is small enough then Cǫ < 1 (the size of the admissible ǫ thus depend on b and, as we will see later, also on a). Iterating this inequality shows that the sequence (θ k ) is bounded in Y b . Combining estimates (5.7) with (5.13) we get Assuming Cǫ < 1, we deduce from the boundness of (θ k ) in Y 3 that (u k ) is bounded in X a . Thus, the solution (u, θ) belongs to X a × Y b and the first part of Proposition 2.5 follows.
Proof of Part (b) of Proposition 2.5. The proof of the second part of Proposition (2.5) is quite similar but relies on the use of slightly different function spaces. So, let a ≥ 2. We define X a as the Banach space of divergence vector fields u = u(x, t) such that, for some C > 0, Such spaces are equipped with their natural norms. Notice that the spaces X a and Y b differ from the their counterparts X a and Y b only by the fact that the time decay conditions are slightly more stringent in the former case.
where I 1 , I 2 and I 3 are as in Lemma 5.5 and and the terms contributing to I 1 are defined below. First, by the zero-mean assumption on θ. Next where we have used the Taylor formula to write the difference K j,3 (x−y, t−s)−K j,3 (x, t−s). From the bounds |K(x, t)| ≤ C|x| −3 and |∇K(x, t)| ≤ C|x| −4 we get where we used the continuous embedding of L 1 1 into Y 4 that follows from the definition of such space. The above pointwise estimate is even better, in our parabolic region |x| ≥ 2 √ t, than what we actually need. Next, can be bounded as follows for all (x, t) such that |x| ≥ 2 √ t (recall that 2 ≤ a < 4).
So far, we proved that for all (x, t) such that |x| ≥ 2 √ t. Our previous L ∞ -bound on L(θ) implies the validity of such estimate in the region |x| ≤ 2 √ t. We thus conclude that L(θ) ∈ X a and Lemma 5.8 follows.
Next Lemma is a simple variant of Lemma 5.6. Lemma 5.9. Let 2 ≤ a < 4 and b ≥ 4. Then, for some constant C > 0, Proof. This Lemma can be easily proved following the steps of estimates (5.13) and (5.21). We thus skip the details.
In the same way, one proves that by our assumptions e t∆ θ 0 ∈ Y b for b ≥ 4. Therefore, going back to the approximation scheme (5.3) and arguing as in the proof of Part (a) of Proposition (2.5) we see that the sequence (θ k ) is bounded in Y b and (u k ) is bounded in X a . Part (b) of Proposition (2.5) follows.

Asymptotic profiles and decay of strong solutions
We denote by E(x) = c |x| the fundamental solution of the Laplacian in R 3 and by E x j ,x k (x) its second order derivatives for x = 0. Notice that E xx,x 3 is a homogeneous function of degree −3. Next lemma describes the asymptotic profile for L(θ)(x, t) as |x| → ∞, by establishing that Lemma 6.1. Let θ = θ(x, t) be any function satisfying the pointwise estimates (2.14), for some 3 < b < 4 and such that θ(t) = θ 0 for all t ≥ 0. Then the j-component L(θ) j (t) of L(θ) can be decomposed as where the remainder function R ′ satisfies, In particular, in the region |x| > > √ t, one has along almost all directions.
Remark 6.2. This idea of obtaining informations on the large time behavior of solutions by first studying their behavior in the parabolic region |x| > > √ t comes from [6].
We now define R ′ (x, t) through the relation and all the previous estimates imply In the case θ 0 = 0, we can use the following variant of Lemma 6.1.
Lemma 6.3. Let θ = θ(x, t) be any function satisfying the second of (2.16), for some 4 < b < 5 and such that θ(t) = θ 0 = 0 for all t ≥ 0. Then the j-component L(θ) j (t) of L(θ) can be decomposed as Proof. We only have to reproduce the proof of the previous Lemma with slight modification. Using the estimate |θ( Next, by the vanishing mean condition I 1,1 = 0. It remains to treat I 1,3 . We can decompose I 1,3 , whose j-component we recall is into the sum of three more terms I 1,3 = I 1,3,1 + I 1,3,2 + I 1,3,3 . Such decomposition is performed exactly in the way we did in the proof of Lemma 6.1. The j-component of the first term is thus (6.6) (I 1,3,1 ) j (x, t) = − t 0 ∇K j,3 (x, t − s) · y θ(y, s) dy ds.
The second term, yθ(y, s) dy ds, can be bounded by the right-hand side of (6.5) using |∇K(x, t)| ≤ C|x| −4 and |θ(x, t)| ≤ can be treated with the Taylor formula. The simple estimate |∇ 2 x K(x, t)| ≤ C|x| −5 allows us to see that also I 1,3,3 is bounded by the right-hand side of (6.5). Therefore, both I 1,3,2 and I 1,3,3 can be included into the remainder term R ′′ (x, t).
We can now establish our main results as simple corollaries: Proof of Theorem 2.6, part (a). Let (u, θ) be a mild solution of the system (5.1), satisfying the pointwise decay estimates (2.13)-(2.14), with a > 3 2 and b > 3. Recall that the spaces X a and Y b decrease as a and b grow. Without restriction we can then assume 3 2 < a < 3 and 3 < b < 4 in our calculations. According to our notations, we can write (6.8) u(x, t) = e t∆ u 0 (x) + B(u, u)(x, t) + L(θ)(x, t).
We now deduce from Theorem 2.6 sharp upper and lower bound estimates in L p -spaces.
with 1 + 1 p = 1 α + a−r 3 (6.10) In this computation, L p r denotes as usual the weighted L p space, whereas L α,p is a Lorentz space. Here we made use of Young convolution inequality, generalised to Lorentz spaces (see [27,Prop. 2.4]).
By comparing the large time behavior of the RHS in expressions (6.9)-(6.10), we deduce the lower bound where t 0 > 0 is some constant depending on all the parameters and the initial data, but independent on t.
Let us now estimate u(t) L p r from below. The proof is based on the asymptotic expansion (2.19). Computing the third order derivatives outside the origin of the fundamental and for all t > 0 we have (6.14) u(t) L p r = ∞. Proof. Indeed, E x j ,x 3 (x) is a homogeneous function of degree 3, smooth outside the origin. Then we can find an open conic set Γ ⊂ R 3 such that |E x j ,x 3 (x)| ≥ C|x| −3 > 0 for all x ∈ Γ, x = 0, for some C > 0. Indeed, by Eq. (2.17), for |x| ≥ A √ t with A > 0 large enough and x ∈ Γ, and t > 0, we have This implies that u j has an infinite · L p r (D A ∩Γ) -norm.

Additional remarks and comments
In this section we collect a few technical remarks on the main results. These are essentially small variants of our statements that can be easily proved with minor modifications to the proofs. Under the same assumptions on θ 0 and replacing the assumption u 0 ∈ L 2 σ ∩ L 3/2 with the weaker condition e t∆ u 0 2 ≤ C(1 + t) −s , for some s ≥ 0, (when u 0 ∈ L 2 σ ∩ L 3/2 , this holds with s = 1/2), we have (7.1) u(t) ≤ C(1 + t) −s * , with s * = min(s, 1/2).
The above decay condition on e t∆ u 0 could also be restated in terms of Besov spaces.
Remark 7.2 (on Theorem 2.6). The second term in the RHS of (2.19) is bounded by C|x| −4 t. If one is interested in studying the asymptotic behavior of u only as |x| √ t → ∞ (with t > 0 not necessarily large), then an additional term in the RHS of (2.19) should be added: − ∇E x h ,x k : t 0 (u h u k )(y, s) dy ds (the : notation means that the h,k symbol has been omitted). Such term is bounded by C|x| −4 t 1/2 , thus justifying its inclusion inside the remainder when |x| > > √ t > > 1. With this additional term, the condition of the remainder can be simplified into lim |x| √ t →∞ R(x,t) t|x| −4 = 0.
Remark 7.3 (on Corollary 2.8). The restriction r + 3 p < a in Part (a) of Corollary 2.8 is natural beacause the decay assumption on u 0 guarantees that e t∆ u 0 ∈ L p r exactly for those r, p satisfying such restriction. However, estimates (2.23) remain valid in the whole range 0 ≤ r + 3 p < 3 (and under the conditions of Part (b) even for 0 ≤ r < 3 p < 4 if we want to estimate u(t) − e t∆ u 0 L p r instead of u(t) L p r in that expression.