1 Introduction↩︎

In this work, we consider the stability of \(N\)-soliton solutions for the modified Camassa–Holm (mCH) equation [1]–[3] \[\begin{align} \label{mch} {m}_{t} + {\left( ( {u}^{2} - {u}_{x}^{2}) m\right) }_{x} + \gamma {u}_{x} = 0,\quad t > 0,\quad x \in \mathbb{R}, \end{align}\tag{1}\] where \(u(t,x)\) is a real-valued function of the time-space variables \((t,x)\), the momentum density is defined as \(m(t,x)=(1-\partial_x^2)u(t,x)\), and the constant \(\gamma>0\) acts as the linear dispersion parameter.

The mCH equation 1 was originally identified as a completely integrable system by Fuchssteiner via the method of recursion operators [3], [4]. It was independently recovered by Olver and Rosenau through the application of tri-Hamiltonian duality within the bi-Hamiltonian framework of the modified Korteweg-de Vries equation [5]. From a physical standpoint, the mCH equation 1 serves as a model for the unidirectional propagation of moderate-amplitude shallow-water waves over a flat bottom [6]. Here, \(u(t,x)\) denotes the horizontal fluid velocity at a specific depth, and the linear dispersion parameter \(\gamma > 0\) corresponds to the critical shallow-water wave speed.

Moreover, the mCH equation 1 is deeply intertwined with other prominent nonlinear dispersive models. Notably, under a characteristic coordinate transformation, it reduces to the integrable short-pulse equation [7], expressed as \[\begin{align} \label{spe} u_{xt} = u + \frac{1}{6}(u^3)_{xx}. \end{align}\tag{2}\] Given that equation 2 is a paradigmatic model for the propagation of ultra-short optical pulses in silica fibers, this connection further underscores the broad physical relevance and intricate mathematical structure of the mCH equation.

A hallmark of the mCH equation 1 is its complete integrability, which is intrinsically tied to its bi-Hamiltonian structure [8], [9]. Specifically, the equation admits two compatible Hamiltonian operators, enabling the formulation \[\begin{align} \label{shmdjg} {m}_{t} = \mathcal{K}_{0}\frac{\delta E}{\delta m} = \mathcal{K}_{1}\frac{\delta H}{\delta m}, \end{align}\tag{3}\] where \(m = u - {u}_{xx}\). The two compatible Hamiltonian operators \(\mathcal{K}_0\) and \(\mathcal{K}_1\) are given by \[\begin{align} \mathcal{K}_{0}=-{\partial}_{x}m{\partial}_{x}^{-1}m{\partial}_{x}-\frac{1}{2}\gamma{\partial}_{x}, \quad \mathcal{K}_{1} = {\partial }_{x} - {\partial }_{x}^{3}, \end{align}\] which correspond to the Hamiltonian functionals \[\begin{align} \tag{4} &E(u)=\int_{\mathbb{R}}\left({{u}^{2} + {u}_{x}^{2}}\right) dx,\\ &H(u)= \frac{1}{4}\int_{\mathbb{R}}\left( {-{u}^{4} - 2{u}^{2}{u}_{x}^{2} + \frac{1}{3}{u}_{x}^{4} - {2\gamma }{u}^{2}}\right) dx.\tag{5} \end{align}\]

Identifying the natural domains for these conserved quantities is essential. The functional \(E(u)\) corresponds to the squared norm of the standard Sobolev space \(H^1(\mathbb{R})\), ensuring it is well-defined for any velocity profile \(u \in H^1(\mathbb{R})\). In contrast, the Hamiltonian \(H(u)\) contains a higher-degree nonlinear term, namely \(\frac{1}{3}u_x^4\). As a result, \(H(u)\) requires higher regularity and is rigorously defined only within \(W^{1,4}(\mathbb{R})\). Consequently, the natural admissible space for the full bi-Hamiltonian system is restricted to the intersection \(H^1(\mathbb{R}) \cap W^{1,4}(\mathbb{R})\) [10].

Moreover, by employing the operator \(\mathcal{K}_1\) and the identity \(m = (1-\partial_x^2)u\), the bi-Hamiltonian system 3 admits an equivalent, more direct representation in terms of the primary velocity variable \(u\), given by \[\begin{align} \label{hmdjgu} \frac{\partial u}{\partial t} = \mathcal{J}\frac{\delta H}{\delta u}, \end{align}\tag{6}\] where the operator \(\mathcal{J}={\partial}_{x}{\left(1-{\partial}_{x}^{2}\right) }^{-1}\) is skew-symmetric and is a bounded operator on \({L}^{2}(\mathbb{R})\). This specific formulation is instrumental for analyzing the underlying dynamics and variational properties of the solutions [9], [11]–[13].

Beyond its elegant algebraic framework, the mCH equation 1 exhibits a rich spectrum of nonlinear wave phenomena, which notably include the existence of solitary waves and multi-soliton solutions [10], [14], [15]. Rigorous analysis establishes that the mCH equation admits smooth solitary wave solutions if and only if the wave propagation speed \(c\) satisfies the condition \(0 < \gamma < c \le 2\gamma\). These traveling wave states assume the standard ansatz \(u(t,x) = \phi_c(x - ct + x_0)\), where \(\phi_c\) denotes the localized wave profile and \(x_0 \in \mathbb{R}\) is an arbitrary phase shift. Substituting this traveling wave ansatz into the mCH equation 1 yields a nonlinear ordinary differential equation governing the profile function: \[\begin{align} \label{profile95eq1} \left(\phi_c^2 - (\partial_\xi\phi_c)^2 - c\right)\left(\phi_c - \partial_\xi^2\phi_c\right) = -\gamma\phi_c, \end{align}\tag{7}\] where \(\xi = x - ct + x_0\). Multiplying 7 by \(\partial_\xi\phi_c\) and integrating once with respect to \(\xi\) yields the corresponding first-order integrated form \[\begin{align} \label{profile95eq2} \left(\phi_c^2 - (\partial_\xi\phi_c)^2 - c\right)^2 = c^2 - 2\gamma\phi_c^2, \end{align}\tag{8}\] which characterizes the homoclinic orbit associated with the smooth solitary wave dynamics.

In addition to single solitary waves, the complete integrability of the mCH equation ensures the existence of complex multisoliton solutions \(U^{(N)}(t,x)\), which describe the nonlinear interaction of \(N\) distinct solitons [16]. A key dynamical feature of these configurations is their purely elastic interaction [17]. As time approaches infinity, the composite solution asymptotically decouples into a linear superposition of \(N\) independent and spatially isolated solitary waves traveling at distinct speeds \(c_i\) for \(i = 1, 2, \dots, N\), namely \[\begin{align} \label{asymptotic95N} U^{(N)}(t,x) \sim \sum_{i=1}^{N} \phi_{c_i}(x - c_i t - x_i), \quad \text{as } t \to \infty, \end{align}\tag{9}\] where the phase parameters \(x_i \in \mathbb{R}\) depend on the corresponding wave speeds \(c_i\). This asymptotic decoupling demonstrates the absence of dispersive effects at spatial infinity, establishing a rigorous analytical foundation for investigating the invariant quantities and stability properties of multisoliton dynamics [18], [19].

The Cauchy problem for the mCH equation has been extensively investigated concerning its well-posedness and wave-breaking phenomena. Initial studies established local well-posedness for initial data in the Sobolev space \(H^s(\mathbb{R})\) with \(s > \frac{7}{2}\) [20], a regularity threshold that was subsequently relaxed to \(s > \frac{5}{2}\) [9], [21]. Analogous to the classical CH model [22], [23], the mCH equation exhibits finite-time wave-breaking, characterized by a uniformly bounded solution with an unbounded spatial derivative [9]. Conversely, global existence is guaranteed for sufficiently regular initial data (\(s \ge 3\)) that satisfy specific momentum threshold conditions, effectively precluding wave-breaking [14].

The stability analysis of solitary waves and multisolitons remains an active area of research in nonlinear dispersive equations. Comparing the mCH equation 1 with the classical Korteweg-de Vries (KdV) [24], [25] and Camassa–Holm (CH) [11], [26] models provides an essential methodological benchmark. These fundamental shallow water equations capture different nonlinear dispersive regimes. Extensive investigations into KdV multisoliton stability culminated in the seminal constrained variational principle developed by Maddocks and Sachs [27] in the natural energy space, a result recently extended to the negative Sobolev space \(H^{-1}(\mathbb{R})\) by Killip and Vişan [28]. Concurrently, CH theory established the orbital stability of single solitons and peakons in \(H^1(\mathbb{R})\) through spectral methods [29]–[31], while multisoliton configurations were subsequently proven stable in \(H^n(\mathbb{R})\) by adapting localized spectral techniques [12], [32], [33]. Furthermore, Li and Zhang [34] utilized nonlocal conservation laws to prove that the \(N\)-solitons of the CH equation are orbitally stable. Additionally, Wu and Tian [35] demonstrated that the \(N\)-solitons of the two-component CH equation are dynamically stable and that its 2-soliton solutions are orbitally stable. These fundamental advancements highlight the pivotal role of higher-order conserved quantities in rigorously establishing the stability of multisolitons.

For the mCH equation, initial investigations established the orbital stability of the smooth single soliton through classical constrained minimization [10]. Utilizing the foundational conserved quantities \(E(m)\) and \(H(m)\), stability was proven within the intersection space \(W^{1,4}(\mathbb{R}) \cap H^1(\mathbb{R})\) for the momentum density \(m = u - u_{xx}\). To circumvent analytical obstacles introduced by the higher-degree nonlinearities in \(H(u)\), subsequent studies shifted to the momentum density formulation \(m\). Leveraging the locally conserved functionals \[\begin{align} &F_1(m) = \int_{\mathbb{R}}\left(\frac{-2}{\gamma}M+\sqrt{2}\gamma^{-\frac{1}{2}}\right)dx,\\ &F_2(m) = \int_{\mathbb{R}}\left(\frac{m_x^2}{2M^5}-\frac{1}{\gamma M}+\sqrt{2}\gamma^{-\frac{3}{2}}\right)dx, \end{align}\] where \(M=\sqrt{m^2+\gamma/2}\), researchers successfully demonstrated the orbital stability of the single smooth solitary wave in \(H^3(\mathbb{R})\) [14].

Recently, this momentum-driven framework was extended to investigate double soliton solutions [15]. Utilizing the multiplier method, a higher-order conserved integral \[\begin{align} F_3(m) = \int_{\mathbb{R}}\left(\frac{m_{xx}^2}{2M^7}+\frac{5m_x^2}{4M^7}-\frac{7m_x^4}{2M^9}+\frac{35\gamma m_x^4}{16M^{11}}-\frac{1}{4\gamma M^3}-\frac{2}{\gamma^2 M}+\frac{5\sqrt{2}}{2}\gamma^{-\frac{5}{2}}\right)dx, \end{align}\] was rigorously derived in [8], [16], [21]. The subsequent stability analysis for the \(2\)-soliton solutions is extended to the higher-regularity Sobolev space \(H^4(\mathbb{R})\) for the primary velocity variable \(u\).

Despite these significant advancements, the study of the mCH equation requires further investigation to relax existing regularity constraints and to extend the stability framework to multi-soliton configurations. The primary objective of this work is to establish the dynamical and orbital stability of \(N\)-soliton solutions in \(H^{N+1}(\mathbb{R})\). Driven by the unique cubic nonlinearities and bi-Hamiltonian structure of the mCH equation, this study addresses three central issues: utilizing its bi-Hamiltonian hierarchy to generate conservation laws for stability analysis in appropriate functional spaces, establishing a rigorous mathematical framework to determine the spectral inertia of generalized \(N\)-solitons, and applying spectral decomposition and constrained variational principles to verify the dynamical and orbital stability of mCH \(N\)-solitons under complex geometric constraints.

Our strategy is to adapt and generalize the framework of Maddocks and Sachs [27], which originally addressed the stability of multi-solitons for the KdV equation via a constrained variational principle. The proof of stability is thus reduced to analyzing the spectral inertia of the second variation of a specially constructed Lyapunov functional. This approach has been extended to various integrable models [12], [35]–[41]. However, the application of this methodology to the mCH \(N\)-solitons presents several significant technical challenges.

Our analysis addresses three principal technical challenges. The first involves the domain characterization of the conservation hierarchy, as higher-order conserved quantities \(\{F_n\}\) often involve higher-order spatial derivatives such as \(m_{xx}\), necessitating restrictive functional topologies like \(H^4(\mathbb{R})\). Identifying a recursive set of conserved quantities that balances mathematical tractability with minimal regularity requirements is therefore essential. The second challenge lies in the spectral analysis of the linearized operator, which for general mCH \(N\)-solitons escalates to a differential operator of order \(2(N+1)\). The presence of these high-order derivative terms and non-constant coefficients exhibiting complex nonlinear interactions precludes the specific algebraic factorizations used in two-soliton scenarios, thereby demanding a fully systematic framework to determine the spectral inertia. Finally, the multidimensional stability constraints induced by \(N\)-soliton interactions pose a formidable structural hurdle. Due to the inherently nonlocal nature of the mCH recursion operator, the constrained variational framework is highly sensitive to the generalized conservation laws and their corresponding Lagrange multipliers, which rigorously requires the explicit evaluation of a highly coupled \(N \times N\) Hessian matrix governed by generalized symmetric polynomials.

In order to systematically overcome the aforementioned difficulties, our analysis employs the following strategies.

To address the initial challenge concerning the domain of conservation laws and stringent regularity constraints, we leverage the intrinsic bi-Hamiltonian structure, circumventing the limitations of the conventional \(\{F_n\}\) sequence. By combining the compatible Hamiltonian operators \(\mathcal{K}_0\) and \(\mathcal{K}_1\), we define the recursion operator

\[\begin{align} \mathcal{K}[m] = \mathcal{K}_1^{-1}\mathcal{K}_0 = -(1-\partial_x^2)^{-1}\left(m\partial_x^{-1}m\partial_x + \frac{\gamma}{2}\right), \end{align}\]

and the explicit inverse of this operator is given by

\[\begin{align} \label{hdgszm} \widehat{\mathcal{K}}[m] = -\frac{2}{\gamma}\left(1 - \frac{m}{M}\partial_x^{-1}\frac{m}{M}\partial_x\right)\left(1-\partial_x^2\right), \end{align}\tag{10}\]

where \(M = \sqrt{m^2 + \frac{\gamma}{2}}\). Applying \(\widehat{\mathcal{K}}[m]\) recursively from the foundational energy \(E_0(u) = E(u)\) generates a novel hierarchy of conserved quantities \(\{E_n\}_{n \ge 0}\) satisfying the recursive relation

\[\begin{align} \label{dggxm} \frac{\delta E_n}{\delta m} = \widehat{\mathcal{K}}[m]\frac{\delta E_{n-1}}{\delta m}, \quad n = 1, 2, \dots, N. \end{align}\tag{11}\]

This recursive scheme yields explicit higher-order functionals, for instance,

\[\begin{align} E_1(m) &= \int_{\mathbb{R}}\left(-2\sqrt{\frac{2}{\gamma}}M+2\right)dx,\\ E_2(m) &= \int_{\mathbb{R}}\left(\sqrt{\frac{\gamma}{2}}\frac{m_x^2}{M^5}+\frac{2}{\gamma}\sqrt{\frac{2}{\gamma}}\frac{m^2}{M}\right)dx. \end{align}\]

To perform the variational analysis directly on the primary velocity profile \(u\), the recursion relations 11 are reformulated as

\[\begin{align} \frac{\delta E_n}{\delta u} = \widehat{\mathcal{R}}(u)\frac{\delta E_{n-1}}{\delta u}, \quad n = 1, 2, \dots, N. \end{align}\]

Here, the transformed recursion operator \(\widehat{\mathcal{R}}(u)\) is defined as \[\begin{align} \widehat{\mathcal{R}}(u) = -\frac{2}{\gamma}\left(1-\partial_x^2\right)\left(1 - \frac{m}{M}\partial_x^{-1}\frac{m}{M}\partial_x\right), \end{align}\] with its inverse \[\begin{align} \mathcal{R}(u) = (1-\partial_x^2)\mathcal{K}[m](1-\partial_x^2)^{-1} = -\left(m\partial_x^{-1}m\partial_x + \frac{\gamma}{2}\right)(1-\partial_x^2)^{-1}. \end{align}\]

Notably, the natural domains for the generated sequence \(\{E_n\}_{n=0}^N\) correspond precisely to the Sobolev spaces \(H^{n+1}(\mathbb{R})\) (as proved in Section 2). This property directly fulfills our objective to lower the regularity threshold and identifies the optimal energy space for the stability analysis.

Equipped with this functional hierarchy, we apply the constrained variational principle by constructing the augmented Lyapunov functional \(\mathcal{S}_{N}\) for the mCH \(N\)-solitons as \[\begin{align} \mathcal{S}_{N}(u) = (-1)^N\left(E_{N}(u) + \sum_{m=0}^{N-1} \mu_m E_m(u)\right),\label{lagrange95N} \end{align}\tag{12}\]

where the Lagrange multipliers \(\mu_{m}\) (\(m=0, 1, \dots, N-1\)) correspond to elementary symmetric functions of the wave speeds \(c_{1}, \dots, c_{N}\) (as detailed in Section 2). We demonstrate that \(U^{(N)}\) is a critical point of \(\mathcal{S}_{N}\), satisfying the generalized Euler-Lagrange equation

\[\begin{align} \frac{\delta E_{N}(u)}{\delta u} + \sum_{m=0}^{N-1} \mu_{m} \frac{\delta E_{m}(u)}{\delta u} = 0,\quad \text{at } u = U^{(N)}.\label{1462095N} \end{align}\tag{13}\]

Establishing the dynamical stability of \(U^{(N)}\) requires showing that it minimizes \(E_{N}\) subject to the constraints \(E_{m}(u) = E_{m}\left(U^{(N)}\right)\) for \(m = 0, 1, \dots, N-1\). This fundamentally requires the self-adjoint second variation operator,

\[\begin{align} \mathcal{L}_N := \mathcal{S}''_{N}(U^{(N)}), \label{VL95N} \end{align}\tag{14}\]

to be positive definite under modulations along the constrained directions.

To address the second challenge regarding the spectral analysis of the high-order differential operator \(\mathcal{L}_N\), we overcome the inapplicability of classical Sturm-Liouville theory by bridging single-soliton spectral data with the composite \(N\)-soliton structure. We focus our study on the linear operator \(L_n = E''_{n}(\phi_c) + \omega E''_{n-1}(\phi_c)\) for the single soliton \(\phi_c\) and determine its spectral properties by leveraging the Inverse Scattering Transform (IST) framework [13], [19], [42]–[45]. By constructing squared eigenfunctions \(F^\pm(x,k)\) from the Jost solutions of the Lax pair, we formulate a generalized basis that facilitates a rigorous analysis of \(L_n\) and its conjugate forms \(\mathcal{J}L_n\) and \(L_n\mathcal{J}\). Specifically, we exploit the fundamental operator identities establishing that \(\mathcal{J}L_n\) and \(L_n\mathcal{J}\) commute with the adjoint recursion operator \(\widehat{\mathcal{R}}^*(\phi_c)\) and the recursion operator \(\widehat{\mathcal{R}}(\phi_c)\), respectively. As \(t \to \infty\), the spectrum of the \(N\)-soliton operator \(\mathcal{L}_N\) converges to the union of the individual single-soliton spectra, allowing us to precisely calculate the inertia index of \(\mathcal{L}_N\) to verify the critical spectral condition \(\mathrm{in}(\mathcal{L}_N) = p(D_N)\). In this formulation, \(\mathrm{in}(\mathcal{L}_N)\) denotes the number of negative eigenvalues of \(\mathcal{L}_N\) and \(p(D_N)\) represents the number of positive eigenvalues of the \(N \times N\) Hessian matrix \(D_N := \left\{ \frac{\partial^{2} \mathcal{S}_{N}}{\partial \mu_{i} \partial \mu_{j}} \right\}\).

Finally, to resolve the structural complexity of the linearized \(N\)-soliton operator, we establish dynamical stability via a refined coercivity estimate. Spectral analysis reveals that \(\mathcal{L}_N\) possesses exactly \(\left[\frac{N+1}{2}\right]\) negative eigenvalues and an \(N\)-dimensional generalized kernel, corresponding to the potential instability directions associated with spatial translations and scalings. Given that \(\mathcal{L}_N\) is a strongly elliptic operator of order \(2(N+1)\), we apply a generalized Gårding’s inequality to obtain a local coercivity bound near the \(N\)-soliton manifold in \(H^{N+1}(\mathbb{R})\). By introducing \(N\) dynamic temporal modulation parameters to continuously project out the neutral translational directions, we extend this local coercivity to a global bound on the \(H^{N+1}\) distance. A standard continuity argument then completes the proof of orbital stability.

Next, we will present our main results as follows. An essential prerequisite for establishing the stability of solitary waves is the existence of global solutions, as demonstrated in [9].

Building upon this well-posedness result, we state our main theorems concerning the stability of the \(N\)-solitons for the mCH equation 1 . A central objective of this work is not only to establish the dynamical and orbital stability of the mCH \(N\)-soliton solutions, but also to further reduce the required regularity of the stability space in [15]. This reduction relies on the fact that the novel hierarchy of conservation laws \(\{E_k\}_{k=0}^{N}\), derived from the bi-Hamiltonian system of the mCH equation, possesses well-defined and optimal regularity domains. Therefore, we first establish the following foundational theorem.

Equipped with the aforementioned hierarchy of conservation laws, we are able to construct a suitable augmented Lyapunov functional framed within an optimal energy space. Utilizing the generalized constrained variational framework developed by Maddocks and Sachs [27], we establish the following theorem regarding the dynamical stability of the mCH \(N\)-solitons.

While Theorem 2 establishes the stability of the entire \(N\)-soliton manifold \(G_{\mathbf{c}}\) under the mCH flow, it is fundamentally important to characterize the stability of the multi-soliton profile up to dynamic spatial shifts. By introducing appropriate temporal modulation parameters to continuously project the perturbed solution onto the translation directions, we can translate the dynamical stability of the manifold into the orbital stability of the wave profile. This leads to our second main result.

Following the statement of the main theorems, we outline the organization of the subsequent proofs. Our general strategy hinges on a detailed spectral analysis of the recursion operators associated with the \(N\)-soliton solutions \(U^{(N)}\). By verifying that the constructed augmented Lyapunov functional \(\mathcal{S}_N\) satisfies the Maddocks-Sachs stability criterion, we employ a modulation argument to conclude both dynamical and orbital stability.

The proof transforms the nonlinear stability problem into a spectral analysis of a specific linearized operator, exploiting the mCH equation’s bi-Hamiltonian recursion operators to construct a suitable Lyapunov functional. The paper is organized as follows:

In Section 2, we review the IST framework and the Hamiltonian formulation. We construct the recursion operator \(\widehat{\mathcal{R}}\) to generate the conservation laws \(\{E_n\}_{n=0}^N\), variationally characterize the \(N\)-solitons as constrained critical points of \(\mathcal{S}_N\), and compute the positive inertia index \(p(D_N)\) of the associated Hessian matrix.

Section 3 focuses on the spectral analysis of the high-order linearized operator \(\mathcal{L}_N = \mathcal{S}_N''(U^{(N)})\). Utilizing IST squared eigenfunctions and the asymptotic decoupling of \(N\)-solitons as \(t \to \infty\), we determine the spectrum of \(\mathcal{L}_N\). This rigorously verifies the critical spectral condition \(\mathrm{in}(\mathcal{L}_N) = p(D_N) = \left[ \frac{N+1}{2}\right]\).

In Section 4, we finalize the stability proofs via constrained variational principles. We establish local coercivity in \(H^{N+1}(\mathbb{R})\) by applying a generalized Gårding’s inequality to \(\mathcal{L}_N\), and employ a dynamic modulation argument to handle translation invariance, which extends the local bound to a global \(H^{N+1}\) estimate and concludes both dynamical and orbital stability.

2 Preliminaries↩︎

In this section, we present fundamental results from the inverse scattering transform theory for the mCH equation (1 ). The section is organized as follows.

2.1 Spectral problem↩︎

As a significant approach for obtaining the orbital stability of the mCH equation, we first introduce the spectral problem of the mCH equation and its inverse scattering transform.

The completely integrable mCH equation 1 possesses the following Lax pair \[\begin{align} \label{lax} \Phi_x(x,t,k)=U(x,t,k)\Phi(x,t,k),\\ \Phi_t(x,t,k)=V(x,t,k)\Phi(x,t,k), \end{align}\tag{15}\] where \(k\) is the spectral parameters, the matrix function \(U\) and \(V\) satisfy the zero- curvature equation \(U_t-V_x+UV-VU=0\) and their can be expressed respectively as \[\begin{align} &U(x,t,\lambda)=\begin{pmatrix} -i\mu & \frac{1}{2}\lambda m(x,t) \\ -\frac{1}{2}\lambda m(x,t) & i\mu \end{pmatrix},\\ &V(x,t,\lambda)=\begin{pmatrix} \frac{2i\mu}{\lambda^2}+i\mu\left(u^2-u_x^2\right) & -\frac{u-2i\mu u_x}{\lambda}-\frac{1}{2}\lambda\left(u^2-u_x^2\right)m \\ \frac{u+2i\mu u_x}{\lambda}+\frac{1}{2}\lambda\left(u^2-u_x^2\right)m & -\frac{2i\mu}{\lambda^2}-i\mu\left(u^2-u_x^2\right) \end{pmatrix}, \end{align}\] and \[\begin{align} \mu(k)=\frac{1}{4}\left(k-\frac{1}{k}\right),\quad \lambda(k)=\frac{1}{\sqrt{2\gamma}}\left(k+\frac{1}{k}\right). \end{align}\] Furthermore, we can easily find that \[\begin{align} \mu=\mu(\lambda,\gamma)=-\frac{i}{2}\sqrt{1-\frac{1}{2}\gamma\lambda^2}. \end{align}\]

We consider solutions \(u(t,x)\) with decaying boundary conditions as \(|x|\to \infty\) and for the continuous spectrum \(\mu(k)\in\mathbb{R}\), the spatial spectral problem 15 admits fundamental matrix solutions \(\Phi(x,k)\) and \(\Psi(x,k)\), satisfying \[\begin{align} \label{L2} \lim_{x\to-\infty}e^{i\mu y\sigma_3}\Phi(x,k)=\mathbb{I},\quad \lim_{x\to+\infty}e^{i\mu y\sigma_3}\Psi(x,k)=\mathbb{I}, \end{align}\tag{16}\] where \(y=x+\int_{+\infty}^{x}(\sqrt{2m^2/\gamma+1}-1)ds\) and \(\sigma_3=\operatorname{diag}\left(1,-1\right)\). These are the so-called Jost solutions. Denote the columns of the Jost solutions by \[\begin{align} \Phi=[\phi^+,\phi^-]=\begin{pmatrix} \phi^+_1 & \phi^-_1 \\ \phi^+_2 & \phi^-_2 \end{pmatrix},\quad \Psi=[\psi^-,\psi^+]=\begin{pmatrix} \psi^-_1 & \psi^+_1 \\ \psi^-_2 & \psi^+_2 \end{pmatrix}. \end{align}\] Using the relation \(\mathrm{Im} \mu=\frac{1+|k|^2}{4|k|^2}\mathrm{Im}k\), standard arguments show that \(\phi^+(x,\lambda)\) and \(\psi^+(x,\lambda)\) are analytic for \(\mathrm{Im} k>0\), while \(\phi^-(x,\lambda)\) and \(\psi^-(x,\lambda)\) are analytic for \(\mathrm{Im}k<0\).

Moreover, the smooth solitary wave solution \(\phi_c\) of the mCH equation 1 was derived by Matsuno [16] via the bilinear method. For our purposes, it can be expressed as \[\begin{align} &u=\frac{16i\mu_1}{\gamma\lambda_1^3}\operatorname{sgn}(\lambda_1)\frac{\cosh(\xi)}{\cosh(2\xi)+\frac{1-4\mu_1^2}{1+4\mu_1^2}},\\ &\xi=2i\mu_1\left(x-\frac{2t}{\lambda_1^2}-\ln\left(\frac{1+\frac{1+2i\mu_1}{1-2i\mu_1}e^{2\xi}}{1+\frac{1-2i\mu_1}{1+2i\mu_1}e^{2\xi}}\right)\right), \end{align}\] where \(\mu_1=\mu(k_1)\), \(\lambda_1=\lambda(k_1)\), and the discrete spectral point \(k_1\) satisfies \(|k_1|=1\) and \(a(k_1)=0\). Consequently, one has \(c=\frac{2}{\lambda_1^2}\).

Furthermore, from our other work [46], the smooth solitary-wave solution \(\phi_c\) admits the complete set \[\begin{align} \label{complete} \Bigg\{ \left(1-\partial_x^2\right)\left((\psi^-_{1})^2 + (\psi^-_{2})^2\right)(x,k),\;\text{for } k\in \mathbb{R}^+; \quad &\left(1-\partial_x^2\right)\left((\psi^-_{1})^2 + (\psi^-_{2})^2\right)(x,k_1^-); \notag \\ & \left(1-\partial_x^2\right)\frac{d}{dk}\left((\psi^-_{1})^2 + (\psi^-_{2})^2\right)(x,k_1^-) \Bigg\} \end{align}\tag{17}\] and the other complete set \[\begin{align} \bigg\lbrace & \left((\psi^-_{1})^2 + (\psi^-_{2})^2\right)_x(x,k),\;\text{for } k\in \mathbb{R}^+; \notag \\ & \left((\psi^-_{1})^2 + (\psi^-_{2})^2\right)_x(x,k_1^-);\; \frac{d}{dk}\left((\psi^-_{1})^2 + (\psi^-_{2})^2\right)_x(x,k_1^-) \bigg\rbrace . \end{align}\] where \(k_1^-\) is the discrete spectral point on the unit circle in the fourth quadrant corresponding to \(\phi_c\). This relation will be helpful for our subsequent research.

2.2 Hamiltonian formation↩︎

As the foundation for analyzing the multi-soliton stability of the mCH equation (1 ), it is crucial to determine the domain of the conservation laws so as to define the space of soliton stability. Hence, we will proceed to establish the proof of Theorem 1

To rigorously establish this theorem, we divide the proof into the following lemmas. First, we determine the exact highest order of spatial derivatives in the variational derivative \(\frac{\delta E_k}{\delta u}\) and verify its non-degeneracy.

Next, we will show that the mCH equation, through the application of the recursion operator, has an expression in terms of the wave speed \(c\), which will be of great help for our subsequent spectral analysis. We denote \(E(u):=E_0(u)\), then the recursion relations 11 can be equivalently expressed in terms of \(u\) as \[\begin{align} \label{hamilton} \frac{\delta E_k}{\delta u}=\widehat{\mathcal{R}}(u)\frac{\delta E_{k-1}}{\delta u}=\cdots=\widehat{\mathcal{R}}^k(u)\frac{\delta E_0}{\delta u}, \end{align}\tag{21}\] where \[\begin{align} \label{inverse} \widehat{\mathcal{R}}(u)=\left(1-\partial_x^2\right)\widehat{\mathcal{K}}[m]\left(1-\partial_x^2\right)^{-1}=-\frac{2}{\gamma}\left(1-\partial_x^2\right)\left(1-\frac{m}{M}\partial_x^{-1}\frac{m}{M}\partial_x\right), \end{align}\tag{22}\] and the inverse of \(\widehat{\mathcal{R}}(u)\) can be written as \[\begin{align} \mathcal{R}(u)=\left(1-\partial_x^2\right)\mathcal{K}[m]\left(1-\partial_x^2\right)^{-1}=-\left(m\partial_x^{-1}m\partial_x+\frac{\gamma}{2}\right)\left(1-\partial_x^2\right)^{-1}. \end{align}\] Near the smooth solitary wave solution \(\phi_c\), the conserved quantities \(E_k, k\in 1, 2, \cdots,N\) will satisfy the following first-variation identities \[\begin{align} \label{24619} \left. \frac{\delta E_k}{\delta u} \right|_{u=\phi_c}=-\omega\left. \frac{\delta E_{k-1}}{\delta u} \right|_{u=\phi_c}=\cdots=(-\omega)^k\left. \frac{\delta E_0}{\delta u} \right|_{u=\phi_c}, \end{align}\tag{23}\] where \(\omega = \frac{2}{c}\). Now multiply (23 ) with \(\frac{\mathrm{d}\phi_c}{\mathrm{d}c}\), one has \[\begin{align} \frac{\mathrm{d}E_k(\phi_c)}{\mathrm{d}c}=-\omega \frac{\mathrm{d}E_{k-1}(\phi_c)}{\mathrm{d}c}=\cdots=(-\omega)^k\frac{\mathrm{d}E_0(\phi_c)}{\mathrm{d}c}, \end{align}\] and therefore \[\begin{align} E_k(\phi_c)=(-1)^k\int_{0}^{c}y^k\frac{\mathrm{d}E_0(\phi_y)}{\mathrm{d}y}\mathrm{d}y, \quad k=1,2, \cdots, N. \end{align}\] Furthermore, we could easily find that \[\begin{align} E_0(\phi_c)=\int_\mathbb{R}\left(\phi_c^2+\phi_{cx}^2\right)dx=\int_\mathbb{R}\phi_c\left(\phi_c-\phi_{cxx}\right)dx, \end{align}\] and one can obtain that \[\begin{align} \frac{dE_0(\phi_c)}{dc}=2\sqrt{\frac{c}{c-\gamma}}>0. \end{align}\] Thus, one has \[\begin{align} E_k(\phi_{c})=(-1)^k 2\int_{0}^{c}y^k\sqrt{\frac{y}{y-\gamma}}\mathrm{d}y, \quad k=1,2, \cdots, N. \end{align}\]

2.3 Variational characterization of the \(N\)-solitons↩︎

In this subsection, we would demonstrate that the mCH \(N\)-solitons satisfy (13 ) if one prescribes the Lagrange multipliers \(\mu_k, k=0,1,\cdots, N-1\) appropriately. This provides a variational characterization of \(U^{(N)}\). The main result in this subsection is as follows.

We consider now the following mCH \(N\)-solitons \(U^{(N)}(t,x;\mathbf{c},\mathbf{x})\) variational principle. The variational expression corresponds to the functional derivative of (13 ) evaluated at \(u = U^{(N)}\). Generally, the \(N\)-soliton solution \(U^{(N)}(t,x)\) does not minimize \(\mathcal{S}_N\) globally but constitutes a constrained, non-isolated critical point for the optimization problem \[\min E_{N}(u(t)) \quad \text{subject to} \quad E_j(u(t)) = E_j(U^{(N)}(t)), \quad j =0,1, \cdots, N-1.\]

Consider the self-adjoint operator \(\mathcal{L}_N(t)\) defined by (14 ), whose second variation properties are essential. Denote by \(n(\mathcal{L}_N(t))\) the count of its negative eigenvalues. These quantities exhibit potential time-dependence. The \(N \times N\) Hessian matrix is defined as \[\begin{align} \label{4461} D_N(t) := \left\{ \frac{\partial^2 \mathcal{S}_N(U^{(N)}(t))}{\partial \mu_{i-1} \partial \mu_{j-1}} \right\}, \end{align}\tag{24}\] with \(p(D(t))\) representing its positive eigenvalue count. Base upon the facts above and approach in [27], the proof of Theorem2 reduces to calculate the exact value of \(n(\mathcal{L}_N)\) and \(p(D_N)\). Concerning \(p(D_N)\), one has the following results.

2.4 Recursion operators around the smooth solitons↩︎

Let us recall that the soliton \(\phi_c(x-ct)\) is a solution of the mCH equation. For simplicity, we denote \(\phi_{c}\) by \(\phi\). Then by (21 ) we have \[\begin{align} E'_k(\phi)=\widehat{R}(\phi)E'_{k-1}(\phi), \end{align}\] where \(\widehat{\mathcal{R}}[\phi]\) is the associated operator \[\begin{align} \widehat{\mathcal{R}}[\phi]=-\frac{2}{\gamma}\left(1-\partial_x^2\right)\left(1-\frac{m_{\phi}}{M_{\phi}}\partial_x^{-1}\frac{m_{\phi}}{M_{\phi}}\partial_x\right). \end{align}\]

To analyze the second variation of the actions, we linearize the equation (21 ) to let \(u=\phi+\epsilon z\), and obtain a relation between linearized Hamiltonian \(E''_{k}(\phi)\) and \(E''_{k-1}(\phi)\) for all \(k\geq 1\). One has the following results.

3 Spectral analysis↩︎

Let \(U^{(N)} = U^{(N)}_{\mathbf{c,x}}\) be the \(N\)-soliton solution to the mCH equation 1 , characterized by the wave speed vector \(\mathbf{c} = (c_1, \dots, c_N)\) and the initial shift vector \(\mathbf{x} = (x_1, \dots, x_N)\). This section is devoted to a comprehensive spectral analysis of the linearized operator \(\mathcal{L}_N = \mathcal{S}_N''(U^{(N)})\), where the augmented Lyapunov functional \(\mathcal{S}_N\) is constructed through the recursion operator \(\widehat{\mathcal{R}}(\phi)\) and the bi-Hamiltonian hierarchy discussed in Section 2.

To determine the spectral properties of \(\mathcal{L}_N\), we adopt the iso-inertial framework proposed by Maddocks and Sachs [27].

As \(U^{(N)}\) constitutes a critical point of the functional \(\mathcal{S}_N\) and fits within the general framework of Hamiltonian systems with symmetries [27], [47], [48], the inertia \(in(\mathcal{L}_N)\) is an invariant of the mCH flow. Consequently, we can evaluate the inertia in the asymptotic limit as \(t \to \infty\). In this limit, the \(N\)-soliton solution \(U^{(N)}\) decomposes into a superposition of \(N\) well-separated single solitons \(\phi_{c_j}\). Thus, the spectrum \(\sigma(\mathcal{L}_N)\) of the linearized operator converges to the union of the spectra of the individual operators \(L_{N,j}\), namely, \[\begin{align} \sigma(\mathcal{L}_N) \longrightarrow \bigcup_{j=1}^{N} \sigma(L_{N,j}), \quad \text{as } t \to \infty, \end{align}\] where \(L_{N,j} := \mathcal{S}_N''(\phi_{c_j})\) is the second variation of the action \(\mathcal{S}_N\) evaluated at the \(j\)-th soliton component. Our objective is to prove that \(in(\mathcal{L}_N) = (\left[ \frac{N+1}{2}\right], N)\).

By virtue of the Euler-Lagrange equation 12 and the symmetric properties of the Lagrange multipliers \(\mu_m\), the operator \(L_{N,j}\) can be expressed as a linear combination of the second variations of the conservation laws \(\{E_k\}_{k=0}^N\), which is \[\begin{align} L_{N,j} = (-1)^N \left( E_N''(\phi_{c_j}) + \sum_{m=0}^{N-1} \mu_m E_m''(\phi_{c_j}) \right). \end{align}\] Utilizing the iterative operator identity established in Proposition 3, the operator \(L_{N,j}\) admits a remarkable factorization in terms of the recursion operator \(\widehat{\mathcal{R}}(\phi_{c_j})\): \[\begin{align} \label{eq:factorization95LNj} L_{N,j} = (-1)^N \left( \prod_{k \neq j}^{N} (\widehat{\mathcal{R}}(\phi_{c_j}) + \omega_k) \right) \left( E_1''(\phi_{c_j}) + \omega_j E_0''(\phi_{c_j}) \right), \end{align}\tag{28}\] where \(\omega_k = 2/c_k\). To proceed, we first analyze the fundamental linearized operator around a single soliton \[\begin{align} L_1 = E_1''(\phi) + \omega E_0''(\phi) = -\sqrt{2\gamma}(1-\partial_x^2)\frac{1}{M^3}(1-\partial_x^2) + \frac{4}{c}(1-\partial_x^2). \end{align}\] The determination of the spectrum for the general operator \(L_{N,j}\) thus relies on analyzing the action of the recursion operator products in 28 on the eigenspaces of \(L_1\). This motivates the subsequent spectral analysis of the recursion operator \(\widehat{\mathcal{R}}(\phi)\) and its adjoint \(\widehat{\mathcal{R}}^*(\phi)\).

A significant consequence of the commutation relations ?? and ?? is that the (adjoint) recursion operator \(\widehat{\mathcal{R}}(\phi)\) (\(\widehat{\mathcal{R}}^*(\phi)\)) and the operator \(L_n \mathcal{J}\) (\(\mathcal{J} L_n\)) share a common set of eigenfunctions. Next, we will investigate the spectrum of the operator \(\mathcal{J} L_n\), which is analytically more accessible than that of \(L_n\). By leveraging ?? , the problem is reduced to the spectral analysis of the adjoint recursion operator \(\widehat{\mathcal{R}}^*(\phi)\). Then we demonstrate that the eigenfunctions of \(\widehat{\mathcal{R}}^*(\phi)\), augmented by the generalized kernel of \(\mathcal{J} L_n\), constitute an orthogonal basis in \(L^2(\mathbb{R})\), establishing a completeness relation. Finally, analyze the spectral information of \(\mathcal{L}_N\) based on the analysis results.

3.1 The spectrum of the recursion operator around the mCH single solitons↩︎

The spectrum of the recursion operator \[\begin{align} \widehat{\mathcal{R}}(\phi)=-\frac{2}{\gamma}\left(1-\partial_x^2\right)\left(1-\frac{m_{\phi}}{M_{\phi}}\partial_x^{-1}\frac{m_{\phi}}{M_{\phi}}\partial_x\right), \end{align}\] and its adjoint operator \(\widehat{\mathcal{R}}^*(\phi)\) are essential for analyzing the linearized Hamiltonian \(L_n\) in 28 . Although these recursion operators are nonlocal and difficult to investigate directly, we can characterize their spectral properties by employing the operator identities in Proposition 4 and the properties of the squared eigenfunctions \(F^\pm (x,k):=(\psi_1^\pm)^2(x,k)+(\psi_2^\pm)^2(x,k)\).

We now consider the adjoint recursion operator \(\widehat{\mathcal{R}}^*(\phi)\). In view of the commutation relation ?? , \(\widehat{\mathcal{R}}^*(\phi)\) shares the same eigenfunctions as \(\mathcal{J}L_n\), which is particularly relevant for the spectral stability analysis of solitons. Recall that: \[\begin{align} \widehat{\mathcal{R}}^*(\phi)=-\frac{2}{\gamma}\left(1-\frac{m_{\phi}}{M_{\phi}}\partial_x^{-1}\frac{m_{\phi}}{M_{\phi}}\partial_x\right)\left(1-\partial_x^2\right). \end{align}\] From 29 , it follows that: \[\begin{align} \label{self} \widehat{\mathcal{R}}^*(\phi)\mathcal{J}F^\pm(x,k)=-\lambda^2\mathcal{J}F^\pm(x,k). \end{align}\tag{31}\] The spectral properties of the adjoint operator \(\widehat{\mathcal{R}}^*(\phi)\) are summarized as follows.

3.2 The spectra of \(\mathcal{J}L_n\) and \(L_n\mathcal{J}\)↩︎

In this subsection, our attention is focused on the spectral analysis of the operators \(\mathcal{J}L_n\) and \(L_n\mathcal{J}\). The main ingredients are the commutation relations ?? and the observation that the eigenfunctions of the adjoint recursion operator \(\widehat{\mathcal{R}}^*(\phi)\) form an orthogonal basis in \(L^2(\mathbb{R})\). It follows that the spectrum of \(\mathcal{J}L_n\) lies on the imaginary axis, which implies directly the spectral stability of solitons.

We now consider the operator \(\mathcal{J}L_n\). Since \(L_n=\widehat{\mathcal{R}}^{n-1}(\phi)L_1\), we first examine the asymptotic constant-coefficient form (i.e., the principal part as \(|x| \to \infty\)) of \(L_1\), which is given by \(-\frac{4}{\gamma}\left(1-\partial_x^2\right)^2+\frac{4}{c}\left(1-\partial_x^2\right)\). Hence, the principal symbol of \(L_1\) is \[\begin{align} \rho_1(\xi)=-\frac{4}{\gamma}\left(1+\xi^2\right)^2+\frac{4}{c}\left(1+\xi^2\right). \end{align}\] Due to the relation \(L_n=(\widehat{\mathcal{R}}(\phi))^{n-1}L_1\), and considering that the symbol of \(\mathcal{J}=\partial_x(1-\partial_x^2)^{-1}\) is \(\frac{i\xi}{1+\xi^2}\) while the asymptotic symbol of \(\widehat{\mathcal{R}}\) is \(-\frac{2}{\gamma}(1+\xi^2)\), we can obtain the principal symbol of \(\mathcal{J}L_n\) as \[\begin{align} \label{34618} \varrho_n(\xi)=(-1)^{n-1}i\xi\left(\frac{2}{\gamma}\right)^{n-1}\left(1+\xi^2\right)^{n-2}\rho_1(\xi). \end{align}\tag{34}\] Based upon the facts above, we have the following statement related to the spectrum of the operator \(\mathcal{J}L_n\) for \(1\leq n\leq N\).

For the adjoint operator of \(\mathcal{J}L_n\), namely, the operator \(L_n\mathcal{J}\) which commutes with the recursion operator \(\widehat{\mathcal{R}}(\phi)\), we have the following result.

On account of Propositions 7 and 8, we now establish the following two function sets. The first set \[\begin{align} \label{set1} \left\lbrace \partial_x F^\pm(x,k),\quad \text{for }k\in\mathbb{R};\quad \partial_x \phi;\quad \frac{\partial\phi}{\partial c}\right\rbrace \end{align}\tag{41}\] consists of linearly independent eigenfunctions and the generalized kernel of the operator \(\mathcal{J}L_n\). Moreover, they are essentially orthogonal under the \(L^2\)-inner product. The second set \[\begin{align} \label{set2} \left\lbrace (1-\partial^2_x) F^\pm(x,k),\quad \text{for }k\in\mathbb{R};\quad m_\phi;\quad \partial_x^{-1}\left(\frac{\partial m_\phi}{\partial c}\right) \right\rbrace \end{align}\tag{42}\] consists of linearly independent eigenfunctions and the generalized kernel of the operator \(L_n\mathcal{J}\).

Given the even symmetry of \(\frac{\partial \phi}{\partial c}\), for \(k,l \in \mathbb{R}\), we can employ the asymptotic behaviors of the Jost solutions 16 to compute their inner products: \[\begin{align} &\int_\mathbb{R}\partial_xF^\pm(x,k) (1-\partial_x^2)\overline{F^\pm(x,l)}dx=\pm 2\pi i \mu(1+k^2)\left|a(k)\right|^2\delta(k-l),\label{34628}\\ &\int_\mathbb{R}\partial_x\phi \partial_x^{-1}\left(\frac{\partial m_{\phi}}{\partial c}\right)dx=-\int_\mathbb{R}\phi\frac{\partial m_{\phi}}{\partial c}dx=-\frac{1}{2}\frac{\partial E_0(\phi)}{\partial c}=-\sqrt{\frac{c}{c-\gamma}}<0,\\ &\int_\mathbb{R}\frac{\partial\phi}{\partial c}m_{\phi}dx=\frac{1}{2}\frac{\partial E_0(\phi)}{\partial c}=\sqrt{\frac{c}{c-\gamma}}>0. \end{align}\tag{43}\] The corresponding closure relation is thus given by \[\begin{align} \label{close} &\pm \int_\mathbb{R} \frac{1}{2\pi i\mu(1+k^2) |a(k)|^2} \partial_xF^\pm(x,k)(1-\partial_y^2)\overline{F^\pm(y,k)}\text{d}k \nonumber\\ &+ \sqrt{\frac{c-\gamma}{c}} \biggl[ -\partial_x\phi(x) \partial_y^{-1}\left(\frac{\partial m_{\phi}}{\partial c}\right)(y)+\frac{\partial \phi}{\partial c}(x)m_\phi(y) \biggr] = \delta(x-y). \end{align}\tag{44}\] This implies any function \(z(x)\) vanishing at infinity can be expanded over the two bases 41 and 42 . In particular, we have the following decomposition of the function \(z\): \[\begin{align} \label{decom} z(x) = \int_\mathbb{R} \partial_xF^\pm(x,k) \alpha^\pm (k)\text{d}k + \beta \partial_x \phi + \eta \frac{\partial \phi}{\partial c}, \end{align}\tag{45}\] where the coefficients \(\alpha^\pm(k), \beta\) and \(\eta\) can be computed explicitly via the closure relation 44 and the completeness relation (17 ). Similarly, the function \(z(x)\) decomposes over the second set 42 by multiplying 44 by \(z(x)\) and integrating with respect to \(x\).

3.3 The spectra of linearized operator around the mCH \(N\)-solitons↩︎

In order to prove Theorem 2, we need to know the spectral information of the operator \(\mathcal{L}_n\) with \(1\leq n\leq N\). More precisely, the inertia of \(\mathcal{L}_n\) has to be determined. The aim of this subsection is to show the following result.

Recall that \(L_{n,j}=S_{n}^{\prime\prime}(\phi_{c_{j}})\) is defined in (28 ). The spectrum of \(\mathcal{L}_{n}\) tends to the unions of \(\sigma(L_{n,j})\), that is \[\sigma(\mathcal{L}_{n}) \to \bigcup_{j=1}^{n}\sigma(L_{n,j}) \quad \text{as} \quad t\to+\infty.\] The result (46 ) follows directly from the following statement which concerning the inertia of the operators \(L_{n,j}\).

4 Proof of the stability results↩︎

This section is devoted to the proof of the stability theorems of this work. Specifically, by verifying the critical Maddocks-Sachs spectral condition \(\mathrm{in}(\mathcal{L}_N) = p(D_N)\), we establish the dynamical stability of the mCH \(N\)-soliton solutions. Furthermore, by employing a generalized Gårding inequality [49] and a refined modulation analysis, we rigorously prove the orbital stability in the optimal regularity space \(H^{N+1}(\mathbb{R})\).

4.1 Dynamical stability of mCH \(N\)-soliton solutions↩︎

In this subsection, we complete the proof of Theorem 2. The proof depends on the following fundamental result originally formulated by Maddocks and Sachs [27].

4.2 Orbital stability of mCH \(N\)-soliton solutions↩︎

We now establish the orbital stability of the mCH \(N\)-solitons (Theorem 3). Under the spectral condition \(n(\mathcal{L}_N)=p(D_N)=\left[\frac{N+1}{2}\right]\), the generalized Maddocks-Sachs framework guarantees the existence of a sufficiently large constant \(C>0\) such that the augmented Lagrangian \[\label{Delta} \Delta_N(u) = \mathcal{S}_N(u)+\frac{C}{2}\sum_{j=0}^{N-1}\bigl(E_j(u)-E_j(U^{(N)})\bigr)^2\tag{50}\] attains a non-degenerate unconstrained minimum at every translate of \(U^{(N)}\). Orbital stability requires establishing three elements: the temporal conservation of \(\Delta_N\), local coercivity near the \(N\)-soliton manifold, and a global continuation argument.

Conservation of \(\Delta_N (u)\). Since the Hamiltonians \(\{E_j(u)\}_{j=0}^N\) are conserved along the flow of 1 , the linear combination \(\mathcal{S}_N(u)\) and consequently \(\Delta_N(u(t))\) remain constant in time, serving as a valid Lyapunov functional.

Coercivity of \(\Delta_N (u)\) around the \(N\)-soliton manifold. Let \(\mathcal{M}_N = \bigl\{U^{(N)}(\cdot\,;\mathbf{y}) \mid \mathbf{y} \in\mathbb{R}^N\bigr\}\) be the \(N\)-dimensional soliton manifold. For any \(u\) sufficiently close to \(\mathcal{M}_N\) in \(H^{N+1}(\mathbb{R})\), the implicit function theorem guarantees a unique spatial shift \(\mathbf{y}=(y_1,\dots,y_N)\) such that the perturbation \(w(x) = u(x)-U^{(N)}(x;\mathbf{y})\) satisfies the orthogonality conditions: \[\label{ortho} \langle w,\partial_{y_j}U^{(N)}\rangle_{L^2} = 0,\qquad j=1,\dots,N.\tag{51}\] Furthermore, \(\|w\|_{H^{N+1}} \le K\,\mathrm{dist}_{H^{N+1}}(u,\mathcal{M}_N)\) for some constant \(K>0\). Expanding \(\Delta_N\) around \(U^{(N)}\) yields: \[\label{Delta95expand} \Delta_N(u)-\Delta_N(U^{(N)}) = \frac{1}{2}\langle\mathcal{L}_N w,w\rangle+R(w) +\frac{C}{2}\sum_{j=0}^{N-1}\Bigl(\langle E_j'(U^{(N)}),w\rangle+Q_j(w)\Bigr)^2,\tag{52}\] where \(\mathcal{L}_N=\mathcal{S}_N''(U^{(N)})\). The higher-order remainders satisfy \(|R(w)| \le C_1\|w\|_{H^{N+1}}^3\) and \(|Q_j(w)| \le C_2\|w\|_{H^{N+1}}^2\).

Gårding inequality for \(\mathcal{L}_N\). The \(N\) iterations of the recursion operator \(\widehat{\mathcal{R}}\) yield a leading differential part for \(E_N''(U^{(N)})\) proportional to \(-\partial_x^{2N+2}\). Incorporating the global \((-1)^N\) parity correction precisely adjusts the leading term of \(\mathcal{L}_N\) to \((-1)^{N+1}\partial_x^{2N+2}\), ensuring strict ellipticity. Explicitly, \(\mathcal{L}_N\) forms a strongly elliptic \(2(N+1)\)-th order self-adjoint operator. Gårding’s inequality provides the lower bound: \[\langle\mathcal{L}_N w,w\rangle \ge \delta\|\partial_x^{N+1}w\|_{L^2}^2 - C\|w\|_{H^N}^2.\] Absorbing intermediate derivatives via interpolation \(\|w\|_{H^N}^2 \le \eta\|\partial_x^{N+1}w\|_{L^2}^2 + C_\eta\|w\|_{L^2}^2\) with sufficiently small \(\eta\) yields the final Gårding inequality: \[\label{Garding95final} \langle\mathcal{L}_N w,w\rangle \ge \gamma_1\|w\|_{H^{N+1}}^2 - \gamma_2\|w\|_{L^2}^2,\tag{53}\] with positive constants \(\gamma_1, \gamma_2\). Combined with the exact spectral structure (i.e., \([\frac{N+1}{2}]\) negative eigenvalues, an essential spectrum bounded away from zero, and an \(N\)-dimensional kernel \(X_0\)), standard compactness arguments (e.g., [27]) upgrade 53 to a full coercivity estimate on the positive subspace: \[\label{coercive95HN} \langle\mathcal{L}_N w,w\rangle \ge \gamma\|w\|_{H^{N+1}}^2 \qquad\text{for all } w \perp X_0 \text{ and } w \perp X_{neg}.\tag{54}\]

Estimate of the constraint terms. Let \(\ell_j(w)=\langle E_j'(U^{(N)}),w\rangle\). Since the gradients \(\{E_j'(U^{(N)})\}\) uniquely span a subspace covering \(X_{neg}\), there exists \(c_0>0\) such that \[\label{ell95estimate} \sum_{j=0}^{N-1}|\ell_j(w)|^2 \ge c_0\,\|P_{neg} w\|_{L^2}^2.\tag{55}\] For sufficiently small \(\|w\|_{H^{N+1}}\), binomial expansion gives: \[\label{constraint95expand} \bigl(\ell_j(w)+Q_j(w)\bigr)^2 \ge \frac{1}{2}|\ell_j(w)|^2-C_3\|w\|_{H^{N+1}}^4.\tag{56}\]

Completion of the coercivity estimate. Inserting 54 , 55 , and 56 into 52 , and utilizing \(w \perp X_0\) from 51 , we obtain for small \(\|w\|_{H^{N+1}}\): \[\Delta_N(u)-\Delta_N(U^{(N)}) \ge \frac{\gamma}{2}\|w\|_{H^{N+1}}^2+\frac{Cc_0}{4}\|P_{neg}w\|_{L^2}^2 -C_4\|w\|_{H^{N+1}}^3.\] Choosing \(\|w\|_{H^{N+1}}\) small enough allows the quadratic terms to absorb the cubic remainder, yielding a uniform constant \(\alpha>0\) such that: \[\label{coercive95HN95final} \Delta_N(u)-\Delta_N(U^{(N)}) \ge \alpha\|w\|_{H^{N+1}}^2 \ge \frac{\alpha}{K^2}\,\mathrm{dist}_{H^{N+1}}(u,\mathcal{M}_N)^2.\tag{57}\]

Orbital stability in \(H^{N+1}\). For initial data satisfying \(\|u_0-U^{(N)}(0)\|_{H^{N+1}} < \delta\), we dynamically choose the modulation parameters \(\mathbf{y}(t)\) to continuously satisfy the orthogonality conditions 51 . By virtue of the temporal conservation of \(\Delta_N\) and the coercivity estimate 57 , we obtain for all \(t \ge 0\): \[\begin{align} \alpha\|w(t)\|_{H^{N+1}}^2 &\le \Delta_N(u(t))-\Delta_N(U^{(N)})\le L\,\|u_0-U^{(N)}(0)\|_{H^{N+1}}^2, \end{align}\] where \(L\) is the local Lipschitz constant for \(\Delta_N\). This yields the uniform bound \(\|w(t)\|_{H^{N+1}} \le C_0\delta\) with \(C_0=\sqrt{L/\alpha}\).

For any \(\varepsilon > 0\), by choosing \(\delta < \varepsilon / C_0\), we ensure that the perturbation remains strictly bounded by \(\varepsilon\) globally in time. Therefore, the solution \(u(t)\) permanently stays within the \(\varepsilon\)-tube of the \(N\)-soliton manifold \(\mathcal{M}_N\), which directly leads to the orbital stability: \[\inf_{\mathbf{r}\in\mathbb{R}^N}\|u(t)-U^{(N)}(t;\cdot+r_1,\dots,\cdot+r_N)\|_{H^{N+1}} < \varepsilon , \qquad\forall t\ge 0.\] This completes the proof of Theorem 3. \(\square\)

This work was supported by the National Natural Science Foundation of China under Grant No. 12371255, the Fundamental Research Funds for the Central Universities of CUMT under Grant No. 2024ZDPYJQ1003, and the Postgraduate Research & Practice Program of Education & Teaching Reform of CUMT under Grant No. 2025YJSJG031.

No data was used for the research described in the article.

The authors declare no conflict of interest.

References↩︎