format: text

text:
# Transverse instability and its long-term development for solitary waves of the (2+1)-Boussinesq equation

## I Introduction

One of the fundamental ways that a solitary wave traveling in one space dimension generates a two space dimensional pattern is through transverse instability. A transverse instability of a line solitary wave is associated with a class of perturbations traveling in a direction transverse to the basic direction of propagation. In addition to establishing the existence of transverse instability, a major question is what implications this instability means for the long-term behaviour of the system: does it settle into a new two-space-dimensional pattern, or collapse ? In this paper we study this sequence of questions for the canonical Boussinesq equation in two space dimensions

$$u_{tt}=(f(u)+\epsilon u_{xx})_{xx}+\sigma u_{yy},$$
(1)
where $`\epsilon =\pm 1`$ and $`\sigma =\pm 1`$. In general, $`f(u)`$ can be any smooth function, but the canonical form of the Boussinesq equation has the form

$$f(u)=D(u^2u)\text{with}D=\pm 1.$$

When $`D=1`$, $`\epsilon =1`$ and $`\sigma =1`$ this equation was derived by Johnson J  to describe the propagation of gravity waves on the surface of water, in particular the head-on collision of oblique waves, and it was derived by Breizman and Malkin BM  in the context of Langmuir waves.

In the absence of the transverse variation (i.e, $`u_y=0`$) and for $`\epsilon =1`$, $`D=1`$ this equation reduces to the so-called ”good” Boussinesq equation, which is well-posed, and for which $`\text{sech}^2`$-solutions exist for any $`c`$ with $`|c|<1`$. These waves are stable when $`\frac{1}{2}<|c|<1`$ BS . For the case $`|c|<\frac{1}{2}`$ it was shown by computer-assisted simulation of the leading term in the Taylor expansion of the Evans function that there is an unstable eigenvalue AS . This result was generalized to include solitary waves with nonzero tails, and rigorously proved using the symplectic Evans matrix in BD1 .

Transverse instability of solitary waves has been widely studied since the seminal work of Zakharov VEZ  on the nonlinear Schrödinger equation and the work of Kadomtsev & Petviashvili KP  on transverse instability of the Korteweg-de Vries soliton. Since then, transverse instability of solitary waves has been investigated for a wide range of models; examples include the nonlinear Schrödinger (NLS) equation and related equations KRZ ; LS ; WAS , Kadomtsev-Petviashvili equation APS ; ARa ; IRS ; KPe , the Zakharov-Kuznetsov equation AR ; B1 ; LS , and water waves B3 . A review of transverse instability for NLS and other related models can be found in Kivshar & Pelinovsky KPe .

In this paper, we will first use a geometric condition as derived in B1  to get an explicit criterion for small transverse wavenumber instability. For this we use the multi-symplectic formulation of (1) in an essential way. To get detailed information for all transverse wavenumbers we compute explicitly the Evans function for the (2+1)-dimensional Boussinesq model linearized about a larger family of line solitary waves (allowing the state at infinity to be nonzero). Plots of the dependence of the growth rate on the transverse wavenumber are presented.

The post-instability behaviour of the nonlinear problem is studied using direct numerical simulation. The numerical evidence confirms the analytic results and suggests that the post-instability in the nonlinear system leads to collapse in all cases. A multi-symplectic pseudospectral discretization BR  is used as a basis for the numerical simulations. The numerical scheme is applied to the full two-dimensional PDE and we observe transverse modulation and further development of the longitudinal and transverse instabilities, resulting in the collapse of the initial line solitary waves. In the parameter region where the analytic criterion indicates that the solitary wave state is longitudinally stable but transversely unstable, simulations support the analytic results and provide insight into the long-term development of this instability.

## II Multi-symplectifying the equations

The Boussinesq system has a range of geometric structures. Firstly, we record the Lagrangian and Hamiltonian structures. Let $`u=\varphi _{xx}`$, then the system is Lagrangian with

$$L=\left[\frac{1}{2}\varphi _{xt}^2+F(\varphi _{xx})+\frac{1}{2}\epsilon \varphi _{xxx}^2+\frac{1}{2}\sigma \varphi _{xy}^2\right]𝑑x𝑑y𝑑t,$$

where $`F()`$ is any function satisfying $`F^{}()=f()`$.

The Boussinesq equation can be represented as a Hamiltonian system in a number of ways (e.g. T ). For example, let

$$H=\left[F(u)\frac{1}{2}\epsilon u_x^2+\frac{1}{2}\mathrm{\Phi }_x^2+\frac{1}{2}\sigma w_y^2+\gamma (uw_x)\right]𝑑x𝑑y,$$

where $`\gamma `$ is a Lagrange multiplier associated with the constraint $`u=w_x`$. With Hamiltonian variables $`(\mathrm{\Phi },u,w,\gamma )`$ the governing equations take the form

$$\begin{array}{ccc}u_t\hfill & =\hfill & \frac{\delta H}{\delta \mathrm{\Phi }}=\mathrm{\Phi }_{xx}\hfill \\ \mathrm{\Phi }_t\hfill & =\hfill & \frac{\delta H}{\delta u}=f(u)+\epsilon u_{xx}+\gamma ,\hfill \\ 0\hfill & =\hfill & \frac{\delta H}{\delta w}=\gamma _x\sigma w_{yy},\hfill \\ 0\hfill & =\hfill & \frac{\delta H}{\delta \gamma }=uw_x.\hfill \end{array}$$
(2)
However, the most interesting form of (1) for the present purposes is the multi-symplectic formulation which can be represented in the canonical form BD3

$$𝐌Z_t+𝐊Z_x+𝐋Z_y=S(Z)Z^6,$$
(3)
where

$$Z=\left(\begin{array}{c}q_1\\ q_2\\ q_3\\ p_1\\ p_2\\ p_3\end{array}\right),\text{ }𝐌=\left(\begin{array}{cccccc}0& 1& 0& 0& 0& 0\\ 1& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0\end{array}\right),\text{with}u(x,y,t)=q_1(x,y,t),$$

$$𝐊=\left(\begin{array}{cccccc}0& 0& 0& 1& 0& 0\\ 0& 0& 0& 0& 1& 0\\ 0& 0& 0& 0& 0& 1\\ 1& 0& 0& 0& 0& 0\\ 0& 1& 0& 0& 0& 0\\ 0& 0& 1& 0& 0& 0\end{array}\right),\text{ }𝐋=\left(\begin{array}{cccccc}0& 0& 1& 0& 0& 0\\ 0& 0& 0& 0& 0& 0\\ 1& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0\end{array}\right),$$

$$S(Z)=F(q_1)\frac{1}{2\epsilon }p_1^2+\frac{1}{2}p_2^2\frac{\sigma }{2}p_3^2.$$

Using $`q_1=u`$ it is straightforward to show that this system is a reformulation of (1).

## III Geometric criterion for transverse instability

An advantage of the multi-symplectic formulation is that there is a geometric condition which is easy to verify for transverse instability of line solitary waves B1 .

Consider the well-known basic family of solitary waves of (1) of the form

$$Z(x,y,t)=\widehat{Z}(\theta ;c,l),\text{ }\theta =xct+ly+\theta _0,$$
(4)
obtained by taking the first component to be a $`\mathrm{sech}^2`$ wave,

$$u(\theta ;c,l)=𝐞_1,\widehat{Z}(\theta ;c,l)=A(c,l)\text{sech}^2(B(c,l)\theta ),$$
(5)
with

$$B(c,l)=\frac{1}{2}\sqrt{\epsilon (D+c^2\sigma l^2)},A(c,l)=6\frac{\epsilon }{D}B^2.$$

Existence of the solitary wave clearly requires $`\epsilon (D+c^2\sigma l^2)>0`$. The other components of $`\widehat{Z}`$ are easily obtained from (5) and the multi-symplectic equations (3).

For the linear stability analysis, let $`Z(x,y,t)=\widehat{Z}(\theta ;c,l)+\mathrm{}[U(\theta ;\lambda ,k)e^{\lambda t+iky}]`$, substitute this into (3) and linearize. Then, if the resulting linear equation has square-integrable solutions $`U(\theta ;\lambda ,k)`$ with $`\mathrm{}(\lambda )>0`$ and $`k`$, we call the basic solitary wave state $`\widehat{Z}(\theta ;c,l)`$ transversely unstable. Assuming that $`\widehat{Z}_\theta `$ is the only square integrable element in the kernel of the linearization operator $`=D^2S(Z)\left[𝐊c𝐌+l𝐋\right]\frac{d}{d\theta }`$, we have the following geometric condition of transverse instability for small $`\lambda `$ and $`k`$. Suppose

$$\mathrm{\Delta }=\left|\begin{array}{cc}𝒜_c& 𝒜_l\\ _c& _l\end{array}\right|>0,\text{where}\{\begin{array}{c}𝒜=\frac{1}{2}_{\mathrm{}}^{\mathrm{}}𝐌\widehat{Z}_\theta ,\widehat{Z}𝑑\theta ,\hfill \\ =\frac{1}{2}_{\mathrm{}}^{\mathrm{}}𝐋\widehat{Z}_\theta ,\widehat{Z}𝑑\theta ,\hfill \end{array}$$
(6)
Then the basic solitary wave $`\widehat{Z}(\theta ;c,l)`$ of (1) is linearly transverse unstable B1 ; B3 .

Using the above definitions of the multi-symplectic matrices $`𝐌`$ and $`𝐋`$, we obtain

$$𝒜=\frac{1}{2}_{\mathrm{}}^{\mathrm{}}\left(q_1\frac{d}{d\theta }q_2q_2\frac{d}{d\theta }q_1\right)𝑑\theta =c_{\mathrm{}}^{\mathrm{}}q_1^2𝑑\theta =cK,$$
(7)
$$=\frac{1}{2}_{\mathrm{}}^{\mathrm{}}\left(q_1\frac{d}{d\theta }q_3q_3\frac{d}{d\theta }q_1\right)𝑑\theta =\sigma l_{\mathrm{}}^{\mathrm{}}q_1^2𝑑\theta =\sigma lK,$$
(8)
where

$$K=_{\mathrm{}}^{\mathrm{}}u^2𝑑\theta =\frac{4}{3}\frac{A^2}{B}=\frac{6\epsilon }{D^2}\left(\sigma l^2c^2D\right)\sqrt{\frac{\sigma l^2c^2D}{\epsilon }}.$$

Substitution of (7) and (8) in (6) yields:

$$\text{sgn}\mathrm{\Delta }=\text{sgn}\left[\frac{\sigma }{c}𝒜\left(𝒜_c+\frac{l}{c}𝒜_l\right)\right]=\text{sgn}\left[\sigma \left(𝒜_c+\frac{l}{c}𝒜_l\right)\right]=$$

$$=\text{sgn}\left[\sigma \left(K+c\frac{}{c}K+l\frac{}{l}K\right)\right]=\text{sgn}\left[\sigma \left(\sigma l^2c^2D\right)\left(4\sigma l^24c^2D\right)\right].$$
(9)
Since the condition for transverse instability requires $`\mathrm{\Delta }>0`$, we have the following result: Suppose

$$\epsilon \sigma \left(4\sigma l^24c^2D\right)>0,$$
(10)
then the basic solitary wave $`\widehat{Z}(\theta ;c,l)`$ is linearly transversely unstable.

The multi-symplectic formulation also provides an expression for the linear growth rate of the instability $`\lambda `$ as a function of the transverse wavenumber $`k`$ for long-wave perturbations B1 :

$$\lambda =\frac{\sqrt{𝒜_c_l𝒜_l_c}}{|𝒜_c|}k+𝒪(k^2)=\frac{\sqrt{\sigma (4\sigma l^24c^2+1)(\sigma l^2c^2+1)}}{4c^21\sigma l^2}k+𝒪(k^2).$$
(11)
This provides the growth rate for $`k`$ small. In the next section, the Evans function will be constructed in order to determine the growth rate for all transverse wavenumbers $`k`$.

In the remainder of this section, we apply the condition (10) for various parameter values.

For the ”good” Boussinesq equation from BS  with $`\epsilon =1`$ and $`D=1`$ the existence and transverse instability requirements are

$$\sigma l^2c^2+1>0\text{and}\sigma \left(\sigma l^2c^2+\frac{1}{4}\right)>0.$$
(12)
respectively. Combining these conditions leads to the following system of inequalities for $`c`$ and $`l`$ when $`\sigma >0`$

$$\frac{1}{4}+\sigma l^2<c^2<1+\sigma l^2,$$
(13)
and for $`\sigma <0`$

$$c^2<\frac{1}{4}+\sigma l^2$$
(14)
These inequalities define the regions in $`(c,l)`$ parameter plane, where the basic solitary wave exits and is linearly transversely unstable, and these regions are presented in Figure 1.

One can do a similar analysis for Johnson’s equation J , where $`\sigma =1`$, $`\epsilon =1`$ and $`D=1`$. The existence requirement is $`l^2<c^21`$ and the instability condition is $`l^2>c^2\frac{1}{4}`$. This result is inconclusive for two reasons. First, the two regions do not overlap so the geometric condition does not predict instability for any parameter values. Secondly, when $`\epsilon =+1`$ the equation is ill-posed as an evolution equation (this can be seen at the linear level where the dispersion relation predicts instability as the wavenumber goes to infinity), and so the question of long time stability is irrelevant.

## IV Evans function analysis of the transverse instability

In this section we use the Evans function formalism in order to analyze the linear transverse stability problem for the Boussinesq model (1) for all values of the transverse wavenumber. We restrict attention to the parameter values of most interest: $`\epsilon =1`$ and $`D=1`$ associated with the “good” Boussinesq, although we put no restriction on $`\sigma `$ (but keeping in mind that $`\sigma =+1`$ is the most interesting case).

However, the class of solitary waves will be enlarged. Namely, we include solitary waves bi-asymptotic to a nontrivial state at infinity, specifically,

$$U(\theta )=U_{\mathrm{}}+6\delta ^2\text{sech}^2\left(\delta \theta \right),\theta =xct+ly,$$
(15)
where

$$\delta =\frac{1}{2}\sqrt{\sqrt{1+4a}c^2+\sigma l^2}\text{and}U_{\mathrm{}}=\frac{1}{2}(1c^2+\sigma l^2)2\delta ^2=\frac{2a}{1+\sqrt{1+4a}}.$$
(16)
The value of the parameter $`a`$ is constrained only by existence of the square root: $`1+4a(c^2\sigma l^2)^2`$.

Here we will not use any geometric structure (although it might be interesting to look more closely in this direction) and so work directly with (1). Let

$$u(x,y,t)=U(\theta )+\mathrm{}\left(\stackrel{~}{u}(\theta )\mathrm{exp}\left[iky+\lambda t\right]\right).$$
(17)
By substituting this expression in (1) and linearising, one obtains the following equation for the complex function $`\stackrel{~}{u}(\theta )`$

$$\stackrel{~}{u}_{\theta \theta \theta \theta }+2(U\stackrel{~}{u})_{\theta \theta }(1c^2+\sigma l^2)\stackrel{~}{u}_{\theta \theta }2(c\lambda +i\sigma kl)\stackrel{~}{u}_\theta +(\lambda ^2+\sigma k^2)\stackrel{~}{u}=0.$$
(18)
After the change of variable $`\stackrel{~}{x}=\delta \theta `$, substitution of the explicit expression for $`U`$ from (15), and dropping the tildes, equation (18) reduces to

$$u_{xxxx}4\left[(13\text{sech}^2x)u\right]_{xx}\gamma u_x+\beta u=0,$$
(19)
where

$$\gamma =\frac{2(c\lambda +i\sigma kl)}{\delta ^3}\text{ and }\beta =\frac{\lambda ^2+\sigma k^2}{\delta ^4}.$$
(20)
To obtain explicit solutions of this equation, we note that by taking $`u=\varphi _{xx}`$ and $`v=13\text{sech}^2x`$ in (19), and integrating twice the equation simplifies to

$$\varphi _{xxxx}4v\varphi _{xx}\gamma \varphi _x+\beta \varphi =0.$$
(21)
Solutions of this equation can be readily found in a manner similar to that in B  (see also APS ). First we note that in the limit $`x\pm \mathrm{}`$, equation (21) reduces to

$$\varphi _{xxxx}4\varphi _{xx}\gamma \varphi _x+\beta \varphi =0.$$
(22)
Substituting now $`\varphi =e^{\mu x}\widehat{\varphi }`$, one can see that $`\mu `$ satisfies the quartic equation

$$\mu ^44\mu ^2\gamma \mu +\beta =0.$$
(23)
Quartics of this form have been analyzed in BD1  (see equation (10.9) there), and when $`\mathrm{}(\beta )>0`$ there are two roots with positive real part and two roots with negative real part. Therefore, the space of solutions decaying as $`x+\mathrm{}`$ is two-dimensional, as is the the space of solutions decaying as $`x\mathrm{}`$.

If the four roots $`\mu _j,\text{ }j=1,..,4`$ of the equation (23) are distinct, the corresponding solutions of (21) are given by

$$\varphi _j(x)=e^{\mu _jx}h_j(x),$$
(24)
with

$$h_j(x)=(4\mu _j^3+8\mu _j\gamma )12\mu _j^2\text{tanh}x.$$
(25)
The case of multiple roots can be handled similarly B . Solutions of the original equation (19) are found by substituting $`u(x)=\varphi (x)_{xx}`$, and the other components of the vector $`𝐯(x)`$ can be obtained by the differentiating the expression for $`u(x)`$.

Localised solutions of the linearised problem exist if one can match the solutions decaying as $`x\mathrm{}`$ with the solutions decaying as $`x\mathrm{}`$. This can be determined by finding the zeros of the so-called Evans function which correspond to the eigenvalues of the linearised problem. To define the Evans function, we write the equation (19) as a first-order system

$$𝐯_x=A(x)𝐯,\text{ }𝐯=\left(\begin{array}{c}u\hfill \\ u_x\hfill \\ u_{xx}\hfill \\ u_{xxx}\hfill \end{array}\right),\text{ }A(x)=\left(\begin{array}{cccc}0& 1& 0& 0\\ 0& 0& 1& 0\\ 0& 0& 0& 1\\ \beta +4v_{xx}& \gamma +8v_x& 4v& 0\end{array}\right)$$
(26)
with $`v=13\text{sech}^2x`$.

Since the trace of the matrix $`A(x)`$ vanishes, the Evans function can be defined as $`E(\lambda ,k)=𝐯_1(x)𝐯_2(x)𝐯_3(x)𝐯_4(x)`$ AGJ . An alternative expression for the Evans function can be derived by using the adjoint system as shown in BD2 . The adjoint system of (26) has the form:

$$𝐰_x=A(x)^{}𝐰,\text{ }𝐰=\left(\begin{array}{c}w_1\hfill \\ w_2\hfill \\ w_3\hfill \\ w_4\hfill \end{array}\right),$$
(27)
where $`A(x)^{}`$ denotes the Hermitian conjugate of $`A`$ $`(A(x)^{}=\overline{A(x)}^T)`$. The equation for $`w_4`$ turns out to be

$$(w_4)_{xxxx}4v(w_4)_{xx}+\overline{\gamma }(w_4)_x+\overline{\beta }(w_4)=0.$$
(28)
This equation is equivalent to (21) up to the change of variables: $`xx`$, $`\gamma \overline{\gamma }`$, $`\beta \overline{\beta }`$, and therefore its solutions can be obtained from (24) by changing $`x`$ for $`x`$ and conjugating them:

$$(w_4)_j=e^{\mu _j^{}x}\overline{h_j(x)},$$
(29)
with $`h_j(x)`$ defined in (25). Other components of the vector $`𝐰(x)`$ can be obtained from (27).

Let $`\mu _1`$ and $`\mu _2`$ be the two roots of the equation (23) with negative real part, and let $`𝐯_j(x)`$ and $`𝐰_j(x)`$, $`j=1,2`$ be the corresponding solution vectors of the linearised (respectively, adjoint) system. Since the matrix $`A(x)`$ in (26) is traceless, we can define the Evans function for the system (26) as follows BD2 :

$$E(\lambda ,k)=\left|\begin{array}{cc}𝐰_1(0),𝐯_1(0)& 𝐰_1(0),𝐯_2(0)\\ 𝐰_2(0),𝐯_1(0)& 𝐰_2(0),𝐯_2(0)\end{array}\right|,$$
(30)
where $`,`$ denotes the complex inner product in $`^4`$. To obtain a unique definition of the Evans function, the scaling $`lim_x\mathrm{}e^{2\mu _jx}𝐰_j(x),𝐯_j(x)=1`$ is used. This normalises the eigenvectors and the adjoint eigenvectors of $`A^{\mathrm{}}=lim_{x\pm \mathrm{}}A(x)`$.

After some lengthy algebra and introducing the scaling, which enforces the asymptotic limit $`E(\lambda ,k)1`$ as $`\lambda \mathrm{}`$, the final expression for the Evans function can be obtained, which we do not present here since it is lengthy (the expression for the Evans function as well as the calculations of the instability growth rate can be downloaded as a Maple-file from the website Web ).

Zeros of the Evans function $`E(\lambda ,k)`$ correspond to the bounded solutions of the linearised stability problem with the wavenumber $`k`$ and the growth rate $`\mathrm{}(\lambda )`$. The leading order terms (in $`k`$ and $`\lambda `$) in the Evans function are in complete agreement with the results of the geometric condition of §3. Note that, since the construction here is based on a basic solitary wave with a nontrivial state at infinity, it is suggestive that the geometric condition B1  extends to such waves.

We illustrate the dependence of the Evans function on the wavespeed and transverse wavenumber in Figure 2. In the left graph, the transverse wavenumber is set to zero, to compare with known results on longitudinal instability. The graph is in complete agreement with known results (e.g. BS ; BD1 ) that the solitary wave is stable for $`\frac{1}{2}<c1`$ and unstable for $`0c<\frac{1}{2}`$. In the right-hand graph in Figure 2 we present the plot of the growth rate $`\mathrm{}(\lambda )`$ as a function of the transverse wavenumber. Note that waves of the good Boussinesq which are longitudinally stable are transverse unstable. Note also that there is a cut-off wavenumber, similar to other cases of transverse instability, such as in the Zakharov-Kuznetsov equation AR .

## V Post-instability simulations

In this section we perform a simulation of the PDE (1) using the multi-symplectic spectral discretization proposed in BR  and applied there to Zakharov-Kuznetsov and shallow-water equations.

The (2+1)-dimensional Boussinesq equation is considered with $`\epsilon =D=1`$ on a finite domain $`(x,y)=[0,L]\times [0,L]`$ with $`L>0`$ some constant, and periodic boundary conditions on both spatial variables. We choose a spatial mesh-size as $`\mathrm{\Delta }x=\mathrm{\Delta }y\mathrm{\Delta }m=L/2N`$ and introduce the discrete two-dimensional Fourier transform defined as

$$U_{kl}=\frac{1}{\sqrt{2N}}\underset{i,j=1}{\overset{2N}{}}u_{ij}e^{\theta _k(i1)\mathrm{\Delta }m\theta _l(l1)\mathrm{\Delta }m},$$

where

$$\theta _k=i\frac{2\pi (k1)}{L}$$

and $`u_{ij}u(m_{ij}),m_{ij}=(i1)\mathrm{\Delta }x+(j1)\mathrm{\Delta }y`$ (cf. F ). Fourier spectral discretization of the (2+1)-dimensional Boussinesq equation yields

$$_{tt}U_{kl}=\overline{\theta }_k^2\left[\epsilon \overline{\theta }_k^2U_{kl}+_{kl}\overline{F}(𝐔)\right]+\sigma \overline{\theta }_l^2U_{kl},$$
(31)
where $`\overline{\theta }_k`$ are the entries of the diagonal matrix defined by the relations

$$\overline{\theta }_k=\theta _k,\text{ for }k=1,\mathrm{},N,$$

$`\overline{\theta }_{N+1}=0`$, and

$$\overline{\theta }_k=\theta _{2Nk+2},\text{ for }k=N+2,\mathrm{},2N$$

which follow from the periodicity of the discrete Fourier transform F , and $`\overline{F}(𝐔)`$ denotes the Fourier transform of the anti-derivative of the function $`f(u)`$ in (1). The same result would be obtained if one applied the spectral discretization to the multi-symplectic formulation (3), as it was done for the Zakharov-Kuznetsov equation in BR .

For the second-order time derivative we used the central difference approximation (time step was chosen to be $`\mathrm{\Delta }t=0.01`$ in all the simulations):

$$_{tt}U_{kl}=\frac{U_{kl}^{n+1}2U_{kl}^n+U_{kl}^{n1}}{\mathrm{\Delta }t^2}.$$
(32)
One should note that the only valid test of this scheme can be done for the “good” Boussinesq equation with $`\sigma >0`$. For $`\sigma <0`$ in the case of the “good” Boussinesq equation, an initial profile independent of $`x`$ would result in a solution which could grow “faster than exponential” because for large transverse wavenumbers, growth rate of the initial data has no upper bound (ill-posedness).

To test the algorithm, we first used it to confirm the results for the dynamics of the one-dimensional solitary waves. The initial profile was taken to be of the form

$$u(x)=\frac{3}{2}\left(1c^2\right)\text{sech}^2\left[\frac{1}{2}\left(1c^2\right)\left(x\frac{L}{2}\right)\right]+\xi (x),$$
(33)
where $`\xi (x)`$ is a small random perturbation. The results are presented in Figures 3 and 4. For $`c=\frac{1}{4}`$ the solitary wave solution is linearly unstable as reported in AS ; BD1 , and the development of this linear instability is shown in Figure 3a). In the case $`c=\frac{3}{4}`$ the numerical results confirm the stability of the solitary wave (see Figure 3b) ). The simulations were run on an interval of the length $`L=256`$ with $`2N=512`$. As a numerical check, the total energy was monitored, and it was found to be well behaved till near the collapse when the significant errors occur, as illustrated in Figure 4.

For the two-dimensional simulations we took an initial profile in the form of the line solitary wave uniform in $`y`$

$$u(x,y,0)=\frac{3}{2}\left(1c^2\right)\text{sech}^2\left[\frac{1}{2}\left(1c^2\right)\left(x\frac{L}{2}\right)\right]+\xi (x,y),$$
(34)
where $`\xi (x,y)`$ is a small random perturbation (in this case $`l=0`$). The length of the square box was chosen to be $`L=64`$ with the number of Fourier modes $`2N=128`$. In the case $`c=\frac{1}{4}`$ the solitary wave (34) is linearly unstable in longitudinal direction as is known from the stability analysis of the 1D equation. In Figure 5b) we can see this instability developing in a similar way as in the 1D case. One can also note in this Figure the development of the stable transverse modulation. Wave collapse in this case is shown in Figure 6a), with the plot of energy as a function of time in Figure 6b). When $`c=\frac{3}{4}`$, the solitary wave is longitudinally stable but transversely unstable, and the development of this instability is presented in Figures 7 and 8. We note that at the initial stage of the evolution there is a transverse modulation developing while the amplitude of the wave is gradually growing Figure 7b), then the instability prevails leading finally to the collapse of the wave Figure 8a). Note that this collapse is clearly a two-dimensional effect, since it does not happen uniformly in the $`y`$-direction. The energy proves to be conserved rather well during the simulations (see Figure 8b) ), although the energy deviates substantially as the wave approaches the stage of collapse.

## VI Concluding remarks

We have considered the transverse instability of line solitary wave solutions of the (2+1)-dimensional Boussinesq equation. Using the multi-symplectic formulation of the system, we derived a geometric condition for this instability for small transverse wavenumbers. With an Evans function approach, the linearised stability equation was analyzed, and this allowed to obtain the dependence of the instability growth rate for all transverse wavenumbers. Numerical simulations support the analytical results about transverse and longitudinal instabilities and demonstrate the development of those instabilities and subsequent wave collapse.

We conclude with an open problem. While analytic theories for collapse of solitary waves for the Boussinesq equation in one space dimension exist T , it is an interesting open problem to develop an analytical technique for predicting collapse for the case of two space dimensions, e.g. a generalization of the virial theorem or the result of T  for example, and moreover, to determine if transverse instability for (1) leads to collapse for all parameter values.

## VII Acknowledgements

The authors would like to thank Sebastian Reich for advice on the numerics and Georg Gottwald for helpful discussions.