2.1. Exact ground state solutions

In 1d the equation (1.1) becomes

(2.1)\begin{equation} i \partial_{t} u -\partial_{x}^{4}u - 2a \partial_{x}^{2}+|u|^{\alpha}u=0. \end{equation}

Letting $u(x,t) = e^{i bt} Q(x)$ with $Q$ real (for this work we take $b \gt 0$), one gets that $Q$ satisfies

(2.2)\begin{equation} Q^{(4)} +2a\,Q^{\prime\prime}+b\,Q-Q^{\alpha+1}=0. \end{equation}

For the parameters as listed below, the following explicit solutions are known, see [Reference Albert2], [Reference Wazwaz50], [Reference Natali and Pastor37], [Reference Posukhovskyi and Stefanov44], which we use later to test our numerical solutions:

• • $\underline{\alpha = 2}$ (subcritical):  for $a \lt 0$ and $b=\frac{16}{25} a^2$,

  (2.3)\begin{equation} Q(x) = \sqrt{\tfrac{6}{5}} \, |a| \, \text{sech}^2 \Big(\sqrt{\tfrac{|a|}{10}} \, x \Big). \end{equation}

• • $\underline{\alpha=8}$ (critical):  for $a \lt 0$ and $b=25 (\frac{a}{13})^2$

  (2.4)\begin{equation} Q(x) = \left(\sqrt{105} \,\tfrac{|a|}{13} \right)^{1/4}\, \text{sech}^{1/2} \Big(2 \sqrt{\tfrac{|a|}{13}} \, x \Big). \end{equation}

• • $\underline{\alpha=10}$ (supercritical):   for $a \lt 0$ and $b=\big(\tfrac{12\, a}{37}\big)^2 $

  (2.5)\begin{equation} Q(x) = \left( \sqrt{714} \, \tfrac{|a|}{37} \right)^{1/5}\, \text{sech}^{2/5} \Big(5\,\sqrt{\tfrac{|a|}{74}}x \Big) . \end{equation}

2.2. Pokhozhaev identities

We record the Pokhozhaev identities in the 1d case, which are useful later. For the first one the equation (2.2) is multiplied by $Q$ and integrated to obtain

(2.6)\begin{equation} \|\partial_x^2 Q\|^2_{L^2} - 2a \|\partial_x Q\|^2_{L^2} + b\|Q\|^2_{L^2} - \|Q\|_{L^{\alpha+2}}^{\alpha+2} = 0. \end{equation}

For the second one, equation (2.2) is multiplied by $x \partial_xQ$ and integrated to get

(2.7)\begin{equation} 3 \|\partial_x^2 Q\|^2_{L^2} - 2a \|\partial_x Q\|^2_{L^2} - {b} \|Q\|^2_{L^2} + \tfrac{2}{\alpha+2} \|Q\|_{L^{\alpha+2}}^{\alpha+2} = 0. \end{equation}

Solving for $\|\partial_x^2 Q\|_{L^2}^2$ and $\|\partial_x Q\|^2_{L^2}$ from (2.6) and (2.7), we obtain

(2.8)\begin{equation} \|\partial_x^2 Q\|^2_{L^2} = {b} \, M[Q] - \tfrac{\alpha+4}{2(\alpha+2)} \|Q\|_{L^{\alpha+2}}^{\alpha+2}, \end{equation}

(2.9)\begin{equation} \|\partial_x Q\|^2_{L^2} = \tfrac{b}{a} \, M[Q] - \tfrac{3\alpha+8}{4a(\alpha+2)} \|Q\|_{L^{\alpha+2}}^{\alpha+2}. \end{equation}

Recalling the energy and the mass from (1.4) and (1.3), we obtain the following relation between the mass, energy and the potential term:

(2.10)\begin{equation} E[Q] = - \frac{b}2 \, M[Q] + \frac{\alpha}{2(\alpha+2)} \|Q\|_{L^{\alpha+2}}^{\alpha+2}, \end{equation}

or the energy in terms of the mass and first derivative,

(2.11)\begin{equation} E[Q] = \frac{\alpha-8}{3\alpha+8} \, \frac{b}{2} \, M[Q] - \frac{2a \,\alpha}{3\alpha+8} \, \| \partial_x Q\|_{L^{2}}^{2}. \end{equation}

From the last expression, one can observe that the following holds:

• • in the subcritical and critical cases, $\alpha \leq 8$: if the lower order dispersion coefficient is positive, $a \gt 0$, then the energy of ground states $Q^{(a)}$ is negative, $E[Q^{(a)}] \lt 0$ (here, the superscript indicates the dependence on $a$);

• • in the critical case ( $\alpha=8$) (2.11) becomes $E[Q^{(a)}] = -\frac{2a \alpha}{3\alpha +8} \|\partial_x Q\|^2_{L^2}$, and hence,

  • – in the pure quartic case $a=0$: the energy is zero, $E[Q^{(0)}]=0$,

  • – when $a \lt 0$: energy is positive, $E[Q^{(a)}] \gt 0$,

  • – when $a \gt 0$: energy is negative, $E[Q^{(a)}] \lt 0$.

• • in the supercritical case $\alpha \gt 8$: the energy is positive when $a \leq 0$.

Furthermore, for the pure quartic case $a=0$, the equations (2.6)-(2.7) yield:

• • the energy is directly proportional to the ground state mass:

  \begin{equation*} E[Q^{(0)}] = \frac{\alpha-8}{2(3\alpha+8)} b \, M[Q^{(0)}]. \end{equation*}

We confirm some of these observations in our numerical computations in Section § 2.4.

2.3. Numerical construction of ground states

In this part, we numerically construct stationary solutions to the equation (2.2) for different parameters. First, we outline the numerical approach and test it for the explicit example in the subcritical case with $\alpha=2$ in (2.3), then we consider examples for various powers $\alpha$ and values of the parameter $a$.

2.3.1. Numerical approach

We are interested in smooth solutions $Q$ of (2.2) that are critical points of the energy and vanish at infinity. These solutions decay exponentially (see, e.g., [Reference Fibich, Ilan and Papanicolaou25]), thus, Fourier spectral methods are very efficient in this case. Concretely, we apply the same approach as in [Reference Klein and Saut34], a Fourier spectral approach with a Newton–Krylov iteration. The solution can be chosen to be real (e.g., see [Reference Bugiera, Lenzmann and Sok15]) and having a positive global maximum at the origin. This is enforced during the iteration.

This means we consider equation (1.7) in the Fourier domain. We define the Fourier transform $\hat{u}$ of a function $u\in L^{2}(\mathbb{R})$ as

(2.12)\begin{align}\hat{u}(k) &= \int_{\mathbb{R}}^{}u(x)e^{-ikx}dx, \nonumber\\ u(x)& = \frac{1}{2\pi}\int_{\mathbb{R}}^{}\hat{u}(k)e^{ikx}dk,\end{align}

and write (2.2) in the form

(2.13)\begin{equation} \mathcal{F}(\hat{Q}):=\hat{Q}- \frac{\widehat{Q^{\alpha+1}}}{k^{4}-2ak^{2}+b} = 0. \end{equation}

Note that the parameter $b$ is always chosen such that $b \gt a^{2}$ implying $k^{4}-2ak^{2}+b \gt 0$. For the purpose of this section we take $b=2$.

To numerically solve the equation (2.13), we approximate the Fourier transform by a discrete Fourier transform (DFT), which is conveniently computed with a fast Fourier transform (FFT). This means that the problem is treated as a periodic problem on $L[-\pi,\pi]$, where $L \gt 0$ is chosen large enough that the considered functions and their relevant derivatives vanish at the domain boundaries to machine precision (we work here with double precision, which is on the order of $10^{-16}$). We introduce the standard discretization of the FFT for $N$ Fourier modes, $h = 2\pi L/N$, $x_{j}=-L\pi+hj$, $j=1,\ldots,N$. In an abuse of notation, we denote the DFT with the same symbol as the Fourier transform. Then the equation (2.13) becomes a system, call it $\mathcal F$, of $N$ nonlinear equations $\mathcal{F}(\hat{Q})=0$ that we solve with a standard Newton iteration

\begin{equation*} \hat{Q}^{(n+1)} = \hat{Q}^{(n)} - \left[\mbox{Jac} \big[ \mathcal{F} (\hat{Q})\big]{\big|_{\hat{Q}^{(n)}}} \right]^{-1} \, \mathcal{F}(\hat{Q}^{(n)}), \end{equation*}

where $\hat{Q}^{(n)}$ is the $n$th iterate, and where $\mbox{Jac} (\mathcal{F})$ denotes the Jacobian of $\mathcal{F}$. The action of the Jacobian on $\mathcal{F}$ is computed iteratively with the Krylov subspace method GMRES [Reference Saad and Schultz45]. Note that a further advantage of Fourier spectral methods is that the DFT $\hat{u}_k$, $k=(0,\ldots,N)/L$ decreases for analytic functions exponentially for large $|k|$. This allows to control the accuracy of the approximation of a function via the highest terms of the DFT, which are of the order of the numerical error, see for instance the discussion in [Reference Trefethen49]. In this paper we always control the spatial resolution in this way.

2.3.2. Test against an exact solution

We now test this method for the explicitly known example (2.3), for which we set, for instance, $b=2$. We choose $L=20$, $N=2^{8}$, and the initial iterate $Q^{(0)}(x)=0.5\,e^{-x^{2}}$. The convergence of Newton iterations obviously depends strongly on the choice of the initial iterate, which is here similar in size to the exact solution (2.3), but with a considerably faster decay for large $x$ (that is, $e^{-x^2}$ vs. $e^{-a|x|}$). Nevertheless, the computation converges after 10 iterations (to be precise, it is stopped once the residual $\|\mathcal{F}\|_{\infty} \, \lt 10^{-10}$).

The difference between numerical and exact solution is on the order of $10^{-13}$, i.e., roughly on the order of the rounding error, which we illustrate in Figure 3. We note that the residual $\|\mathcal{F} \|_{\infty}$ for the exact solution is on the order of $10^{-15}$.

Figure 3. The difference between the numerically computed $Q$ from (2.13) and the exact solution from (2.3).

2.3.3. Examples

Below we show several examples of solutions to (1.7) with a fixed value $b=2$, varying the nonlinearity power $\alpha$ and the coefficient of the lower (2nd order) dispersion $a$. Here, we take $L=10$, $N=2^{10}$ and the initial iterate $Q^{(0)}(x) = 1.5e^{-x^{2}}$ unless noted otherwise.

On the left of Figure 4 one can see the profiles of the ground states for a fixed $a=1$ with varying nonlinearity power $\alpha$: between $2$ (subcritical case) and $10$ (supercritical case), noting that $\alpha=8$ corresponds to the critical case. All solutions are non-monotonic and all have a depression into negative values around the central hump and then continuing out with damped oscillations; this is due to the higher order dispersion (and mixed dispersion), which breaks positivity: in this equation the fourth order dispersive term is coupled with the second-order term, and unless the second order dispersion is ‘helping’ the higher order dispersion with the very negative coefficient ( $a\ll 0$), the profile will have oscillations; similar phenomena are seen in other equations, for instance, in the case of the Benjamin equation [Reference Albert, Bona and Restrepo3, Reference Benjamin9]. There is no decisive effect of different nonlinearities on the overall shape of the solutions for a fixed $a$, just the overall height is slightly decreasing and the larger values of the nonlinearity $\alpha$ lead to slightly smaller amplitudes.

On the right of Figure 4, one can observe changes in the profiles for a (fixed) cubic nonlinearity ( $\alpha=2$) with varying $a$, the coefficient of the lower (second order) dispersion. In particular, the more negative $a$ becomes, the less and less oscillations can be seen. We discuss this in further details in the next figure.

Figure 4. Profiles of ground state solutions to (2.2) with $b=2$. Left: $a=1$, $2 \leq \alpha \leq 10$. Right: cubic nonlinearity ( $\alpha=2$), coefficient of lower dispersion $a= -1,0,1$.

As far as the numerical computations of the profiles, we point out that some relaxation is needed for negative values of $a$ when computing the profile: instead of $Q^{(n+1)}$ of the Newton iteration, the value $\mu Q^{(n+1)}+(1-\mu)Q^{(n)}$ with $\mu\in(0,1)$ is chosen as the new iterate (for instance, with $\mu=0.1$).

To make further clarification about the oscillatory vs. monotone nature of the profiles, we plot the solution for $a=-\sqrt{2}$ on the left of Figure 5 (for this example we chose the initial iterate $Q^{(0)}(x)=2e^{-x^{2}}$), noting that it is positive and monotone (as a function of $|x|$). This property of positivity and monotonicity (in $|x|$) will be shared by all profiles with $a \leq -\sqrt 2$, which is a confirmation of the results on positivity and radiality of ground states in the regime $a \leq -\sqrt{b}$ in a general case from [Reference Bonheure and Nascimento13] and [Reference Bonheure, Casteras, dos Santos and Nascimento12]. As we increase $a$ from the value $-\sqrt 2$, we start observing more and more oscillatory behaviour of the profiles while decreasing in height as shown in Figure 4, and as $a\to \sqrt{2}$, solutions become very oscillatory, with height eventually diminishing to zero. We show the (almost limiting) case $a=1.4$ on the right of Figure 5 (the profiles exist for $a \lt \sqrt b = \sqrt 2$). For this computational example we chose $L=40$, $N=2^{11}$ and the initial iterate $Q^{(0)}(x)=0.5e^{-x^{2}}$ with the relaxation discussed above. Note that the height of the solution has decreased to 0.35 and the number of oscillations has significantly increased compared to the cases $a=0, 1$ as in Figure 4.

Figure 5. Profiles of ground state solutions to (2.2) with $b=2$ and cubic nonlinearity ( $\alpha=2$). Left: $a=-\sqrt{2}$. Right: $a=1.4$.

2.4. Bifurcation of ground states

We investigate the dependence of ground state solutions $Q$ on the parameters $a$ and $b$ more carefully.

2.4.1. Dependence on $a$

First, we fix $b=2$ and study how the properties of ground state solutions to (2.2) change with the parameter $a$, the coefficient of the lower order dispersion in (1.1). For a given value of $a$, we denote the ground states as $Q^{(a)}$ and recall that ground states only exist for $a \lt \sqrt b \equiv \sqrt 2$.

We consider the critical case $\alpha=8$ and in Table 1 give the values of the ground states mass $M[Q^{(a)}]$ for several values of $a$. One can observe that the mass is not monotonic in this dependence.

To study this further we investigate the dependence of $a$ on several quantities of $Q^{(a)}$, such as the mass, energy, and the $L^\infty$ norm, see Figure 6.

Figure 6. $\alpha=8$. Dependence of the ground state mass $M[Q^{(a)}]$ on the parameter $a$ for a fixed $b=2$. Mass (left), energy (middle), $L^{\infty}$ (right).

Table 1. Mass $M[Q^{(a)}]$ for different $a$ with fixed $b=2$, $\alpha=8$

One notices that the mass has a minimum around $a \approx 1$, while the energy and the sup norm are decreasing as $a \to \sqrt 2$. The plot of the $L^\infty$ norm shows that the ground states decrease in their height as $a$ increases, being consistent with Figure 5 and the right plot of Figure 4, and since the mass is increasing, they gain more and more oscillatory behaviour as $a \to \sqrt 2$ as shown in Figure 5. Note from the middle plot of Figure 6 that the energy is zero when $a=0$ (the scaling invariant case), $E[Q^{(0)}] = 0$, it is positive for $a \lt 0$ and negative for $a \gt 0$, as it was proved at the end of Section § 2.2. This will play an important role in the investigation of blow-up in Section § 6.

We investigate the subcritical case $\alpha = 2$ in Figure  7 and observe that, opposite to the critical case, the mass is decreasing as $a$ increases to $\sqrt 2$ value (left plot). The energy is increasing, having a small dip around $a\approx 1.3$ (in the middle plot), and the sup norm decreases as before (right plot), indicating that the oscillatory envelope of the ground state is decreasing in height, and since the mass is decreasing, the decay of the envelope of oscillations is much faster than in the critical case.

Figure 7. $\alpha=2$. Dependence of the ground state mass $M[Q^{(a)}]$ on the parameter $a$ for a fixed $b=2$. Mass (left), energy (middle), $L^{\infty}$ (right).

2.4.2. Dependence on $b$

We next fix $a=1$, the lower order dispersion coefficient (note that in this case both dispersions work against each other), and track the dependence of the ground state $Q^{(1)}$ quantities on $b$. We investigate cases from subcritical to supercritical: $\alpha = 2,4,6,8,10$, to show a new phenomenon about the ground states.

We plot cases $\alpha = 2,4,6$ in Figure 8 and $\alpha=8,10$ in Figure 9.

Figure 8. Dependence in sub-critical cases $\alpha=2$ (top), $\alpha=4$ (middle), $\alpha=6$ (bottom), of $M[Q^{(1)}]$ and $E[Q^{(1)}]$ on the parameter $b$ for a fixed $a=1$ (left and middle columns). Dependence of energy as a function of mass, $E = E(M)$ (right column).

Figure 9. Dependence in the critical $\alpha=8$ (top) and super-critical $\alpha=10$ (bottom) cases, of $M[Q^{(1)}]$ and $E[Q^{(1)}]$ on the parameter $b$ for a fixed $a=1$ (left and middle columns). Dependence of energy as a function of mass, $E = E(M)$ (right column).

First, note that while in the sub-critical case $\alpha=2$ all graphs look to be monotone, in the cases $\alpha=4,6$, there is a different behaviour: initially, the mass is decreasing as $b$ increases and then it starts increasing, while the energy does exactly the opposite; the reverse behaviour is consistent with the dependence shown in (2.10). This change in monotonicity, when plotted as a function, where the energy dependence on the mass (in the computed range of $b$ between $1$ and $4$), i.e., $E = E(M)$, shows the appearance of two branches in the energy, see the right plots in Figure 8. This means that for the same mass there are two solutions for the ground state. This, in some sense, resembles the behaviour of the NLS equation with a combined nonlinearity or a double well potential, for example, see [Reference Carles, Klein and Sparber17]. We, therefore, call the ground state from the upper branch an unstable branch, and the lower one – a stable branch. As we show later, when the ground state from the upper branch is perturbed such that it has a larger mass than the unperturbed ground state, it will jump to the lower branch (with the lower energy) and will try to approach asymptotically that ground state. On the other hand if the perturbation leads to a situation with less mass than the unperturbed state, the initial data are simply dispersed. This is investigated in more details in Section 4.2.

In Figure 9, we plotted the dependencies of the same conserved quantities (mass and energy) on the parameter $b$ in the critical and supercritical cases. Observe that in the critical case the behaviour of the quantities is similar to the sub-critical cases with $\alpha=4$ and $6$, producing a bifurcation in the energy vs. mass plot, $E = E(M)$, see top row in Figure 9.

In the supercritical case, $\alpha=10$, the dependence of mass and energy becomes monotone (as in the case $\alpha=2$), and thus, no more bifurcation is present in this case, see the bottom row of Figure 9 (at least in the range of $b$ that we computed).

For completeness in the critical case, we also provide the dependence on $b$ of the potential term, $L^{10}$-norm, in Figure 10, which shows a linear dependence on $b$.

Figure 10. Dependence of $\|Q^{(1)}\|_{L^{10}}^{10}$ on $b$ in the critical case $\alpha=8$.

Having examined the plethora of ground states, we proceed in the next section onto studying the time evolution of ground states and other solutions to the equation (1.1).

4.1. Near soliton dynamics, non-branching case

We start with considering the sub-critical case of (2.1) with $\alpha=2$, since this nonlinearity did not show any branching in our investigation of the energy vs. mass $E(M)$ dependence of solutions to (2.2), in Figure 8.

We fix $b=2$ and take initial data as the ground state perturbations of the form

(4.1)\begin{equation} u(x,0) = A \, Q(x),\quad A \sim 1. \end{equation}

$\scriptstyle\blacklozenge$ ( $a=0$) We begin with the pure quartic case $a=0$, since it is a scaling-invariant case, and hence, the ground state (1.8) scales as well.

The $L^{\infty}$ norm of the solution with the initial data (4.1) and $A=1.01$ is shown in plot (A) on the left top of Figure 11. It can be seen that the $L^{\infty}$ norm of the solution grows slightly and then oscillates, approaching in amplitude a level of the value slightly higher than the unperturbed ground state. This is due to the fact that while analytically the perturbations are infinitesimally small, numerically we have to consider a finitely small perturbation in order to observe a visible effect of the perturbation (in finite time). The perturbation leads to a ground state, but of a slightly higher $L^\infty$ norm. A fit to the ground state for a parameter $b$ according to (1.8), namely, $Q_{b} = (b/2)^{1/\alpha}Q((b/2)^{1/4}x)$, is done in the following way: $(b/2)^{1/\alpha}$ is given by the maximum of the modulus of the solution at the final time divided by the maximum of $Q$. The difference between the modulus of the solution at the final time and the fitted ground state $|Q_b|$ is shown on the right of the same figure. It is on the order of $10^{-3}$, since the final state is not yet reached. Note there will always be oscillations around the final state, since there is no dissipation in the systems (see for instance [Reference Carles, Klein and Sparber17]), and since we are working on the torus, no radiation can escape to infinity. In fact, the size of the torus has to be chosen large enough in order to limit the effects of radiation reentering the computational domain, see also [Reference Arbunich, Klein and Sparber4, Reference Klein, Sparber and Markowich35] for similar perturbations of the NLS ground states. The jumps in the $L^{\infty}$ norm are due to the fact that it is determined at the collocation points of the FFT, and thus, the result depends on whether the actual maximum of $|u|$ is on a grid point or not.

Figure 11. Solution to (2.1), $\alpha=2$, $a=0$, $b=2$ with $u_0=1.01 Q$ (top) and $u_0=0.99 Q$ (bottom). Left: time dependence of the $L^{\infty}$ norm. Right: the difference of the modulus of the solution at $t=10$ and the rescaled ground state $Q_b$.

The situation is similar for a factor $A=0.99$. As can be seen on the bottom left of Figure 11, or plot (C), the $L^{\infty}$ norm of the solution decreases slightly at the beginning and then decreasingly oscillates in amplitude around a slightly smaller (shorter in height) ground state. The difference of the modulus of the solution at the final time with the fitted ground state $Q_b$ (for a fitted $b$ as discussed above) can be seen on the right of the same figure, plot (D). The difference is again on the order of $10^{-3}$, which provides a good numerical confirmation about the asymptotic approach to a soliton state, and thus, in a sense a soliton resolution.

$\scriptstyle\blacklozenge$ If $a\neq0$, there is no simple scaling in $b$ of the form (1.8) as it exists when $a=0$. Thus, we can only consider the $L^{\infty}$ norm in these cases. In the top row of Figure 12 we show the $L^{\infty}$ norms of the solution $u(t)$ with initial data (4.1), $u_0 = A\, Q^{(a)}$, $A=1.01$, for $a=\pm1$. The behaviour is similar (converging with oscillations) to what is shown on the left of Figure 11.

Figure 12. Time dependence of the $L^\infty$ norms of solutions to (2.1), $\alpha=2$, $b=2$, with $u_0=A Q^{(a)}$ with $A=1.01$ (top), $A=0.99$ (bottom). Left: $a=1$. Right: $a=-1$.

On the bottom row of Figure 12 we show solutions with perturbations by $A=0.99$ and $a=\pm1$. The behaviour is similar to the case of $a=0$ shown on the left of Figure 11.

We conclude that the ground states in the case of subcritical nonlinearity $\alpha=2$ exhibit a stable behaviour for different values of parameter $a$. As no branching of energy occurs (at least in the range of parameters that we considered), it also shows how a stable ground state behaves under small perturbations in the case of combined dispersions, that is, small amplitude changes in $A$ force the solution to slightly increase or decrease in order to ‘find’ the rescaled version of itself and asymptotically (and in oscillatory manner with very small amplitude oscillations) approach it. We refer to this behaviour as the stable perturbation of the ground state, when considering different branches of ground states below.

4.2. Near soliton dynamics, branching case

Here we consider the nonlinearities where the branching of the energy was observed, recall Figures 8 and 9.

$\scriptstyle\blacklozenge$ Subcritical case. We take $\alpha=6$, fix $a=1$, and consider the initial data $u_0=A\, Q^{(1)}$ with different values of $b$ and $A \sim 1$. In our first example, we consider $b=1.1$, noting that the mass is relatively large, see Figure 8(G), therefore, we expect this ground state solution to be unstable, which is indeed the case, shown in Figure 13. We point out that the behaviour of these perturbations is qualitatively different from those in Figure 11 and 12, the perturbations with $A \lt 1$ scatter to zero, those with $A \gt 1$ create large oscillations, of different nature, than those shown in Figure 11(A), 12(A)&(B), and seem to approach a different state (a stable ground state).

Figure 13. Time dependence of the $L^\infty$ norm of solutions to (2.1), $\alpha=6$, $a=1$ with $u_0=.99Q^{(1)}$ (left) and $u_0=1.01Q^{(1)}$ (right) for $b=1.1$.

To further confirm the unstable behaviour and the unstable branch, we take the initial data, of approximately the same mass, and compare two solutions. For example, fixing the mass around $2.8$, we find that two values of $b$ can produce the ground state solutions, namely, $b\sim$1.3 and $b\sim$3.5. We fix these values, compute their energies, and then run the simulations to compare the behaviour of perturbations. We have for $b=1.3$, $M[Q^{(1)}] = 2.7$, $E[Q^{(1)}]=-1.23$ and for $b=3.5$, $M[Q^{(1)}] = 2.9$, $E[Q^{(1)}]=-1.59$. Since the energy is larger for smaller $b$, we conclude that the ground state for a smaller $b$ value should be unstable. We confirm this in Figure 14.

Figure 14. Solution to (2.1), $\alpha=6$, $a=1$ with $b=1.3$ (top), $b=3.5$ (bottom), and $u_0=0.99 Q^{(1)}$ (left), $u_0=1.01 Q^{(1)}$ (right). The initial data are chosen such that the mass is approximately the same, but the ground state in the top row has higher energy, thus, unstable.

$\scriptstyle\blacklozenge$ Critical case. We also study the critical case, $\alpha=8$, to confirm existence of stable and unstable branches of ground states. We continue with $a=1$ and consider three cases of the parameter $b = 1.1, 1.3, 4$. The first two values produce the ground state with the energies on the lower branch in Figure 9(C), and thus, are expected to be unstable. We confirm that in Figure 15.

Figure 15. Solution to (2.1), $\alpha=8$, $a=1$, with initial data $u_0$ and the parameter $b$ as indicated.

One can observe a stable behaviour for $b=4$ in the left column of Figure 15, with small oscillations (eventually converging to some final state); however, the other two values of $b$ show the unstable behaviour with dispersing down to zero for amplitude perturbations $A \lt 1$ and having large oscillatory behaviour for $A \gt 1$. In particular, one can notice the convergence in the right top of Figure 15 to a higher level, indicating the convergence to a stable ground state. We also note that in the critical case, we do not see formation of blow-up right after the value $A=1$, which is different from the classical case. Thus, this is a qualitatively different behaviour of ground states from the classical (as in the NLS, e.g., see [Reference Holmer and Roudenko29]) ground state sharp threshold for blow-up vs. scattering behaviour. This happens in the non-scaling invariant cases ( $a \neq 0$), since the lower order dispersion (with either sign) cannot be scaled away.

5.1. Behaviour of solutions for Gaussian initial data

We start with the following Gaussian initial data

(5.1)\begin{equation} u_0(x)=A e^{-x^{2}}, \quad A \gt 0. \end{equation}

(Note that $M[u_0]=A^2 \sqrt{\frac{\pi}{2}}$.)

The computation is done on large tori in order to avoid emission of radiation towards one boundary of the computational domain and its reappearance on the other side (due to the periodicity of the setting), which may have a strong effect on the results. But since, as can be seen below, there is a large amount of radiation, this cannot be completely avoided even on large tori. We consider $x\in[-100\pi,100\pi]$ with $N^{12}$ Fourier modes and $N_{t}=10^{4}$ time steps for $t\in[0,10]$. The Fourier coefficients decrease to machine precision during the whole computation, and the discretized energy is conserved in relative error to the order of $10^{-9}$.

$\bullet$ Gaussian, $a=0$. Fixing $a=0$ and $\alpha=2$ in (2.1), we show its solution with the initial condition (5.1) and $A=2$ in Figure 16. The initial maximal peak height drops down from $2$ to about $1.2$ and an asymptotic profile appears to emerge, although on a non-negligible radiation background. The $L^{\infty}$ norm of the solution is shown in the middle of Figure 16.

Figure 16. Left: time evolution of the solution to (2.1), $\alpha=2$, $a=0$, with $u_0=2\,e^{-x^2}$. Middle: time dependence of the $L^{\infty}$ norm of the solution. Right: solution at the final computational time (blue) and a fitted rescaled ground state (green).

The norm appears to oscillate asymptotically approaching some final state, though the radiation in the background becomes visible after some time (see left and middle plots after about $t=5$). The interpretation of the asymptotic final state as a ground state is supported by the right plot of the same figure, where we show a fit to a rescaled soliton according to (1.8). This is one of the biggest advantages of the pure quartic case, having the scaling. This also confirms the soliton resolution, where localized smooth generic data splits into a rescaled and shifted soliton(s) and radiation.

$\bullet$ Gaussian, $a=1$. Considering the same Gaussian initial data (5.1) for the equation (2.1) but with $a=1$, we get the solution shown on the left of Figure 17. The height of the main peak drops down from $2$ to around $1.5$ and then oscillates to the asymptotic final state while emitting non-trivial amount of radiation. This interpretation is in accordance with the $L^{\infty}$ norm of the solution on the left of Figure 17. The norm appears to oscillate around the height of the final state, but since the computation is done on a large torus, the radiation is clearly visible for large times (see middle plot after about $t=4.5$). The solution at the final computational time also suggests that it can be considered as a ground state, though no scaling as in (1.8) is available, and thus, we do not perform any fitting.

Figure 17. Left: time evolution of solution to (2.1), $\alpha=2$, $a=1$, with $u_0=2\,e^{-x^2}$. Middle: Time dependence of $L^{\infty}$ norm. Right: solution at the final computational time.

$\bullet$ Gaussian, $a=-1$. We next consider the opposite sign of the parameter $a$, namely, $a=-1$, where the ground state $Q^{(-1)}$ has a larger maximum in amplitude than in the cases $a=0$ or $a=1$, as can be seen on the right of Figure 4. The solution to the bi-harmonic NLS equation (2.1) with the same Gaussian initial condition (5.1) and $A=2$, appears to radiate away (or scatter), no stable structures emerge.

However, if we consider the initial condition (5.1) with a larger amplitude, for example, $A=3$, an asymptotic final state plus radiation is observed as in previous examples, see Figure 18. The solution approaching an asymptotic final state is confirmed by the $L^{\infty}$ norm of the solution in the middle of Figure 18, whereas before, the solution is oscillatory approaching the asymptotic profile; we plot the solution at the final computational time on the right. Since there is no scaling in this case either, there is no matching to the ground state provided.

Figure 18. Left: time evolution of solution to (2.1), $\alpha=2$, $a=-1$, with $u_0=3\,e^{-x^2}$. Middle: Time dependence of $L^{\infty}$ norm. Right: solution at the final computational time.

We conjecture that the scattering, or dispersion of the solution down to zero, happens if the mass of the initial data is below some threshold that is less than $\min \{Q^{(0)}, Q^{(a)} \}$ for a given power $\alpha \lt 8$ (and $a \lt \sqrt{b}$), which would be interesting to investigate further. Above this threshold, localized smooth solutions confirm the soliton resolution.

6.1. Behaviour above the mass of ground states

We look at different solutions in the critical case, recalling from Section 4.2 (Figure 15, top row, and Figure 20, right plot) that perturbations of ground states $Q^{(a)}$ have a somewhat different behaviour than in the standard NLS case: for very small perturbations of $A Q^{(a)}$, $a \neq 0$, with $1 \lt A \lt 1+\epsilon$, the solutions do not blow up, but asymptotically (in oscillatory manner) approach another stable state. For the first several examples of ground state perturbations we fix $b=2$, later we vary this value.

We first consider cases $a \leq 0$ and then $a \gt 0$, since for small $\epsilon \gt 0$ the mass $M[Q^{(0+\epsilon)}] \lt M[Q^{(0)}] \lt M[Q^{(0-\epsilon)}]$ as shown in Figure 6.

$\scriptstyle\blacklozenge$ Case $\underline{a=0}$. We consider $u_0(x)=A\, Q^{(0)}(x)$, with $A \approx 1$. The time evolution of the $L^\infty$ norm of solutions with $A=0.9$ or $1.1$ is shown in Figure 19 on the left and right, correspondingly. In the first case the solution is dispersed, whereas in the second case it blows up in finite time. We also recall that the energy $E[Q^{(0)}] = 0$ and the value of the ground state mass does not depend on $b$ in this pure quartic dispersion case, as the solution can be rescaled.

Figure 19. Time dependence of the $L^\infty$ norm of solutions to (2.1), $\alpha=8$, $b=2$, $a=0$, with $u_0 = A\, Q^{(0)}$, $A=0.9$ (left) and $A= 1.1$ (right).

Furthermore, since the scaling invariance holds in the $a=0$ case, we are able to fit the final computational state with the (numerical) rescaled ground state profile of $Q^{(0)}$, see details in Section 6.2.2.

$\scriptstyle\blacklozenge$ Case $\underline{a \lt 0}$. We fix $a=-1$ and observe that the perturbation of $u_0 = A\, Q^{(-1)}$ with $A = 0.95$ produces a bounded (oscillatory) behaviour and with $A=1.1$ the solution blows up in finite time.

Next, we consider $a=-2$ and compute the corresponding ground state solution $Q^{(-2)}$ from (2.2), which has no oscillations in this case (since $a \lt -\sqrt 2$). We consider perturbations $u_0 = A Q^{(a)}$ with $A=1.1$. This initial data produces a blow up in finite time, see Figure 25 (we discuss this case in more details below in terms of the blow-up rate and profile).

$\scriptstyle\blacklozenge$ Case $\underline{a \gt 0}$. We fix $a=1$ and consider $u_0=A\, Q^{(1)}$, recalling that for $A=0.99$ the solution dispersed and for $A=1.01$ it asymptotically (in oscillatory manner) approached some stable state, see Figure 15. This behaviour was observed for small perturbations with $1 \lt A \lesssim 1.1$. Only when $A \gt rsim 1.2$ we observe that the perturbations blow up in finite time, see Figure 20. We note that $M[Q^{(1)}] \lt M[Q^{(0)}]$ (see Table 1), and when $A=1.10055$, we get $M[Q^{(0)}] = M[A\, Q^{(1)}]$. Thus, in this case, we observe that

• • scattering happens below the mass of $Q^{(1)}$ (and in this case, $M[Q^{(1)}] \lt M[Q^{(0)}]$),

• • there is a gap between $M[Q^{(1)}]$ and $M[Q^{(0)}]$, where solutions with perturbed $A\, Q^{(1)}$ data do not scatter but instead approach, in oscillatory manner, some final state,

• • when the mass of the initial data gets above $M[Q^{(0)}]$, solutions blow up in finite time.

Figure 20. Time dependence of the $L^\infty$ norm of the solution $u(t)$ for $b=4$, $\alpha=8$, $a=1$ with $u_0=AQ^{(1)}$, $A=0.9$ (left) and $A=1.2Q^{(1)}$ (right).

Fixing $a = 0.5$ and considering the ground state $Q^{(0.5)}$ of (2.2) and its perturbations $u_0 = A\, Q^{(0.5)}$, we observe that $A = 1.05$ produces a bounded (oscillatory) behaviour and with $A=1.1$ the solution blows up in finite time. Noting again that $M[1.1 Q^{(0.5)}] = 3.23 \gt M[Q^{(0)}]$ and $M[1.05 Q^{(0.5)}] = 2.95 \lt M[Q^{(0)}]$ gives a positive confirmation to the statement that the behaviour of solutions is determined by the size compared to either the ground state $M[Q^{(1)}]$ or the scaling-invariant ground state $Q^{(0)}$ and its mass.

Fixing $a=1.35$, we first note that $M[Q^{(1.35)}] \approx 3.41 \gt M[Q^{(0)}]$, and recalling the left plot of Figure 6, which indicates that this ground state is on the unstable branch, we investigate the perturbations of $u_0 = A\, Q^{(1.35)}$. For $A = 1.2$ the solution produces a bounded (oscillatory) behaviour and with $A=1.5$ the solution blows up in finite time. While $M[1.2Q^{(1.35)}]$ is greater than that of the scaling invariant ground state mass, $M[Q^{(0)}]$, we point out that the observed behaviour could be an indication that we are dealing with the unstable branch of ground states. At larger magnitudes of $A$ we do observe blow-up in finite time.

6.1.1. Behaviour of other types of data above the threshold

We consider different types of initial data to investigate further the behaviour of solutions above the corresponding ground state’s mass in the critical case.

$\scriptstyle\blacklozenge$ We take super-Gaussian initial data, compute its mass, and investigate its evolution,

\begin{equation*} u_0 = A\, e^{-x^4}, \qquad M[u_0] = 2^{\frac34} \Gamma(\tfrac54) A^2. \end{equation*}

We note that the mass of super-Gaussian data is smaller than the mass of the pure quartic ground state $Q^{(0)}$, that is, $M[u_0] \lt M[Q^{(0)}]$, when (approximately) $A \lt 1.4$.

In Fig 21 we show the $L^\infty$ norm of solutions to (1.1), $\alpha=8$, $a=1$, for the super-Gaussian initial data with $A=1.3$ on the left and with $A=1.4$ on the right. For values of $A$ larger than $1.4$, we observe finite time blow-up. Thus, the blow-up occurs for the initial data with the mass above that of the scaling-invariant ground state $M[Q^{(0)}]$ (though the equation has the second-order dispersion term with $a=1$), then there is a gap where the solution does not seem to neither blow-up nor scatter down to zero, but instead approach some final state, and below certain value the solutions scatters to zero.

Figure 21. Left column: time dependence of the $L^\infty$ norm of solution to (2.1), $\alpha=8$, $a=1$, with $u_0 = A e^{-x^4}$, $A=1.3$ (left) and $A=1.4$ (right).

It seems to be plausible to conjecture that if the mass of the non-scaling-invariant ground state is less than the mass of $Q^{(0)}$, i.e., if $M[Q^{(a)}] \lt M[Q^{(0)}]$, then the solution scatters to zero in that case. There seems to be a small gap where solutions do not disperse to zero, nor blow-up, but rather approach another non-zero state. Eventually, the perturbations with larger $A$ blow-up in finite time, starting from a certain threshold above that mass of $Q^{(0)}$ (or above $M[Q^{(a)}]$, if that value is larger).

$\scriptstyle\blacklozenge$ Another example we look at is of slower decay than Gaussian or super-Gaussian, namely,

\begin{equation*} u_0=A\, \text{sech}\, x, \quad M[u_0] = 2 A^2. \end{equation*}

The condition $M[u_0] \lt M[Q^{(0)}]$ for such data holds when $A \lt 1.23$. We show the dynamics of solutions for such data in Figure 22 with $A=1.2$ on the left, which disperses, and with $A=1.3$ on the right, which oscillates toward some final state. For values much larger than 1.3, $A \gg 1.3$, we observe blow-up.

To summarize the results of this part, we find that in the critical case, the thresholds for scattering (down to zero) or blow-up are not easily identified for different parameters of lower dispersion $a$ or initial magnitudes and profiles (and $b$); the values of the mass of corresponding ground states, and especially the scaling-invariant value $M[Q^{(0)}]$ should be taken into account, as well as the values (and possibly other properties) of the ground states $Q^{(a)}$. Furthermore, the stable and unstable branches of the ground states play a significant role in the solutions behaviour, possibly identifying thresholds for scattering down to zero or approaching some final states, or blow-up in finite time. We do show a new behaviour in the critical case, which is significantly different from classical dichotomy of scattering vs. blow-up, as we observe a ‘gap’ when solutions neither scatter nor blow-up but rather approach some final state.

Figure 22. Left column: time dependence of the $L^\infty$ norm of solution to (2.1), $\alpha=8$, $a=1$, with $u_0 = A \, \text{sech} x$, $A=1.2$ (left) and $A=1.3$ (right).

We next investigate self-similarity in blow up solutions.

6.2. Self-similar blow-up

A typical approach to study the self-similar blow-up solutions in equations with scaling symmetry, for example, as (1.5), is to consider the method of dynamical rescaling (see [Reference Sulem and Sulem46, § 6.1.2]) with a time dependent factor $L(t)$. To be precise, set

(6.1)\begin{equation} \xi=\frac{x}{\lambda(t)},\quad \frac{d\tau}{dt}=\frac{1}{\lambda(t)^{4}},\quad u(t,x) = \frac1{\lambda(t)^{4/\alpha}} \,U (\tau,\xi), \end{equation}

with the factor $\lambda(t)$ to be chosen such that it tends to zero as $t\to t^{*}$ (finite blow-up time $t^{*} \lt \infty$), whereas the rescaled time variable $\tau$ tends to infinity in this limit. This time dependent change of variables replaces the equation (2.1) with

(6.2)\begin{equation} iU_{\tau}-i(\ln \lambda)_{\tau}(\tfrac4{\alpha}\, U+\xi U_{\xi})-U_{\xi\xi\xi\xi} -2a\lambda^{2}U_{\xi\xi}+|U|^{\alpha}U=0. \end{equation}

In order to look for singular self-similar type solutions near the blow-up, the following ansatz

\begin{equation*} U(\tau,\xi) = e^{ib\tau} \mathcal R^{(\tau)}(\xi), \end{equation*}

is used to obtain the profile equation (we drop the superscript $\tau$ in the profile $\mathcal R^{(\tau)}$ for brevity):

(6.3)\begin{equation} -b\, \mathcal R -i(\ln \lambda)_{\tau}(\tfrac4{\alpha}\, \mathcal R + \xi \mathcal R_{\xi}) - \mathcal R_{\xi\xi\xi\xi} -2 a \lambda^{2} \mathcal R_{\xi\xi}+|\mathcal R|^{\alpha} \mathcal R=0. \end{equation}

In the limit $\tau\to\infty$ (i.e., $t \to t^*$), the blow-up profile $\mathcal R^{(\tau)}$ is expected to become $\tau$-independent, denote it by $\mathcal R^\infty$. Setting $A^{\infty}:=\lim_{\tau\to\infty}(\ln \lambda)_{\tau}$ and observing that, since $\lambda(t) \to 0$ as $\tau \to \infty$ ( $t \to t^*$), the term proportional to $a$ vanishes and we obtain the profile equation from (6.3) as follows

(6.4)\begin{equation} -b \, \mathcal R^{\infty}-i\,A^{\infty}(\tfrac4{\alpha} \, \mathcal R^{\infty}+\xi \mathcal R^{\infty}_{\xi}) - \mathcal R^{\infty}_{\xi\xi\xi\xi} +|\mathcal R^{\infty}|^{\alpha} \mathcal R^{\infty}=0. \end{equation}

In the critical case ( $\alpha=8$), the coefficient $A^{\infty}$ becomes $0$ (recall that in the NLS equation the rate of how fast this coefficient decreases to zero determines the correction in the stable blow-up, so called ‘log-log’ correction of the square-root rate), and the profile equation (6.4) reduces to the equation $b\, \mathcal R^{\infty} + \mathcal R^{\infty}_{\xi\xi\xi\xi} -|\mathcal R^{\infty}|^{8} \mathcal R^{\infty}=0$, the same as in (1.7) with $a=0$, and hence, instead of $\mathcal R^\infty$ from now on, we can simply write $Q$ in the critical case as in (1.7), namely, in the one dimensional case:

(6.5)\begin{equation} b\, Q + \partial^4_{\xi}Q -|Q|^{8} Q=0. \end{equation}

Solving for $Q$ as we did in Section 2.3.1, we obtain profiles for different values of $b$. Observe that since the critical equation preserves the $L^2$ norm, the mass of the ground states for different values of $b$ is the same, and thus, we simply denote it by $M[Q]$. We obtain (computed for different values of $b$ and also listed in Table 1)

(6.6)\begin{equation} M[Q] = 2.986792978326142, \end{equation}

which is consistent, for example, with the results in [Reference Fibich, Ilan and Papanicolaou25, (41)].

To investigate self-similar blow-up numerically, we fix $a=0$ in (2.1) and take the Gaussian initial data (5.1), $u_0 = A\, e^{-x^2}$ with $A=1.7$. The reason for this choice of the amplitude $A$ is that the mass of the initial data $M[u_0] = A^2 \, \sqrt{\tfrac{\pi}{2}} $ is about $1.1$ times the mass of the solitary wave $M[Q^{(0)}]$ (computed in (6.6)), and thus, blow-up is expected. Since the mass ratio is sufficiently close to 1, it allows us to track the numerical solution for some non-trivial finite time before the amplitude gets too large to be handled numerically, see Figure 23.

Figure 23. Solution to (2.1), $\alpha=8$, $a = 0$, and $u_0=1.7\, e^{-x^2}$. Left: time dependence of the $L^\infty$ norm. Middle: energy conservation $\Delta = \log_{10}|E - E_0|$. Right: Blow up rate and linear fitting on a log scale, the blow-up occurs at $t_* = 0.0210375$, the fitting on a log scale up to the time step $1.5 \times 10^{-10}$ is too delicate to give a conclusive value for the rate parameter.

A note about the blow-up time: an estimated value for the blow-up time is determined experimentally. Then the time interval up to the (numerical) blow-up time is subdivided into several intervals with increasingly smaller time steps as it gets closer to the blow-up. For the last interval the time step $h\approx 10^{-10}$, which is a reasonable implementation in a double precision setting. Such small time steps are needed because of the dependence of the dynamically rescaled time $\tau$ on $t$ in (6.1). The $L^\infty$ norm of the solution is shown on the left of Figure 23. Fitting the $L^\infty$ norm of the numerical solution close to the blow-up core as discussed above, we obtain the numerical blow-up time $t^* = 0.2102004$. In the middle of Figure 23, we track the conservation of the energy during this simulation, namely, the quantity $\delta = \log_{10}|E - E_0|$ (the computation is terminated once the relative energy conservation drops below $10^{-3}$).

6.2.1. Rate of self-similar blow-up

Using the fact that (1.1) is locally well-posed in $H^2$, a similar argument from the classical existence theory as in the critical NLS case (see [Reference Sulem and Sulem46, Thm 5.9]) gives the lower bound on $\Vert \Delta u(t)\Vert_{L^2}$ close to the blow-up time (in the scaling invariance case with $a=0$), see details in [Reference Baruch, Fibich and Mandelbaum7, Theorem 5],

(6.7)\begin{equation} \Vert \Delta u(t)\Vert_{L^2(\mathbb R)}\geq\frac{C(u_0)}{\sqrt{t^*-t}}, \end{equation}

where $C(u_0) \gt 0$ is some constant that depends on the initial data. A further conjecture on the blow-up rate in the critical case (pure quartic dispersion) is given in [Reference Baruch, Fibich and Mandelbaum7, Conjecture 8].

To investigate the blow-up rate, we numerically track the $L^\infty$ norm of the solution, see Figure 23. For that we examine the numerical time dependence of the factor $\lambda(t)$, from the ansatz (6.1), via the time dependence of $\|u(t)\|_{L^\infty}$. As discussed above, the decay of $A^\infty \equiv (\ln L)_\tau \searrow 0$ controls the blow-up rate (as well as the convergence to the blow-up profile). Specifically, in the $L^{2}$-critical case for a stable blow-up an algebraic decrease $L\propto 1/\tau$ is expected, which implies that near the blow-up time, the time dependence is

(6.8)\begin{equation} \lambda(t) \propto (t^{*}-t)^{1/3}, \end{equation}

whereas an exponential decrease $\lambda\propto e^{-\beta \tau}$ is expected for a stable blow-up in the $L^2$-supercritical case, leading to

(6.9)\begin{equation} \lambda(t)\propto (t^{*}-t)^{1/4}. \end{equation}

Using the rescaling (6.1) and $\alpha=8$, we obtain

(6.10)\begin{equation} \|u(t)\|_{L^\infty} = \frac1{\lambda(t)^{1/2}} \|U(\tau)\|_{L^\infty} \cong \frac{1}{\big( (t^*-t)^{1/3}\big)^{1/2}} = \frac1{(t^*-t)^{1/6}}, \end{equation}

and

(6.11)\begin{equation} \|\partial_x^2 u(t)\|_{L^2} = \frac{c}{\lambda(t)^{2}} \cong \frac{1}{\big( (t^*-t)^{1/3}\big)^{2}} = \frac1{(t^*-t)^{2/3}}. \end{equation}

Confirmation towards the rate in (6.10) can be seen on the right of Figure 25 instead, whereas Figure 23 does not provide a conclusive answer. This rate appears to be different from the conjectured rates in [Reference Baruch, Fibich and Mandelbaum7].

6.2.2. Asymptotic profile of self-similar blow-up

Recalling the example of Gaussian data discussed in Figure 23, we show the asymptotic blow-up profile on the left of Figure 24. In our simulations we use the profile $Q=Q^{(0)}$ from (6.5), which is the same as (1.7). We rescale it and fit with the final computational state to check the matching of the asymptotic profile in the case of the finite time blow-up solution. The fitting is done as follows: the maximum of the modulus of the solution to equation (1.1) is divided by the maximum of $Q$ to give according to (1.8) the value of $b^{1/\alpha}$. With this value of $b$ relation (1.8) gives $Q_{b}$ for the computed $Q$.

Figure 24. Blow-up profile and the fitting error for the solution of (2.1), $\alpha=8$, $a=0$, with $u_0=1.7 e^{-x^2}$, evolution of which is shown in Figure 23. Left: profile $\|u(t)\|_{L^\infty}$ (blue) at time $t_m: t^*-t_m = 1.8028\times 10^{-5}$, fitted to a rescaled ground state $Q$ (red) from (6.5). Right: difference on a log scale between the solution and the fitted ground state.

$\scriptstyle\blacklozenge$ ( $a=0$) When $a=0$ (hence, scaling invariance holds), we are able to fit the final computational state with the (numerical) rescaled ground state profile of $Q^{(0)}$, for the initial datum $u_0 = 1.1Q^{(0)}$, we omit the figure for brevity, as it is similar to the next example.

$\scriptstyle\blacklozenge$ ( $a\neq 0$) We take $a=-2$ and compute the corresponding ground state solution $Q^{(-2)}$ (which has no oscillations in this case, since $a \lt -\sqrt 2$) from (2.2). The mass of this ground state is given in Table 1. Then we consider perturbations $u_0 = A\, Q^{(a)}$ for different $A$. In Figure 25 we show the blow-up solution with the initial condition $u_0 = 1.1 Q^{(-2)}$, where we plot the solution and the fitting by the rescaled ground state $Q^{(0)}_{b'}$ at the blow-up time (namely, the last computable time before the blow-up with the step time difference on the order of $10^{-6}$) on the left plot and the error in the log scale on the right. The bottom row shows our fitting of the blow-up rate on the log scale, where we fit $(t^*-t)$ vs. $\|u_{xx}\|_{L^2}$ or $\|u\|_{L^\infty}$ norms close to the blow-up time. We obtain the power for the rate $\|u_{xx}\|_{L^2}$ to be around 0.7, and for $\|u\|_{L^\infty}$ around 0.139, which provides some positive confirmation towards (6.10) and (6.11). (We do stress that it is computationally very challenging to obtain more refined and more accurate rates in a double precision approach.)

Figure 25. Blow-up profile for the solution of (1.1), $\alpha=8$, $a = -2$, with $u_0=1.1\, Q^{(-2)}$ and fitting it with the rescaled ground state $Q^{(0)}$. The solution blows up at $t^* = 0.19177$. Top left: Profile $|u(t)|$ at time $t_m: t^*-t_m = -4.1132\times 10^{-6}$ (blue) fitted to a re-scaled ground state $Q^{(0)}$ (red). Top right: Difference on a log scale between the solution and the fitted ground state. Bottom left: Blow-up of $\|u_{xx}\|_{L^2}$ experimentally fitted to the slope $0.7$. Bottom right: Blow-up of $\|u\|_{L^\infty}$ experimentally fitted to the slope of $0.139$.

We next consider the equation (2.1) with $a=1$ (thus, the ground state $Q^{(1)}$ from (2.2) has more oscillation than the scaling-invariant ground state $Q^{(0)}$ in the pure quartic case, $a=0$) and we study the blow-up behaviour in that case. Note that the mass of the ground state $Q^{(1)}$ is smaller than that in the pure quartic case, $M[Q^{(1)}] \lt M[Q^{(0)}]$, see Table 1.

In Figure 26 instead of a slightly perturbed ground state data (which has similar results), we show the behaviour of the Gaussian initial condition $u_0 = 1.7 e^{-x^2}$, which has mass above $M[Q^{(0)}]$. We show that it blows up in finite time with the profile converging to the rescaled ground state $Q^{(0)}$, see top left plot in Figure 26. The difference between the solution (at the final computed time) and the rescaled soliton $Q^{(0)}$ is shown in the top right plot. We also show the fitted rates in the bottom row, as in the case of $a=-2$, we get similar values.

Figure 26. Blow-up profile and fitting with the rescaled ground state for the solution of (1.1) with $a = 1$ and Gaussian initial data $u(x,0)=1.7 e^{-x^2}$. The solution blows up at $t^* = 0.021933$. Snapshot of the profile at $t_m: t^*-t_m = 1.3\times 10^{-6}$. Top left: Profile $|u(t)|$ at time $t_m$ (blue) fitted to a rescaled ground state $Q^{(0)}$ (red). Top right: Difference on a log scale between the solution and the fitted ground state. Bottom left: Blow-up of $\|u_{xx}\|_{L^2}$ experimentally fitted to the slope $0.65$. Bottom right: Blow-up of $\|u\|_{L^\infty}$ experimentally fitted to the slope of $0.13$.

In Figure 27 we consider hyperbolic secant initial condition $u_0 = A \, \text{sech} \, x$ with the magnitude $A = 1.45$, so that the mass is above $M[Q^{(0)}]$. The plot on the left shows that the solution blows up in finite time with the profile converging to the rescaled ground state $Q^{(0)}$.

Figure 27. Blow up profile and fitting with the rescaled ground state for the solution of (1.1) with $a = 1$, $\alpha=8,$ and initial data $u(x,0)=1.45 \,\text{sech}(x)$. Top left: Profile $|u(t)|$ at time $t_m$ (blue) fitted to a rescaled ground state $Q^{(0)}$ (red). Top right: Difference on a log scale between the solution and the fitted ground state.

7.1. Rate and Profile

To confirm the blow up rate, we use the rescaling (6.1) to deduce

(7.1)\begin{equation} \|u(t)\|_{L^\infty} = \frac1{\lambda(t)^{4/\alpha}} \|U(\tau)\|_{L^\infty} \cong \frac{1}{\big( (t^*-t)^{1/4}\big)^{4/\alpha}} = \frac1{(t^*-t)^{1/\alpha}}. \end{equation}

In particular, in the power case $\alpha=10$, the rate is $1/10$, which we confirm in our fittings in Figure 29. We note that in this supercritical case, it is easier to check the rate numerically, as it converges to a specific profile very fast (unlike the critical case, so these features are similar to the standard NLS equation).

As far as the profile is concerned, one can see in the top left plot of Figure 29 that the profile is different from the ground state in this case $Q^{(0)}$ (with $\alpha=10$). This is consistent with the standard NLS equation.

Figure 29. Supercritical case $\alpha = 10$. Blow up profile and fitting with the rescaled ground state for the solution of (1.1) with $a = 1$ and Gaussian initial data $u(x,0)=1.7 e^{-x^2}$. The solution blows up at $t^* = 0.00375322$. Snapshot at $t_m: t^*-t_m = 2.0\times 10^{-8}$. Clockwise from top left: Profile $|u(t)|$ at time $t$ (blue) fitted to a re-scaled ground state (red, given by the asymptotic solution for $a = 0$); Difference on a log scale between the solution and the fitted ground state; Blow-up of $\|u_{xx}\|_{L^2}$ experimentally fitted to rate $0.59$; Blow-up of $\|u\|_{L^\infty}$ experimentally fitted to rate of $0.1$. The case with negative $a$ does not provide qualitatively new results.

To check further on rates and profiles of the solution, we considered super-Gaussian and sech initial data, and obtained similar results on the rate and profiles.

The principal difference of the supercritical case to the critical case is that the asymptotic profile and the rescaled ground state $Q^{(0)}$ do not coincide, while the rate is easier to track and fit, and it is consistent with the theoretical prediction as in (7.1).