1School of Mathematics and Physics, North China Electric Power University, Beijing 102206, China2College of Science, China University of Petroleum, Beijing 102249, China3State Key Laboratory of Heavy Oil Processing, China University of Petroleum, Beijing 102249, China
Received July 29, 2019; Accepted March 10, 2020; Published April 16, 2020
Based on the Darboux transformation and N-soliton solutions, we obtain the explicit formulas of arbitrary-order multi-pole (MP) solutions of the Hirota equation via some limit technique. Then, by an improved asymptotic analysis method relying on the balance between exponential and algebraic terms, we derive the accurate expressions of all asymptotic solitons in the double- and triple-pole solutions. Moreover, we study the soliton interactions in the MP solutions and especially emphasize some unusual properties, like the interacting solitons separate from each other in a logarithmical law, the strengths of interaction forces decrease exponentially with their relative distances, and the constant–velocity asymptotic soliton does not experience a phase shift upon an interaction.
It is known that the N-soliton solution of an integrable nonlinear partial differential equation (NPDE) is associated to that the reflection coefficient admits N simple poles in the terminology of the inverse scattering transform (IST).1,2) When those simple poles merge (i.e., the multiplicity of discrete eigenvalues is higher than one), one can obtain the multi-pole (MP) solutions,1,3–5) in which the interacting solitons form weakly bound states. Early in 1972, Zakharov and Shabat1) first obtained the double-pole solution of the focusing nonlinear Schrödinger equation (NLSE), as the degeneration of the two-soliton solution when two distinct poles coalesce into one. Later on, Olmedilla3) derived the higher-order MP solutions of the NLSE and studied the asymptotic behavior of the double- and triple-pole solutions. Lately, Schiebold4) gave a rigorous and complete asymptotic description of an arbitrary MP solution based on some operator-theoretic approach. Meanwhile, the MP solutions (mostly in the double-pole case) have also been identified in the modified Korteweg–de Vries (MKdV) equation,6) sine-Gordon equation7) and many other integrable NPDEs.8–18) However, it does not say that the MP degeneration can occur for all types of N-soliton solutions, e.g., the MP solutions are not available for the multi-dark-soliton solutions of the defocusing NLSE.19)
Different from the usual N-soliton interactions, the MP solutions can exhibit a strong interaction in the near-field region, and the relative distance between two interacting solitons grows logarithmically with \(|t|\).1,3–5) This peculiar property has drawn certain attention to the MP solutions both in mathematics and physics. In addition, the existence of two-soliton solution with the logarithmic relative distance has also been shown in the nonintegrable subcritical and supercritical NLSEs.20) In the context of optical fibers, the MP solutions can describe the interactions of multiple chirped pulses with the same amplitudes and group velocities when they are initially input with no phase difference.21,22) Upon such an interaction, the force between two solitons is attractive and its strength decreases exponentially with their initial distance.21,22)
In this paper, we study the MP solution dynamics of the Hirota equation23) \begin{equation} i q_{t} + \frac{1}{2}q_{xx}+ |q|^{2}q + i \varepsilon q_{xxx} + 6i \varepsilon |q|^{2}q_{x} = 0, \end{equation} (1) which is an integrable higher-order generalization of NLSE. Equation (1) has primarily arisen from optical fibers to describe the femtosecond soliton pulse propagation in a single-mode fiber for a certain parameter regime,24,25) where x and t are, respectively, the retarded time and transverse spatial variable, the terms \(q_{xx}\), \(|q|^{2} q\), \(q_{xxx}\), \(|q|^{2} q_{x}\) represent the group velocity dispersion, self-phase modulation, third-order dispersion and self-steepening effects, respectively.26) Also, this equation is relevant for the vortex string motion for a three-dimensional Euler incompressible fluid.27,28) As the second member of the NLSE hierarchy, Eq. (1) is solvable by the IST and its integrability can be seen from the Painlevé property and Lax pair,29) an infinite number of conservation laws,30) Darboux transformation (DT),31,32) etc. In the past decades, the initial-value problem of Eq. (1) has been studied from different aspects, like the global existence and uniqueness of smooth solutions,33) long-time asymptotic behavior via the nonlinear steepest descent method,34) and global attractor with zero order dissipation.35) Meanwhile, researchers have constructed its exact solutions of physical importance (including the soliton solutions,29,32,36) breather solutions,36–38) and rogue-wave solutions36,37,39)), and have detailed their dynamical evolution behaviors as well.
We notice that although the double-pole solution of Eq. (1) was obtained in Refs. 13 and 28, the explicit formulas for the higher-order MP solutions have not been reported in the existing literature. In this work, we first employ the DT to construct the N-soliton solutions in the double-Wronskian form. On this basis, we obtain the N-pole solutions by some limit technique (which is usually combined with the DT in constructing the rogue-wave solutions36,37,40,41) and algebraic soliton solutions42,43)) when all the spectral parameters coalesce into the same one. Second, by an improved asymptotic analysis method44,45) which relies on the balance between the exponential and algebraic terms, we derive the accurate expressions of all asymptotic solitons in the double- and triple-pole solutions. Third, we study the soliton interactions in the MP solutions, emphasizing the difference from the multi-soliton interactions. Similarly to the NLSE,3,4,22) the MP solutions of Eq. (1) are found to display such properties: all the asymptotic solitons have the same amplitudes, the interacting solitons separate from each other in a logarithmical law, the strengths of interaction forces decrease exponentially with their relative distances, the constant–velocity asymptotic solitons experience no phase shift upon an interaction (see the asymptotic solitons \(u^{\pm}_{1}\) in the triple-pole solution).
It should be mentioned that very recently Ref. 46 constructed the MP solutions of Eq. (1) from Hirota's direct method as well as Darboux–Crum transformations, and detailed their asymptotic and scattering behavior. In comparison, the difference of our work lies in that the N-pole solutions are obtained explicitly in the determinant form, and the improved asymptotic analysis allows us to derive the asymptotic expressions of the MP solutions with arbitrary order. Moreover, we give the explicit formulas to measure the pairwise-soliton interaction forces, which can quantitatively explain why the MP solutions exhibit the absorb-emit scattering process.46)
Equivalently, Eq. (1) can be transformed into the complex MKdV equation \begin{equation} u_{\eta} + u_{\xi\xi\xi} + 6 |u|^{2} u_{\xi} = 0, \end{equation} (2) by the variable transformations (e.g., see Refs. 24 and 29) \begin{equation} \begin{split} &q(x,t) = u(\xi,\eta)\mathrm{e}^{i \left(\frac{x}{6\varepsilon}-\frac{t}{108 \varepsilon^{2}}\right)},\\ &\eta = \varepsilon t,\quad \xi = x-\frac{t}{12\varepsilon}. \end{split} \end{equation} (3) Thus, we turn to construct the exact MP solutions of Eq. (2) and study their described dynamical behavior. The structure of this paper is arranged as follows: In Sect. 2, we review the DT for the complex MKdV equation (2), and then obtain the N-pole solutions as the degeneration of N-soliton solutions. In Sects. 3 and 4, by an improved asymptotic analysis method, we derive the asymptotic solitons of the double- and triple-pole solutions, and reveal the relevant soliton interaction properties. In Sect. 5, we address the conclusions and discussions of this paper.
2. Darboux Transformation and N-pole Solutions
To begin with, we write the Lax pair of Eq. (2) in the form \begin{align} \Psi_{\xi} &= (\lambda U_{0} +U_{1})\Psi, \end{align} (4a) \begin{align} \Psi_{\eta} &= (\lambda^{3} V_{0} + \lambda^{2} V_{1} +\lambda V_{2} + V_{3})\Psi, \end{align} (4b) with \begin{align*} &U_{0} = \begin{pmatrix} 1 &\mathbf{0}\\ \mathbf{0}&-1 \end{pmatrix},\quad U_{1} = \begin{pmatrix} 0 &u\\ -\bar{u} &\mathbf{0}\end{pmatrix},\\ &V_{0} = - 4U_{0},\quad V_{1} = - 4U_{1},\notag\\ &V_{2} = - 2 \begin{pmatrix} |u|^{2} &u_{\xi}\\ \bar{u}_{\xi} &|u|^{2} \end{pmatrix},\\ &V_{3} = \begin{pmatrix} u \bar{u}_{\xi} -u_{\xi}\bar{u} &-u_{\xi\xi}-2|u|^{2} u\\ \bar{u}_{\xi\xi}+2|u|^{2}\bar{u} &\bar{u} u_{\xi} - \bar{u}_{\xi} u \end{pmatrix}, \end{align*} where \(\Psi=(f, g)^{T}\) (the superscript T signifies the vector transpose) is the vector eigenfunction, λ is the spectral parameter, and \(\bar{u}\) represents the complex conjugate of u. It can be checked that the compatibility condition \(\Psi_{\xi\eta}=\Psi_{\eta\xi}\) exactly yields Eq. (2).
The symmetry of Lax pair (4) shows that if \((f, g)^{T}\) is a nonzero solution of (4) with \(\lambda=\lambda_{1}\), then \((-\bar{g},\bar{f})^{T}\) also solves the same Lax pair at \(\lambda = -\bar{\lambda}_{1}\). Then, based on the work in Ref. 32, we present the N-fold DT of Eq. (2) [or equivalently Eq. (1)] as follows: \begin{align} &\Psi^{[N]} = T^{[N]}\Psi,\quad T^{[N]} = \begin{pmatrix} \lambda^{N} - \displaystyle\sum_{n=1}^{N}a_{n}(\xi,\eta) \lambda^{n-1} &- \displaystyle\sum_{n=1}^{N}b_{n}(\xi,\eta) (-\lambda)^{n-1}\\ - \displaystyle\sum_{n=1}^{N}c_{n}(\xi,\eta) \lambda^{n-1} &\lambda^{N} - \displaystyle\sum_{n=1}^{N}d_{n}(\xi,\eta)(-\lambda)^{n-1} \end{pmatrix}, \end{align} (5a) \begin{align} &u^{[N]} = u +2(-1)^{N-1} b_{N},\quad \bar{u}^{[N]} = \bar{u} + 2 c_{N}, \end{align} (5b) where the functions \(a_{n}\), \(b_{n}\), \(c_{n}\), and \(d_{n}\) (\(1\leq n\leq N\)) are uniquely determined by \begin{equation} \begin{split} &T^{[N]}|_{\lambda=\lambda_{k}}\Psi_{k} = \mathbf{0},\\ &T^{[N]}|_{\lambda=-\bar{\lambda}_{k}}\Phi_{k} = \mathbf{0}, \end{split} \end{equation} (6) with \(1\leq k\leq N\), \(\Psi_{k}=(f_{k}, g_{k})^{T}\) and \(\Phi_{k}=(-\bar{g}_{k},\bar{f}_{k})^{T}\) (\(1\leq k\leq N\)) being \(2N\) linearly-independent solutions of Lax pair (4) at different spectral parameters \(\lambda_{k}\)'s and \(-\bar{\lambda}_{k}\)'s. In particular, via Cramer's rule we can obtain \(b_{N}\) and \(c_{N}\) from Eq. (6) in the determinant form \begin{equation} b_{N} = (-1)^{N-1}\frac{\tau_{N+1,N-1}}{\tau_{N,N}},\quad c_{N} = \frac{\tau_{N-1,N+1}}{\tau_{N,N}}, \end{equation} (7) with \begin{equation} \tau_{M, L} = \begin{vmatrix} F_{N\times M} &G_{N\times L}\\ -\bar{G}_{N\times M} &\bar{F}_{N\times L} \end{vmatrix}\quad (M+L = 2N), \end{equation} (8) where \begin{align*} F_{N\times M}&=(\lambda_{k}^{m-1}f_{k})_{\substack{1{\leqslant} k{\leqslant} N,\\1{\leqslant} m{\leqslant} M}}, \\ {G_{N\times L}}&=[(-\lambda_{k})^{m-1}g_{k}]_{\substack{1{\leqslant} k{\leqslant} N,\\1{\leqslant} m{\leqslant} L}}, \\ {\bar{G}_{N\times M}}&=[(-\bar{\lambda}_{k})^{m-1}\bar{g}_{k}]_{\substack{1{\leqslant} k{\leqslant} N,\\1{\leqslant} m{\leqslant} M}}, \\ {\bar{F}_{N\times L}}&{{}=(\bar{\lambda}_{k}^{m-1}\bar{f}_{k})_{\substack{1{\leqslant} k{\leqslant} N,\\1{\leqslant} m{\leqslant} L}}.} \end{align*} It can be verified that the transformed eigenfunction \(\Psi^{[N]}=(f^{[N]}, g^{[N]})^{T}\) exactly satisfies the Lax pair (4) with the potential given by \(u^{[N]}\) in Eq. (5b), which implies that \(u^{[N]}\) solves the complex MKdV equation (2). In addition, we require \(\lambda_{k}\not\in\mathbb{R}\) for all \(1\leq k\leq N\) to avoid the trivial iteration of the DT.
Taking the trivial solution \(u=0\) as a seed, we implement the above DT-iterated algorithm and obtain the bright N-soliton solutions in the double-Wronskian form32) \begin{equation} u^{[N]} = 2\frac{\tau_{N+1,N-1}}{\tau_{N,N}},\quad \bar{u}^{[N]} = 2\frac{\tau_{N-1,N+1}}{\tau_{N,N}}, \end{equation} (9) in which \(\tau_{M, L}\) (\(M=N-1, N, N+1\)) are defined in Eq. (8) with \((f_{k}, g_{k})^{T}\) as the solution of Lax pair (4) at the zero potential: \begin{equation} \begin{pmatrix} f_{k}\\ g_{k} \end{pmatrix}= \begin{pmatrix} \alpha_{k} e^{\lambda_{k} \xi-4 \lambda_{k}^{3}\eta}\\ \beta_{k} e^{-\lambda_{k} \xi+4 \lambda_{k}^{3}\eta} \end{pmatrix}\quad (1\leq k\leq N), \end{equation} (10) where \(\alpha_{k},\beta_{k}\in \mathbb{C}\backslash\{0\}\) and \(\lambda_{k}\in \mathbb{C}\backslash\mathbb{R}\) are arbitrary parameters. In general, solutions (9) describe the standard elastic collisions among N solitons when \(\lambda_{k}\)'s are different from one another.
Considering the degeneracy \(\lambda_{1},\lambda_{2},\ldots,\lambda_{N}\to \chi\) with \(\chi\in \mathbb{C}\backslash\mathbb{R}\), the N-soliton solutions (9) can yield a chain of MP solutions which are expressible in the mixed exponential-algebraic form. According to the way in Ref. 40, we assume \(\lambda_{k}=\chi(1+\delta_{k})\) with \(\delta_{k}\) being small parameter, and introduce \(\delta_{k}\) into the constants \(\alpha_{k}\) and \(\beta_{k}\) in solution (10) in the form \begin{equation} \alpha_{k} = \alpha' e^{\sum_{l=1}^{\infty}s_{l}\delta_{k}^{l-1}},\quad \beta_{k} = \beta' e^{- \sum_{l=1}^{\infty}s_{l}\delta_{k}^{l-1}}, \end{equation} (11) where \(1\leq k\leq N\), \(\alpha'\), \(\beta'\) and \(s_{l}\)'s are parameters in \(\mathbb{C}\). Next, we expand all the elements of \(\tau_{M,L}\) in the Taylor series of \(\delta_{k}\) as follows: \begin{equation} \begin{split} \lambda_{k}^{m-1}f_{k} = \sum_{j=0}^{\infty}f^{(j,m-1)}(\xi,\eta)\delta_{k}^{j},\\ \lambda_{k}^{m-1}g_{k} = \sum_{j=0}^{\infty}g^{(j,m-1)}(\xi,\eta)\delta_{k}^{j}, \end{split} \end{equation} (12) where \(1\leq k\leq N\), \(1\leq m\leq\max\{M, L\}\), and the expansion coefficients \(f^{(j,m-1)}\) and \(g^{(j,m-1)}\) are defined by \begin{equation} \begin{split} f^{(j,m-1)} &= \frac{1}{j !}\frac{\partial^{j}(\lambda_{k}^{m-1} f_{k})}{\partial\lambda^{j}_{k}}|_{\delta_{k}=0},\\ g^{(j,m-1)} &= \frac{1}{j !}\frac{\partial^{j}(\lambda_{k}^{m-1} g_{k})}{\partial\lambda_{k}^{j}}|_{\delta_{k}=0}. \end{split} \end{equation} (13) Then, by substituting Eq. (12) into \(\tau_{M,L}\) and using the determinant properties, we have \begin{align} \tau_{M,L} &= \sum_{\substack{i_{1},\ldots, i_{N} {\in} \mathbb{Z}^{+};\\ j_{1}, \ldots, j_{N} {\in} \mathbb{Z}^{+}}}\prod_{k=1}^{N}\delta_{k}^{i_{k}-1}\bar{\delta}_{k}^{j_{k}-1}\notag\\ &\quad \times \begin{vmatrix} f^{(i_{1}-1,0)} &\cdots &f^{(i_{1}-1,M-1)} &g^{(i_{1}-1,0)} &\cdots &(-1)^{L-1}g^{(i_{1}-1,L-1)}\\ \vdots &\ddots &\vdots &\vdots &\vdots &\vdots\\ f^{(i_{N}-1,0)} &\cdots &f^{(i_{N}-1,M-1)} &g^{(i_{N}-1,0)} &\cdots &(-1)^{L-1}g^{(i_{N}-1,L-1)}\\ - \bar{g}^{(j_{1}-1,0)} &\cdots &(-1)^{M} \bar{g}^{(j_{1}-1,M-1)} &\bar{f}^{(j_{1}-1,0)} &\cdots &\bar{f}^{(j_{1}-1,L-1)}\\ \vdots &\ddots &\vdots &\vdots &\ddots &\vdots\\ - \bar{g}^{(j_{N}-1,0)} &\cdots &(-1)^{M} \bar{g}^{(j_{N}-1,M-1)} &\bar{f}^{(j_{N}-1,0)} &\cdots &\bar{f}^{(j_{N}-1,L-1)} \end{vmatrix}\notag\\ &= \Delta_{N} \begin{vmatrix} f^{(0,0)} &\cdots &f^{(0,M-1)} &g^{(0,0)} &\cdots &(-1)^{L-1}g^{(0,L-1)}\\ \vdots &\ddots &\vdots &\vdots &\vdots &\vdots\\ f^{(N-1,0)} &\cdots &f^{(N-1,M-1)} &g^{(N-1,0)} &\cdots &(-1)^{L-1}g^{(N-1,L-1)}\\ - \bar{g}^{(0,0)} &\cdots &(-1)^{M} \bar{g}^{(0,M-1)} &\bar{f}^{(0,0)} &\cdots &\bar{f}^{(0,L-1)}\\ \vdots &\ddots &\vdots &\vdots &\ddots &\vdots\\ - \bar{g}^{(N-1,0)} &\cdots &(-1)^{M} \bar{g}^{(N-1,M-1)} &\bar{f}^{(N-1,0)} &\cdots &\bar{f}^{(N-1,L-1)} \end{vmatrix}+ \text{small terms}, \end{align} (14) where \(\Delta_{N}\) is given by \begin{equation} \Delta_{N} = \sum_{\substack{i_{1}, \ldots, i_{N}{\in} [N];\\ j_{1}, \ldots, j_{N}{\in} [N]}}(-1)^{p_{i_{1}, \ldots, i_{N}}+p_{j_{1}, \ldots, j_{N}}}\prod_{k=1}^{N}\delta_{k}^{i_{k}-1} \bar{\delta}_{k}^{j_{k}-1}, \end{equation} (15) with \([N]:=\{1,\ldots,N\}\), \(i_{1}\neq\cdots\neq i_{N}\), \(j_{1}\neq\cdots\neq j_{N}\), and \begin{align*} p_{i_{1},\ldots,i_{N}} &= \begin{cases} 1, &\text{$(i_{1},\ldots,i_{N})$ is an odd permutation},\\ 0, &\text{$(i_{1},\ldots,i_{N})$ is an even permutation}, \end{cases}\\ p_{j_{1},\ldots,j_{N}} &= \begin{cases} 1, &\text{$(j_{1},\ldots,j_{N})$ is an odd permutation},\\ 0, &\text{$(j_{1},\ldots,j_{N})$ is an even permutation}. \end{cases}\end{align*}
Inserting (14) into Eq. (9) and taking the limit at \(\delta_{k}\to 0\), we finally derive the N-pole solutions as follows: \begin{equation} u^{[N]} = 2\frac{\tau'_{N+1,N-1}}{\tau'_{N,N}},\quad \bar{u}^{[N]} = 2\frac{\tau'_{N-1,N+1}}{\tau'_{N,N}}, \end{equation} (16) with \begin{equation} \tau'_{M, L} = \begin{vmatrix} F'_{N\times M} &G'_{N\times L}\\ -\bar{G}'_{N\times M} &\bar{F}'_{N\times L} \end{vmatrix}\quad (M+L = 2N), \end{equation} (17) where the block matrices \begin{align*} F'_{N\times M}&=(f^{(j,m-1)})_{\substack{0{\leqslant} j{\leqslant} {N{-}1},\\ 1{\leqslant} m{\leqslant} M}}, \\ {G'_{N\times L}}&=[(-1)^{m-1}g^{(j,m-1)}]_{\substack{0{\leqslant} j{\leqslant} {N{-}1},\\ 1{\leqslant} m{\leqslant} L}}, \\ \bar{G}'_{N\times M}&=[(-1)^{m-1}\bar{g}^{(j,m-1)}]_{\substack{0{\leqslant} j{\leqslant} {N{-}1},\\ 1{\leqslant} m{\leqslant} M}}, \\ {\bar{F}'_{N\times L}}&{{}= (\bar{f}^{(j,m-1)})_{\substack{0{\leqslant} j{\leqslant} {N{-}1},\\1{\leqslant} m{\leqslant} L}}.} \end{align*} Specifically, by calculating the Taylor expansions in (12) with induction, the coefficients \(f^{(j,m-1)}\) and \(g^{(j,m-1)}\) can be determined in the recursive form \begin{align} &f^{(0,0)} = \alpha' e^{\chi \xi-4 \chi^{3}\eta},\quad g^{(0,0)} = \beta' e^{-\chi \xi+4\chi^{3}\eta}, \end{align} (18a) \begin{align} &f^{(0,m-1)} = \chi^{m-1}f^{(0,0)},\quad g^{(0,m-1)} = \chi^{m-1}g^{(0,0)}\quad (2\leq m\leq N+1), \end{align} (18b) \begin{align} &f^{(j,0)} = \frac{1}{j}(\chi\xi -12\chi^{3}\eta+s_{1})f^{(j-1,0)}\\ &\quad\qquad-\frac{1}{j}(24\chi^{2}\eta-2s_{2})f^{(j-2,0)}- \frac{3}{j}(4\chi^{2}\eta-s_{3})f^{(j-3,0)}\notag\\ &\quad\qquad +\sum_{l=4}^{j} \frac{l}{j}s_{l} f^{(j-l,0)}\quad (1\leq j\leq N-1), \end{align} (18c) \begin{align} &g^{(j,0)} = - \frac{1}{j}(\chi\xi -12\chi^{3}\eta+s_{1})g^{(j-1,0)}\\&\quad\qquad+\frac{1}{j}(24\chi^{2}\eta-2s_{2}) g^{(j-2,0)} + \frac{3}{j}(4\chi^{2}\eta-s_{3})g^{(j-3,0)}\notag\\ &\quad\qquad-\sum_{l=4}^{j} \frac{l}{j}s_{l} g^{(j-l,0)} \quad (1\leq j\leq N-1), \end{align} (18d) \begin{align} &f^{(j,m-1)} = \chi^{m-1}f^{(j,0)}+\sum_{l=0}^{j-1}(-1)^{j-l+1}\frac{(m+j-l-2)!}{(j-l)!(m-2)!}f^{(l,m-1)}\notag\\ &\quad (1\leq j\leq N-1,\ 2\leq m\leq N+1), \end{align} (18e) \begin{align} &g^{(j,m-1)} = \chi^{m-1}g^{(j,0)}+\sum_{l=0}^{j-1}(-1)^{j-l+1}\frac{(m+j-l-2)!}{(j-l)!(m-2)!}g^{(l,m-1)}\notag\\ &\quad (1\leq j\leq N-1,\ 2\leq m\leq N+1). \end{align} (18f) Therefore, with substitution of Eqs. (18) into (16), one can obtain the explicit formula for an MP solution with arbitrary order \(N\geq 2\).
3. Double-pole Solution
Taking \(N=2\) in Eq. (16), we present the double-pole solution as follows: \begin{equation} u = \frac{-8\bar{\gamma}_{1}\mu e^{2\theta} [\chi (\bar{\zeta} -\bar{\chi}) e^{2\theta+2\bar{\theta}}-|\gamma_{1}|^{2}\bar{\chi} (\zeta +\chi)]}{|\chi |^{2} e^{4\theta+4\bar{\theta}} +2|\gamma_{1}|^{2}(2| \zeta |^{2}+ |\chi |^{2})e^{2\theta+2\bar{\theta}}+| \gamma_{1}|^{4} |\chi |^{2}}, \end{equation} (19) with \begin{align*} \zeta &= 2 \mu (\theta-8\chi^{3}\eta + s_{1}),\quad \bar{\zeta} = 2 \mu (\bar{\theta}-8\bar{\chi}^{3}\eta + \bar{s}_{1}),\\ \theta &= \chi \xi - 4 \chi^{3} \eta,\quad \bar{\theta} = \bar{\chi}\xi - 4 \bar{\chi}^{3} \eta,\quad \gamma_{1} = \frac{\beta'}{\alpha'}, \end{align*} where \(\mu=\text{Re}(\chi)\neq 0\), \(\nu=\text{Im}(\chi)\), \(s_{1}\) is an arbitrary constant in \(\mathbb{C}\). Next, via an improved asymptotic analysis method,44) we try to derive all the asymptotic solitons in solution (19), and then study its described two-soliton interaction dynamics.
Asymptotic analysis
In the following, we denote the real and imaginary parts of θ by \(\theta_{R}\) and \(\theta_{I}\), respectively. In view that \(\theta_{I}\) is expressible in terms of \(\theta_{R}\) by the relation \begin{equation} \theta_{I} = \frac{\nu}{\mu}\theta_{R}-8\nu(\mu^{2}+\nu^{2})\eta, \end{equation} (20) one can regard that solution (19) is explicitly dependent on \(\theta_{R}\) and η. This fact provides an important basis for us to get the accurate asymptotic expressions of solution (19).
First, we say that the asymptotic solitons of solution (19) cannot be located in any straight line \(\mathcal{L}:\mu\xi- c\tilde{\eta}=\text{const}\) with \(\tilde{\eta}=\text{sgn}(\mu)\eta\). From the difference \(\theta_{R}-(\mu\xi- c\tilde{\eta})=(c- 4 |\mu|\mu^{2} +12|\mu|\nu^{2})\tilde{\eta}\), we immediately have the asymptotic behavior of \(\theta_{R}\) along \(\mathcal{L}\) when \(\tilde{\eta}\rightarrow\pm \infty\): \begin{equation} \theta_{R} = \begin{cases} \pm \infty, &\text{$c > 4|\mu|(\mu^{2}-3\nu^{2})$},\\ \mathcal{O}(1), &\text{$c = 4|\mu|(\mu^{2}-3\nu^{2})$},\\ \mp \infty, &\text{$c < 4|\mu|(\mu^{2}-3\nu^{2})$}, \end{cases} \end{equation} (21) where \(\theta_{R}=\mathcal{O}(\eta)\) when \(c\neq 4|\mu|(\mu^{2}-3\nu^{2})\). For the three cases of \(\theta_{R}\) in (21), we obtain that solution (19) is dominated by \begin{equation} u \sim \begin{cases} - \dfrac{16\bar{\gamma}_{1} \mu^{2}(\bar{\theta}-8\bar{\chi}^{3}\eta)}{\bar{\chi} e^{2\bar{\theta}}}, &\text{$\theta_{R}\rightarrow +\infty$},\\ \dfrac{16\mu^{2}(\theta-8\chi^{3}\eta)}{\chi\gamma_{1} e^{-2\theta}}, &\text{$\theta_{R}\rightarrow -\infty$},\\ \dfrac{\bar{\gamma}_{1}\bar{\chi} e^{-2\bar{\theta}}}{\bar{\theta}-8\bar{\chi}^{3}\eta}-\dfrac{\chi e^{2\theta}}{\gamma_{1}(\theta-8\chi^{3}\eta)}, &\text{$\theta_{R} = \mathcal{O}(1)$}. \end{cases} \end{equation} (22) Noticing that \(|e^{2\bar{\theta}}|\gg|\bar{\theta}-8\bar{\chi}^{3}\eta|\) (\(\theta_{R}\rightarrow +\infty\)), \(|e^{-2\theta}|\gg|\theta-8\chi^{3}\eta|\) (\(\theta_{R}\rightarrow -\infty\)), and \(|\theta-8\chi^{3}\eta|\gg e^{\pm\theta_{R}}\) [\(\theta_{R} =\mathcal{O}(1)\)], thus all the limits in (22) are 0 as \(\tilde{\eta}\rightarrow\pm \infty\). That is, solution (19) has no asymptotic soliton lying in any straight line of the \(\xi\tilde{\eta}\) plane.
Next, we consider that the asymptotic solitons of solution (19) are located in some curves \(\Gamma: F(\xi,\tilde{\eta})=0\). Because the asymptotic limit of solution (19) is 0 when \(\theta_{R}=0\) or \(\mathcal{O}(1)\) [see the third case in (22)], we know that there must be \(|\theta_{R}|\to \infty\) along Γ as \(|\tilde{\eta}|\rightarrow \infty\). Meanwhile, we recall that solution (19) depends only on \(\tilde{\eta}\) and \(\theta_{R}\) with substitution of (20) for \(\theta_{I}\). Then, according to the inequality relations \(e^{\theta_{R}}\gg \theta_{R}\) (\(\theta_{R}\rightarrow +\infty\)) and \(e^{\theta_{R}}\ll\theta_{R}\) (\(\theta_{R}\rightarrow -\infty\)), there must exist some balance between \(\tilde{\eta}\) and \(e^{\theta_{R}}\) for an asymptotic soliton located in Γ. To obtain such balance relation explicitly, we assume that \begin{equation} \tilde{\eta} \sim W_{0} e^{2\alpha\theta_{R}}\quad (|\tilde{\eta}|\to \infty), \end{equation} (23) where \(W_{0}\) and α are two constants to be determined. With replacement of \(\theta_{I}\) in the algebraic terms by \(\theta_{R}\) via (20), the dominant behavior of solution (19) for different α as \(|\theta_{R}|\to \infty\) can be given by \begin{align} u \sim \begin{cases} \dfrac{\mu e^{2\theta}}{8|\mu|\gamma_{1}\kappa \tilde{\eta}}, &\text{$\alpha > 1$},\\ \dfrac{128|\mu| \mu\bar{\gamma}_{1}\bar{\kappa}\tilde{\eta}e^{2\theta}}{e^{4\theta_{R}} + 1024|\gamma_{1}|^{2}\mu^{2}|\kappa|^{2}\tilde{\eta}^{2}}, &\text{$\alpha = 1$},\\ \dfrac{128|\mu|\mu\bar{\gamma}_{1}\bar{\kappa}\tilde{\eta}}{e^{2\bar{\theta}}}, &\text{$0 < \alpha < 1$}, \end{cases}&\notag\\ (\theta_{R}\rightarrow +\infty), & \end{align} (24) and \begin{align} u \sim \begin{cases} - \dfrac{128|\mu|\mu\kappa\tilde{\eta}}{\gamma_{1} e^{-2\theta}}, &\text{$-1 < \alpha < 0$},\\ - \dfrac{128|\mu|\mu\bar{\gamma}_{1}\kappa\tilde{\eta} e^{2\theta}}{|\gamma_{1}|^{2} + 1024\mu^{2}|\kappa|^{2}\tilde{\eta}^{2}e^{4\theta_{R}}}, &\text{$\alpha = -1$},\\ - \dfrac{\mu\bar{\gamma}_{1} e^{-2\bar{\theta}}}{8|\mu|\bar{\kappa}\tilde{\eta}}, &\text{$\alpha < -1$}, \end{cases}&\notag\\ (\theta_{R}\rightarrow -\infty), & \end{align} (25) where \(\kappa=\chi^{2}+i\nu\bar{\chi}\). Apparently, it follows from (24) and (25) that u goes to 0 as \(\theta_{R}\rightarrow\pm \infty\) when \(\alpha\neq\pm 1\), in which no asymptotic soliton is available. As a result, the asymptotic solitons of solution (19) are formed only when the balance \(\tilde{\eta} e^{2\theta_{R}}\) or \(\tilde{\eta} e^{-2\theta_{R}}=\mathcal{O}(1)\) is met. Considering that such a balance may correspond to \(\tilde{\eta}\rightarrow +\infty\) and \(\tilde{\eta}\rightarrow -\infty\), we obtain that there are four asymptotic solitons in solution (19) and their expressions can be given as follows: \begin{align} u_{1}^{+} &= \frac{2\mu \bar{\gamma}_{1} \bar{\kappa} e^{2 i \theta_{I}}}{|\gamma_{1}| |\kappa|}\text{sech} [2 \theta_{R}-\ln(32 |\mu| | \gamma_{1}| | \kappa| \tilde{\eta})]\notag\\ &\quad (\tilde{\eta} e^{-2\theta_{R}} = \mathcal{O}(1),\ \tilde{\eta}\rightarrow + \infty), \end{align} (26a) \begin{align} u_{1}^{-} &= \frac{2 \mu \bar{\gamma}_{1} \kappa e^{2 i \theta_{I}}}{|\gamma_{1}| |\kappa|}\text{sech} \left[2 \theta_{R}+\ln\left(-\frac{32 |\mu| | \kappa | \tilde{\eta}}{| \gamma_{1}|}\right) \right]\notag\\ &\quad (\tilde{\eta} e^{2\theta_{R}} = \mathcal{O}(1),\ \tilde{\eta}\rightarrow -\infty), \end{align} (26b) \begin{align} u_{2}^{+} &= - \frac{2 \mu \bar{\gamma}_{1} \kappa e^{2 i \theta_{I}}}{|\gamma_{1}| |\kappa|}\text{sech} \left[2 \theta_{R}+\ln\left(\frac{32 |\mu| | \kappa | \tilde{\eta}}{| \gamma_{1}|}\right) \right]\notag\\ &\quad (\tilde{\eta} e^{2\theta_{R}} = \mathcal{O}(1),\ \tilde{\eta}\rightarrow +\infty), \end{align} (26c) \begin{align} u_{2}^{-} &= - \frac{2 \mu \bar{\gamma}_{1} \bar{\kappa} e^{2 i \theta_{I}}}{|\gamma_{1}| |\kappa|}\text{sech} [2 \theta_{R} - \ln(-32|\mu||\gamma_{1}| |\kappa|\tilde{\eta})]\notag\\ &\quad (\tilde{\eta} e^{-2\theta_{R}} = \mathcal{O}(1),\ \tilde{\eta}\rightarrow -\infty), \end{align} (26d) where we have labeled the four asymptotic solitons based on that the slope of center trajectory of \(u^{+}_{i}\) at \(\tilde{\eta}>0\) is the same as that of \(u^{-}_{i}\) at \(\tilde{\eta}<0\) (\(i=1,2\)), as given in Eq. (28).
As seen from Eqs. (26a)–(26d), \(u_{i}^{\pm}\) (\(i=1,2\)) represent four bright asymptotic solitons with their center trajectories given by the curves \begin{align} &\Gamma^{+}_{1}{:}\ \tilde{\eta} = \frac{e^{2\theta_{R}}}{32 |\mu| | \gamma_{1}| |\kappa|}\quad (\tilde{\eta} > 0), \end{align} (27a) \begin{align} &\Gamma^{-}_{1}{:}\ \tilde{\eta} = - \frac{| \gamma_{1}|e^{-2\theta_{R}}}{32 |\mu||\kappa|}\quad (\tilde{\eta} < 0), \end{align} (27b) \begin{align} &\Gamma^{+}_{2}{:}\ \tilde{\eta} = \frac{| \gamma_{1}|e^{-2\theta_{R}}}{32 |\mu ||\kappa|}\quad (\tilde{\eta} > 0), \end{align} (27c) \begin{align} &\Gamma^{-}_{2}{:}\ \tilde{\eta} = - \frac{e^{2\theta_{R}}}{32 |\mu| | \gamma_{1}| |\kappa|}\quad (\tilde{\eta} < 0), \end{align} (27d) along which \(|u_{i}^{\pm}|\) (\(i=1,2\)) reach their respective maxima. Correspondingly, the four curves' slopes in the \(\xi\tilde{\eta}\) plane can be obtained by \begin{align} K^{+}_{1} &= \frac{2\mu}{8|\mu|(\mu^{2}-3\nu^{2})+1/\tilde{\eta}}\quad (\tilde{\eta} > 0), \end{align} (28a) \begin{align} K^{-}_{1} &= \frac{2\mu}{8|\mu|(\mu^{2}-3\nu^{2})-1/\tilde{\eta}}\quad (\tilde{\eta} < 0), \end{align} (28b) \begin{align} K^{+}_{2} &= \frac{2\mu}{8|\mu|(\mu^{2}-3\nu^{2})-1/\tilde{\eta}}\quad (\tilde{\eta} > 0), \end{align} (28c) \begin{align} K^{-}_{2} &= \frac{2\mu}{8|\mu|(\mu^{2}-3\nu^{2})+1/\tilde{\eta}}\quad (\tilde{\eta} < 0). \end{align} (28d) It is shown from (28) that \(K^{+}_{1}\lessgtr\frac{1}{4\text{sgn}(\mu)(\mu^{2}-3\nu^{2})}\) when \(\tilde{\eta}>\max\left\{0,\frac{1}{8|\mu|(3\nu^{2}-\mu^{2})}\right\}\) for \(\mu\gtrless 0\), \(K^{-}_{1}\lessgtr \frac{1}{4\text{sgn}(\mu)(\mu^{2}-3\nu^{2})}\) when \(\tilde{\eta}<\min\left\{0,\frac{1}{8|\mu|(\mu^{2}-3\nu^{2})}\right\}\) for \(\mu\gtrless 0\), \(K^{+}_{2}\gtrless \frac{1}{4\text{sgn}(\mu)(\mu^{2}-3\nu^{2})}\) when \(\tilde{\eta}>\max\left\{0,\frac{1}{8|\mu|(\mu^{2}-3\nu^{2})}\right\}\) for \(\mu\gtrless 0\), \(K^{-}_{2}\gtrless \frac{1}{4\text{sgn}(\mu)(\mu^{2}-3\nu^{2})}\) when \(\tilde{\eta}<\min\left\{0,\frac{1}{8|\mu|(3\nu^{2}-\mu^{2})}\right\}\) for \(\mu\gtrless 0\). Letting \(l_{1}\) represent the direction \(\xi-4\text{sgn}(\mu)(\mu^{2}-3\nu^{2})\tilde{\eta}=0\), we can acquire the distribution of four asymptotic solitons of solution (19) in the \(\xi\tilde{\eta}\) plane as follows:
(i)
If \(\mu>0\) and \(\mu^{2}-3\nu^{2}>0\), the asymptotic solitons \(u^{\pm}_{1}\) is situated between the direction \(l_{1}\) and \(\pm\xi\) axis, \(u^{\pm}_{2}\) is between the direction \(l_{1}\) and \(\pm\tilde{\eta}\) axis.
(ii)
If \(\mu>0\) and \(\mu^{2}-3\nu^{2}<0\), the asymptotic solitons \(u^{\pm}_{1}\) is situated between the direction \(l_{1}\) and \(\pm\tilde{\eta}\) axis, \(u^{\pm}_{2}\) is between the direction \(l_{1}\) and \(\mp\xi\) axis.
(iii)
If \(\mu<0\) and \(\mu^{2}-3\nu^{2}>0\), the asymptotic solitons \(u^{\pm}_{1}\) is situated between the direction \(l_{1}\) and \(\mp\xi\) axis, \(u^{\pm}_{2}\) is between the direction \(l_{1}\) and \(\pm\tilde{\eta}\) axis.
(iv)
If \(\mu<0\) and \(\mu^{2}-3\nu^{2}<0\), the asymptotic solitons \(u^{\pm}_{1}\) is situated between the direction \(l_{1}\) and \(\pm\tilde{\eta}\) axis, \(u^{\pm}_{2}\) is situated between the direction \(l_{1}\) and \(\pm\xi\) axis.
(v)
If \(\mu>0\) and \(\mu^{2}-3\nu^{2}= 0\), the asymptotic solitons \(u^{+}_{1}\), \(u^{+}_{2}\), \(u^{-}_{1}\), and \(u^{-}_{2}\) are distributed counterclockwise in the four quadrants of the \(\xi\tilde{\eta}\) plane.
(vi)
If \(\mu<0\) and \(\mu^{2}-3\nu^{2}= 0\), the asymptotic solitons \(u^{+}_{2}\), \(u^{+}_{1}\), \(u^{-}_{2}\), and \(u^{-}_{1}\) are distributed counterclockwise in the four quadrants of the \(\xi\tilde{\eta}\) plane.
Since the slopes \(K_{i}^{\pm}\) (\(i=1,2\)) all tend to the same value \(\frac{1}{4\text{sgn}(\mu)(\mu^{2}-3\nu^{2})}\) asymptotically as \(|\tilde{\eta}|\to\infty\), all the center trajectories of \(u_{i}^{\pm}\) (\(i=1,2\)) are approximately parallel to the line \(l_{1}\) when \(|\tilde{\eta}|\gg 1\). Particularly for \(\mu^{2}-3\nu^{2} = 0\), the soliton center trajectories tend to be asymptotically parallel to the \(\tilde{\eta}\) axis. To illustrate, we plot the double-pole solution (19) and mark the four asymptotic solitons in Figs. 1 and 2, which are respectively associated with the cases \(\mu>0\), \(\mu^{2}-3\nu^{2}>0\) and \(\mu>0\), \(\mu^{2}-3\nu^{2}<0\).
Figure 1. (Color online) (a) 3D plot of the double-pole solution (19) with \(\chi=\frac{1}{2}-\frac{i}{4}\), \(\alpha'=1\), \(\beta'=1+i\), and \(s_{1}=1-\frac{i}{2}\). (b) Contour plot of solution (19) with the same parameters.
Figure 2. (Color online) (a) 3D plot of the double-pole solution (19) with \(\chi=\frac{1}{2}-\frac{i}{3}\), \(\alpha'=1\), \(\beta'=1+i\), and \(s_{1}=1-\frac{i}{2}\). (b) Contour plot of solution (19) with the same parameters.
To test the validity of our asymptotic analysis, we compare the exact solution (19) with the asymptotic expressions (26) at different values of η. As shown in Fig. 3, the asymptotic solitons \(u_{i}^{\pm}\) (\(i=1,2\)) have a good agreement with the exact solution in the far-field region of the \(\xi\tilde{\eta}\) plane, which indicates that Eqs. (26a)–(26d) are valid in describing the asymptotic solitons of solution (19).
Figure 3. (Color online) (a) Comparison of the asymptotic solitons \(u_{1}^{+}\) (blue dashed) and \(u_{2}^{+}\) (red dashed) with the exact solution (19) (black solid). (b) Comparison of the asymptotic solitons \(u_{1}^{-}\) (blue dashed) and \(u_{2}^{-}\) (red dashed) with the exact solution (19) (black solid). The relevant parameters take the same values as those in Fig. 1.
Properties of two-soliton interactions
Based on the asymptotic expressions in Eqs. (26a)–(26d), we discuss the two-soliton interactions in solution (19) from the following aspects:
(i)
The maximum values that \(|u_{1,2}^{\pm}|\) reach along the curves \(\Gamma^{\pm}_{1,2}\) show that the four asymptotic solitons have the same amplitudes \(A^{\pm}_{1,2} = 2|\mu|\).
(ii)
Taking the inverse of the slopes in Eq. (28), we obtain that the asymptotic solitons \(u^{\pm}_{i}\) (\(1\leq i\leq 2\)) have the \(\tilde{\eta}\)-dependent velocities: \begin{align} v_{1}^{+} &= \frac{8|\mu|(\mu^{2}-3\nu^{2}) +1/\tilde{\eta}}{2\mu}, \end{align} (29a) \begin{align} v_{1}^{-} &= \frac{8|\mu|(\mu^{2}-3\nu^{2}) -1/\tilde{\eta}}{2 \mu}, \end{align} (29b) \begin{align} v_{2}^{+} &= \frac{8|\mu|(\mu^{2}-3\nu^{2}) -1/\tilde{\eta}}{2\mu}, \end{align} (29c) \begin{align} v_{2}^{-} &= \frac{8|\mu|(\mu^{2}-3\nu^{2}) +1/\tilde{\eta}}{2\mu}. \end{align} (29d) It can be seen that the velocities \(v^{+}_{1,2}\) at \(\tilde{\eta}>0\) are respectively equal to \(v^{-}_{1,2}\) at \(\tilde{\eta}<0\), and that all the velocities slowly approach a common limit \(4\text{sgn}(\mu)(\mu^{2}-3\nu^{2})\) as \(|\tilde{\eta}|\to\infty\).
(iii)
The phase shifts between the envelopes of \(u_{i}^{+}\) and \(u_{i}^{-}\) (\(i=1, 2\)) are also \(\tilde{\eta}\)-dependent and can be explicitly given by \begin{equation} \delta\phi_{1} = - \delta\phi_{2} = - 2\ln(32 |\mu||\kappa||\tilde{\eta}|), \end{equation} (30) which exhibits the monotonic variation with \(|\tilde{\eta}|\) in a logarithmical manner (see Fig. 4).
(iv)
By calculating the position difference between \(u_{1}^{+}\) and \(u_{2}^{+}\) (or between \(u_{1}^{-}\) and \(u_{2}^{-}\)), we obtain the two-soliton relative distance as follows: \begin{equation} d = \frac{\ln(32|\mu| |\kappa| |\tilde{\eta}|)}{|\mu|}\quad \left(|\tilde{\eta}| \geq \frac{1}{32|\mu| |\kappa|}\right), \end{equation} (31) which increases logarithmically with \(|\tilde{\eta}|\). Meanwhile, we use the second derivative of d with respect to \(|\tilde{\eta}|\) (i.e., the acceleration) \begin{equation} a = - 1024|\mu||\kappa|^{2} e^{-2|\mu|d}, \end{equation} (32) to measure the interaction force between two solitons. In view of \(a<0\), the force is attractive and its strength decays exponentially to 0 as the relative distance increases (see Fig. 5).
Figure 4. (Color online) The phase shifts between the envelopes of \(u_{i}^{+}\) and \(u_{i}^{-}\) (\(i=1,2\)) versus \(|\tilde{\eta}|\), where the parameters are the same as those in Fig. 1.
Figure 5. (Color online) The two-soliton relative distance (red) and interaction force (blue) versus \(|\tilde{\eta}|\), where the parameters are the same as those in Fig. 1.
4. Triple-pole Solution
By taking \(N=3\) in Eq. (16), we can obtain the explicit triple-pole solution (which is omitted here for saving space). In a similar way, we will make an asymptotic analysis of the triple-pole solution, and discuss the relevant three-soliton interaction dynamics.
Asymptotic analysis
Like the procedure in Sect. 3.1, we first consider the possibility of asymptotic solitons lying in some straight lines. By using Eq. (21), we obtain the asymptotic limits of solution (16) with \(N=3\) along \(\mu\xi- c\tilde{\eta}=\text{const}\) as \(\tilde{\eta}\rightarrow\pm \infty\) as follows: \begin{equation} u \sim \begin{cases} 0, &\text{$c \neq 4|\mu|(\mu^{2}-3\nu^{2})$},\\ \dfrac{4 \mu \bar{\gamma}_{1} e^{2 \theta}}{e^{4\theta_{R}}+|\gamma_{1}|^{2}}, &\text{$c = 4|\mu|(\mu^{2}-3\nu^{2})$}. \end{cases} \end{equation} (33) Here, the second case shows that there are two asymptotic solitons of the same form \begin{align} u_{1}^{+} = u_{1}^{-} = \frac{2\mu \bar{\gamma}_{1} e^{2i\theta_{I}}}{|\gamma_{1}|}\text{sech} [2\theta_{R}-\ln(|\gamma_{1}|)] &\notag\\ (\tilde{\eta}\to \pm \infty), & \end{align} (34) which are localized in the straight line \(\mu\xi-4|\mu|(\mu^{2}-3\nu^{2})\tilde{\eta}=\text{const}\).
Next, we continue to work out the asymptotic solitons along some curves in the triple-pole solution. By virtue of (20), we replace \(\theta_{I}\) with \(\theta_{R}\) in the algebraic terms of the solution, and then compare the dominant balance between the terms \(\tilde{\eta}\) and \(e^{\theta_{R}}\). It turns out that the asymptotic solitons occur as \(|\tilde{\eta}|\rightarrow\infty\) only for the dominant balance \(\tilde{\eta}\sim e^{\theta_{R}}\) or \(\tilde{\eta}\sim e^{-\theta_{R}}\). Therefore, we derive the expressions of other four asymptotic solitons as follows: \begin{align} u_{2}^{+} &= \frac{2 \mu\bar{\gamma}_{1}{\bar{\kappa}}^{2} e^{2 i \theta_{I}}}{|\gamma_{1}| |\kappa|^{2}}\text{sech} [2 \theta_{R}-2\ln(16 \sqrt{2|\gamma_{1}|}|\mu||\kappa| \tilde{\eta})],\notag\\ &\quad (\tilde{\eta} e^{-\theta_{R}} = \mathcal{O}(1),\ \tilde{\eta}\rightarrow +\infty), \end{align} (35a) \begin{align} u_{2}^{-} &= \frac{2 \mu\bar{\gamma}_{1}\kappa^{2} e^{2 i \theta_{I}}}{|\gamma_{1}||\kappa|^{2}}\text{sech} \left[2 \theta_{R}+2\ln\left(-\frac{16\sqrt{2}|\mu| |\kappa| \tilde{\eta}}{\sqrt{| \gamma_{1}|}}\right) \right],\notag\\ &\quad (\tilde{\eta} e^{\theta_{R}} = \mathcal{O}(1),\ \tilde{\eta}\rightarrow -\infty), \end{align} (35b) \begin{align} u_{3}^{+} &= \frac{2 \mu\bar{\gamma}_{1}\kappa^{2} e^{2 i \theta_{I}}}{|\gamma_{1}| |\kappa|^{2}}\text{sech} \left[2 \theta_{R}+2\ln\left(\frac{16\sqrt{2}|\mu| |\kappa|\tilde{\eta}}{\sqrt{| \gamma_{1}|}}\right) \right],\notag\\ &\quad (\tilde{\eta} e^{\theta_{R}} = \mathcal{O}(1),\ \tilde{\eta}\rightarrow +\infty), \end{align} (35c) \begin{align} u_{3}^{-} &= \frac{2 \mu\bar{\gamma}_{1}\bar{\kappa}^{2} e^{2 i \theta_{I}}}{|\gamma_{1}| |\kappa|^{2}}\text{sech} [2 \theta_{R}-2\ln(- 16\sqrt{2|\gamma_{1}|}|\mu||\kappa|\tilde{\eta})],\notag\\ &\quad (\tilde{\eta} e^{-\theta_{R}} = \mathcal{O}(1),\ \tilde{\eta}\rightarrow -\infty), \end{align} (35d) where we have also labeled the four asymptotic solitons in view that the slope of center trajectory of \(u^{+}_{i}\) at \(\tilde{\eta}>0\) is the same as that of \(u^{-}_{i}\) at \(\tilde{\eta}<0\) (\(i=2,3\)), as given in Eqs. (37c)–(37e).
Based on Eqs. (34) and (35), we know that all the asymptotic solitons \(u_{i}^{\pm}\) (\(1\leq i\leq 3\)) are of the bright type, and that the center trajectories of those asymptotic solitons can be given by \begin{align} &\Gamma^{+}_{1}{:}\ e^{\theta_{R}} = \sqrt{|\gamma_{1}|} &\quad &(\tilde{\eta} > 0), \end{align} (36a) \begin{align} &\Gamma^{-}_{1}{:}\ e^{\theta_{R}} = \sqrt{|\gamma_{1}|} &\quad &(\tilde{\eta} < 0), \end{align} (36b) \begin{align} &\Gamma^{+}_{2}{:}\ \tilde{\eta} = \frac{e^{\theta_{R}}}{16\sqrt{2|\gamma_{1}|}|\mu||\kappa|} &\quad &(\tilde{\eta} > 0), \end{align} (36c) \begin{align} &\Gamma^{-}_{2}{:}\ \tilde{\eta} = - \frac{\sqrt{|\gamma_{1}|}e^{-\theta_{R}}}{16\sqrt{2}|\mu ||\kappa|} &\quad &(\tilde{\eta} < 0), \end{align} (36d) \begin{align} &\Gamma^{+}_{3}{:}\ \tilde{\eta} = \frac{\sqrt{|\gamma_{1}|}e^{-\theta_{R}}}{16\sqrt{2} |\mu ||\kappa|} &\quad &(\tilde{\eta} > 0), \end{align} (36e) \begin{align} &\Gamma^{-}_{3}{:}\ \tilde{\eta} = - \frac{e^{\theta_{R}}}{16\sqrt{2|\gamma_{1}|}|\mu||\kappa|} &\quad &(\tilde{\eta} < 0), \end{align} (36f) whose slopes are respectively as follows: \begin{align} &K^{+}_{1} = \frac{\mu}{4|\mu|(\mu^{2}-3\nu^{2})} &\quad &(\tilde{\eta} > 0), \end{align} (37a) \begin{align} &K^{-}_{1} = \frac{\mu}{4|\mu|(\mu^{2}-3\nu^{2})} &\quad &(\tilde{\eta} < 0), \end{align} (37b) \begin{align} &K^{+}_{2} = \frac{\mu}{4|\mu|(\mu^{2} -3\nu^{2}) +1/\tilde{\eta}} &\quad &(\tilde{\eta} > 0), \end{align} (37c) \begin{align} &K^{-}_{2} = \frac{\mu}{4|\mu|(\mu^{2} -3\nu^{2})-1/\tilde{\eta}} &\quad &(\tilde{\eta} < 0), \end{align} (37d) \begin{align} &K^{+}_{3} = \frac{\mu}{4|\mu|(\mu^{2} -3\nu^{2}) -1/\tilde{\eta}} &\quad &(\tilde{\eta} > 0), \end{align} (37e) \begin{align} &K^{-}_{3} = \frac{\mu}{4|\mu|(\mu^{2} -3\nu^{2})+1/\tilde{\eta}} &\quad &(\tilde{\eta} < 0). \end{align} (37f) It is suggested by Eqs. (37a)–(37f) that \(K^{+}_{2}\lessgtr K^{+}_{1}\) when \(\tilde{\eta}>\max\left\{0,\frac{1}{4|\mu|(3\nu^{2}-\mu^{2})}\right\}\) for \(\mu\gtrless 0\), \(K^{-}_{2}\lessgtr K^{-}_{1}\) when \(\tilde{\eta}<\min\left\{0,\frac{1}{4|\mu|(\mu^{2}-3\nu^{2})}\right\}\) for \(\mu\gtrless 0\), \(K^{+}_{3}\gtrless K^{+}_{1}\) when \(\tilde{\eta}>\max\left\{0,\frac{1}{4|\mu|(\mu^{2}-3\nu^{2})}\right\}\) for \(\mu\gtrless 0\), \(K^{-}_{3}\gtrless K^{-}_{1}\) when \(\tilde{\eta}<\min\left\{0,\frac{1}{4|\mu|(3\nu^{2}-\mu^{2})}\right\}\) for \(\mu\gtrless 0\). On this basis, we can determine the distribution of six asymptotic solitons of the triple-pole solution in the \(\xi\tilde{\eta}\) plane as follows:
(i)
If \(\mu>0\) and \(\mu^{2}-3\nu^{2}>0\), the asymptotic solitons \(u^{+}_{2}\), \(u^{+}_{1}\), and \(u^{+}_{3}\) are distributed counterclockwise in the first quadrant, \(u^{-}_{2}\), \(u^{-}_{1}\), and \(u^{-}_{3}\) are distributed counterclockwise in the third quadrant.
(ii)
If \(\mu>0\) and \(\mu^{2}-3\nu^{2}<0\), the asymptotic solitons \(u^{+}_{2}\), \(u^{+}_{1}\), and \(u^{+}_{3}\) are distributed counterclockwise in the second quadrant, \(u^{-}_{2}\), \(u^{-}_{1}\), and \(u^{-}_{3}\) are distributed counterclockwise in the fourth quadrant.
(iii)
If \(\mu<0\) and \(\mu^{2}-3\nu^{2}>0\), the asymptotic solitons \(u^{+}_{3}\), \(u^{+}_{1}\), and \(u^{+}_{2}\) are distributed counterclockwise in the second quadrant, \(u^{-}_{3}\), \(u^{-}_{1}\), and \(u^{-}_{2}\) are distributed counterclockwise in the fourth quadrant.
(iv)
If \(\mu<0\) and \(\mu^{2}-3\nu^{2}<0\), the asymptotic solitons \(u^{+}_{3}\), \(u^{+}_{1}\), and \(u^{+}_{2}\) are distributed counterclockwise in the first quadrant, \(u^{-}_{3}\), \(u^{-}_{1}\), and \(u^{-}_{2}\) are distributed counterclockwise in the third quadrant.
(v)
If \(\mu>0\) and \(\mu^{2}-3\nu^{2}= 0\), the asymptotic solitons \(u^{+}_{2}\), \(u^{+}_{3}\), \(u^{-}_{2}\), and \(u^{-}_{3}\) are distributed counterclockwise in the four quadrants, whereas the asymptotic solitons \(u^{\pm}_{1}\) are along the \(\pm\tilde{\eta}\) axis.
(vi)
If \(\mu<0\) and \(\mu^{2}-3\nu^{2}= 0\), the asymptotic solitons \(u^{+}_{3}\), \(u^{+}_{2}\), \(u^{-}_{3}\), and \(u^{-}_{2}\) are distributed counterclockwise in the four quadrants, whereas the asymptotic solitons \(u^{\pm}_{1}\) are along the \(\pm\tilde{\eta}\) axis.
Asymptotically as \(|\tilde{\eta}|\to\infty\) the slopes \(K_{2,3}^{\pm}\) tend to \(K_{1}^{\pm} =\frac{1}{4\text{sgn}(\mu)(\mu^{2}-3\nu^{2})}\), so that the center trajectories of \(u_{2,3}^{\pm}\) are approximately parallel to those of \(u_{1}^{\pm}\) when \(|\tilde{\eta}|\gg 1\). Particularly for \(\mu^{2}-3\nu^{2} = 0\), all the soliton center trajectories are asymptotically parallel to the \(\tilde{\eta}\) axis. In Figs. 6 and 7, we present two examples for the distribution of six asymptotic solitons in the triple-pole solution, which are respectively associated with the cases \(\mu>0\), \(\mu^{2}-3\nu^{2}>0\) and \(\mu>0\), \(\mu^{2}-3\nu^{2}<0\).
Figure 6. (Color online) (a) 3D plot of the triple-pole solution with \(\chi=\frac{1}{2}-\frac{i}{4}\), \(\alpha'=1\), \(\beta'=1+i\), \(s_{1}=s_{2}=1-\frac{i}{2}\). (b) Contour plot of the triple-pole solution with the same parameters.
Figure 7. (Color online) (a) 3D plot of the triple-pole solution with \(\chi=\frac{1}{2}-\frac{i}{3}\), \(\alpha'=1\), \(\beta'=1+i\), \(s_{1}=s_{2}=1-\frac{i}{2}\). (b) Contour plot of the triple-pole solution with the same parameters.
To further show the validity of the above asymptotic analysis, we also make a comparison of the exact triple-pole solution with the asymptotic expressions (34) and (35) in Fig. 8. It can be seen that the asymptotic expressions give a good approximation to the exact solution at large values of \(|\eta|\).
Figure 8. (Color online) (a) Comparison of the asymptotic solitons \(u_{1}^{+}\) (red dashed), \(u_{2}^{+}\) (blue dashed), and \(u_{3}^{+}\) (purple dashed) with the exact triple-pole solution (black solid). (b) Comparison of the asymptotic solitons \(u_{1}^{-}\) (red dashed), \(u_{2}^{-}\) (blue dashed), and \(u_{3}^{-}\) (purple dashed) with the exact triple-pole solution (black solid). The relevant parameters take the same values as those in Fig. 6.
Properties of three-soliton interactions
Similarly to Sect. 3.2, based on the asymptotic expressions in Eqs. (34) and (35), we also discuss the three-soliton interactions in the triple-pole solution from the following four aspects:
(i)
Like the double-pole case, all the asymptotic solitons \(u_{i}^{\pm}\) (\(1\leq i\leq 3\)) have the same amplitudes \(2|\mu|\).
(ii)
By taking the inverse of the slopes in Eqs. (37a)–(37f), the velocities of six asymptotic solitons are presented as follows: \begin{align} v_{1}^{+} &= v_{1}^{-} = \frac{4|\mu|(\mu^{2} -3\nu^{2})}{\mu}, \end{align} (38a) \begin{align} v_{2}^{+} &= \frac{4|\mu|(\mu^{2} -3\nu^{2}) +1/\tilde{\eta}}{\mu}, \end{align} (38b) \begin{align} v_{2}^{-} &= \frac{4|\mu|(\mu^{2} -3\nu^{2}) -1/\tilde{\eta}}{\mu}, \end{align} (38c) \begin{align} v_{3}^{+} &= \frac{4|\mu|(\mu^{2} -3\nu^{2}) -1/\tilde{\eta}}{\mu}, \end{align} (38d) \begin{align} v_{3}^{-} &= \frac{4|\mu|(\mu^{2} -3\nu^{2}) +1/\tilde{\eta}}{\mu}, \end{align} (38e) where \(v^{+}_{1}=v^{-}_{1}\) are constants, whereas \(v^{\pm}_{2,3}\) slowly approaches \(v^{\pm}_{1}\) as \(|\tilde{\eta}|\to\infty\). It should be noted that the velocities \(v^{+}_{2,3}\) at \(\tilde{\eta}>0\) are, respectively, equal to \(v^{-}_{2,3}\) at \(\tilde{\eta}<0\) although they are \(\tilde{\eta}\)-dependent.
(iii)
There is no phase shift between the envelopes of \(u_{1}^{+}\) and \(u_{1}^{-}\) since they are represented by the same asymptotic expression (34). However, the asymptotic solitons \(u^{\pm}_{2,3}\) still experience the \(\tilde{\eta}\)-dependent phase shifts, namely, \begin{equation} \delta\phi_{2} = - \delta\phi_{3} = - 4\ln(16\sqrt{2}|\mu||\kappa||\tilde{\eta}|), \end{equation} (39) which also change with \(|\tilde{\eta}|\) in a logarithmical manner (see Fig. 9).
(iv)
Likewise, the relative distances and interaction forces between two of asymptotic solitons \(u_{i}^{\pm}\) (\(1\leq i\leq 3\)) can be obtained by \begin{align} d_{\text{I,II}} &= \frac{1}{2}d_{\text{III}} = \frac{\ln(16\sqrt{2}|\mu| |\kappa| |\tilde{\eta}|)}{|\mu|}, \end{align} (40) \begin{align} a_{\text{I,II}} &= \frac{1}{2}a_{\text{III}} = - 512|\mu||\kappa|^{2}e^{-2|\mu|d_{\text{I,II}}}, \end{align} (41) where \(|\tilde{\eta}|\geq\frac{1}{16\sqrt{2}|\mu| |\kappa|}\), the subscripts “I, II, III” respectively represent the distances and interaction forces between the solitons \(u^{\pm}_{1}\) and \(u^{\pm}_{2}\), between \(u^{\pm}_{1}\) and \(u^{\pm}_{3}\), and between \(u^{\pm}_{2}\) and \(u^{\pm}_{3}\). Because of \(a_{\text{I,II,III}}<0\), all the pairwise-soliton interaction forces are of the attractive type. As shown in Fig. 10, the solitons' relative distances grow logarithmically with \(|\tilde{\eta}|\), while the force strengths decrease exponentially with the distances and finally tend to 0 as \(|\tilde{\eta}|\to\infty\).
Figure 9. (Color online) The phase shifts between the envelopes of \(u_{i}^{+}\) and \(u_{i}^{-}\) (\(1\leq i\leq 3\)) versus \(|\tilde{\eta}|\), where the parameters are the same as those in Fig. 6.
Figure 10. (Color online) The solitons' relative distances \(d_{\text{I,II}}\) (black solid), \(d_{\text{III}}\) (red solid) and two-soliton interaction forces \(a_{\text{I,II}}\) (black dashed), \(a_{\text{III}}\) (red dashed) versus \(|\tilde{\eta}|\), where the parameters are the same as those in Fig. 6.
5. Conclusions and Discussions
In this paper, we have constructed the exact N-pole solutions for the Hirota equation, and have studied their asymptotic behavior and soliton interactions. Finally, we address the conclusions and discussions of this paper from the following three aspects:
First, based on the DT and N-soliton solutions, we have obtained the explicit formulas of MP solutions with arbitrary order. Note that the MP solutions are just the degenerate cases of N-soliton solutions when all the spectral parameters coincide with one another. In our derivation, we have employed some limit technique to deal with the coalescence of multiple spectral parameters. In analogy to the double-Wronskian representation of the N-soliton solutions, the MP solutions (16) can also be expressed as the ratio of two determinants in which all the elements are recursively determined in Eq. (18).
Second, we have studied the asymptotic behavior of the double- and triple-pole solutions via an improved asymptotic analysis method. The key point of such method is to find the dominant balances between the exponential and algebraic terms, which in turn determine the center trajectories of asymptotic solitons. Following the same procedure, we can in principle perform an asymptotic analysis of the MP solutions with arbitrary order. It should be noted that one must replace \(\theta_{I}\) in the algebraic terms of the solutions by \(\theta_{R}\) via (20), so that the correct dominant relations can be found between \(\tilde{\eta}\) and \(e^{\theta_{R}}\).
Third, we have discussed the soliton interactions in the double- and triple-pole solutions based on their asymptotic expressions. It is known that there is no interaction force among asymptotic solitons in the usual N-soliton solutions since all of them have the constant velocities. Quite differently, the MP solutions exhibit the attractive force between two solitons, and its strength decays exponentially to 0 with their relative distance. Such an attraction causes that the soliton velocities and phase differences vary slowly with \(\tilde{\eta}\). However, the MP solutions still describe the elastic interactions among solitons with equal amplitudes because the physical quantities of all interacting solitons are retained upon interactions.
Acknowledgments
This work was partially supported by the National Natural Science Foundation of China (Grant Nos. 11705284 and 61505054), and by the program of China Scholarship Council (Grant No. 201806445009). T.X. appreciates the hospitality of the Department of Mathematics & Statistics at McMaster University during his visit in 2019.
References
1 V. E. Zakharov and A. B. Shabat, Sov. Phys. JETP 34, 62 (1972). Google Scholar
2 M. J. Ablowitz and P. A. Clarkson, Solitons, Nonlinear Evolution Equations and Inverse Scattering (Cambridge University Press, Cambridge, U.K., 1992). Google Scholar
3 E. Olmedilla, Physica D 25, 330 (1987). 10.1016/0167-2789(87)90107-2 Crossref, Google Scholar
4 C. Schiebold, Nonlinearity 30, 2930 (2017). 10.1088/1361-6544/aa6d9a Crossref, Google Scholar
5 M. Pichler and G. Biondini, IMA J. Appl. Math. 82, 131 (2017). 10.1093/imamat/hxw009 Crossref, Google Scholar
6 M. Wadati and K. Ohkuma, J. Phys. Soc. Jpn. 51, 2029 (1982). 10.1143/JPSJ.51.2029 Link, Google Scholar
7 H. Tsuru and M. Wadati, J. Phys. Soc. Jpn. 53, 2908 (1984). 10.1143/JPSJ.53.2908 Link, Google Scholar
8 M. Li, J. H. Xiao, and W. J. Liu, Phys. Lett. A 375, 549 (2011). 10.1016/j.physleta.2010.12.031 Crossref, Google Scholar
9 W. R. Sun, B. Tian, H. Zhong, and R. X. Liu, Wave Random Complex 26, 168 (2016). 10.1080/17455030.2015.1125039 Crossref, Google Scholar
10 K. W. Chow, R. H. J. Grimshaw, and E. Ding, Wave Motion 43, 158 (2005). 10.1016/j.wavemoti.2005.09.005 Crossref, Google Scholar
11 D. J. Zhang, S. L. Zhao, and Y. Y. Sun, Rev. Math. Phys. 26, 1430006 (2014). 10.1142/S0129055X14300064 Crossref, Google Scholar
12 B. Yang and Y. Chen, Nonlinear Anal. 45, 918 (2019). 10.1016/j.nonrwa.2018.08.004 Crossref, Google Scholar
13 D. W. C. Lai, K. W. Chow, and K. Nakkeeran, J. Mod. Opt. 51, 455 (2004). 10.1080/09500340408235537 Crossref, Google Scholar
14 H. Takahashi and K. Konno, J. Phys. Soc. Jpn. 58, 3085 (1989). 10.1143/JPSJ.58.3085 Link, Google Scholar
15 K. W. Chow, C. C. Mark, and C. Rogers, J. Comput. Appl. Math. 190, 114 (2006). 10.1016/j.cam.2004.12.042 Crossref, Google Scholar
16 Y. Y. Wang, L. Chen, C. Q. Dai, J. Zheng, and Y. Fan, Nonlinear Dyn. 90, 1269 (2017). 10.1007/s11071-017-3725-5 Crossref, Google Scholar
17 Y.-F. Zheng, G.-Q. Huang, and J. Lin, Acta Phys. Sin. 67, 214207 (2018). 10.7498/aps.67.20180786 Crossref, Google Scholar
18 T. Xu, Y. Chen, M. Li, and D. X. Meng, Chaos 29, 123124 (2019). 10.1063/1.5121776 Crossref, Google Scholar
19 L. D. Faddeev and L. A. Takhtajan, Hamiltonian Methods in the Theory of Solitons (Springer, Berlin, 1987). Crossref, Google Scholar
20 T. V. Nguyen, C. R. Math. 357, 13 (2019). 10.1016/j.crma.2018.11.012 Crossref, Google Scholar
24 J. K. Yang, Phys. Rev. Lett. 91, 143903 (2003). 10.1103/PhysRevLett.91.143903 Crossref, Google Scholar
25 R. F. Rodriguez, J. A. Reyes, A. Espinosa-Ceron, J. Fujioka, and B. A. Malomed, Phys. Rev. E 68, 036606 (2003). 10.1103/PhysRevE.68.036606 Crossref, Google Scholar
26 G. P. Agrawal, Nonlinear Fiber Optics (Academic Press, San Diego, CA, 2012) 5th ed. Google Scholar
27 Y. Fukumoto and T. Miyazaki, J. Fluid Mech. 222, 369 (1991). 10.1017/S0022112091001143 Crossref, Google Scholar
28 F. Demontis, G. Ortenzi, and C. van der Mee, Physica D 313, 61 (2015). 10.1016/j.physd.2015.09.009 Crossref, Google Scholar
29 A. Mahalingam and K. Porsezian, Phys. Rev. E 64, 046608 (2001). 10.1103/PhysRevE.64.046608 Crossref, Google Scholar
30 J. Kim, Q. H. Park, and H. J. Shin, Phys. Rev. E 58, 6746 (1998). 10.1103/PhysRevE.58.6746 Crossref, Google Scholar
31 H. Q. Zhang and S. S. Yuan, Nonlinear Dyn. 89, 531 (2017). 10.1007/s11071-017-3469-2 Crossref, Google Scholar
32 T. Xu, B. Tian, and F.-H. Qi, Z. Naturforsch. A 67, 39 (2012). 10.5560/zna.2011-0055 Crossref, Google Scholar
33 B. L. Guo and S. B. Tan, Sci. China, Ser. A 35, 1425 (1992). Google Scholar
34 L. Huang, J. Xu, and E. G. Fan, Nonlinear Anal. 26, 229 (2015). 10.1016/j.nonrwa.2015.05.011 Crossref, Google Scholar
35 R. F. Zhang and B. L. Guo, Appl. Math. 23, 57 (2008). 10.1007/s11766-008-0108-1 Crossref, Google Scholar
36 Y. S. Tao and J. S. He, Phys. Rev. E 85, 026601 (2012). 10.1103/PhysRevE.85.026601 Crossref, Google Scholar
37 A. Ankiewicz, J. M. Soto-Crespo, and N. Akhmediev, Phys. Rev. E 81, 046602 (2010). 10.1103/PhysRevE.81.046602 Crossref, Google Scholar
38 C. Liu, Y. Ren, Z. Y. Yang, and W. L. Yang, Chaos 27, 083120 (2017). 10.1063/1.4999916 Crossref, Google Scholar
39 S. Y. Chen and Z. Y. Yan, Appl. Math. Lett. 95, 65 (2019). 10.1016/j.aml.2019.03.020 Crossref, Google Scholar
40 B. L. Guo, L. M. Ling, and Q. P. Liu, Phys. Rev. E 85, 026607 (2012). 10.1103/PhysRevE.85.026607 Crossref, Google Scholar
41 L. C. Zhao, B. L. Guo, and L. M. Ling, J. Math. Phys. 57, 043508 (2016). 10.1063/1.4947113 Crossref, Google Scholar
42 M. Li, T. Xu, and D. X. Meng, J. Phys. Soc. Jpn. 85, 124001 (2016). 10.7566/JPSJ.85.124001 Link, Google Scholar
43 T. Xu and D. E. Pelinovsky, Phys. Lett. A 383, 125948 (2019). 10.1016/j.physleta.2019.125948 Crossref, Google Scholar
44 T. Xu, S. Lan, M. Li, L. L. Li, and G. W. Zhang, Physica D 390, 47 (2019). 10.1016/j.physd.2018.11.001 Crossref, Google Scholar
45 M. Li, X. L. Yue, and T. Xu, Phys. Scr. 95, 055222 (2020). 10.1088/1402-4896/ab4503 Crossref, Google Scholar
46 J. Cen and A. Fring, Physica D 397, 17 (2019). 10.1016/j.physd.2019.05.005 Crossref, Google Scholar