Multi-breathers and high order rogue waves for the nonlinear Schr\"odinger equation on the elliptic function background

TL;DR

This paper constructs and analyzes multi-breather and rogue wave solutions of the focusing nonlinear Schrödinger equation on elliptic function backgrounds, revealing their complex dynamics and stability properties.

Contribution

It introduces explicit multi-breather and rogue wave solutions on elliptic backgrounds using Darboux transformation and theta functions, advancing understanding of their behaviors and instabilities.

Findings

Multi-breather solutions expressed via theta functions.

Asymptotic analysis of breather dynamics at large times.

Exact rogue wave solutions derived from breather limits.

Abstract

We construct the multi-breather solutions of the focusing nonlinear Schr\"odinger equation (NLSE) on the background of elliptic functions by the Darboux transformation, and express them in terms of the determinant of theta functions. The dynamics of the breathers in the presence of various kinds of backgrounds such as dn, cn, and non-trivial phase-modulating elliptic solutions are presented, and their behaviors dependent on the effect of backgrounds are elucidated. We also determine the asymptotic behaviors for the multi-breather solutions with different velocities in the limit $t\to\pm\infty$, where the solution in the neighborhood of each breather tends to the simple one-breather solution. Furthermore, we exactly solve the linearized NLSE using the squared eigenfunction and determine the unstable spectra for elliptic function background. By using them, the Akhmediev breathers arising from these modulational instabilities are plotted and their dynamics are revealed. Finally, we provide the rogue-wave and higher-order rogue-wave solutions by taking the special limit of the breather solutions at branch points and the generalized Darboux transformation. The resulting dynamics of the rogue waves involves rich phenomena: depending on the choice of the background and possessing different velocities relative to the background. We also provide an example of the multi- and higher-order rogue wave solution.

Equations

$$\mathrm{i}{\psi }_t+\frac{1}{2}{\psi }_{xx}+|\psi {|}^2\psi =0,$$

$$\begin{array}{rl}{\Phi }_x& =\mathbf{U}(\mathbf{Q};\lambda )\Phi ,\mathbf{U}(\mathbf{Q};\lambda )=-\mathrm{i}\lambda {\sigma }_3+\mathrm{i}\mathbf{Q},\\ {\Phi }_t& =\mathbf{V}(\mathbf{Q};\lambda )\Phi ,\mathbf{V}(\mathbf{Q};\lambda )=-\mathrm{i}{\lambda }^2{\sigma }_3+\lambda \mathrm{i}\mathbf{Q}+{\mathbf{V}}_0,{\mathbf{V}}_0=\frac{1}{2}{\sigma }_3(\mathrm{i}{\mathbf{Q}}^2-{\mathbf{Q}}_x),\end{array}$$

$${\sigma }_3=\left[\begin{array}{cc}1& 0\\ 0& -1\end{array}\right],\mathbf{Q}=\left[\begin{array}{cc}0& \psi \\ {\psi }^{\ast }& 0\end{array}\right].$$

$${\sigma }_2=\left[\begin{array}{cc}0& -\mathrm{i}\\ \mathrm{i}& 0\end{array}\right].$$

$$\begin{array}{rl}{\mathbf{L}}_x& =[\mathbf{U},\mathbf{L}],\\ {\mathbf{L}}_t& =[\mathbf{V},\mathbf{L}].\end{array}$$

$$\mathrm{i}{p}_t+\frac{1}{2}{p}_{xx}+2|\psi {|}^{2}p+{\psi }^2{p}^{\ast }=0,$$

$$\mathbf{L}(x,t;\lambda )=\left[\begin{array}{cc}-\mathrm{i}f(x,t;\lambda )& g(x,t;\lambda )\\ h(x,t;\lambda )& \mathrm{i}f(x,t;\lambda )\end{array}\right],$$

$$\begin{array}{rl}f(x,t;\lambda )& =-\frac{\mathrm{i}}{2}\left[{\varphi }_1(x,t;\lambda ){\phi }_2(x,t;\lambda )+{\phi }_1(x,t;\lambda ){\varphi }_2(x,t;\lambda )\right],\\ g(x,t;\lambda )& ={\varphi }_1(x,t;\lambda ){\varphi }_2(x,t;\lambda ),h(x,t;\lambda )=-{\phi }_1(x,t;\lambda ){\phi }_2(x,t;\lambda ),\end{array}$$

$$\Phi (x,t;\lambda )=\left[\begin{array}{cc}{\varphi }_1(x,t;\lambda )& {\varphi }_2(x,t;\lambda )\\ {\phi }_1(x,t;\lambda )& {\phi }_2(x,t;\lambda )\end{array}\right],$$

$$\begin{array}{rl}{f}_x=& {\psi }^{\ast }g-\psi h,g_x=-2(\psi f+\mathrm{i}\lambda g),h_x=2({\psi }^{\ast }f+\mathrm{i}\lambda h),\\ f_t=& \left(\lambda {\psi }^{\ast }-\frac{\mathrm{i}}{2}{\psi }_x^{\ast }\right)g-\left(\psi \lambda +\frac{\mathrm{i}}{2}{\psi }_x\right)h,g_t=-2\left[\left(\lambda \psi +\frac{\mathrm{i}}{2}{\psi }_x\right)f+\left(\mathrm{i}{\lambda }^2-\frac{\mathrm{i}}{2}|\psi {|}^2\right)g\right],\\ h_t=& 2\left[\left(\lambda {\psi }^{\ast }-\frac{\mathrm{i}}{2}{\psi }_x^{\ast }\right)f+\left(\mathrm{i}{\lambda }^2-\frac{\mathrm{i}}{2}|\psi {|}^2\right)h\right].\end{array}$$

$$-\frac{\mathrm{i}}{2}[\frac{1}{2}{g}_{xx}+\psi ({\psi }^{\ast }g-\psi h)]=\frac{1}{2}(\mathrm{i}{\psi }_x+2\lambda \psi )f+\mathrm{i}{\lambda }^{2}g,$$

$$\mathrm{i}{g}_t+\frac{1}{2}{g}_{xx}+2|\psi {|}^{2}g-{\psi }^{2}h=0.$$

$$-\mathrm{i}{h}_t+\frac{1}{2}{h}_{xx}+2|\psi {|}^{2}h-{\psi }^{\ast 2}g=0.$$

$$\mathbf{L}={\mathbf{L}}_0(x,t){\lambda }^2+{\mathbf{L}}_1(x,t)\lambda +{\mathbf{L}}_2(x,t).$$

$$\begin{array}{rl}& {\sigma }_3{\mathbf{L}}_0^{\mathrm{off}}=0,\\ & {\mathbf{L}}_{0,x}+2\mathrm{i}{\sigma }_3{\mathbf{L}}_1^{\mathrm{off}}-\mathrm{i}[\mathbf{Q},{\mathbf{L}}_0]=0,\\ & {\mathbf{L}}_{1,x}+2\mathrm{i}{\sigma }_3{\mathbf{L}}_2^{\mathrm{off}}-\mathrm{i}[\mathbf{Q},{\mathbf{L}}_1]=0,\\ & {\mathbf{L}}_{2,x}-\mathrm{i}[\mathbf{Q},{\mathbf{L}}_2]=0,\end{array}$$

$$\begin{array}{rl}& {\sigma }_3{\mathbf{L}}_0^{\mathrm{off}}=0,\\ & 2\mathrm{i}{\sigma }_3{\mathbf{L}}_1^{\mathrm{off}}-\mathrm{i}[\mathbf{Q},{\mathbf{L}}_0]=0,\\ & {\mathbf{L}}_{0,t}+2\mathrm{i}{\sigma }_3{\mathbf{L}}_2^{\mathrm{off}}-\mathrm{i}[\mathbf{Q},{\mathbf{L}}_1]-[{\mathbf{V}}_0,{\mathbf{L}}_0]=0,\\ & {\mathbf{L}}_{1,t}-\mathrm{i}[\mathbf{Q},{\mathbf{L}}_2]-[{\mathbf{V}}_0,{\mathbf{L}}_1]=0,\\ & {\mathbf{L}}_{2,t}-[{\mathbf{V}}_0,{\mathbf{L}}_2]=0,\end{array}$$

$${\mathbf{Q}}_t-2\mathrm{i}{\alpha }_2{\sigma }_3\mathbf{Q}+{\alpha }_1{\mathbf{Q}}_x=0.$$

$$\begin{array}{rl}\mathbf{L}(x,t;\lambda )=& \exp \left[{\int }_{(x_0,t_0)}^{(x,t_0)}\mathbf{U}(x^{\prime },t^{\prime };\lambda )\mathrm{d}{x}^{\prime }+{\int }_{(x,t_0)}^{(x,t)}\mathbf{V}(x^{\prime },t^{\prime };\lambda )\mathrm{d}{t}^{\prime }\right]\mathbf{M}(\lambda )\\ & \exp \left[-{\int }_{(x_0,t_0)}^{(x,t_0)}\mathbf{U}(x^{\prime },t^{\prime };\lambda )\mathrm{d}{x}^{\prime }-{\int }_{(x,t_0)}^{(x,t)}\mathbf{V}(x^{\prime },t^{\prime };\lambda )\mathrm{d}{t}^{\prime }\right]\end{array}$$

$$\begin{array}{r}\det (\mathbf{L}(x,t;\lambda ))=f(x,t;\lambda {)}^2-g(x,t;\lambda )h(x,t;\lambda )=P(\lambda ),\\ P(\lambda ):=\prod _{i=1}^4(\lambda -{\lambda }_i)={\lambda }^4-s_1{\lambda }^3+s_2{\lambda }^2-s_3\lambda +s_4,\end{array}$$

$$f(x,t;\lambda )={\lambda }^2+{\alpha }_1\lambda +({\alpha }_2-\frac{1}{2}|\psi {|}^2),g(x,t;\lambda )=-h^{\ast }(x,t;{\lambda }^{\ast })=\mathrm{i}\lambda \psi -\frac{1}{2}{\psi }_x+\mathrm{i}{\alpha }_1\psi =\mathrm{i}\psi (\lambda -\mu )$$

$$\begin{array}{rl}2{\alpha }_1& =-s_1,{\alpha }_1^2+2{\alpha }_2=s_2,\\ 2{\alpha }_1({\alpha }_2-\frac{1}{2}\nu )-\nu (\mu +{\mu }^{\ast })& =-s_3,({\alpha }_2-\frac{1}{2}\nu {)}^2+\nu \mu {\mu }^{\ast }=s_4.\end{array}$$

$$\mu +{\mu }^{\ast }=\frac{s_1}{2}-\frac{q}{\nu },\mu {\mu }^{\ast }=-\frac{1}{\nu }({\nu }^2-2p\nu +p^2-4s_4),$$

$$\mu ,{\mu }^{\ast }=\frac{s_1}{4}-\frac{q\pm \mathrm{i}\sqrt{-R(\nu )}}{2\nu },$$

$$R(\nu )={\nu }^3+(s_1^2/4-2p){\nu }^2-(s_{1}q+4s_4)\nu +q^2.$$

$$[{\lambda }^2+({\alpha }_1+\mathrm{i}\sqrt{{\nu }_1})\lambda +({\alpha }_2-{\nu }_1/2)-\mathrm{i}\sqrt{{\nu }_1}\mu ][{\lambda }^2+({\alpha }_1-\mathrm{i}\sqrt{{\nu }_1})\lambda +({\alpha }_2-{\nu }_1/2)+\mathrm{i}\sqrt{{\nu }_1}\mu ].$$

$${\alpha }_1+\mathrm{i}\sqrt{{\nu }_1}=-{\lambda }_1-{\lambda }_3,{\alpha }_1-\mathrm{i}\sqrt{{\nu }_1}=-{\lambda }_2-{\lambda }_4.$$

$${\nu }_1=-\frac{1}{4}({\lambda }_1+{\lambda }_3-{\lambda }_2-{\lambda }_4{)}^2.$$

$${\nu }_2=-\frac{1}{4}({\lambda }_1+{\lambda }_4-{\lambda }_2-{\lambda }_3{)}^2,{\nu }_3=-\frac{1}{4}({\lambda }_1+{\lambda }_2-{\lambda }_3-{\lambda }_4{)}^2.$$

$${\mu }_x=-2\mathrm{i}f(\mu )=-2\mathrm{i}\sqrt{P(\mu )},$$

$${\mu }_t=-2\mathrm{i}\left(\mu +\frac{\mathrm{i}}{2}(\ln \psi {)}_x\right)\sqrt{P(\mu )}=-\mathrm{i}{s}_1\sqrt{P(\mu )}.$$

Full text

1. Introduction

As a universal model in the description of the propagation of a quasi-monochromatic wave in a weakly nonlinear medium, the nonlinear Schrödinger equation (NLSE) has played a central role in mathematical physics and integrable partial differential equations [18]. It is relevant in numerous applications of natural sciences such as water waves, nonlinear optics, plasmas, Bose-Einstein condensates. For the zero background, a variety of methods either for the analysis of the Cauchy problem or the systematic construction of the multi-soliton solutions has been well established for a long time. A nonzero with plane wave background can be unstable under the perturbation, and this modulation instability (MI) is considered as the main cause for the formation of rogue waves (RWs) in nature. The integrable nature of the NLSE allows one to construct a family of exact solutions corresponding to this nonzero plane wave background by different methods such as the generalized Darboux transformation method, the dressing method, the inverse scattering method and Hirota bilinear method. Since the 70s, there are lots of works on the integrable NLSE model; giving a complete account of the reference on the relative subject is out of the scope of this introduction.

Recent study by Grinevich and Santini [19] (equations (48)-(49) ) shows finite-gap soliton solutions of genius-1 or genius-2 may be developed from a Cauchy problem of the NLSE with a very small perturbation. For the genius-1 case, a nonzero background expressed by elliptic function is formed before the rogue wave solution or the Akhmediev breather solution is formed. In this sense, the rogue wave or Akhmediev breather solutions with complicated backgrounds such as elliptic functions are more natural and are expected to occur in experiments. Actually, Akhmediev et al. have studied this kind of solution in a simple setting [24]. Due to the mathematical difficulty, the study of constructing soliton solutions with elliptic function background is much less and nontrivial. In this paper, one aim is to construct general breathers and rogue waves of the NLSE with the elliptic function background systematically.

From the theoretical point of view, the systematic construction of the multi-soliton solutions on the elliptic function background can be broadly classified into two categories. The first method is to take the limit in the higher-genus solutions constructed by the algebro-geometric method [3] or the Riemann-Hilbert method. The other way is using the inverse scattering method or the Darboux transformation [7, 9, 15, 20, 21, 30] to construct the breathers by the linear algebra. However, both methods have different difficulties. For the first way, it is hard to analyze the different limitation and deal with lots of elliptic integration, while for the second way, the obtained formula by the linear algebra is easy, but it is difficult to analyze the properties for the determinant.

Firstly, we review several known results on this aspect. In the section 4.5 of monograph [3], the authors considered the degenerate finite-gap solution—multi-phase modulations of the cnoidal wave for the defocusing NLSE by taking the limit and used the complex coordinate transformation $x\to\mathrm{i}x$ to derive multi Akhmediev breathers on the ${\rm dn}$ function background [3]. Recently, Shin considered the single soliton dynamics in phase-modulated lattices by Darboux transformation by modified squared wavefunction approach [34]. A compact formula for the single breather solution was not, however, given explicitly, and hence the shift of the crest and phase was not analyzed. The exact solitons on the elliptic function background are constructed by the eigenfunction symmetry method [13, 32]. Breathers, rogue waves and higher order rogue waves on the dn and cn background are obtained by numerically solving the Lax pair equation and Darboux transformation [2, 24]. The solitons for the defocusing NLSE on the elliptic function background was obtained by the inverse scattering method [36], where the seed eigenfunction of the Lax pair was obtained using the result of [37], and the solution was expressed by the determinant of theta function. The rogue waves for the focusing NLSE on the cn and dn background are obtained by combining the Darboux transformation with nonlinearization of the Lax pair [12]. Meanwhile, there are some relevant works on the rogue waves on the higher genus solutions [4], in which the analytical criterion to generalized rogue waves are obtained by the Riemann-Hilbert method. The bounded ultra-elliptic solutions for the defocusing NLSE are considered by the effective integration method [38]. The various initial data for the semiclassical NLSE are related with the rogue-wave phenomena [33, 5, 17].

Even though there are some contributions on this subject, there is not systematic theory to construct the breathers and rogue waves on the elliptic function background. What is more important is that there is no method to analyze the solutions generated by the Darboux transformation. In this work, we would like to address these issues. The key step to generate the solution by the Darboux transformation is solving the Lax pair equation with the seed solution. Once it can be obtained, plugging the vector solutions of the Lax pair into the Bäcklund transformation induced by the Darboux transformation, we get the new solutions. Since the formula of the Bäcklund transformation is algebraic, the new solutions are obtained exactly, and furthermore, this operation can be iterated as many times as desired. Therefore, the family of the multi-fold and higher-order solutions can be obtained recursively. The closed-form multi-breather solutions can be easily plotted by just setting the parameters, and no numerical computation is necessary.

However, since the solutions generated by the Darboux transformation involve lots of free non-trivial parameters, if we want to understand their effect precisely, the formulas found from the Bäcklund transformation are not always useful for observation of the general laws of the dynamics. Thus we need to reduce the formula in a proper form, so that the dynamics of solutions can be analyzed. In this work, we use the addition formulas of theta functions to rewrite the expressions of the multi-soliton solutions. In this way, the general dynamics can be analyzed directly through the asymptotic analysis of the determinant formula.

It is well known that the breather solutions in the focusing NLSE involve the special and important limiting cases, which have a close relation to the modulational instability. The elliptic-function solutions for the focusing NLSE are also known to be unstable [16]. It is therefore natural to construct the exact solutions to describe the unstable process. Thus, the second aim in this work is to construct the Akhmediev breathers on the elliptic function background, which is helpful to understand the mechanism of modulational instability.

Finally, for the rogue waves, it is also related with the modulational instability with infinite period. The comprehensive construction for higher order rogue waves with the closed expressions has been open so far. In this work, we prove that the formulas of higher order rogue waves can be represented as a rational functions of $x$, $t$, the Jacobi elliptic functions ${\rm sn}$, ${\rm cn}$, ${\rm dn}$ and the second type of elliptic integral.

The rest of this paper is organized as the follows. In section 2 we solve the Lax pair with the elliptic function by combining the algebro-geometric method and the modified squared wavefunction approach. As a byproduct of the modified squared wavefunction approach, we also solve the linearized NLSE directly, which is useful for the analysis of the modulational instability. In section 3, we construct the breathers and multi-breathers on the elliptic function background. Furthermore, with the aid of the addition formula for the theta function, we rewrite them in a determinant form. Their dynamics are also investigated for various backgrounds. In section 4, we determine the unstable spectra for the linearized NLSE, and analyze the corresponding Akhmediev breathers on the elliptic-function background, which is related with the modulational instability. In section 5, we provide the rogue waves and high order rogue waves. We further discuss the effect of the choice of backgrounds to the dynamics of rogue waves. Section 6 is devoted to the conclusion and discussion.

2. The elliptic function solution of NLSE and its Lax pair

In this section, we derive the elliptic function solution for NLSE firstly and then the corresponding solutions for the Lax pair are solved with the aid of algebraic geometry method. The scheme is given in details in the monograph [3]. Here we just focus on the genus one case, in which we can use the elliptic functions and integrals to represent these solutions instead of the Abel map and Riemann theta function. In the following, we give the details to construct the solutions.

The NLSE reads

[TABLE]

which admits the following Lax pair:

[TABLE]

where

[TABLE]

In other words, the compatibility condition $\Phi_{xt}=\Phi_{tx}$ or the zero curvature condition $\mathbf{U}_{t}-\mathbf{V}_{x}+[\mathbf{U},\mathbf{V}]=0$ yields the NLSE (2.1). It is readily to see the symmetry for the potential function $\mathbf{U}$ and $\mathbf{V}$: $\mathbf{U}^{{\dagger}}(\mathbf{Q};\lambda^{*})=-\mathbf{U}(\mathbf{Q};\lambda)$ and $\mathbf{V}^{{\dagger}}(\mathbf{Q};\lambda^{*})=-\mathbf{V}(\mathbf{Q};\lambda)$, where † denotes the hermitian conjugate. Based on the above symmetry, if $\Phi(x,t;\lambda)$ is a matrix solution for Lax pair (2.2), then $\Phi^{{\dagger}}(x,t;\lambda^{*})$ is a matrix solution for the adjoint Lax pair of (2.2). As a general symmetry, if $\Phi(x,t;\lambda)$ satisfies the Lax pair (2.2), then $({\rm i}\sigma_{2})\Phi(x,t;\lambda)^{\mathrm{T}}(-{\rm i}\sigma_{2})$ satisfies the adjoint Lax pair of (2.2), where $\sigma_{2}$ is a Pauli matrix:

[TABLE]

We define a new matrix function $\mathbf{L}(x,t;\lambda)=-\frac{1}{2}\Phi(x,t;\lambda)\sigma_{3}({\rm i}\sigma_{2})\Phi(x,t;\lambda)^{\mathrm{T}}(-{\rm i}\sigma_{2})$ which satisfies the following stationary zero curvature equations:

[TABLE]

It should be stressed that the compatibility condition for the above stationary zero curvature conditions also give the NLSE (2.1) by the compatibility condition of the Lax pair and Jacobi identity.

2.1. Linearized NLSE and its solution

The above stationary zero curvature equation (2.3) can be used to construct the solutions of linearized nonlinear Schrödinger equation:

[TABLE]

where $p(x,t)$ is a perturbed function. To give the exact solution from the stationary zero curvature equations, we rewrite the matrix function $\mathbf{L}(x,t;\lambda)$ with in the form:

[TABLE]

where

[TABLE]

and

[TABLE]

$[\phi_{i},\varphi_{i}]^{\rm T}$ are two vector solutions for Lax pair (2.2). Then we write the stationary zero curvature equation in more details as

[TABLE]

From the above equations, we obtain

[TABLE]

then we find that

[TABLE]

Similarly, we have

[TABLE]

If $h=-g^{*}$, then the above equation turns into the linearized NLSE (2.4). But in general, $h$ does not equals to $-g^{*}$ in the above setting. To satisfy this requirement, we use the symmetry to construct the combination of solution. Since $[\phi_{i}(\lambda),\varphi_{i}(\lambda)]^{\rm T}$, $i=1,2$ is a solution for Lax pair (2.2), then $[-\varphi_{i}^{*}(\lambda),\phi_{i}^{*}(\lambda)]^{\rm T}$ is a solution for Lax pair (2.2) by replacing $\lambda$ with $\lambda^{*}$. Thus, we find that the function $g(\lambda^{*})=\varphi_{1}^{*}(x,t;\lambda)\varphi_{2}^{*}(x,t;\lambda)$ and $h(\lambda^{*})=-\phi_{1}^{*}(x,t;\lambda)\phi_{2}^{*}(x,t;\lambda)$ satisfy the linearized equation (2.8) and (2.9) respectively. Then we conclude that the exact solution for the linearized NLSE (2.4) is $p=\phi_{1}(x,t;\lambda)\phi_{2}(x,t;\lambda)+\varphi_{1}^{*}(x,t;\lambda)\varphi_{2}^{*}(x,t;\lambda).$ This solution is useful to analyze the spectral stability of the background solution.

2.2. Elliptic function solution for NLSE and the corresponding solution of Lax pair

Starting from the stationary zero curvature equations, we set the solution $\mathbf{L}$ as the polynomial of $\lambda$. For the genus one solution, we assume that the matrix function $\mathbf{L}$ is a quadratic polynomials of $\lambda$:

[TABLE]

Plugging the above equation (2.10) into equation (2.3) and comparing the coefficient of $\lambda$, we obtain that: for the $x$-equation

[TABLE]

which infers that $\mathbf{L}=-{\rm i}\sigma_{3}\alpha_{0}\lambda^{2}+(-{\rm i}\sigma_{3}\alpha_{1}+{\rm i}\alpha_{0}\mathbf{Q})\lambda+\alpha_{0}\mathbf{V}_{0}+{\rm i}\alpha_{1}\mathbf{Q}+\alpha_{2}(-{\rm i}\sigma_{3})$; for the $t$-equation

[TABLE]

from which we obtain that $\alpha_{0,t}=\alpha_{1,t}=\alpha_{2,t}=0$ and

[TABLE]

On the other hand, suppose that the matrix solution $\mathbf{L}(x,t;\lambda)$ can be solved formally with the form:

[TABLE]

which is non-degenerate if the constant matrix $\mathbf{M}(\lambda)$ is non-singular. Thus, combining with the above ansatz (2.10), the determinant of the matrix $\mathbf{L}(x,t;\lambda)$ can take the following form

[TABLE]

where $\lambda_{i}$s are zeros of the polynomial which characterize the genus one solution. Here we see that $\alpha_{0}=1$. By the symmetry relationship and the existence and uniqueness of ordinary differential equation, we know that the solution satisfies $\mathbf{L}(x,t;\lambda)=-\mathbf{L}^{{\dagger}}(x,t;\lambda^{*}).$ Moreover, through the determinant relation (2.14), we obtain that $\prod_{i=1}^{4}(\lambda-\lambda_{i})=\prod_{i=1}^{4}(\lambda-\lambda_{i}^{*}).$ In this way, we can suppose $\lambda_{j}=\lambda_{j+2}^{*}$, $j=1,2.$ In summary, the functions $f(x,t;\lambda)$, $g(x,t;\lambda)$ and $h(x,t;\lambda)$ are given by:

[TABLE]

where we introduced $\mu(x,t):=-\frac{{\rm i}}{2}(\ln\psi)_{x}-\alpha_{1}$ in the last expression. Now we give the relationship between the coefficients $s_{1\leq i\leq 4}$ and functions $\psi(x,t)$, $\mu(x,t)$, $\nu=|\psi|^{2}$, $\alpha_{2}.$ By comparing the equation (2.14), we find that

[TABLE]

Moreover, from the latter two equations of (2.16), we obtain the algebraic relation between $\mu$, $\mu^{*}$ and $\nu$:

[TABLE]

where $q=\frac{s_{1}}{2}(s_{2}-\frac{s_{1}^{2}}{4})-s_{3}$, $p=s_{2}-\frac{s_{1}^{2}}{4}.$ Then $\mu$, $\mu^{*}$ can be solved exactly:

[TABLE]

where

[TABLE]

If there exists a root $\nu_{1}$ such that $R(\nu_{1})=0$, then the polynomials (2.14) $\det(\mathbf{L}(x,t;\lambda))$ can be decomposed into

[TABLE]

Based on the relation between the roots and coefficients in equations (2.14) and (2.20), we know that

[TABLE]

Furthermore, we have the relation

[TABLE]

In a similar procedure, for the other two roots $\nu_{2}$ and $\nu_{3}$, we have

[TABLE]

Here, we take advantage of the effective integration technique in [23]. In what follows, we will see that the parameters $\nu_{1}$, $\nu_{2}$ and $\nu_{3}$ will determine the dynamics of solution $\psi$. Because of the above relationship between $\nu_{i}$ $(i=1,2,3)$ and $\lambda_{j}$ $(j=1,2,3,4)$, the parameters $\lambda_{j}$ also determine the dynamics of solution.

Next we derive the Dubrovin type equation for the root $\mu(x,t)$. Since $g_{x}=-2(\psi f+{\rm i}\lambda g)$ and $g={\rm i}\psi(\lambda-\mu)$ [equations (2.7) and (2.15)], by taking $\lambda=\mu$, it follows that

[TABLE]

and

[TABLE]

Furthermore, from equations (2.13) and (2.18) and the relation $\mu-\mu^{*}=-\frac{{\rm i}}{2}\frac{\nu_{x}}{\nu}$, we know that

[TABLE]

Introduce the notations $\lambda_{1,3}=\alpha\pm{\rm i}\gamma$, $\lambda_{2,4}=\beta\pm{\rm i}\delta$, then it follows from formulas (2.21) and (2.22) that

[TABLE]

We can see that $\nu_{1}\leq 0\leq\nu_{2}<\nu_{3}$. By simple mathematical analysis on the differential equation (2.25), we find that $\nu$ is a periodic function on the interval $[\nu_{2},\nu_{3}]$, which is indeed written by an elliptic function:

[TABLE]

where $m=(\nu_{3}-\nu_{2})/(\nu_{3}-\nu_{1}).$ Moreover, from equations (2.13), (2.18), and (2.25), we know that

[TABLE]

When $\nu_{2}>0$, the function $\psi(x,t)$ can integrate exactly

[TABLE]

When $\nu_{2}=0$, since there is a singularity located at $\sqrt{\nu_{3}-\nu_{1}}(x+s_{1}t)=0\mod 2K$ in equations (2.28), ($K\equiv K(m)$ is the complete elliptic integral of the first kind), the function $\psi(x,t)$ can be represented as

[TABLE]

Remark 1**.**

We consider the Galilean transformation to conceal the parameter $s_{1}$, using the symmetry of the NLSE $\psi(x,t)\rightarrow\psi(x-\frac{s_{1}}{2}t,t){\rm e}^{{\rm i}\frac{\theta}{2}},$ $\theta=s_{1}(x-\frac{s_{1}}{4}t)$. We also apply a similar transformation on the Lax pair to achieve the symmetry $(\lambda,\Phi)\rightarrow(\widehat{\lambda},\widehat{\Phi})$, where $\widehat{\lambda}=\lambda-\frac{s_{1}}{4}$, $\widehat{\Phi}={\rm e}^{\frac{\rm i}{2}\theta\sigma_{3}}\Phi(x-\frac{s_{1}}{2}t,t)$. Thus we have the correpsponce between the spectral curve for $\lambda$ and $\widehat{\lambda}$, where the new spectral curve is given by $\widehat{P}(\widehat{\lambda}):=\widehat{\lambda}^{4}-\hat{s}_{1}\widehat{\lambda}^{3}+\hat{s}_{2}\widehat{\lambda}^{2}-\hat{s}_{3}\widehat{\lambda}+\hat{s}_{4}$ with $\hat{s}_{1}=0$, $\hat{s}_{2}=s_{2}-\frac{3}{8}s_{1}^{2}$, $\hat{s}_{3}=\frac{s_{1}^{3}}{8}-\frac{s_{1}s_{2}}{2}+s_{3}$, and $\hat{s}_{4}=-\frac{3s_{1}^{4}}{256}+\frac{s_{1}^{2}s_{2}}{16}-\frac{s_{1}s_{3}}{4}+s_{4}$. The variables $\nu_{1\leq i\leq 3}$ are invariant under the transformation. So is the parameter $q=\sqrt{-\nu_{1}\nu_{2}\nu_{3}}$. Thus, under the Galilean transformation, the solution (2.29) can convert into a stationary solution which is better to analyze the dynamics without loss the generality:

[TABLE]

So in the following, we set the parameter $s_{1}=0$ based on the above analysis.

Remark 2**.**

We give the Weierstrass elliptic function representation. Set $\nu(x)=4(\frac{s_{2}}{6}-\wp(2x+\omega_{3}))$, then $\wp(2x+\omega_{3})$ satisfies $[\frac{{\rm d}\wp(2x+\omega_{3})}{{\rm d}x}]^{2}=16\times\left[4\wp(2x+\omega_{3})-g_{2}\wp(2x+\omega_{3})-g_{3}\right]=64(\wp(2x+\omega_{3})-e_{1})(\wp(2x+\omega_{3})-e_{2})(\wp(2x+\omega_{3})-e_{3})$, where $g_{2}=\frac{s_{2}^{2}}{12}+s_{4}$ and $g_{3}=\frac{s_{2}^{3}}{216}+\frac{s_{3}^{2}}{16}-\frac{s_{2}s_{4}}{6}$, $e_{1}=\frac{1}{12}(\nu_{2}+\nu_{3}-2\nu_{1})$, $e_{2}=\frac{1}{12}(\nu_{1}+\nu_{3}-2\nu_{2})$, $e_{3}=\frac{1}{12}(\nu_{1}+\nu_{2}-2\nu_{3})$. And $\wp(\omega_{i})=e_{i}$, $i=1,2,3,$ $\omega_{1}$ and $\omega_{3}$ are the half-period of the $\wp$ function. This is consistent with the Jacobi elliptic function representation since it takes the maximum at $x=0$.

We now proceed to find the solutions of the corresponding Lax pair. Let us introduce the notation $y^{2}=\det(\mathbf{L}(x,t;\lambda))$, i.e. $\det({\rm i}y\mathbb{I}-\mathbf{L}(x,t;\lambda))=0$. So the kernel of matrix ${\rm i}y\mathbb{I}-\mathbf{L}(x,t;\lambda)$ are $K_{1}(x,t;\lambda)=(1,r_{1}(x,t;\lambda))^{\rm T}$ and $K_{2}(x,t;\lambda)=(1,r_{2}(x,t;\lambda))^{\rm T},$ where

[TABLE]

Suppose there exists a fundamental matrix solution for the Lax pair (2.2): $\widehat{\Phi}(x,t;\lambda)$ with the initial data $\widehat{\Phi}(0,0;\lambda)=\mathbb{I}$, and the solution matrix $\widehat{\Phi}(x,t;\lambda)\mathbf{L}(0,0;\lambda)$ is also a matrix solution for the Lax pair (2.2). On the other hand, we know that $\mathbf{L}(x,t;\lambda)\widehat{\Phi}(x,t;\lambda)$ is also a matrix solution for the Lax pair (2.2). By the uniqueness and existence of ordinary differential equation, we arrive at $\widehat{\Phi}(x,t;\lambda)\mathbf{L}(0,0;\lambda)=\mathbf{L}(x,t;\lambda)\widehat{\Phi}(x,t;\lambda).$ Now we again consider the solution $\Phi(x,t;\lambda)\equiv(\Psi_{1}(x,t;\lambda),\Psi_{2}(x,t;\lambda))$. Since the Lax pair is linear, we can normalized the first component of the vector solutions $\Psi_{1}(0,0;\lambda),$ $\Psi_{2}(0,0;\lambda)$ to be unity, and set their second components as $\psi_{1}(0,0;\lambda)$ and $\psi_{2}(0,0;\lambda)$, respectively. In other words,

[TABLE]

Based on this, we have $[{\rm i}y\mathbb{I}-\mathbf{L}(x,t;\lambda)]\Phi(x,t;\lambda)=0,$ which yields $r_{i}(x,t;\lambda)=\frac{\varphi_{i}(x,t;\lambda)}{\phi_{i}(x,t;\lambda)}$, $i=1,2.$

Since $f(x,t;\lambda)$ and $g(x,t;\lambda)$ are known, then the functions $\psi_{i}(x,t;\lambda)$, $i=1,2$ can be exactly obtained. The left unknown variables for Lax pair are $\phi_{i}(x,t;\lambda)$, $i=1,2.$ Indeed, based on the previous known variables, $\phi_{i}(x,t;\lambda)$, $i=1,2$ can be derived from the integration:

[TABLE]

Inserting equation (2.15) into above equation (2.32), we have

[TABLE]

where $\beta_{1}=2(\lambda^{2}+\frac{s_{2}}{2}-y)$, $C_{1}=2\lambda\beta_{1}-s_{3}.$ By using a similar integration, we obtain

[TABLE]

where $\beta_{2}=2(\lambda^{2}+\frac{s_{2}}{2}+y)$ and $C_{2}=2\lambda\beta_{2}-s_{3}.$ We also have $\varphi_{i}=\phi_{i}r_{i},$ $i=1,2.$ Thus, the solution for the Lax pair with the elliptic potential function are determined uniquely.

Remark 3**.**

Though the expressions (2.33) and (2.34) seem to have singularities on the points $\nu(x)-\beta_{i}=0$, $i=1,2$, they in fact have no singularity. Here we check it for, e.g., $\phi_{1}$. Let $\wp(\omega_{\pm})=-\lambda^{2}/2-s_{2}/12\pm\sqrt{P}/2$, then $C_{1}=\mp 4{\rm i}\wp^{\prime}(\omega_{+})$ and $C_{2}=\mp 4{\rm i}\wp^{\prime}(\omega_{-})$ since $[\wp^{\prime}(\omega_{\pm})]^{2}=4\wp(\omega_{\pm})-g_{2}\wp(\omega_{\pm})-g_{3}.$ Now we can expand at the neighborhood of $(\omega_{+}-\omega_{3})/2$: $\nu(x)-\beta_{1}=4\left[\wp(\omega_{+})-\wp(2x+\omega_{3})\right]=\mp 8\wp^{\prime}(\omega_{+})(x-(\omega_{+}-\omega_{3})/2)+O((x-(\omega_{+}-\omega_{3})/2)^{2})$. It follows that $\exp(\theta_{1})=O([(x-(\omega_{+}-\omega_{3})/2]^{\mp 1/2})$. This proves that $\phi_{1}$ has no singularity.

To look for a fundamental matrix solution with better symmetry, we give the integral representation for $\varphi_{i}$, $i=1,2$:

[TABLE]

where

[TABLE]

and $D_{2}={\rm i}\sqrt{\frac{\nu(0)-\beta_{1}}{\nu(0)-\beta_{2}}}.$ By ignoring the constant factors of vector solutions, we obtain the following theorem:

Theorem 1**.**

One fundamental solution for Lax pair (2.2) with the elliptic functions (2.29) or (2.30) can be represented as

[TABLE]

where $\theta_{1}$, $\theta_{2}$ are given in equations (2.33), (2.34).

The above formulas involve the incomplete elliptic integral of the third kind. A direct way is using the integral representation. In what follows, we use the theta function to represent the integration.

The fundamental solutions (2.36) can be normalized at $(x,t)=(0,0)$:

[TABLE]

with $\Phi^{\rm N}(0,0;\lambda)=\mathbb{I},$ which can be proved by the fact that it is an entire function of $\lambda$ by the Proposition 2.1 in [8].

2.3. Theta function representation for solutions

Let $z_{0}=K+2{\rm i}l$ such that $l\in[-K^{\prime}/2,K^{\prime}/2]$ and $\alpha=\sqrt{\nu_{3}-\nu_{1}}$, $m=\frac{\nu_{2}-\nu_{1}}{\nu_{3}-\nu_{1}}$, and we introduce the parameterization:

[TABLE]

Thus $q={\rm i}m\alpha^{3}{\rm sn}(K+2{\rm i}l|m){\rm cn}(K+2{\rm i}l|m){\rm dn}(K+2{\rm i}l|m)$, $\nu(x)=m\alpha^{2}({\rm sn}^{2}(K+2{\rm i}l|m)-{\rm sn}^{2}(\alpha x|m))$. Hereafter, for simplicity we omit the modulus parameter $m$. Rewriting the formula (B.4),

[TABLE]

and using it, we obtain that

[TABLE]

On the other hand, on account of addition formula for theta function (see Appendix B or [26]), we obtain

[TABLE]

Then we obtain the theta function representation:

[TABLE]

From the formula $Z(u)=Z(u+2K)$ and $Z(-u)=-Z(u)$, we know that $Z(K+2{\rm i}l)\in{\rm i}\mathbb{R}.$

Furthermore, we derive the theta-function representation for the solutions of the Lax pair. Firstly, we need to parameterize the spectral curve $y^{2}=\prod_{i=1}^{4}(\lambda-\lambda_{i})$ with the aid of equation (2.23): $\lambda(z)=\mu({\rm i}(z-l)/\alpha)$ and $y(z)=\frac{\alpha}{2}\frac{{\rm d}\mu({\rm i}(z-l)/\alpha)}{{\rm d}z}$, where

[TABLE]

and $k=\sqrt{m},$ $k^{\prime}=\sqrt{1-m}.$ So $s_{2}=\lambda_{1I}^{2}+\lambda_{2I}^{2}-2\lambda_{R}^{2}=\frac{1}{2}(\nu_{1}+\nu_{2}+\nu_{3}).$ We use the addition formula for the Jacobi elliptic function to reduce $\mu({\rm i}(z-l)/\alpha)$ with a better representation:

[TABLE]

where we used the addition formula for Jacobi elliptic function:

[TABLE]

If $l=\pm K^{\prime}/2$, the above formula can be replaced by $\lambda(z)=\frac{\alpha\mathrm{i}}{2}k^{2}{\rm sn}({\rm i}(z-l))\,{\rm cn}({\rm i}(z-l))/{\rm dn}({\rm i}(z-l)).$ For the parameter $y(z)$, we can rewrite it as

[TABLE]

where $C=K+\mathrm{i}K^{\prime}.$ Moreover, by addition formula of the Jacobi elliptic function, we obtain a simple formula of $\lambda^{2}(z)$:

[TABLE]

Furthermore, combining the above result, we have

[TABLE]

It follows that

[TABLE]

Moreover, we have

[TABLE]

and

[TABLE]

On account of the above equations, through the similar calculation as (2.40) and (2.41), we obtain

[TABLE]

By the addition formula for theta functions (see Appendix B), we obtain

[TABLE]

Similarly, we have

[TABLE]

Thus

[TABLE]

Meanwhile, we have

[TABLE]

and

[TABLE]

Moreover, we have

[TABLE]

where we used the following shift formula in the second equality:

[TABLE]

In a similar procedure as above, we have

[TABLE]

Together with the equations (2.50), (2.56) and (2.57), we conclude that

Theorem 2**.**

The vector solutions can be represented with the theta function form:

[TABLE]

and

[TABLE]

In what follows, we would like use the different versions of vector solutions to construct the breathers and rogue waves on the elliptic function background by the Darboux transformation. For the breather solution, we use the theta function, which is helpful to give asymptotic analysis on that. For the rogue wave and high order rogue waves, it is convenient to construct the exact solutions by the Jacobi cn, sn and dn function to represent them.

3. Multi-breather solutions

In this section, we consider the multi-breather solution by the multi-fold Darboux matrix (see Appendix A). This step is just by inserting the solution into the formula of Bäcklund transformation. However, if we would like to analyze the dynamics of the solution, the expressions should be simplified into a proper form.

We consider the general $N$- breather solution on the elliptic function background. Let us choose a special solution using (2.58) and (2.59)

[TABLE]

where $\alpha_{j}$ is a complex parameter, z_{j}^{\prime}=-\mathrm{i}(K+\mathrm{i}K^{\prime})-z_{j},\ $\omega(x,t)=\alpha\,Z(K+2\mathrm{i}l)\,x-\mathrm{i}s_{2}\,t$, and

[TABLE]

By the symmetry relationship, we arrive at

[TABLE]

where $z_{i}^{*}\,{}^{\prime}=-\mathrm{i}(K+\mathrm{i}K^{\prime})-z_{i}^{*}$. Consequently the general breather solution can be expressed with the Darboux transformation:

[TABLE]

where $\psi$ is given by (2.42) and

[TABLE]

Let us rewrite (3.2) in a simpler form using the Sherman-Morrison-Woodbury-type matrix identity

[TABLE]

Firstly, the numerator of (3.3) is

[TABLE]

By the addition formula (B.8) and shift formula (B.9), we obtain that

[TABLE]

Using the addition formula (B.8), we have

[TABLE]

Thus the $(i,j)$-component of $\mathbf{M}$ is given by

[TABLE]

where

[TABLE]

On the other hand, we have

[TABLE]

We can simplify $\Phi_{1}(z_{j})\Phi_{2}^{*}(z_{i})-\psi\frac{1}{\alpha\mathrm{i}}\frac{\vartheta_{3}\vartheta_{4}(\frac{\alpha x}{2K})}{\vartheta_{2}\vartheta_{4}}\Delta_{j,i}$, which is necessary to calculate the numerator in the formula (3.4), using the following expressions found from the addition formula (B.8)

[TABLE]

and

[TABLE]

and

[TABLE]

and

[TABLE]

In summary, we have the following theorem:

Theorem 3**.**

The multi-breather solutions can be represented as the following determinant form:

[TABLE]

where

[TABLE]

and

[TABLE]

In what follows, we analyze the interaction between breathers with different velocity. Introduce the notation $\beta_{i}=\alpha{\rm Re}(Z(\mathrm{i}(z_{i}^{\prime}-l))-Z(\mathrm{i}(z_{i}-l)))$ and $v_{i}=-\frac{2{\rm Im}(y(z_{i}))}{\beta_{i}}.$ Assume that $\beta_{i}>0$ and $v_{1}<v_{2}<\cdots<v_{k}<\cdots<v_{n}.$ Firstly, we split the $(x,t)$ plane into $2n$ different pieces of regions $R_{i}^{\pm}$ exterior to the region of interaction (we approximate region as $\{(x,t)|x^{2}+t^{2}>r\}$, $r$ is big enough) by the $n$ lines $l_{i}:=\{x=v_{i}\,t+\gamma_{i,\pm}\}$, $i=1,2,\cdots,n$ (Fig. 1).

Rewrite the expression of $\Delta_{i,j}$:

[TABLE]

and

[TABLE]

Then the determinants $\det(\Delta)$ and $\det(\mathbf{H})$ in breather solution can be rewritten as the following matrix form:

[TABLE]

where

[TABLE]

and

[TABLE]

with $\xi_{k}=\beta_{k}(x-v_{k}t)+\mathrm{i}(\alpha{\rm Im}(Z(\mathrm{i}(z_{k}^{\prime}-l))-Z(\mathrm{i}(z_{k}-l)))\,x-2{\rm Re}(y(z_{k}))\,t)+\ln(\alpha_{k}).$ Fixing the line $x-v_{k}t={\rm const}$, and taking the limit $t\to-\infty$, we can obtain that

[TABLE]

where

[TABLE]

and

[TABLE]

We list the concrete sequences of $\xi_{i}^{[a,b]}$ and $\eta_{j}^{[a,b]}$:

[TABLE]

and

[TABLE]

and

[TABLE]

and

[TABLE]

To calculate the determinant of above matrix. We use the results from the Appendix D of Ref. [36], which may be regarded as a higher-order version of Fay’s identity. The matrices $\mathbf{M}^{[a,b]}$ and $\mathbf{K}^{[a,b]}$s are the theta function version Cauchy matrices, whose determinants can be represented as

[TABLE]

We observe that above four pairs of vectors have different $k$-th components. The other components are the same. Furthermore, the solution can be reduced as

[TABLE]

where

[TABLE]

and

[TABLE]

Then the asymptotic behavior for the multi-breather solution with different velocity can be concluded as the following theorem:

Theorem 4**.**

When $t\to-\infty$, the breather solution can be asymptotically expressed as

[TABLE]

where $\beta_{k}^{[-\infty]}=\prod_{i=1}^{k-1}\frac{\gamma_{i}}{\gamma_{i}^{*}}\prod_{i=k+1}^{n}\frac{\gamma_{i}^{*}}{\gamma_{i}}$

[TABLE]

and $\alpha^{[-\infty]}_{k}=\alpha_{k}\Delta_{k}^{[-\infty]}$, $s_{k}^{[-\infty]}=\sum_{i=1}^{k-1}\mathrm{i}(z_{i}^{*}-z_{i})-\sum_{i=k+1}^{n}\mathrm{i}(z_{i}^{*}-z_{i}),$

[TABLE]

with

[TABLE]

In a similar procedure, we can consider the asymptotic analysis when $t\to+\infty$, and a similar expression for the asymptotic breather solution $\psi_{k}^{[+\infty]}$ can be obtained by replacing $\alpha^{[-\infty]}_{k}$, $s_{k}^{[-\infty]}$ with $\alpha^{[+\infty]}_{k}$, $s_{k}^{[+\infty]}$, where $\alpha^{[+\infty]}_{k}=\alpha_{k}\Delta_{k}^{[+\infty]}$, $s_{k}^{[+\infty]}=-\sum_{i=1}^{k-1}\mathrm{i}(z_{i}^{*}-z_{i})+\sum_{i=k+1}^{n}\mathrm{i}(z_{i}^{*}-z_{i})$, and $\beta_{k}^{[+\infty]}=\prod_{i=1}^{k-1}\frac{\gamma_{i}^{*}}{\gamma_{i}}\prod_{i=k+1}^{n}\frac{\gamma_{i}}{\gamma_{i}^{*}}$

[TABLE]

Actually, this proves the interaction between breathers is hard to say elastic or inelastic. Through the above asymptotic analysis, we find the expression of breathers before and after the interaction just change some parameters. Actually, these parameters will affect the dynamics of breathers, since the background solution is not uniform as the zero or plane wave background.

In the following, we give some examples to illustrate the dynamics for the single breather and their interaction. Firstly we show the single breather solution, i.e. the formula (3.7) when $n=1$. From this formula, we show the asymptotic behavior for the single breather solution. We have the asymptotic analysis:

[TABLE]

as $\beta_{1}(x-v_{1}t+\gamma_{1})\to-\infty$, and

[TABLE]

as $\beta_{1}(x-v_{1}t+\gamma_{1})\to+\infty.$ From the above asymptotic expression, we know the phase differences is $2\ln\left(\frac{\vartheta_{4}(\frac{\mathrm{i}(z_{1}^{*}-l)}{2K})\vartheta_{2}(\frac{\mathrm{i}(z_{1}+l)}{2K})}{\vartheta_{4}(\frac{\mathrm{i}(z_{1}-l)}{2K})\vartheta_{2}(\frac{\mathrm{i}(z_{1}^{*}+l)}{2K})}\right)$ and the shift of crest is $4\mathrm{Im}(z_{1})/\alpha$. The parameter $\alpha_{1}$ will affect the maximum of amplitude. If ${\rm Im}(\lambda(z_{1}))>0$, the maximum value $\nu_{3}+2{\rm Im}(\lambda(z_{1}))$ ($\nu_{3}$ refers to equation (2.38)) appears at $(x,t)=(0,0)$ by choosing $\alpha_{1}=-1$ and the locally minimum value $\nu_{3}-2{\rm Im}(\lambda(z_{1}))$ appears at $(x,t)=(0,0)$ by choosing $\alpha_{1}=1$. Instead, if ${\rm Im}(\lambda(z_{1}))<0$, the maximum value emerges when $\alpha_{1}=1$ and the locally minimum value emerges for $\alpha_{1}=-1.$

To show the dynamics for the single breather solution, we choose the exact parameters and plot the figure. We discuss them in the following four cases.

Example 1: The breather solutions on the dn solution background.

In this case, we choose the parameters $l=\frac{K^{\prime}}{2}$. By the shift formula (see Appendix B), the single breather solution can be rewritten as the following form:

[TABLE]

where

[TABLE]

and $C=K+\mathrm{i}K^{\prime}$,

[TABLE]

Even though the breather solution above presents a compact form, it is not readily to obtain the properties from the above expression directly. The maximum of the modulus $|\psi_{1}|$ attains to $\max(|\psi_{0}|)+2|\textrm{Im}(\lambda(z_{1}))|$, where

[TABLE]

Through the expression, we find that the value of $|\psi_{0}|$ attains to the maximum at $x=2mK$, $m\in\mathbb{Z}.$ Thus the maximum of the breather solution can merely attain to the maximum value at these lines $x=2mK$, $m\in\mathbb{Z}.$ In what follows, we give two exact examples on the breather solutions. Firstly, we consider the stationary breather by choosing the parameters of background solution $\alpha=1$, $k=\sqrt{m}=\frac{\sqrt{2}}{2}$. It follows that $\nu_{1}=0$, $\nu_{2}=\frac{1}{2}$, $\nu_{3}=1$, $s_{2}=\frac{1}{2}(\nu_{1}+\nu_{2}+\nu_{3})=\frac{3}{4}$. By taking the parameters of breathers $z_{1}=-\frac{K^{\prime}(1/2)}{2}-\frac{K(1/2)\mathrm{i}}{3}$ and $\alpha_{1}=-1$, we plot Fig. 2. For these values, we find the parameters $\lambda_{1,3}\approx\pm\mathrm{i}0.854$ and $\lambda_{2,4}\approx\pm\mathrm{i}0.146$. Furthermore, the parameter $y_{1}\approx 0.476$ infers that the period of breather in the direction of $t$ is $T=\frac{\pi}{|{\rm Im}(y_{1})|}\approx 6.602.$ The parameters $\lambda(z_{1})\approx 0.984\mathrm{i}$ and $\alpha_{1}=-1$ can deduce that the maximum of the breather is $1+2{\rm Im}(\lambda(z_{1}))=2.967$ at $x=0$ and $t=mT$, $m\in\mathbb{Z}.$ The parameter $\alpha(Z(C-\mathrm{i}(z_{1}+l))-Z(\mathrm{i}(z_{1}-l)))\approx-0.978$ shows that $|E_{1}|\to\infty$ as $x\to-\infty$ and $|E_{1}|\to 0$ as $x\to+\infty.$ It follows that

[TABLE]

and

[TABLE]

where

[TABLE]

The above asymptotic behaviors can be seen and verified in Fig. 3.

Now we consider the non-stationary breather solution. In this case, we just illustrate the solutions by graphs. The properties are similar as the stationary breathers (Fig. 2 (c)). So we skip the detailed analysis of their dynamics.

Example 2: The breather solutions on the cn solution background.

Choosing the parameter $l=0$ and $n=1$ in (3.7), we obtain the one-breather solution with ${\rm cn}$ background. Fixing the parameters of background solution $\alpha=1$, $k=\sqrt{m}=\frac{1}{\sqrt{2}}$, we have $\nu_{1}=-\frac{1}{2}$, $\nu_{2}=0$, $\nu_{3}=\frac{1}{2}$, $s_{2}=0$. By taking the parameters $\alpha_{1}=1$, $z_{1}=-K(1/2)+\mathrm{i}\frac{K^{\prime}(1/2)}{4}$, we obtain the plot Fig. 4(a). For these parameters, we obtain $\lambda(z_{1})=-1.07\mathrm{i}$, $y(z_{1})\approx 1.19$, and the maximum value is about $2.85$ at the origin. Changing the parameters $\alpha_{1}=1/20$, we obtain the breather with different shape (Fig. 4(c)). If we choose the parameters $z_{1}=-\frac{9}{10}K(1/2)+\frac{K^{\prime}(1/2)\mathrm{i}}{4}$ and $\alpha_{1}=1$, the non stationary breather is obtained (figure 4(b)). The maximum value is about $2.56$ at the origin. Comparing Figs. 2 and 4, we see that the two kinds of breathers with different backgrounds possess evidently different dynamics.

Example 3: The breather solutions on the nontrivial-phase solution background.

Here we consider the case where the elliptic-function background has phase modulation. In particular, we choose the parameters $k=\sqrt{m}=\frac{1}{\sqrt{2}}$, $l=\frac{K(1/2)}{4}$, $\alpha=1$, which can deduce that $\nu_{1}=-\alpha^{2}{\rm dn}^{2}(K(1/2)+2\mathrm{i}l)$, $\nu_{2}=-k^{2}\alpha^{2}{\rm cn}^{2}(K(1/2)+2\mathrm{i}l)$, $\nu_{3}=k^{2}\alpha^{2}{\rm sn}^{2}(K(1/2)+2\mathrm{i}l)$, $s_{2}=\frac{1}{2}(\nu_{1}+\nu_{2}+\nu_{3})$. By choosing the parameter $z_{1}=-0.5+1.223\mathrm{i}$, $\alpha_{1}=1$, we can plot the figure of breather(Fig.  5 (a)). Moreover, $\lambda(z_{1})=0.087-0.852\mathrm{i}$, $y_{1}\approx-0.586$. The period of solution is $T\approx 5.35$. The maximum value is about $2.546$ at the origin.

For the second case, by choosing the parameters of breather $z_{1}=-\frac{K(1/2)}{2}+\frac{K^{\prime}(1/2)\mathrm{i}}{8}$, $\alpha_{1}=1$, which infers that $\lambda(z_{1})\approx 0.797-0.480\mathrm{i}$, $y_{1}=-0.665+0.719\mathrm{i}$, we can plot the breathers(Fig. 5 (b)). The maximum value of breather is about $1.801$.

Example 4: Two-breather solutions.

Now we consider the two-breather solutions and its analysis. Meanwhile, we use the numeric way to verify the above asymptotic analysis.

We consider the two-breather solution with ${\rm cn}$ background in detail. We choose the parameters of the background solution $l=0$, $\alpha=1$, $k=\frac{1}{\sqrt{2}}$, which determine the parameter $s_{2}=0$ and the background solution. Taking the parameters of breathers as $z_{1}=-\frac{9K(1/2)}{10}+\frac{K^{\prime}(1/2)\mathrm{i}}{4}$, $z_{2}=\frac{9K(1/2)}{10}+\frac{K^{\prime}(1/2)\mathrm{i}}{4}$, $\alpha_{1}=\alpha_{2}=1$, we can obtain the two breathers (figure 6 (a)). The maximum value is about $4.416$ at the origin.

We illustrate the asymptotic analysis on the two cn breathers, i.e two breather solution with cn elliptic function background. It is seen that the asymptotic expression for $t\to-\infty$, which are obtained in Theorem 4, is consistent with two cn breathers [Fig. 7(a) and (b)]. The asymptotic behaviors for $t\to+\infty$ also show good agreement [Fig. 7(c) and (d)].

As another case, we consider the two-breather solution on the non-trivial phase background. By taking the parameters $k=\frac{1}{\sqrt{2}}$, $l=\frac{K(1/2)}{4}$ and $\alpha=1$, we obtain $\nu_{1}=-{\rm dn}^{2}(K(1/2)+2\mathrm{i}l)$, $\nu_{2}=-\frac{1}{2}{\rm cn}^{2}(K(1/2)+2\mathrm{i}l)$, $\nu_{3}=\frac{1}{2}{\rm sn}^{2}(K(1/2)+2\mathrm{i}l)$ and the background solution. On the other hand, taking the parameters $\alpha_{1}=\alpha_{2}=1$, $z_{1}=-\frac{9K(1/2)}{10}+\frac{K^{\prime}(1/2)\mathrm{i}}{4}$, $z_{2}=-0.5+1.2335\mathrm{i}$, we can plot two breathers (Fig. 6 (b)). The maximum value is about $4.198$ located at the origin. The interaction can be obtained in the above and be verified as the cn case too.

4. Modulational instability and Akhmediev breathers on the elliptic function background

Based on the linearized NLSE (2.4), we analyze the spectral stability for the elliptic-function solution of NLSE. As shown previously, the elliptic-function solution can be represented as an abstract form $\psi=u(x)\mathrm{e}^{\mathrm{i}s_{2}t}$. Let us assume that the perturbed solution has the form

[TABLE]

then the linearized equation is converted into the following spectral problem:

[TABLE]

With the aid of the stationary Lax pair, we can solve the linear spectral problem with stationary periodic solution exactly. Thus the spectrum for the spectral problem can be determined.

As shown in Sec. 2.1, $p(x,t)=\phi_{1}(x,t;\lambda)\phi_{2}(x,t;\lambda)+\varphi_{1}^{*}(x,t;\lambda)\varphi_{2}^{*}(x,t;\lambda)$ solves the linearized perturbed equation. Together with the vector solution for the Lax pair, we can set them as the form of equation (4.1). It follows that

[TABLE]

with $\Omega=2y;$ or

[TABLE]

with $\Omega=-2y.$ We have four different roots $\lambda=\xi_{i}$, $i=1,2,3,4$ for fixed $\Omega.$ By the transformation $\lambda(z)=\mu(\mathrm{i}(z-l)/\alpha)$, we know that this establishes the corresponding between the Riemann surface and the torus. Thus we have four different points $\mu(\mathrm{i}(z_{i}-l)/\alpha)=\xi_{i}.$ Inserting the four different roots $\xi_{i}$s into equations (4.3) and (4.4), we obtain four independent solutions for the linearized spectral problem (4.2). So the spectral stability can be obtained by the analysis of solutions $p_{1}(x)$, $p_{2}(x)$, the bound properties of which will determine the spectrum of spectral problem (4.2). Moreover, the coefficient of the exponential function in $p_{1}(x)$ and $p_{2}(x)$, which may be interpreted as a crystal momentum, determines the spectrum: i.e. $\textrm{Re}(\alpha Z(\mathrm{i}(z-l))+\mathrm{i}\lambda)=0$ or $\textrm{Re}(\alpha Z(K+\mathrm{i}K^{\prime}-\mathrm{i}(z+l))+\mathrm{i}\lambda)=0$. We know the determinant of solution $\Phi(x,t;\lambda)$ is a constant, so that $\textrm{Re}(\alpha Z(\mathrm{i}(z-l))+\alpha Z(K+\mathrm{i}K^{\prime}-\mathrm{i}(z+l))+2\mathrm{i}\lambda)=0$. Then the spectral condition becomes $\textrm{Re}(Z(\mathrm{i}(z-l))-Z(K+\mathrm{i}K^{\prime}-\mathrm{i}(z+l)))=0.$ Following the method proposed by Deconinck and Segal [16], the spectrum can be obtained numerically by the tangent vector field

[TABLE]

and the starting point $z=\frac{1}{2}(K^{\prime}\pm\mathrm{i}K)$ and $z=\frac{1}{2}(-K^{\prime}\pm\mathrm{i}K)$. By the formula $\lambda(z)=\mu(\mathrm{i}(z-l)/\alpha)$, this determines the spectrum in the $\lambda$ plane. By the theory of Riemann surfaces, we know that the Abel map

[TABLE]

establishes the conformal mapping between the genus $1$ Riemann surface $\mathbb{S}$ and the rectangular region $\{z|-K^{\prime}\leq{\rm Im}(z-l)\leq K^{\prime},\,\,\,\,-K\leq{\rm Re}(z-l)\leq K\}$, by fixed the integral path in the Riemann surface, i.e., the relation (2.43). Through the above analysis, the spectrum with respect to $\lambda$ can be obtained by the relation (2.43). Also, the spectrum with respect to $\Omega=\pm 2y$ can be obtained via (2.44).

In this section, our main goal is not to analyze the stability properties for the NLSE, but to construct a special solution which can describe the modulational instability. In the plane wave background, this type of solution is called Akhemediev breather, so we also use the same name in the case of the elliptic backgrounds.

Firstly, we consider the ${\rm dn}$ case by fixing the parameters $k=\frac{1}{\sqrt{2}}$, $\alpha=1$, $l=\frac{K^{\prime}}{2}$. The resultant unstable spectra are shown in Fig. 8. This figure shows that there are two lines $z=s\pm\frac{K^{\prime}}{2}$, $(-K\leq z-l\leq K)$ in $z$ plane. The unstable spectrum in $\lambda$ plane consists of two line segments located on the imaginary axis. In $\Omega$ plane, it is represented by one line segment on the imaginary axis.

For the second case, we consider the ${\rm cn}$ case ($\nu_{1}=-1/2$, $\nu_{2}=0$ and $\nu_{3}=1/2$, i.e. $l=0$, $\alpha=1$ and $m=k^{2}=1/2$). Here, we demonstrate the above-mentioned method in Ref. [16] to derive the unstable spectrum. Firstly, we have $K(k)$ and $K^{\prime}(k)$. From the four starting points $z=\frac{\mathrm{i}}{2}(K(k)\pm\mathrm{i}K^{\prime}(k)),\frac{\mathrm{i}}{2}(-K(k)\pm\mathrm{i}K^{\prime}(k))$, the curve can be replaced by the tangent line with the vector $\pm v_{T}(z)$ [equation (4.5)] for the sufficiently small segment, where $\pm$ determines the direction. Combining the eight segments, we obtain the unstable spectrum with respect to $z$. Furthermore, using the formulas (2.43) and (2.44), we have the figure (Fig. 9). It shows that there are two waist type curves in $z$ plane. Moreover, there are two almost vertical curves in $\lambda$ plane, and one figure-eight curve in $\Omega$ plane.

In a similar way, we obtain the figure for the non-trivial phase case (Fig. 10). It shows that there are two non-symmetric waist type curves in $z$ plane. Moreover, there are two non-symmetric almost vertical curves in $\lambda$ plane, and two figure-eight curves in $\Omega$ plane. For the general discussion on the unstable spectrum on the elliptic function solution of NLSE, we can refer to the reference [16]. Here we just provide some special cases to illustrate the relationship between the unstable spectrum and the exact breather solutions in the following.

An interesting application of the above unstable spectra is that we can construct the exact solution to illustrate the modulation dynamics to describe the unstable process. We can obtain such solutions by setting the breather parameter $z_{1}$ in one-breather solution [equation (3.7) with $n=1$] to be values on the lines in Figs. 8-10. Below we provide a few examples.

Let the parameters be $k=\frac{1}{\sqrt{2}}$, and $\alpha=1$, $l=\frac{K^{\prime}}{2}$, which infer $s_{2}=\frac{3}{4}$. Choosing the parameter $z_{1}=-0.527-\frac{K(1/2)}{2}\mathrm{i}\approx-0.527-0.927\mathrm{i}$, $\alpha_{1}=-1$, we obtain the Akhmediev dn breather (Fig 11 (a)). The maximum value is about $2.384$ located at the origin. From the figure 11 (a), we see that there are lots of peaks located on the line $t=0$, but the value of peaks is less than $2.384$, which can be obtained by the analysis of formula (A.2) and is distinct with the Akhmediev breather in the plane wave background. In this setting, it is special that there is no shift of crest under the special parameters setting (Fig 11 (b)), by the above analysis (3.11) and (3.12) for the single breather solution. Here we should point out that this phenomenon with no shift on the crest is a general rule in the dn background. Through the figure 11 (a) and the analysis, we show that this exact breather solution describes the unstable dynamical behavior induced by the unstable spectrum.

Let us set $k=\frac{1}{\sqrt{2}}$, $l=0$, and $\alpha=1$, which implies $s_{2}=0$. Choosing the parameters $z_{1}\approx-0.771-0.560\mathrm{i}$, $\alpha_{1}=-1$, we arrive at the Akhmediev cn breather, which admits the phase and crest shift simultaneously (Fig 12 (a)). The maximum value of the peak is about $1.264$ located at origin. As the analysis in the above, we know that the other peaks’ value is less than $1.264$. In this setting, we can obtain the shift of crest and phase through above general analysis. The dynamics of breather describes the unstable process of the cn solution.

Finally, let us take the parameters $k=\frac{1}{\sqrt{2}}$, $\alpha=1$, and $l=\frac{K(1/2)}{8}$, which yield $\nu_{1}=-{\rm dn}^{2}(K(1/2)+2\mathrm{i}l)$, $\nu_{2}=-\frac{1}{2}{\rm cn}^{2}(K(1/2)+2\mathrm{i}l)$, $\nu_{3}=\frac{1}{2}{\rm sn}^{2}(K(1/2)+2\mathrm{i}l)$ and $s_{2}=\frac{1}{2}(\nu_{1}+\nu_{2}+\nu_{3})$. Choosing parameters $z_{1}=-0.663-0.503\mathrm{i}$, $\alpha_{1}=-1$, we plot Fig. 12 (b). The maximum value of breather is about $1.414$ located at origin, which is greater than the other peaks located on the line $t=0$. This breather describes the unstable dynamics for the elliptic-function background with non-trivial phase modulation.

5. Rogue waves on the elliptic function background

In this section, we consider the rogue waves, which are located at the branch points. In this case, above two vector solutions (2.36) merge into one vector solution. Thus we need to look for another linear independent vector solution. Here we would like to utilize the limit technique to deal with this problem. Due to the theta function representations include the parameters by Abel map, we still use the original parameter to derive the formulas of rogue waves. Suppose we expand the parameter $y(\lambda)$ at the branch point $\lambda_{1}$: $\lambda=\lambda_{1}+\epsilon_{1}^{2}$, here $\epsilon_{1}$ is a small parameter, then we have the expansions:

[TABLE]

where

[TABLE]

It follows that

[TABLE]

and

[TABLE]

Moreover, the elements in the vector solution can be expanded in the following:

[TABLE]

where

[TABLE]

We need the combinatorial formulas:

[TABLE]

and

[TABLE]

Let us introduce the notation

[TABLE]

and

[TABLE]

then we arrive at

[TABLE]

where $[\cdot]$ represents the floor function. In a similar fashion, we get the expansion for $1/(\nu(x)-\beta_{1}(\lambda))$:

[TABLE]

where we have defined

[TABLE]

With the help of equations (5.3) and (5.6), it follows that

[TABLE]

where

[TABLE]

Then we obtain the expansion for the phase factor:

[TABLE]

Remark 4**.**

Here we point out that the integrations $\int_{0}^{x}I_{j}(s){\rm d}s,$ $(i\geq 1)$ do not possess integrals of the third kind. In general, the integration $\int_{0}^{x}I_{j}(s){\rm d}s$ can be expressed in terms of elliptic functions and elliptic integrals of the second and third kind via the formulas (B.1),(B.2) and (B.3). It is hard to give a claim that the integrations $\int_{0}^{x}I_{j}(s){\rm d}s,$ $(i\geq 1)$ do not have the elliptic integral of the third by the pure algebraic calculations. It is much easier to give a proof by analyzing the singularity point at the point $x=\xi\mod\Lambda\in\{2mK+n\mathrm{i}\,K^{\prime}|m,n\in\mathbb{Z}\}$ when $\nu(\xi)-(2\lambda_{1}^{2}+s_{2})=0$. By using the conformal transformation (2.43) and (2.44), the integration

[TABLE]

has the logarithmic singularity at $x=\xi\mod\Lambda$. Expanding the above integration at $\lambda=\lambda_{1}$ is equivalent to the expansion at $z=z_{s}\equiv\frac{K^{\prime}+\mathrm{i}K}{2}=-\mathrm{i}(\alpha\,\xi+z_{0}).$ Firstly, we have the following two expansions:

[TABLE]

and

[TABLE]

Thus we have

[TABLE]

which infers that the logarithmic singularity can not appear at the higher-order term of $(z-z_{s}).$ Moreover, the higher-order terms of $(z-z_{s})$ do not involve elliptic integrals of the third kind since the leading order term involves the logarithmic singularity. From the conformal transformation

[TABLE]

we can deduce that, at the critical point $z=z_{s}$, we have the expansion $\lambda(z)-\lambda_{1}=(z-z_{s})^{2}H^{2}(z-z_{s})$, where $H(z-z_{s})$ is an analytic function with $H(0)\neq 0.$ So $\epsilon_{1}=(z-z_{s})H(z-z_{s})$ and

[TABLE]

which implies that there appears no elliptic integral of the third kind in the coefficients of higher-order terms of $\epsilon_{1}$.

By the Schur polynomial, we have

[TABLE]

where

[TABLE]

Summarizing, we obtain the expansion for the first component of vector solution:

[TABLE]

To derive the expansion for the second component, we need to expand the following expression:

[TABLE]

In fact, the second component of the vector solution has the expansion

[TABLE]

where

[TABLE]

In a similar procedure, we obtain the expansion at $\lambda=\lambda_{2}$:

[TABLE]

where $\Theta_{2}^{[i]}$, $\phi_{2}^{[i]}$ and $\varphi_{2}^{[i]}$ are $\Theta_{1}^{[i]}$, $\phi_{1}^{[i]}$ and $\varphi_{1}^{[i]}$ by exchanging the variables $\lambda_{1}$ to $\lambda_{2}.$ By the Remark 4 and the integration formulas (B.1), (B.2) and (B.3), we know that the functions $\phi_{k}^{[i]}$, $\varphi_{k}^{[i]}$ ($k=1,2;i\in\mathbb{Z}^{+}$) are the rational function of $\mathrm{sn}$, $\mathrm{cn}$, $\mathrm{dn}$ and the second type of elliptic integration $E(u)$.

Summarizing the above, we can obtain the following theorem of rogue waves:

Theorem 5**.**

The general rogue wave formulas obtained by the high order Darboux transformation (in the Appendix A) can be represented as

[TABLE]

where $\mathbf{M}_{n}=\nu(x)\mathbf{X}_{n}^{{\dagger}}\mathbf{D}_{n}\mathbf{X}_{n}+\mathbf{Y}_{n}^{{\dagger}}\mathbf{D}_{n}\mathbf{Y}_{n}$, $\mathbf{K}_{n}=\mathbf{M}_{n}-2\mathbf{X}_{n,1}^{{\dagger}}\mathbf{Y}_{n,1}$, $\mathbf{X}_{n,1}=[\mathbf{X}_{1}^{[1]},\mathbf{X}_{1}^{[2]}]$, $\mathbf{Y}_{n,1}=[\mathbf{Y}_{1}^{[1]},\mathbf{Y}_{1}^{[2]}]$ and

[TABLE]

[TABLE]

and $l,k=1,2$, $n_{1}+n_{2}=n,$

[TABLE]

and $\mathbf{X}_{1}^{[l]}$ and $\mathbf{Y}_{1}^{[l]}$ represent the first row of $\mathbf{X}^{[l]}$ and $\mathbf{Y}^{[l]}$, respectively.

The general expressions for the high order rogue waves are very complicated. However, the extremal value at the point $(x,t)=(0,0)$ can be obtained through the determinant formula. Since the NLSE possesses the scaling symmetry, we set the amplitude for the background be $1$: i.e., $\nu_{3}=1$ without loss of generality. Then, the vector fundamental solution at $(x,t)=(0,0)$ can be written as

[TABLE]

At the branch points $\lambda=\lambda_{i}+\epsilon_{i}^{2}$, we know that $y(-\epsilon_{i})=-y(\epsilon_{i})$, $i=1,2$. So we have the expansion

[TABLE]

Moreover, we have

[TABLE]

where $\mathbf{M}_{n}=\mathbf{X}_{n}^{{\dagger}}\mathbf{D}_{n}\mathbf{X}_{n}$, $\mathbf{X}_{1,n}=[\phi_{1}^{[1]}(\lambda_{1}),\phi_{1}^{[3]}(\lambda_{1}),\cdots,\phi_{1}^{[2n_{1}-1]}(\lambda_{1}),\phi_{1}^{[1]}(\lambda_{2}),\phi_{1}^{[3]}(\lambda_{2}),\cdots,\phi_{1}^{[2n_{2}-1]}(\lambda_{2})],$ and

[TABLE]

By the formula of the Cauchy determinant and taking its special limit, we find

[TABLE]

where

[TABLE]

and $1\leq l,k\leq 2$, $\delta_{1,1}=1$ if $i=j=1$, otherwise $\delta_{1,1}=0.$

Here we stress that the high order rogue waves given previously do not possess the free parameters. Actually, since the fundamental solution can involve lots of free parameters, it is easy to bring the free parameters into the formulas of rogue waves as in reference [21].

Now we consider the several explicit examples of the rogue waves. Firstly, we consider the first order rogue wave. We fix the parameters $\lambda_{1}=b+c\mathrm{i}$, $\lambda_{2}=-b+d\mathrm{i}$, $\lambda_{3}=\lambda_{1}^{*}$ and $\lambda_{4}=\lambda_{2}^{*}$. Using the above expansion, we can obtain that

[TABLE]

where $\theta_{0}(x,t)=q\int_{0}^{x}\frac{{d}s}{\nu(s)}-s_{2}t$, $p(x)=\frac{1}{u_{1}(x)}+\mathrm{i}t+J(x)+k_{1}$, $u_{1}(x)=\nu(x)-(2\lambda_{1}^{2}+s_{2})$, $C_{1}=2\lambda_{1}(2\lambda_{1}^{2}+s_{2})-s_{3}$, $k_{1}\in\mathbb{C}$, and the integration

[TABLE]

where $\alpha=\sqrt{4b^{2}+(c+d)^{2}}$, $A_{1}=\frac{-\alpha}{2cA_{0}}$, $A_{2}=\frac{1}{2c}$, $A_{3}=\frac{2\alpha d}{A_{0}}$, $A_{0}=(\lambda_{1}-\lambda_{2})(\lambda_{1}-\lambda_{2}^{*})$. Thus the rogue wave solution can be obtained as

[TABLE]

The parameter $k_{1}$ determines not only the location of the rogue wave but also its shape and dynamics, since it is influenced by the periodic modulation of the background. The rogue waves on cn background and dn background at upper branch point were given in [12]. If we switch the parameter $\lambda_{1}$ to $\lambda_{2}=-b+d\mathrm{i}$, then we obtain another family of parameters $A_{1}=\frac{-\alpha}{2dA_{0}}$, $A_{2}=\frac{1}{2d}$, $A_{3}=\frac{2\alpha c}{A_{0}}$, $A_{0}=(\lambda_{2}-\lambda_{1})(\lambda_{2}-\lambda_{1}^{*})$, which will determine another family of rogue wave solutions arising from the different branch point.

Here if we intend to prove the asymptotic behavior when $x,t\to\infty$, we need to use the theta-function representation mentioned previously. It is easy to prove that the rogue wave solution approaches the background solution by the shift of phase and crest (dn,cn or non-trivial phase function). Since the calculation is similar to Sec. 3, here we omit the detailed asymptotic analysis.

Using Theorem 5, we can also obtain the two-rogue-wave solution:

[TABLE]

where

[TABLE]

and $\phi_{i}$ and $\varphi_{i}$ are given by equations (5.12).

Example 1: The rogue waves and two rogue waves on the dn solution background

Since the NLSE possesses the scaling symmetry, we can fix the amplitude of background. So we choose the parameters $b=0$, $c=\frac{1}{2}(1+\kappa)$ and $d=\frac{1}{2}(1-\kappa)$, $(|\kappa|<1)$, then $\nu_{1}=0$, $\nu_{2}=\kappa^{2}$ and $\nu_{3}=1$. Moreover, $m=1-\kappa^{2}$, $q=-s_{3}=0$, $s_{2}=\frac{1}{2}(1+\kappa^{2})$. It follows that the elliptic function solution is

[TABLE]

By fixing the background solution, two different families of rogue wave solution can be obtained by the formula (5.13) depending on the choice of different branch points. When $\kappa\to 0$, we have $\psi\to\exp(\mathrm{i}t/2)$. So when $\kappa$ is small and approaches zero, the rogue waves would exhibit the similar scenario as the classical rogue waves on the plane wave background. Since there exist two different types of rogue waves, it is shown that while one of the two is similar to the classical one, the other one is different from the classical one, which locates on the lower branch point. When $\kappa$ is small enough, these rogue waves vanish and become the plane waves. On the other hand, when $\kappa\to 1$, then $\psi\to{\rm sech}(x)$. In this case, there still exist two types of rogue waves. It is shown that the rogue waves seem a second-order soliton in the center part. In this case, two types of rogue waves have the similar structure.

Choosing the parameters $\kappa=\frac{1}{3}$, $\lambda_{1}=\frac{2}{3}\mathrm{i}$, $k_{1}=0$ and inserting them into formula (5.13), we obtain the rogue wave solution on the dn background (Fig. 13 (a)). The maximum value of the rogue wave is $\frac{7}{3}$. It is seen that the central part of the dn rogue wave under this parameter looks likes a second-order soliton. For the other part, it is still likely the dn solution. If we change the parameter with $k_{1}=2$, the shape of the rogue wave changes (Fig. 13 (b)). The maximum value of this rogue wave is about $1.89$, which is less than $7/3$. Meanwhile the emergent peak is deviated from the line $x=0.$

Choosing another type parameters $\kappa=\frac{1}{3}$, $\lambda_{2}=\frac{1}{3}\mathrm{i}$, $k_{1}=0$ and plugging them into formula (5.13), we obtain another type of rogue waves. The maximum value is about $\frac{5}{3}$ (Fig. 13(c)). In [12], the authors do not give this case since the integration is singular. In the same way as above, if the parameter $k_{1}$ varies, the shape also changes. Here we do not show the figure.

It is interesting that the two rogue waves can exist simultaneously, which can be achieved by the two-rogue wave formula (5.14). Taking the parameters $\kappa=\frac{1}{3}$, $\lambda_{1}=\frac{2}{3}\mathrm{i}$, $\lambda_{2}=\frac{1}{3}\mathrm{i}$, $k_{1}=k_{2}=0$ and inserting them into the formula of two rogue waves (5.14), we obtain the two rogue waves with the maximum value $3$ at origin (Fig. 13(d)). In this parameter setting, the two rogue waves attain the maximum value. If we change the parameters $k_{1}$ and $k_{2}$, we obtain two types of rogue wave situated at different locations.

Example 2: The rogue waves and two rogue waves on the cn solution background

Since the NLSE possesses the scaling symmetry, we can fix the amplitude of background. So we choose the parameters $c=d=1/2$, then $\nu_{1}=-4b^{2}$, $\nu_{2}=0$ and $\nu_{3}=1$. Moreover, $m=1/(1+4\,b^{2})$, $q=-s_{3}=0$, $s_{2}=\frac{1}{2}-2b^{2}$. It follows that the elliptic-function solution is

[TABLE]

By fixing the background solution, two different families of rogue wave solutions originating from two branch points can be obtained by the formula (5.13). In this case, we find that two kinds of rogue waves move to the left and right. When $\kappa\to 0$, we have $\psi\to\exp(\mathrm{i}t/2)$. So when $\kappa$ is small and is close to zero, the rogue waves exhibit the similar scenario as the classical rogue waves on the plane wave background. Since there exist two different types of rogue waves, it is shown that one of them is similar to the classical one, while the other one is different from the classical one, which locates on the lower branch point. When $\kappa$ is small enough, these rogue waves merge into the mere plane wave solution. On the other hand, when $\kappa\to 1$, i.e., $\psi\to{\rm sech}(x)$, there still exist two families of rogue waves. It is shown that the rogue waves looks like a second-order soliton in the center part, and under this case, two rogue waves have the similar structure.

Choosing the parameters $b=\frac{1}{2}$, $\lambda_{1}=-\frac{1}{2}+\frac{\mathrm{i}}{2}$, $k_{1}=0$, we plot the figure for the rogue wave on the ${\rm cn}$ function background. The maximum value equals to 2 located at origin (Fig. 14 (a)). Also, if we change the parameters $\lambda_{2}=\frac{1}{2}+\frac{\mathrm{i}}{2}$, $k_{2}=0$, we obtain another type of rogue waves. The maximum value equals to 2 located at origin (Fig. 14 (b)).

For the two rogue waves, we set the parameters $b=\frac{1}{2}$, $\lambda_{1}=-\frac{1}{2}+\frac{\mathrm{i}}{2}$, $\lambda_{2}=\frac{1}{2}+\frac{\mathrm{i}}{2}$, $k_{1}=k_{2}=0$, and then we find that the maximum value of two rogue waves equals to 3 located at origin (Fig. 14(c)). If we change the parameters $k_{1}$ and $k_{2}$, two rogue waves change their locations.

Example 3: The rogue waves and two rogue waves on the general elliptic-function solution background

Since the NLSE possesses the scaling symmetry, we can fix the amplitude of background. So we choose the parameters $b\neq 0$, $c=\frac{1}{2}(1+\kappa)$ and $d=\frac{1}{2}(1-\kappa),$ $(\kappa\neq 0)$, then $\nu_{1}=-4b^{2}$, $\nu_{2}=\kappa^{2}$ and $\nu_{3}=1$. Moreover, $m=(1-\kappa^{2})/(1+4b^{2})$, $q=-s_{3}=-2b\kappa$, $s_{2}=-2b^{2}+\frac{1}{2}(1+\kappa^{2})$. It follows that the elliptic function solution is

[TABLE]

By fixing the background solution, the two different families of the rogue wave solutions can be obtained by the formula (5.13) corresponding to the two branch points. In this case, we obtain two rogue waves which have the asymmetry pairs.

Inserting the parameters $b=\frac{\sqrt{5}}{4}$, $\kappa=\frac{1}{3}$, $\lambda_{1}=-\frac{\sqrt{5}}{4}+\frac{\mathrm{i}}{3}$, $k_{1}=0$ into the formula (5.13), we obtain the rogue waves on the non-trivial phase background (Fig. 15(a)). The maximum value is $\frac{5}{3}$ located at origin.

Similarly, inserting the parameters $b=\frac{\sqrt{5}}{4}$, $\kappa=\frac{1}{3}$, $\lambda_{2}=\frac{\sqrt{5}}{4}+\frac{2\mathrm{i}}{3}$, $k_{1}=0$ into the formula (5.13), we arrive at another type of rogue waves (Fig. 15(b)). In this parameter setting, the maximum value is $\frac{7}{3}$.

Plugging the parameters $b=\frac{\sqrt{5}}{4}$, $\kappa=\frac{1}{3}$, $\lambda_{1}=-\frac{\sqrt{5}}{4}+\frac{\mathrm{i}}{3}$, $\lambda_{2}=\frac{\sqrt{5}}{4}+\frac{2\mathrm{i}}{3}$, $k_{1}=k_{2}=0$ into the formula (5.13), we obtain the two rogue waves, which attain the maximum amplitude $3$ (Fig. 15(c)).

Example 4: Multi high order rogue waves

For the general higher-order rogue waves, we have given the general formula above. In practice, since the expression $y_{1}$ involves the roots of complex number, we can absorb the parameter $y_{1}$ into parameter $\epsilon_{i}$ to make the expressions more compact. On the other hand, the general expression for the elements are very complicated, and hence we only provide the second-order and two-second order rogue waves in the ${\rm cn}$ background. For the systematic treatment of the dynamics of general higher-order rogue waves, we would like to exhibit them in the other works.

By the computer plotting, we exhibit the figures for the fundamental second-second order rogue waves (Fig. 16), in which the exact formulas are shown in Appendix C. Both figures show the same solution by different coordinates to display the dynamics clearly. From the left figure, it is seen that four solitons collide at the origin, where they show the maximum value $5$. When $x$ and $t$ become large enough, the rogue waves tend to the background solution ${\rm cn}(2x|\frac{1}{2}){\rm e}^{-\mathrm{i}t},$ which can not be observed from the picture directly. The right figure shows the detail of the dynamics near the origin $(x,t)=(0,0)$, implying the more complicated profile than that of the plane-wave background.

6. Conclusion and Discussion

In this work, we systematically constructed the breathers and rogue waves for the focusing NLSE on the elliptic function background by the Darboux transformation or its generalized version. What is more, we provided the asymptotic analysis of breathers with different velocity by using the addition formula of theta functions. To verify the asymptotic analysis, we plot some figures of interaction of two breathers. Furthermore, we establish the relation between modulational instability and the Akhmediev breather or rogue wave solutions on the elliptic function background by the linear stability analysis, which describe the dynamical destabilization process of modulational instability in the same way as the case of the plane-wave background. The method provided here can be easily extended to the other integrable model, such as the sine-Gordon equation, modified KdV equation, derivative nonlinear Schödinger equation and so on.

For the physical significance, the solutions found in this work will be observed in physical experiments such as nonlinear fibre optics [27] with oscillating background and the rogue waves in the presence of the soliton trains in Bose-Einstein condensates [35, 25, 14]. Generalization of the solutions to the multi-component NLSE with modulated background will also be an interesting future issue, and there are several relevant works from both experiment and theory [22, 31, 29].

The long-time asymptotics for the focusing NLSE with plane wave background without discrete spectrum tends to an elliptic function in a triangle region [10]. It can be expected that the long-time asymptotics with discrete spectrum can be represented by the breathers in this work. Meanwhile, the classical inverse scattering method on plane wave background can not involve the rogue waves and high order rogue waves directly; we need to perform a sensitive limiting procedure to obtain them. Recently, Bilman and Miller [8] proposed a new scheme called robust inverse scattering transform to deal with the spectral singularity with the aid of the Darboux transformation. Meanwhile, it is hopeful to be used in analyzing the long-time asymptotics for the rogue waves and high order rogue wave on the plane wave background. The rogue waves and high order rogue waves are the possible results for their long-time dynamics. For both of the above works, they need to be analyzed in the future work with the aid of Riemann-Hilbert representation [28, 6].

Appendix A Darboux transformation on the elliptic function background

In this appendix, we give the relative results on the Darboux transformation. The theory for the Darboux transformation is very mature.

Elementary Darboux transformation and its multi and higher order version

The elementary Darboux transforation can be stated by the following theorem, which is essential to yield the multi-fold and higher order versions [15, 30, 21].

Theorem 6**.**

Assume that there is an entire matrix function $\Phi(x,t;\lambda)$ for the Lax pair (2.2) with the potential function $\psi(x,t)\in\mathbf{L}^{\infty}(\mathbb{R}^{2})\cup\mathbf{C}^{\infty}(\mathbb{R}^{2})$, which is normalized at $(x,t)=(0,0)$: $\Phi(0,0;\lambda)=\mathbb{I}$, then the Darboux transforation

[TABLE]

where $\Phi_{1}=\Phi(x,t;\lambda_{1})(1,c_{1})^{\rm T}\equiv(\phi_{1},\varphi_{1})^{\rm T}$ is a special solution $\lambda=\lambda_{1}$, converts the Lax pair (2.2) into the Lax pair with same form (2.2) by replacing the potential function $\psi(x,t)$ with a new one $\psi_{1}(x,t)$

[TABLE]

and the normalized matrix solution for the new system is $\Phi[1]=\mathbf{T}_{1}(x,t;\lambda)\Phi(x,t;\lambda)\mathbf{T}_{1}(0,0;\lambda)^{-1}$, which is also an entire matrix function with respect to $\lambda$.

In the background of the elliptic function, we need to the entire fundamental matrix solution $\Phi^{\rm N}(x,t;\lambda)$ given in (2.37).

Now we analyze the properties for the elementary Darboux transformation and Bäcklund transformation. Through the Bäcklund transformation and mean value inequality, we know that the amplitude of the new potential function ranges in the interval $[\min(|\psi|)-2\mathrm{Im}(\lambda_{1}),\max(|\psi|)+2\mathrm{Im}(\lambda_{1})]$. Now if we need to choose the solution $q_{1}$ to attend the maximum value $\max(|\psi|)+2\mathrm{Im}(\lambda_{1}),$ then $\varphi_{1}/\phi_{1}=-\mathrm{i}.$ For simplicity, we assume that $\psi_{1}(0,0)=1$ and $\varphi_{1}(0,0)=-\mathrm{i}.$ For instance, if we choose the seed solution as the plane wave solution $\psi=\psi_{BC}=1$. Meanwhile, the rogue wave solutions is located on the branch point $\lambda=\mathrm{i},$ in which case the vector solution $\Psi_{1}$ is not a periodic solution, but a rational function. To make the solution attend the maximum at the origin, we should impose the initial condition for the vector solution as $\Phi_{1}(0,0)=k(1,\mathrm{i})^{\rm T}.$

Actually, even though the formula for the Bäcklund transformation is simple, the new generated function is not easy to analyze by the classical tools of mathematical analysis for the elliptic-function background. For the breathers on the elliptic function background, we define that the center of breather is located at the maximum value point, which we choose the solution $\Phi_{1}=\Phi(x,t;\lambda_{1})(1,-\mathrm{i})^{T}$. In practice, the local minimum value can also be seen as the center of breather, i.e. $\Phi_{1}=\Phi(x,t;\lambda_{1})(1,\mathrm{i})^{T}$. For the other parameter choices, the determination of the center of breathers need the use of the numerical method.

Based on the above elementary Darboux transformation, we can iterate it to obtain the multi- and higher order ones. If we iterate the Darboux transformation with different spectral parameters $\lambda_{i}$, $i=1,2,\cdots,n$, then we obtain the multi-fold Darboux matrix. Otherwise, if we iterate it with the same spectral parameter, then the higher order Darboux matrix can be obtained. For the multi-fold Darboux matrix, it is purely algebraic procedure. On the other hand, for the higher order Darboux matrix, because the iteration is always performed on the same spectral point which is a moving singularity point for the higher order iteration, we should use the Laurent series expansion to continue the higher order iteration.

Firstly, we consider the multi-fold Darboux matrix and its Bäcklund transformation. Suppose we have $n$ different spectral parameters $\lambda_{i}$, $i=1,2,\cdots,n$ and its corresponding solutions $\Phi_{i}=\Phi(x,t;\lambda_{i})(1,c_{i})^{T}$, then the $n$ times of iteration can reduce as the $n$-fold Darboux matrix:

[TABLE]

where

[TABLE]

The corresponding Bäcklund transformation is

[TABLE]

The higher order Darboux matrix corresponds to the iterated Darboux matrix with the same spectral parameter, which is useful in constructing the higher order soliton and higher order rogue waves. In the elliptic-function background, there are two pairs of conjugated roots. Thus we consider the following higher order Darboux matrix with two different parameters $\lambda_{1}$ and $\lambda_{2}$. To construct the higher order Darboux matrix, we expand the vector solution $\Phi(x,t;\lambda)$ at the points $\lambda=\lambda_{1}$ and $\lambda=\lambda_{2}$:

[TABLE]

Now we construct the higher order Darboux matrix and the corresponding Bäcklund transformation based on $\Phi$. We consider the one from vector solution $\Phi$:

[TABLE]

where

[TABLE]

The corresponding Bäcklund transformation can be represented as

[TABLE]

The trivial Darboux matrix and the translation

We would like to discuss the solutions constructed from the Darboux transformation. Actually, even though the formula for the Darboux transformation itself is simple, it is not very easy to write down the formula in a compact form. In this subsection, we deal with this problem. Firstly, we give the proof that, in this case, the Darboux transformation is equivalent to the translation and phase gauge.

By the symmetry, we know that $(\mathrm{i}\sigma_{2})\Phi^{\rm T}(-\mathrm{i}\sigma_{2})$ satisfies the adjoint Lax pair. Thus the squared wave function can be constructed as the following way:

[TABLE]

On the other hand, we know the Darboux matrix

[TABLE]

Now we apply the Darboux matrix to construct a new squared wave function:

[TABLE]

Here we just need to verify that there is no singularity at point $\lambda_{1}$ and $\lambda_{1}^{*}$. It is seen that

[TABLE]

has no singularity at point $\lambda=\lambda_{1}.$ So is the squared wave function $\widehat{\mathbf{L}}(x,t;\lambda)$. Similarly, we know that

[TABLE]

which infers that the squared wave function $\widehat{\mathbf{L}}(x,t;\lambda)$ has no singularity at point $\lambda=\lambda_{1}^{*}.$ Therefore the squared wave function $\widehat{\mathbf{L}}(x,t;\lambda)$ is an entire function. It is readily to see that the squared wave function $\widehat{\mathbf{L}}(x,t;\lambda)$ has the same asymptotic behavior as $\mathbf{L}(x,t;\lambda).$ Meanwhile, by virtue of properties of the Darboux matrix, we obtain that $\det(\mathbf{L})=\det(\widehat{\mathbf{L}}).$ This is to say, the matrices $\mathbf{L}$ and $\widehat{\mathbf{L}}$ share the same structure. So the coefficients satisfy the same differential relation. Then the solutions transformed by the above special Darboux transformation is equivalent to the translation of coordinates and the phases. By direct calculation, the Darboux matrix can be rewritten as

[TABLE]

So the transformed squared wave function $\widehat{\mathbf{L}}(\lambda)$ can be rewritten as

[TABLE]

where

[TABLE]

and

[TABLE]

The matrix functions $\widehat{\mathbf{L}}_{1}$ and $\widehat{\mathbf{L}}_{0}$ shares the same differential equation as $\mathbf{L}_{1}$ and $\mathbf{L}_{0}$ as in equations (2.11) and (2.12).

The aim for this subsection is to prove the shift of crest and phase, if we choose the special vector solution $\Phi_{1}$ with $c_{1}=0$ or $c_{1}=\infty.$ The new obtained solution is an elliptic function solution with the shift of crest and phase. The explicit shift of crest and phase can be calculated by the theta function given in section 3.

Appendix B The formula of theta function

Here we need to use the integration formulas: (Ref. [11], 336.00-336.03)

[TABLE]

and

[TABLE]

and

[TABLE]

The elliptic integral of the second kind is given by

[TABLE]

where $\vartheta_{4}(u)$ and $Z(u)$ are the theta function and the Jacobi zeta function introduced in (B.6) and (B.7).

The elliptic integral of the third kind is given by (see Appendix B of Ref. [36] and 17.4.35 of Ref. [1] for derivation.)

[TABLE]

Theta functions and related formulas:

We define the theta function by [26]

[TABLE]

and

[TABLE]

We also write $\vartheta_{j}(u,q)=\vartheta_{j}(u|\tau)$ with the nome $q=\mathrm{e}^{\mathrm{i}\pi\tau}$, where $\tau=\mathrm{i}K^{\prime}/K$, $K=K(m)$, $K^{\prime}=K(1-m)$. The notation $\vartheta_{j}=\vartheta_{j}(0|\tau)$ is always used. The Jacobi elliptic functions in terms of theta functions are

[TABLE]

The elliptic parameters are given by $m=\frac{\vartheta_{2}^{4}}{\vartheta_{3}^{4}}$ and $m^{\prime}=1-m=\frac{\vartheta_{4}^{4}}{\vartheta_{3}^{4}}$.

We introduce the Jacobi zeta function:

[TABLE]

The three-term Weierstrass addition identities (or Fay’s identities) are given by [26]

[TABLE]

The shift formulas by half-periods are:

[TABLE]

Appendix C The second-second order rogue waves

In this appendix, we give a detailed calculation on the construction of second-second order rogue waves.

Choosing the parameters $b=\frac{\sqrt{3}}{2}$, $c=d=0.5$, and $\lambda=\lambda_{i}+\frac{\mathrm{i}\epsilon_{i}^{2}}{a_{i}^{2}}$, $a_{1}=\sqrt{-3-\mathrm{i}\sqrt{3}}$, $a_{2}=\sqrt{-3+\mathrm{i}\sqrt{3}}$, then we have

[TABLE]

where

[TABLE]

and

[TABLE]

Then the vector solutions can be expanded as

[TABLE]

and

[TABLE]

We just list the odd order coefficients for the above expansions, which would be used to construct the high order rogue waves:

[TABLE]

and

[TABLE]

and

[TABLE]

The $(2,2)$-rogue waves can be represented as

[TABLE]

where

[TABLE]

and

[TABLE]

and

[TABLE]

Acknowledgments

B.-F. acknowledges the partial support by NSF under Grant No. DMS-1715991, and National Natural Science Foundation of China under Grant No.11728103. Liming Ling is supported by National Natural Science Foundation of China(Contact No. 11771151), Guangdong Natural Science Foundation (Contact No. 2017A030313008), Guangzhou Science and Technology Program(No. 201707010040). The author Liming Ling thanks the Professor Peter D. Miller for his encourage, help, guidance and discussion during the scholar visiting of one year in the University of Michigan. Meanwhile, he thanks Dr. Deniz Bilman for his valuable discussions. The work of D.A.T. is supported by the Ministry of Education, Culture, Sports, Science (MEXT)-Supported Program for the Strategic Research Foundation at Private Universities “Topological Science” (Grant No. S1511006).