High-Frequency Instabilities of a Boussinesq-Whitham System: A Perturbative Approach

We analyze the spectral stability of small-amplitude, periodic, traveling-wave solutions of a Boussinesq-Whitham system. These solutions are shown numerically to exhibit high-frequency instabilities when subject to bounded perturbations on the real line. We use a formal perturbation method to estimate the asymptotic behavior of these instabilities in the small-amplitude regime. We compare these asymptotic results with direct numerical computations. This is the second paper in a series of three that investigates high-frequency instabilities of Stokes waves.


Introduction
We investigate small-amplitude, 2π/κ-periodic, traveling waves of a Boussinesq-Whitham system proposed by Hur and Pandey [15] and Hur and Tao [16]: (1.1) In this model, η(x, t) represents the displacement of a wave profile from its equilibrium depth h 0 , u(x, t) is the horizontal velocity along η, and K is a Fourier multiplier operator defined so that the linearized dispersion relation of (1.1) matches that of the Euler water wave problem (WWP) [27]. For functions f ∈ L 1 per (−π/κ, π/κ), K is defined as where · denotes the Fourier transform of f : Alternatively, K can be defined in physical variables as the pseudo-differential operator where D = −i∂ x . For the remainder of this manuscript, we refer to (1.1) as the Hur-Pandey-Tao-Boussinesq-Whitham system, or HPT-BW for short. The HPT-BW system is Hamiltonian [15] with H = 1 2 WWP of the same name [22,24,25]. In Section 2, we derive a power series expansion for HPT-BW Stokes waves in a small parameter ε that scales with the amplitude of the waves. Perturbing Stokes waves yields a spectral problem after linearizing the governing equations of the perturbations. The spectral elements of this problem define the stability spectrum of Stokes waves; see Section 3. The stability spectrum is purely continuous [7,18,23], but Floquet theory decomposes the spectrum into an uncountable union of point spectra. Each of these point spectra is indexed by the Floquet exponent [11,14,17].
The stability spectrum inherits quadrafold symmetry from the Hamiltonian structure of (1.1), i.e., the spectrum is invariant under conjugation and negation [14,18]. Because of quadrafold symmetry, all elements of the stability spectrum have non-positive real component only if the stability spectrum is a subset of the imaginary axis. Therefore, HPT-BW Stokes waves are spectrally stable only if the spectrum is on the imaginary axis. Otherwise, the Stokes waves are spectrally unstable.
If the aspect ratio κh 0 is sufficiently large, both HPT-BW [15] and WWP Stokes waves [3,4,5] have stability spectra near the origin that leave the imaginary axis for 0 < |ε| ≪ 1, resulting in modulational instability. Using the Floquet-Fourier-Hill (FFH) method [11], recent numerical work by [8] and [12] shows, respectively, that HPT-BW and WWP Stokes waves also have stability spectra away from the origin that leave the imaginary axis, regardless of κh 0 . These spectra give rise to the so-called high-frequency instabilities [13], shown schematically in Figure 1.
High-frequency instabilities arise from the collision of nonzero stability eigenvalues of zero-amplitude (ε = 0) Stokes waves. At these collided spectral elements, a Hamiltonian-Hopf bifurcation occurs, resulting in a locus of spectral elements bounded away from the origin that leave the imaginary axis as |ε| increases. We refer to this locus of spectral elements as a high-frequency isola.
High-frequency isolas are difficult to find using numerical methods like FFH. To capture the isola closest to the origin, for example, the interval of Floquet exponents that parameterizes the isola has width O ε 2 ( Figure 2). For isolas further from the origin, this width appears to decay geometrically in ε. To compound these difficulties, the isolas drift away from their initial collision sites (Figure 2), meaning that the Floquet exponent that gives rise to the collided spectral elements at ε = 0 is not contained in the interval parameterizing the corresponding isola for |ε| sufficiently large. To circumvent these difficulties, one must supply the numerical method with asymptotic expressions for the interval of Floquet exponents corresponding to the desired isolas. We discover these expressions for high-frequency isolas of HPT-BW. The high-frequency isola closest to the origin of HPT-BW Stokes waves with κ = g = h 0 = 1 and ε = 2.5 × 10 −4 (red), 5 × 10 −4 (burgundy), 7.5 × 10 −4 (purple), and 10 −3 (blue). The collision that generates the isola when ε = 0 is indicated by the red star. The imaginary axis is recentered to show the magnitude of the drift of the isola from the collision site λ 0 . (Right) Interval of Floquet exponents that parameterize the isola on the left as a function of ε. The solid black lines indicate the boundaries of the interval, while the dashed black line gives the Floquet exponent of the most unstable spectral element of the isola. The colored dots provide the Floquet data for the correspondingly colored isola in the left plot. The red star indicates the Floquet exponent µ 0 that generates the isola when ε = 0. The Floquet axis is recentered to show the magnitude of the drift of the Floquet exponents from µ 0 .
Our motivation for studying the HPT-BW system, apart from its inherent interest, is that it retains the full dispersion relation (both branches) of the more complicated WWP. Our goal is the application of the perturbation method developed herein to the WWP. The first step towards this goal was the investigation of the stability spectra of Stokes waves of the Kawahara equation [9]. The investigations presented here constitute our second step, before proceeding to the finite-depth WWP next [10].
For a given high-frequency isola of an HPT-BW Stokes wave, we obtain (i) an asymptotic range of Floquet exponents that parameterize the isola, (ii) an asymptotic estimate for the most unstable spectral element of the isola, (iii) expressions of curves that are asymptotic to the isola, (iv) wavenumbers for which the given isola is not present. Our approach is inspired by a perturbation method outlined in [2], but modified appropriately for higher-order calculations. We compare all asymptotic results with numerical results computed by the FFH method.

Small-Amplitude Stokes Waves
In a traveling frame moving with velocity c, x → x − ct and (1.1) becomes √ gh 0 yields the following system: The parameter α is chosen to map 2π/κ-periodic solutions of (2.1) to 2π-periodic solutions of (2.2). Consequently, α = κh 0 > 0, the aspect ratio of the solutions, and or, alternatively, for f ∈ L 1 per (−π, π) with the Fourier transform (1.3) redefined over (−π, π). Stokes wave solutions of (2.2) are independent of time. Equating time derivatives in (2.2) to zero and integrating in x, we find where I j are integration constants. For each α > 0, there exists a three-parameter family of infinitely differentiable, even, small-amplitude, 2π-periodic solutions of (2.5), provided I j are sufficiently small [15]. We call these solutions the HPT-BW Stokes waves, denoted (η S (x; ε, I j ), u S (x; ε, I j )) T , where ε is a small-amplitude parameter defined implicitly in terms of the first Fourier mode of η S (x; ε, I j ): Remark. Redefining c → c − I 1 and u → u − I 1 in (2.5) implies I 1 = 0 without loss of generality. If we also equate I 2 = 0, our Stokes waves reduce to a one-parameter family of solutions to (2.5) such that (2.6) ensures that η S (x; ε) ∼ ε cos(x) as ε → 0. We restrict to this case for simplicitly, but the methodology in Sections 4 and 5 are unchanged if I 2 = 0. For series representations of Stokes waves that include I j , see [15].
The Stokes waves and their velocity may be expanded as power series in ε: where η j (x) and u j (x) are analytic, even, and 2π-periodic for each j. Substituting these expansions into (2.5) (with I j = 0) and following a Poincaré-Lindstedt perturbation method [27], one determines η j (x), u j (x), and c 2j order by order. In Appendix A, we report expansions of η S , u S , and c up to fourth order in ε; this is sufficient for our asymptotic calculations of high-frequency isolas discussed in Sections 4 and 5.

The Stability Spectrum of Stokes Waves
Consider perturbations to (η S , u S ) T of the form where |ρ| ≪ 1 is a parameter independent of ε and H and U are sufficiently smooth, bounded functions of x on R for all t ≥ 0. When (3.1) is substituted into (2.2), terms of O ρ 0 cancel by (2.5) (with I j = 0). Equating terms of O(ρ), the perturbation (H, U ) T solves the linear system ∂ ∂t where primes denote differentiation with respect to x. Formally separating variables, where (H, U) T solves the spectral problem Since the entries of the matrix operator above are 2π-periodic, one can use Floquet theory 1 to solve (3.4) for (H, U) T . These solutions take the form where µ ∈ [−1/2, 1/2] is called the Floquet exponent and h, u ∈ H 1 per (−π, π). Substituting (3.5) into (3.4) results in a spectral problem for w = (h, u) T : For sufficiently small ε, (3.6) has a countable collection of eigenvalues λ ε,µ for each Floquet exponent µ [17]. The union of these eigenvalues over µ ∈ [−1/2, 1/2] recovers the purely continuous spectrum of (3.4) for fixed ε; this is the stability spectrum of HPT-BW Stokes waves. We use the Floquet-Fourier-Hill method [11] to compute the stability spectrum numerically. If there exists a µ such that there is a λ ε,µ with Re (λ ε,µ ) > 0, then there exists a perturbation (3.3) that grows exponentially in time, and Stokes waves of amplitude ε are spectrally unstable. If no such µ is found, then the Stokes waves are spectrally stable. Because of the quadrafold symmetry mentioned in the introduction, Stokes waves are spectrally stable if and only if their stability spectrum is a subset of the imaginary axis.
For some µ = µ 0 , nonzero eigenvalues of L 0,µ0 with double multiplicity may give rise to Hamiltonian-Hopf bifurcations and, thus, to high-frequency instabilities for 0 < |ε| ≪ 1. These eigenvalues exist provided there exists µ 0 , m, and n such that λ (σ1) 0,µ0,n = λ (σ2) 0,µ0,m = 0. (3.11) We view (3.11) as a collision of two simple, nonzero eigenvalues. It can be shown that such a collision occurs only if σ 1 = σ 2 [1,13,15]. Theorem 4 in Appendix B shows that, for any p ∈ Z \ {0, ±1}, there exist unique µ 0 , m, and n that satisfy (3.11) with m − n = p. Thus, there are a countably infinite number of nonzero eigenvalue collisions in the zero-amplitude stability spectrum; each of which has potential to develop a high-frequency instability in the small-amplitude stability spectrum.
Remark. Using results in Appendix B, it can be shown that the Krein signatures [20] of the colliding eigenvalues have opposite signs. This is a second necessary criterion for the occurance of high-frequency instabilities [13,21].
Remark. The WWP shares the same collided eigenvalues with HPT-BW, since (3.9) is also the dispersion relation of the WWP.

High-Frequency Instabilities: p = 2
We use perturbation methods to investigate the high-frequency instability that develops from the collision of λ The p = 2 isola develops from the spectral data where γ j are arbitrary, nonzero constants. As |ε| increases, we assume the spectral data vary analytically [1] with ε: where we suppress functional dependencies for ease of notation. We normalize w so that or, alternatively, so that This normalization ensures that h 0 fully resolves the n th Fourier mode of h, a convenient choice for the perturbation calculations that follow. With this normalization, The arbitrary constant γ 0 will be determined at higher order, leading to a unique expression for w 0 .
Remark. The eigenvalue corrections λ j derived below are independent of the normalization chosen for w.
If λ 0 is a semi-simple, isolated eigenvalue of L 0,µ0 , we may justify (4.2) using analytic perturbation theory [19], provided the Floquet exponent is fixed. For ε sufficiently small, this method of proof gives two spectral elements on the isola. Numerical and asymptotic calculations show that these spectral elements quickly leave the isola as a result of the change in its Floquet parameterization with ε ( Figure 2, Figure 5, Figure 10). To account for this variation, we allow the Floquet exponent to depend on ε as well: Remark. In the calculations that follow, explicit expressions of select quantities are suppressed for ease of readability. The interested reader may consult the supplemental Mathematica file hptbw_isolap2.nb for these expressions.

The O (ε) Problem
Substituting the expansions of the Stokes wave (2.7), spectral data (4.2), and Floquet exponent (4.6) into the spectral problem (3.6) and collecting terms of O(ε), we find . (4.8) The inhomogeneous terms on the RHS of (4.7) can be evaluated using expressions for η 1 , u 1 , and w 0 . Each of these quantities are finite linear combinations of 2π-periodic sinusoids. As a result, the inhomogeneous terms can be rewritten as a finite Fourier series, and (4.7) becomes where T 1,j depend on µ 0 , α, and γ 0 ; see the Mathematica file for details.
With the solvability conditions satisfied, we solve for the particular solution of w 1 in (4.9). Combining with the nullspace of L 0,µ0 − λ 0 , where β 1,j are arbitrary constants and W 1,j are found in the Mathematica file. Enforcing the normalization condition (4.4), one finds β 1,n = 0. For ease of notation, let β 1,m → γ 1 so that Using (4.14), the spectral problem (3.6) at O ε 2 is where L 1 | µ1=0 is the same as above, but evaluated at µ 1 = 0, and One can evaluate the inhomogeneous terms of (4.18) using η j , u j , and w j−1 for j ∈ {1, 2}. These inhomogeneous terms can be expressed as a finite Fourier series, giving It can be shown that T 2,n−1 = 0.
Proceeding similarly to the previous order, solvability conditions for (4.19) are Expressions for S 2,j and P 2,j have no dependence on γ 0 , γ 1 , µ 2 , or λ 2 ; see the attached Mathematica file for details. Conditions (4.20a) and (4.20b) form a nonlinear system for γ 0 and λ 2 . Solving for λ 2 yields A direct calculation shows that where S 2 is given in the attached Mathematica file. Then, and That c g−1 (µ 0 + m) = c g1 (µ 0 + n) follows from Lemma 2 in Appendix B. A plot of S 2 as a function of α suggests that S 2 > 0 for all values of α > 0 ( Figure 3). We conjecture that HPT-BW Stokes waves of any wavenumer experience a p = 2 high-frequency instability at O ε 2 . For µ 2 ∈ (M 2,− , M 2,+ ), a quick calculation shows that (4.24) parameterizes an ellipse asymptotic to the numerically observed p = 2 high-frequency isola (Figure 4). The ellipse has semi-major and -minor axes that scale with ε 2 , and the center of the ellipse drifts along the imaginary axis like ε 2 from λ 0 , the collision point at ε = 0.
The midpoint of (M 2,− , M 2,+ ) maximizes the real part of λ 2 . Thus, the most unstable spectral element of the isola has Floquet exponent and its real and imaginary components are respectively. These expansions agree well with the FFH results, see Figure 5.

High-Frequency Instabilities: p = 3
According to Theorem 4 in Appendix B, the p = 3 high-frequency instability is the second-closest to the origin. As will be seen, this instability arises at O ε 3 . Let µ 0 correspond to the unique Floquet exponent in [−1/2, 1/2] that satisfies the collision condition (3.11) with m − n = 3. Then, the spectral data (4.1) give rise to the p = 3 high-frequency instability. We assume these data and the Floquet exponent vary analytically with ε. For uniqueness, we normalize the eigenfunction w according to (4.4) so that w 0 is given by (4.5). We proceed as in the p = 2 case.
Remark. In the calculations that follow, explicit expressions of select quantities are suppressed for ease of readability. The interested reader may consult the supplemental Mathematica file hptbw_isolap3.nb for these expressions.

The O (ε) Problem
Substituting expansions (2.7), (4.2), and (4.6) into the spectral problem (3.6), equating terms of O(ε), and using expression for η 1 , u 1 , and w 0 to simplify, we find where γ 1 is arbitrary and expressions for W 1,j are found in the supplemental Mathematica file. Because T 1,n−1 and T 1,m+1 are identical to their p = 2 counterparts, W 1,n−1 and W 1,m+1 are as well.
We conjecture that HPT-BW Stokes waves with α = 1.1862... are not succeptible to the p = 3 instability, even beyond O ε 3 . Indeed, in the next subsection, we find that λ 4 is purely imaginary, so Stokes waves with aspect ratio α = 1.1862... do not exhibit p = 3 instabilities to O ε 4 . Assuming α = 1.1862..., µ 3 ∈ (−M 3 , M 3 ) parameterizes an ellipse asymptotic to the p = 3 highfrequeny isola; see Figure 7. The ellipse has semi-major and -minor axes that scale with ε 3 . The center of this ellipse drifts along the imaginary axis like ε 2 due to the purely imaginary correction found at O ε 2 .
The interval of Floquet exponents that parameterizes the p = 3 isola is (5.14) The width of this interval is an order of magnitude smaller than that of the p = 2 isola. Consequently, the p = 3 isola is more challenging to find numerically than the p = 2 isola, at least for methods similar to FFH (Table 1). For α = 1 and |ε| < 5 × 10 −4 , (5.14) provides an excellent approximation to the numerically computed interval of Floquet exponents (Figure 8). Fourth-order corrections are necessary to improve agreement between (5.14) and numerical computations for larger ε, see Section 5.4 below.
Choosing µ 3 = 0 maximizes the real part of λ 3 . Thus, the most unstable spectral element of the p = 3 isola has Floquet exponent Table 1. Intervals of Floquet exponents that parameterize the p = 2 and p = 3 high-frequency isolas with ε = 10 −3 and α = 1/2, 1, and 2. The first digit for which the boundary values disagree is underlined and colored red. If a uniform mesh of Floquet exponents in [−1/2, 1/2] is used for numerical methods like FFH, the spacing of the mesh must be finer than ε 2 to capture the p = 2 instability and ε 3 to capture the p = 3 instability. The intervals vary with α as well, making it difficult to adapt and refine a uniform mesh to find high-frequency isolas.
where µ 2 is as in (5.5b), and its real and imaginary components are   The real (blue) and imaginary (red) components of the most unstable spectral element of the isola as a function of ε (zero-and second-order imaginary and Floquet corrections removed for better visibility). Solid curves illustrate perturbation calculations. Circles illustrate FFH results.
respectively. The expansion for λ r, * is in excellent agreement with numerical results using the FFH method ( Figure 8). As with (5.14), corrections to µ * and λ i, * at O ε 4 improve the agreement between numerical and asymptotic results for these quantities. Before proceeding to O ε 4 , we solve (5.9) for w 3 , assuming solvability conditions (5.4a) and (5.4b) and normalization condition (4.4) are satisfied. We find where γ 3 is arbitrary and W 3,n−2 = 0 (since T 3,n−2 = 0).
We refer to this equality as the regular curve condition: it ensures that the curve asymptotic to the p = 3 isola is continuous near its intersections with the imaginary axis. From the regular curve condition, we get As expected, the Floquet parameterization and imaginary component of the p = 3 isola have a nonzero correction at O ε 4 . These corrections improve the agreement between numerical and asymptotic results observed at the previous order (Figure 9, Figure 10). No corrections to the real component of the isola are found at fourth order.

Conclusions
We have extended a formal perturbation method, first introduced in [2], to obtain asymptotic behavior of the largest (p = 2, 3) high-frequency instabilities of small-amplitude, HPT-BW Stokes waves. In particular, we have computed explicit expressions for (i) the interval of Floquet exponents that asymptotically parameterize the p th isola, (ii) the leading-order behavior of its most unstable spectral elements, (iii) the leading-order curve asymptotic to the isola, and (iv) wavenumbers that do not have a p th isola. Items (i)-(iii) can be extended to higher-order if necessary using the regular curve condition. In all instances, our perturbation calculations are in excellent agreement with numerical results computed by the FFH method [11].
We restrict to p = 2 and p = 3 in this work, but our method can provide asymptotic expressions for p > 3 isolas. We conjecture that this method yields the first real-component correction of the isola at O (ε p ), similar to the cases p = 2 and p = 3. If correct, this conjecture highlights the main difficulty of computing higher-order high-frequency instabilities, both numerically and perturbatively.
The asymptotic expressions derived in this paper are intimidating and cumbersome. Although it is satisfying to have asymptotic expressions for the results previously obtained only numerically, this is not the main point of our work. Rather, (i) the perturbation method demonstrated allows one to approximate an entire isola at once, going beyond standard eigenvalue perturbation theory [19], (ii) the results obtained constitute a first step toward a proof of the presence of the high-frequency instabilities, and (iii) the asymptotic expressions for the range of Floquet exponents allow for a far more  Acknowledgements: This research was funded partially by the ARCS Foundation Fellowship.

Appendix B. Collision Condition
Up to redefining m and n, (3.11) simplifies to With k = µ 0 + n and p = m − n, (8.1) becomes We refer to (8.2) as the collision condition. We prove that, for each p ∈ Z \ {0, ±1}, there exists a unique k(p; α) that satisfies the collision condition. These solutions k(p; α) are distinct from each other (for each α > 0) and result in an infinite number of distinct collision points on the imaginary axis, according to (3.11). First, we establish important monotonicity properties of Ω σ (k), defined in (3.9).
The proof of the third claim is immediate since {Ω 1 (k(p; α))} is strictly increasing.