Finite-dimensional leading order dynamics for the fast diffusion equation near extinction

We study the vanishing behavior of non-negative solutions to the fast diffusion equation on a bounded smooth domain $\Omega$ in ${\mathbf{R}}^{N}$ with Dirichlet boundary condition, that is,

For nonlinearity exponents $m\in(0,1)$, the Dirichlet condition turns the diffusion flux $mw^{m-1}$ singular on the boundary, which has the effect that for bounded solutions the mass escapes entirely in finite time $T>0$ [17, 18, 2].

We will restrict our attention to the Sobolev supercritical range $m\in(\frac{N-2}{N+2},1)$, and thus, for dimensions $N>2$ we penetrate the so-called very fast diffusion regime $m\in(0,\frac{N-2}{N})$. In the latter, certain integrability conditions on the initial data are needed in order to produce bounded solutions. We refer to Chapter 4 of the recent survey paper [4] for a discussion.

Close to the extinction time $T$, the dynamics are essentially of separated-variables type, and we consider the accordingly rescaled solutions

Introducing the reciprocal exponent $p=\frac{1}{m}\in(1,\frac{N+2}{N-2})$, the quasi-linearity in the fast diffusion equation (1.1) is translated from the spatial derivatives to the time derivative, and the resulting equation is reminiscent of the semi-linear heat equation,

It was proved by Berryman and Holland [2] that solutions to this rescaled equation approach in infinite time a stationary solution to the semi-linear elliptic equation

In fact, in spite of a possible non-uniqueness of the elliptic problem itself [12], Feireisl and Simondon showed that the dynamics (1.1) pick a unique solution of (1.3), which is determined by the initial datum $v(t=0)$ (or $w(\tau=0)$, respectively) alone [11].

Later on, Bonforte, Grillo and Vázquez [5] proved that the convergence to stationarity proceeds uniformly in the relative error $h=\frac{v-V}{V}$, that is,

Moreover, they provide sub-optimal estimates on the rate of convergence in the entropy sense for exponents $p$ close to $1$. Notice that the relative error is a natural quantity to consider because it divides out the degeneracy at the boundary, which is enforced by the Dirichlet conditions. Instead, as we will see in (1.5), the degeneracy is transferred into the linear operator.

Establishing the sharp rate of convergence has been a problem of considerable interest for many years, and it was only very recently (partially) successfully settled. Indeed, Bonforte and Figalli [3] compute the sharp rates of convergence — at least for $C^{2,\alpha}$ generic domains including the ball, and the authors of the present paper together with McCann [8] extend the results (using a completely different approach) by formulating a dichotomy: Either the relative error decays at least exponentially in time with rate given by the smallest stable eigenvalue, or the decay is algebraically $O(1/t)$ or slower. In fact, exponential convergence can be established if the limiting solution $V$ is ordinary, that is, if the manifold generated by all (weak) solutions of (1.3) has near $V$ the dimension of the kernel of the linear operator $\mathcal{L}$ associated to the dynamics of $h$ (see below), and if the error relative to $V$ embeds this manifold differentiably into $L^{p+1}(\Omega)$, which is the Lebesgue space into which $W^{1,2}(\Omega)$ embeds compactly precisely if $p<\frac{N+2}{N-2}$. Notice that the critical exponent $\frac{N+2}{N-2}$ distinguishes good from bad behavior in the treatment of the semi-linear elliptic equation (1.3), the rescaled equation (1.2) or the semi-linear heat equation. Again, we refer to [4] for more information.

To sensitize the reader to difficulties that had to one has to face when studying the convergence rates for the large-time dynamics (1.4), we consider the evolution equation for the relative error,

where $\mathcal{L}$ is the linear operator relative to $V$, that is,

and $\mathcal{M}(h)$ is the nonlinearity, given by

Observe that we may write $\mathcal{M}(h)$ in a different form $\mathcal{N}(h)$, which involves $h$ and $\partial_{t}h$, i.e., for any solution $h$ to (1.5), it holds that

Notice that the linear operator $\mathcal{L}$ is self-adjoint in the weighted Lebesgue inner product generating $L^{2}_{p+1}:=L^{2}(\Omega,V^{p+1}dx)$, and it is thus invertible without imposing boundary conditions on $\partial\Omega$, see also Theorem 2.2 below.

It is readily checked that $\mathcal{L}$ admits the negative (unstable) eigenvalue $1-p$, which corresponds to constant eigenfunctions. The occurrence of a negative eigenvalue and thus the presence of potentially unstable modes rules out soft arguments that are nowadays standard in the derivation of equilibration rates. Moreover, if $\mathcal{L}$ has a non-trivial kernel, neutral (zero) modes may potentially effect the decay behavior.

The strategy of Bonforte and Figalli to partially overcome these challenges relies on proving that the unstable modes corresponding to negative eigenvalues cannot be active during the evolution in light of Feireisl and Simondon’s convergence [11] and of Bonforte, Grillo, and Vázquez’ improvement [5]. Focusing on generic domains on which the neutral modes are absent [19], Bonforte and Figalli are able to show that the relative error is uniformly decaying with a rate $\lambda$ not smaller than the first positive eigenvalue,

Their derivation was subsequently simplified by Agaki [1], though the comparison of his bounds in terms of an energy rather than the relative error exploits Jin and Xiong’s recent boundary regularity theory [13]. On arbitrary smooth domains, however, the kernel is not necessarily trivial, hence the limiting stationary solution is in general not isolated, and thus, the proposed method is not applicable.

Finally, Jin and Xiong established unconditionally that the convergence is at least algebraic (though with an non-explicit power) [14].

In [8], in which the authors with McCann bridge the wide gap between upper exponential bounds and algebraic lower bounds on the rate of convergence, a different approach is developed. It reduces the derivation of decay estimates to the analysis of a dynamical system that fits into the setting of Merle and Zaag’s ODE lemma [16]. Exploiting the uniform convergence in (1.4), the lemma implies that the dynamics are eventually dominated either by the neutral or by the stable modes. The main (technical) work is spent in [8] to obtain new smoothing estimates, that allow to uniformly control temporal and (tangential) spatial derivatives of the solution at positive times by a weighted $L^{2}$ norm of the initial datum.

If a solution decays exponentially, apart from providing the leading order estimates on the decay of the relative error (1.9), it is showed in [8] that the large-time behavior is described by the linear dynamics up to an error of the order $e^{-2\lambda t}$ which is generated by the quadratic behavior of $N(h)$ in (1.8) near the limit $h=0$. More precisely, if the solution is decaying at rate $\lambda$ (1.9), it can be approximated by an eigenmode expansion corresponding to the eigenvalues in the interval $(\lambda,2\lambda)$. If a solution decays non-exponentially, it is eventually dominated by the neutral mode and the dynamics of dominating modes becomes non-linear in an essential way. No precise asymptotics, such as the rate of convergence and the dominating profile, is known except of the lower and upper algebraic bounds provided by [9] and [13], respectively. The problem of determining the exact algebraic or near-algebraic decay behavior in general setting is still unsettled for equations with gradient flow structures. In recent work [7], the first named author and Hung address a higher order asymptotics for a class of equations algebraically converging to smooth compact limits without boundary.

The goal of the present paper is to further refine the information on the large time dynamics. Our main result provides an unconditional convergence of the decaying solution (1.4) towards an approximate solution that lives on a finite-dimensional manifold. More precisely, denoting the (discrete [3, Lemma 2.1]) spectrum of $\mathcal{L}$ by $\sigma(\mathcal{L})=\{\lambda_{k}\,:\,k\in\mathbb{N}\}$ and assuming that the eigenvalues $\lambda_{k}$ are strictly increasing, we introduce for a given integer $K\in\mathbb{N}$ the eigenspace

and consider the orthogonal decomposition

so that $E_{s}$ is spanned by the eigenfunctions with $\lambda_{k}>\lambda_{K}$. Let $\lambda_{-}\in(\lambda_{K},\lambda_{K+1})$ be fixed.

In order to understand the meaning of the second estimate in (1.12), we recall that the decay of the stationary solution towards the domain boundary is essentially linear, that is, there is a constant $C$ such that

for any $x\in\Omega$, see, for instance, Theorem 1.1 in [10] or Theorem 5.9 in [6]. Hence, $\tilde{h}$ is a solution to the original equation (1.5) if its amplitude is small and if it is not growing too much towards the boundary. The smallness in amplitude applies (at least, after a waiting time) thanks to the convergence in (1.4). In Lemma 2.1 below, we will recall that also the second estimate is eventually satisfied.

In the construction of $\theta_{\epsilon}$ below we will see that the center manifold

is invariant under the evolution of a suitably truncated version of the dynamics (1.5). The truncation is of the kind that it is inactive for solutions to the original equation after a certain waiting time, provided that the solution is uniformly small in the sense of (1.12), cf. Lemma 2.1 below. The need of the truncation in the formulation of the theorem is a consequence of a lack of suitable $L^{q}$ regularity estimates for the linear version of (1.5) that we believe to hold true but whose proof is beyond the scope of the present paper. Nonetheless, the statement contributes to our understanding of the large-time asymptotics of the relative error. Choosing $\lambda_{K}=0$, we have the following dichotomy according to our earlier findings [8]: If the solution vanishes at an exponential rate, the approximate solution is trivial, $\tilde{h}=0$. Otherwise, if the solution is at most algebraically decaying, there is an approximate solution in the finite-dimensional space that describes the given solution up to an exponentially decaying error. It is the restriction to the finite-dimensional manifold that we consider as the major progress provided in this contribution.

In principle, that choice of $K$ in the statement of the theorem is arbitrary, and we could also choose $\lambda_{K}>0$. In this case, we could gain some information on the higher-order asymptotics similar to, e.g., [20, 21]. Because in the general setting under consideration eigenfunctions are not explicitly known, pursuing this direction is a rather abstract maneuver that we prefer to avoid.

The paper is organized as follows: In Section 2, we introduce the truncated model and discuss its well-posedness. Section 3 then provides the construction of the invariant manifolds and the proof of Theorem 1.1.

The construction of invariant manifolds necessitates the availability of a semi-flow that is defined globally on the underlying Banach space. It is thus customary to suitably truncate the nonlinearity in (1.5). We will do it in such a way that the divergence structure of the leading order term in partially retained.

For a standard cut-off function $\eta:\mathbb{R}\to\mathbb{R}_{\geq 0}$ that has a support on the interval $[-2,2]$ and satisfies $\eta=1$ on the smaller interval $[-1,1]$ and a positive real $\epsilon$, let us introduce the two truncations

and

With that we consider our evolutionary problem with a truncated nonlinearity,

where we have set

The choice of our truncations guarantees that for functions that are small in amplitude, $|h|\leq\epsilon$ and which are at most mildly increasing towards the domain boundary in the sense $V|\nabla h|\leq\epsilon$, it holds that $\mathcal{M}^{\epsilon}[h]=\mathcal{M}[h]$, and thus, the original equation (1.5) and its truncated variants (2.3) are both equivalent. This observation is the origin of (1.12).

Before turning to the analysis of the truncated equation, we remark that our decaying solution (1.9) to the original equation indeed solves the truncated equation after a possible waiting time.

In view of (1.9), the smallness condition on $h$ in the hypothesis of the lemma is always true for sufficiently large times.

We start now by discussing the well-posedness of the truncated problem.

Here, we have used $L^{2}_{2}$ to denote the weighted space $L^{2}(\Omega,V^{2}dx)$, and the Sobolev space $H^{1}_{2}$ is the one for which both the function $h$ and its gradient $\nabla h$ belong to $L^{2}_{2}$. However, the precise weight in the lower order term is arbitrary as can be seen in the following Hardy-type inequality, that will be used during the proof,

and which was established in Lemma 5.2 of [8]. Actually, via an application of Young’s inequality, the weight on the $L^{2}$ part can indeed be chosen arbitrary large.

From now on, we use $\|\cdot\|$ and $\langle\cdot,\cdot\rangle$ to denote the norm and the inner product in $L^{2}_{p+1}$.

In this section, we construct invariant manifold for the truncated equation. Our argument follows a method proposed by Koch [15], which relies on a time discretization. It is thus enough to study the evolution of the time-one maps. Koch’s method was recently used in [20, 21] to construct invariant manifolds in the context of the porous medium equation and the thin film equation both on ${\mathbf{R}}^{N}$.

Let us denote the orthogonal projections onto the eigenspaces $E_{c}$ and $E_{s}$ by $P_{c}$ and $P_{s}$, respectively. Note that $P_{c}+P_{s}=I$. Let us introduce the restricted linear operators

and the associated semi-group operators,

If $\Pi_{k}$ denotes the projection onto $\lambda_{k}$-eigenspace, $h_{k}(t):=\Pi_{k}L^{t}h_{0}$ solves the linear equation $\partial_{t}h_{k}+\lambda_{k}h_{k}=0$ and hence

Note that we have bounds on the operator norms

Indeed, this follows by

and a similar computation for $L_{s}$.

Choose two numbers $\lambda_{+}$ and $\lambda_{-}$ such that

and also denote

Observe we have an order in $\Lambda$

Let us define the time-$t$-map of truncated problem $S^{\epsilon}:L^{2}_{p+1}\to L^{2}_{p+1}$ by

where $h(t)$ is the unique weak solution of truncated problem starting from $h_{0}$ given by Theorem 2.2. In particular, we denote the time-$1$-map by $S^{\epsilon}:=S^{\epsilon}_{1}$.

We may decompose the operator $S^{\epsilon}_{1}$ as

If $\epsilon$ is small enough, the nonlinear part $R^{\epsilon}$ is a contraction. We will omit the super script truncation parameter $\epsilon$ for $S$, $R$ and others whenever there is no confusion.

Let us define the convenient norm

which is equivalent to $\|h\|_{L^{2}_{p+1}}$. For a bi-directed sequence of functions in $L^{2}_{p+1}$, say $\{h_{k}\}_{k\in\mathbb{Z}}$, let us define the norm $\vvvert\cdot\vvvert_{\Lambda_{-},\Lambda_{+}}$ and the space

Similarly, we define the norm $\vvvert\cdot\vvvert_{\Lambda_{-}}$ and the space