INTRODUCTION

The forced generation of waves from a boundary is a fundamental problem in applied mathematics with broad applications in fluid dynamics, 1 material science, 2 seismology, 3 and astronomy 4 just to name a few. The fundamental mathematical problem can be posed as follows. Consider a linear evolutionary equation subject to time-harmonic forcing at a localized source or boundary ∝ 𝑒 -𝑖𝜔 0 𝑡 and decay at infinity. We would expect the solution to exhibit spatial dependence via the wavevector 𝐤 0 consistent with the system's dispersion relation such that 𝑢(𝐱, 𝑡) ∝ 𝑒 𝑖(𝐤 0 ⋅𝐱-𝜔 0 𝑡) near the source or boundary. However, for any evolutionary equation whose spatial order is higher than one, the wavevector 𝐤 0 is not unique. In the context of the time-independent Helmholtz equation (∇ 2 + 𝑘 2 0 )𝑢 = -𝑓(𝐱), Sommerfeld's eponymous radiation condition in three-dimensional space, lim 𝑟→∞ 𝑟( 𝜕𝑢 𝜕𝑟 -𝑖𝑘 0 𝑢) = 0, was formulated to select the unique wave that propagates away from a source or boundary to ensure well-posedness on an infinite domain. 5 Another, closely related approach to this problem is called the limiting absorption principle for the resolvent 𝑅(𝑘 0 ) = (Δ + 𝑘 2 0 ) -1 . 6 The unique solution is selected by taking a limit of the resolvent whose spectrum is perturbed off the real, essential spectrum into the appropriate half of the complex plane. Ultimately, this boils down to selecting the appropriate branch of the square root.

Both the radiation condition and the limiting absorption principle are mathematical workarounds to the fundamental problem that the Helmholtz equation exhibits two linearly independent, decaying solutions (in two or more dimensions). A physically inspired approach to the problem is to consider a wavemaker in which a boundary or source time-harmonically forces an initially quiescent (𝑢(𝐱, 0) = 0) medium. In the context of the forced wave equation 𝑢 𝑡𝑡 -Δ𝑢 = -𝑒 -𝑖𝜔 0 𝑡 𝑓(𝐱), solving this initial-boundary value problem and taking 𝑡 → ∞ is known as the limiting amplitude principle, which selects the same outgoing wave as the Sommerfeld radiation condition and corresponds to the causal solution by perturbing the essential spectrum of the Helmholtz resolvent into the upper half of the complex plane. 7,8 A simple, physically inspired calculation 1 demonstrates the limiting amplitude principle. Consider a linear, first order in time evolutionary equation for 𝑢(𝑥, 𝑡) on the half-line 𝑥 > 0 with real-valued dispersion relation 𝜔 = 𝜔(𝑘), the asymptotic Dirichlet boundary condition 𝑢(0, 𝑡) ∼ 𝑒 -𝑖𝜔 0 𝑡 , 𝑡 → ∞ with 𝜔 0 ∈ ℝ, and decay lim 𝑥→∞ 𝑢(𝑥, 𝑡) = 0 for each 𝑡 ∈ ℝ. Anticipating the wave 𝑢(𝑥, 𝑡) ∼ 𝑒 𝑖(𝑘 0 𝑥-𝜔 0 𝑡) , 𝑡 → ∞, 𝑥 = (1),

(1)

we seek the correct wave number branch 𝑘 0 = 𝑘(𝜔 0 ) selected from the multivalued inverse 𝑘(𝜔 0 ) = 𝜔 -1 (𝜔 0 ). The limiting amplitude principle states that, initially, the solution is zero, 𝑢(𝑥, 0) = 0 for each 𝑥 > 0. Consequently, the causal solution can be modeled by perturbing the frequency into the upper half of the complex plane 𝜔 0 → 𝜔 0 + 𝑖𝜖, 0 < 𝜖 ≪ |𝜔 0 | so that the boundary condition becomes 𝑢(0, 𝑡) = 𝑒 -𝑖𝜔 0 𝑡 𝑒 𝜖𝑡 with initially slow amplitude growth to eventually attain the asymptotically periodic boundary condition. Assuming that 𝑘(𝜔) is analytic in a neighborhood of 𝜔 0 , we expand 𝑘(𝜔 0 + 𝑖𝜖) = 𝑘(𝜔 0 ) + 𝑖𝜖 d𝑘 d𝜔 (𝜔 0 ) + ⋯ = 𝑘 0 + 𝑖 𝜖 𝑐 𝑔 (𝜔 0 ) + ⋯ ,

where 𝑐 𝑔 ≡ d𝜔∕d𝑘 is the group velocity. Then, the anticipated solution takes the form 𝑢(𝑥, 𝑡) ∼ 𝑒 𝜖𝑡 𝑒 𝑖(𝑘(𝜔 0 )𝑥-𝜔 0 𝑡) = 𝑒 𝜖[𝑡-𝑥∕𝑐 𝑔 (𝜔 0 )]+⋯ 𝑒 𝑖(𝑘 0 𝑥-𝜔 0 𝑡) .

Since 𝑢(𝑥, 𝑡) ∼ 𝑒 (𝑖𝑘 0 -𝜖∕𝑐 𝑔 )𝑥-(𝑖𝜔 0 -𝜖)𝑡 , in order to ensure lim 𝑥→∞ 𝑢(𝑥, 𝑡) = 0 in (3) as 𝜖 → 0 + , we require 𝑘 0 ∈ ℝ and 𝑐 𝑔 (𝜔 0 ) > 0, or Im 𝑘 0 > 0.

These inequalities are the radiation condition for the wavemaker problem. The positive group velocity for real wave numbers corresponds to waves being emitted from the boundary into the domain. For complex 𝑘, the positive imaginary part implies spatial decay. The boundary case 𝑐 𝑔 (𝜔 0 ) = 0 is not included in this argument because 𝑘(𝜔) has a branch point at 𝜔 0 ; therefore, it is not analytic. Zero group velocity corresponds to a critical value of the frequency that represents a transition between propagating and spatially decaying waves. We will prove the general radiation condition 𝑘 0 ∈ ℝ and 𝑐 𝑔 (𝜔 0 ) ≥ 0, or Im 𝑘 0 > 0.

Although we have given a physically inspired, plausible argument for the radiation condition (4), its proof is more involved.

In this paper, we prove that the radiation condition (4) uniquely determines the asymptotic form of the solution (1) by constructing the Dirichlet-to-Neumann (D-N) map for the wavemaker initial-boundary value problem

𝑢 𝑡 = 𝑖𝐴 0 𝑢 + 𝐴 1 𝑢 𝑥 + 𝑖𝐴 2 𝑢 𝑥𝑥 + 𝐴 3 𝑢 𝑥𝑥𝑥 , 𝑥,𝑡 > 0,

𝑢(𝑥, 0) = 0, 𝑥 > 0,

𝑢(0, 𝑡) = 𝑔 0 (𝑡), 𝑡 > 0,

where 𝐴 -2 ≥ 0 to avoid instability, 𝐴 0 , 𝐴 1 , 𝐴 2 ∈ ℝ, and 𝐴 3 ≤ 0 so that only one boundary condition is required at 𝑥 = 0. By scaling 𝑡 and 𝑥, we can restrict 𝐴 3 ∈ {-1, 0} and 𝐴 -2 ∈ {0, 1} without loss of generality. For real wave numbers, these restrictions guarantee the real-valued dispersion relation

The choice of parameters {𝐴 -2 = 𝐴 0 = 𝐴 2 = 0, 𝐴 1 = 𝐴 3 = -1} corresponds to the linear Korteweg-de Vries (KdV) equation, also known as the first Stokes equation, whereas for {𝐴 -2 = 1, 𝐴 1 = -1, 𝐴 0 = 𝐴 2 = 𝐴 3 = 0}, Equation (5a) becomes the linear Benjamin-Bona-Mahoney (BBM) equation. We denote 𝜔(𝑘) as the dispersion relation and distinguish it from the input frequency 𝜔 0 at the boundary. In 1997, Fokas introduced a new approach for the analytical solution of initial-boundary value problems. 9 The main motivation for Fokas's method, otherwise known as the unified transform, was twofold:

(i) to provide an analog of the inverse scattering transform method, which is used for solving the initial value problem of completely integrable systems like the KdV and NLS equations, when such systems are considered in domains that involve a boundary, for example, the half-line or the finite interval, thereby giving rise to nonlinear initial-boundary value problems; (ii) to provide a natural analog of the Fourier transform, which yields the solution to linear evolution equations in the initial value problem setting, when these linear equations are instead formulated in the setting of initial-boundary value problems.

While the difficulty associated with task (i) can be easily appreciated due to the presence of nonlinearity, task (ii) comes with significant challenges as well, despite the fact that it concerns linear equations. This becomes evident through the realization that, for general nonzero boundary conditions, no classical spatial transform is available for solving initial-boundary value problems that involve linear evolution equations of spatial order higher than two. Moreover, even for second-order linear equations, classical techniques fail when the problem is nonseparable. The Fokas method effectively addresses both of the above objectives as a unified way of treating initial-boundary value problems for a plethora of linear and nonlinear equations, formulated in various domains and supplemented with different kinds of boundary conditions. Since its inception, the method has been further developed within the completely integrable systems and applied analysis communities in various directions such as the analytical solution of nonlinear equations via the Riemann-Hilbert/dbar problem approach, the solution of linear evolution and elliptic equations, asymptotics, and numerical studies; see the books 10,11 as well as the review articles 12,13 for a detailed list of references on the above topics and beyond.

It should be emphasized that the solution formulae obtained via the unified transform have significant advantages even when the relevant problems can be solved via classical approaches. Importantly, those portions of the formulae that involve the boundary data possess integrals with exponentially decaying integrands, due to the fact that integration takes place along complex contours in the spectral plane. As a result, the solution formulae are uniformly convergent up to the boundary of the domain-a feature that turns out particularly helpful for our purposes here.

In this work, we employ the Fokas method to prove the existence and uniqueness of the D-N map for the wavemaker problem (5) and use it to construct a series representation of the solution for fixed 𝑥 > 0 as 𝑡 → ∞, thereby providing a constructive proof of the radiation condition (4). We also obtain the large 𝑡 asymptotics along rays 𝑥∕𝑡 = constant of the wavemaker problem for the physically important linear KdV and linear BBM equations. At this point, we mention that the works 14-17 also consider the question of asymptotic periodicity for the linear KdV, the linear BBM, and a linear advection-dispersion equation; in particular, the authors employ either a Laplace transform in the temporal variable, a Green's function approach, or the Fokas method. Finally, we report experiments on a two-fluid system that agree quantitatively with solutions of the linear BBM wavemaker problem.

We are now ready to precisely state the above results. We begin with the former result, which is actually valid for generic, nonzero initial data, provided that these have sufficient smoothness and decay (e.g., of Schwartz type).

Theorem 1 (D-N map). Consider the evolution equation (5a) formulated on the half-line {𝑥 > 0} with a sufficiently smooth and decaying initial datum and the Dirichlet boundary datum (5c) that is infinitely smooth and asymptotically periodic in time, that is, 𝑔 0 ∈ 𝐶 ∞ (0, ∞) and 𝑔 0 (𝑡) ∼ ∑ 𝑛∈ℤ 𝑎 𝑛 𝑒 𝑖𝑛𝜔 0 𝑡 , 𝑡 → ∞, 𝜔 0 ≥ 0, 𝑎 𝑛 ∈ ℂ with the infinite sum being absolutely convergent. Assume further that the initial datum and the boundary datum 𝑔 0 (𝑡) satisfy the necessary compatibility conditions at the origin. Suppose that either {𝐴 -2 = 0, 𝐴 3 < 0} (local, third order) or {𝐴 -2 > 0, 𝐴 3 = 0} (nonlocal, second order) with 𝐴 2 ∕(𝜔 0 𝐴 -2 ) ∉ ℤ in the latter case. Assume that the solution 𝑢(𝑥, 𝑡) and its boundary values 𝑢 𝑥 (0, 𝑡), 𝑢 𝑥𝑥 (0, 𝑡) are asymptotically time-periodic in the sense of (21) with the additional assumptions (22) and (23). Then, in the case that the equation 𝐴 3 𝑘 3 + (𝐴 2 + 𝑛𝜔 0 𝐴 -2 )𝑘 2 -𝐴 1 𝑘 -𝐴 0 + 𝑛𝜔 0 = 0 has complex roots, the associated generalized D-N map satisfies to leading order the asymptotic formula

where 𝑘 0 (𝑛) is the unique solution of the equation 𝐴 3 𝑘 3 + (𝐴 2 + 𝑛𝜔 0 𝐴 -2 )𝑘 2 -𝐴 1 𝑘 -𝐴 0 + 𝑛𝜔 0 = 0 that lies along the positively oriented boundary 𝜕𝐷 + of the region 𝐷 + = {𝑘 ∈ ℂ ∶ Im(𝑘) > 0 and Re(

) < 0}.

Conjecture 1. Theorem 1 holds for the case that the equation 𝐴 3 𝑘 3 + (𝐴 2 + 𝑛𝜔 0 𝐴 -2 )𝑘 2 -𝐴 1 𝑘 -𝐴 0 + 𝑛𝜔 0 = 0 has three real roots.

Theorem 1 is established in Section 2, where the overall framework for the general, fiveparameter equation (5a) is first laid out and the analysis specific to the cases {𝐴 -2 = 0, 𝐴 3 < 0} (local, third order) and {𝐴 -2 > 0, 𝐴 3 = 0} (nonlocal, second order) addressed by Theorem 1 is then carried out in detail. Conjecture 1 is proven in the special cases of the linear KdV and linear BBM equations in Sections 2.2 and 2.3.

Corollary 1 (Radiation condition). Under the assumptions of Theorem 1 and Conjecture 1, 𝑘 0 (𝑛) is the unique branch of the dispersion relation (6) satisfying 𝑛𝜔 0 = -𝜔(𝑘 0 (𝑛)) if and only if the radiation condition (4) holds. Furthermore, the asymptotic solution of the wavemaker problem (5) takes the form 𝑢(𝑥, 𝑡) ∼ ∑ 𝑛∈ℤ 𝑎 𝑛 𝑒 𝑖(𝑘 0 (𝑛)𝑥+𝑛𝜔 0 𝑡) , 𝑡 → ∞, 𝑥 = (1).

Finally, there exist real numbers 𝜔 ± 𝑐𝑟 such that Im 𝑘 0 (1) = 0 if and only if 𝜔 - 𝑐𝑟 ≤ 𝜔 0 ≤ 𝜔 + 𝑐𝑟 .

A few remarks are in order.

1. In the nonlocal, second-order case {𝐴 -2 > 0, 𝐴 3 = 0}, the scenario 𝐴 2 ∕(𝜔 0 𝐴 -2 ) ∈ ℤ, which is not covered by Theorem 1, is discussed separately at the end of Section 2.3.

2.

The assumption 𝑔 0 ∈ 𝐶 ∞ (0, ∞) implies that the Fourier series coefficients 𝑎 𝑛 decay faster than any power of 𝑛 ∈ ℤ; thus, the infinite sum in ( 7) is absolutely convergent. To see this, note that 𝑘 0 (𝑛) grows like 𝑛 1 3 and 𝑛 1 2 in the third-order and second-order cases, respectively. 3. The uniqueness of the D-N map (7) follows from the requirement that 𝑘 0 (𝑛) is the unique root of the cubic 𝐴 3 𝑘 3 + (𝐴 2 + 𝑛𝜔 0 𝐴 -2 )𝑘 2 -𝐴 1 𝑘 -𝐴 0 + 𝑛𝜔 0 that lies along the boundary of the domain 𝐷 + . As explained in Section 2, this requirement is easy to justify for 𝑛𝜔 0 outside of the interval [𝜔 - 𝑐𝑟 , 𝜔 + 𝑐𝑟 ] (supercritical regime). However, within the interval (subcritical regime), things become more subtle and uniqueness follows by studying the large 𝑡 limit of the Fokas solution formulae for fixed 𝑥 > 0. This is done in Section 2 for the linear KdV and linear BBM equations, where, by taking advantage of the analyticity and exponential decay properties of the integrands involved in the solution formulae, we perform suitable contour deformations that allow us to quickly derive the leading-order, long-time asymptotics of the corresponding solutions for each fixed 𝑥 > 0. We emphasize that this approach is not valid when 𝑡 → ∞ along rays 𝑥∕𝑡 = constant, which is addressed in Theorems 2 and 3. Importantly, these calculations show that the choice of 𝑘 0 (𝑛) as the unique root of the aforementioned cubic that lies on 𝜕𝐷 + is directly related to the fact that 𝐷 + corresponds to the region in the upper half of the complex 𝑘-plane where the imaginary part of the dispersion relation 𝑘(𝜔) of Equation (5a) is positive, as in the radiation condition (4). Moreover, in the general case covered by Conjecture 1 that goes beyond the linear KdV and BBM equations, the same analysis can be carried out to reach the same criterion for 𝑘 0 (𝑛). This analysis addresses a relevant uniqueness issue that was identified in Ref. 18. 4. In the case of real roots of the aforementioned cubic, the root 𝑘 0 (𝑛) is the only one that corresponds to a nonnegative group velocity. Indeed, the group velocity associated with Equation (5a) is 𝑑 𝑑𝑘

Thus, for 𝑘 ∈ ℝ, the condition that 𝑘 ∈ 𝜕𝐷 + is equivalent to the group velocity being nonnegative at that solution as in the radiation condition (4). In the case of complex roots, the root 𝑘 0 (𝑛) is the unique root with Im 𝑘 0 (𝑛) > 0 that leads to a spatially decaying solution. 5. The calculations performed in Section 2 are under the assumption of smooth data that lead to a smooth solution with sufficient decay as 𝑥 → ∞. Nevertheless, depending on the precise choice of coefficients in (5a), these conditions can be relaxed significantly-see, for example, the works 19,20 where the Fokas method is combined with harmonic analysis techniques for the study of well-posedness of the KdV equation on the half-line and the finite interval with data in appropriate Sobolev spaces. Such an investigation of optimal regularity properties is not our task here. 6. A major assumption in the hypothesis of Theorem 1 is that, for zero initial data and asymptotically periodic Dirichlet data, the solutions of the models corresponding to (5a) with the choices {𝐴 -2 = 0, 𝐴 3 < 0} and {𝐴 -2 > 0, 𝐴 3 = 0}, as well as the various boundary values of these solutions, are asymptotically periodic in time. The explicit formulae obtained via the Fokas method for the wavemaker problems of the linear KdV and the linear BBM equations (which provide the unique solutions to these problems-see the Appendix as well as 21 ) allow us to verify the asymptotic time periodicity assumed in Theorem 1 in the case of these two physical models. More specifically, in Section 2, it is shown that by taking the asymptotic limit 𝑥 fixed and 𝑡 → ∞ in the solution formulae, we obtain the asymptotically time-periodic solution (8), and therefore, verify the periodicity assumption in Theorem 1. Then, in Sections 3 and 4, we apply the steepest descent method to the solution formulae and derive the leadingorder asymptotics of these solutions as 𝑡 → ∞ and 𝑥∕𝑡 = (1), that is, the behavior of these solutions along rays in the 𝑥-𝑡 plane. This type of limit provides detailed information about the propagation of waves away from the origin in the wavemaker problem. We also demonstrate that the asymptotic results from each region 𝑥 = (1) and 𝑥 = (𝑡) match as 𝑡 → ∞.

The precise asymptotic solution statements to the wavemaker problems are as follows.

Theorem 2 (Linear KdV wavemaker problem). Consider the linear KdV wavemaker problem with input frequency 𝜔 0 ,

There exists a critical frequency 𝜔 𝑐𝑟 ∶= 2∕(3 3). Assuming that the solution is sufficiently smooth, its leading-order asymptotic behavior as 𝑡 → ∞ and 𝜉 ∶= 𝑥∕𝑡 = (1) is as follows:

(I) A time-periodic sinusoidal function, when 0 < 𝜔 0 < 𝜔 𝑐𝑟 and 0 < 𝜉 < 𝑐 𝑔 ∶= 1 -3𝑘 2 0 < 1 with 𝑘 0 = 𝑘 3 (𝜔 0 ) given by ( 93), namely,

Moreover, this is an asymptotic approximation that is also valid for 𝜉𝑡 = (1), 𝑡 → ∞. (II) A sinusoidal function in 𝑡 with an algebraically decaying amplitude of (𝑡 -1∕2 ), when 0 < 𝜔 0 < 𝜔 𝑐𝑟 and 𝑐 𝑔 < 𝜉 < 1, or when 𝜔 0 > 𝜔 𝑐𝑟 and 0 < 𝜉 < 1 with the condition 𝑙 1 (𝜉, 𝜔 0 ) < 0 where 𝑙 1 is defined in (105), that is,

(III) A sinusoidal function in 𝑡 with an algebraically decaying amplitude of (𝑡 -1∕2 ) associated with an exponentially decaying sinusoidal function when 𝜔 0 > 𝜔 𝑐𝑟 and 0 < 𝜉 < 1 with the condition 𝑙 1 (𝜉, 𝜔 0 ) > 0 where 𝑙 1 is defined in (105), that is,

with 𝑘 0 = 𝑘 3 (𝜔 0 ) now given by (94). While the first term is beyond all orders in 𝑡 when 𝜉 = (1), this is an asymptotic approximation that is also valid for 𝜉𝑡 = (1), 𝑡 → ∞, where the first term is the dominant contribution. (IV) An exponential, decaying in time function when 𝜉 > 1, that is,

Theorem 3 (Linear BBM wavemaker problem). Consider the linear BBM wavemaker problem with input frequency 𝜔 0 ,

There exists a critical frequency 𝜔 𝑐𝑟 ∶= 1∕2. Assuming that the solution is sufficiently smooth, its leading-order asymptotic behavior as 𝑡 → ∞ and 𝜉 ∶= 𝑥∕𝑡 = (1) is described by: (I) A time-periodic sinusoidal function, when 0 < 𝜔 0 < 𝜔 𝑐𝑟 and 0 < 𝜉 < 𝑐 𝑔 ∶=

Moreover, this is an asymptotic approximation that is also valid for 𝜉𝑡 = (1), 𝑡 → ∞. (II) A sinusoidal function in 𝑡 with an algebraically decaying amplitude of (𝑡 -1∕2 ), when 0 < 𝜔 0 < 𝜔 𝑐𝑟 and 𝑐 𝑔 < 𝜉 < 1, or when 𝜔 0 > 𝜔 𝑐𝑟 and 0 < 𝜉 < 1 with the condition 𝑙 3 (𝜉, 𝜔 0 ) < 0 where 𝑙 3 is defined in (130), that is,

)) .

(

(III) A sinusoidal function in 𝑡 with an algebraically decaying amplitude of (𝑡 -1∕2 ) associated with an exponentially decaying sinusoidal function when 𝜔 0 > 𝜔 𝑐𝑟 and 0 < 𝜉 < 1 with the condition 𝑙 3 (𝜉, 𝜔 0 ) > 0 where 𝑙 3 is defined in (130), that is,

While the first term is beyond all orders in 𝑡 when 𝜉 = (1), this is an asymptotic approximation that is also valid for 𝜉𝑡 = (1), 𝑡 → ∞, where the first term is the dominant contribution. (IV) An exponentially decaying in time function when 𝜉 > 1, that is,

) )1∕2 (

14679590, 2024, 3, Downloaded from https://onlinelibrary.wiley.com/doi/10.1111/sapm.12665 by University Of Colorado Librari, Wiley Online Library on [28/04/2024]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License

The positive group velocity 𝑐 𝑔 in Theorems 2 and 3 is defined to be 𝑐 𝑔 ∶= 𝜔 ′ (𝑘 0 ), where 𝑘 0 is the unique root of the dispersion relation 𝜔(𝑘) that lies on 𝜕𝐷 + defined in Theorem 1. The ray 𝑥∕𝑡 = 𝑐 𝑔 for large time 𝑡 → ∞ represents the leading edge of the wave's transition from temporally periodic to temporally damped. Theorems 2 and 3 are established in Sections 3 and 4, respectively. Moreover, in Section 5, the analytical results of Section 4 are found to be in excellent agreement with experimental observations of a viscous fluid conduit modeled by the linear BBM equation in the long-wavelength and small-amplitude regime. A lighter, less viscous fluid is injected into a tall reservoir of fluid to establish a vertical fluid conduit rising buoyantly. The spontaneous introduction of a sinusoidal injection rate amounts to the wavemaker problem. As in previous experiments, 22 we observe a critical frequency below which waves propagate and above which waves spatially decay. New measurements include spatiotemporal wave patterns along rays 𝑥∕𝑡 = constant that agree quantitatively with our large 𝑡 predictions. These measurements include nondecaying waves, algebraically decaying waves, and exponentially decaying waves. The transition between nondecaying and algebraically decaying waves in the subcritical regime is observed to occur when 𝑥∕𝑡 coincides with the group velocity, in agreement with asymptotic predictions.

The universal nature of the Fokas method allows us to formulate a conjecture about the D-N map of the general nonlocal, third-order model (5a) in the case of asymptotically periodic Dirichlet data on the half-line. Although we do not prove this conjecture in the case of general coefficients 𝐴 -2 , 𝐴 0 , 𝐴 1 , 𝐴 2 , 𝐴 3 , we anticipate that one should be able to establish it for any specific choice of these coefficients by suitably adapting the arguments leading to the proofs of Theorems 1-3. (5) with an asymptotically periodic Dirichlet boundary condition 𝑢(0, 𝑡) = 𝑔 0 (𝑡) ∼ ∑ 𝑛∈ℤ 𝑎 𝑛 𝑒 𝑖𝑛𝜔 0 𝑡 , 𝑡 → ∞, where 𝜔 0 ≥ 0, 𝑎 𝑛 ∈ ℂ. Then, under the assumption that the solution 𝑢(𝑥, 𝑡) along with its boundary values 𝜕 𝑗 𝑥 𝑢(0, 𝑡), 𝑗 ∈ ℕ, is asymptotically time-periodic, the associated generalized D-N map satisfies the asymptotic formula of Theorem 1. Furthermore, the uniqueness of the D-N map is equivalent to the radiation condition (4).

Conjecture. Consider the general initial-boundary value problem

The paper is organized as follows. Theorem 1 is proven, under various assumptions, in Section 2. In Section 2.1, we construct the asymptotic D-N map (7) for the general model (5a), which is used to obtain the series representation (8) of the asymptotic solution. The asymptotic D-N map involves the root 𝑘 0 (𝑛) of the cubic polynomial 𝜔(𝑘 0 (𝑛)) = 𝑛𝜔 0 that is proved to be the unique root lying on 𝜕𝐷 + obtained in the restricted settings when 𝐴 -2 = 0 (Section 2.2) or 𝐴 3 = 0 (Section 2.3), provided that the frequency 𝑛𝜔 0 lies outside an interval. The proof for all real frequencies 𝑛𝜔 0 is obtained for the linear KdV (Section 2.2) and linear BBM (Section 2.3) equations where we also prove the asymptotic periodicity of the solution. Consequently, the radiation condition (4) is proven for the important physical cases of the linear KdV and linear BBM equations. Theorems 2 and 3 are proven in Sections 3 and 4, respectively. Viscous core-annular fluid experiments that quantitatively realize Theorems 1 and 3 are presented in Section 5. We conclude with a discussion of the implications of this work in Section 6.

THE D-N MAP VIA THE FOKAS METHOD

As noted in the introduction, beyond its main purpose as a universal approach for the solution of linear evolution equations in the initial-boundary value problems setting, the Fokas method provides an effective and efficient way of determining the (generalized) D-N map associated with such problems when the boundary conditions are asymptotically periodic in time. This aspect of the method has been the subject of the recent articles by Fokas and van der Weele 18 and Fokas, Pelloni, and Smith 23 for problems on the half-line and the finite interval, respectively.

Here, we consider the Dirichlet problem on the half-line for the general third-order model (5a), which includes the linear KdV equation discussed in Ref. 18 but additionally involves a mixed spatiotemporal derivative like the one present in the linear BBM equation. Through our analysis, we demonstrate that the Fokas method for the characterization of the D-N map in the case of asymptotically time-periodic boundary conditions remains effective for this general third-order model. We note that the analysis of this section is valid for generic initial data and not just for the zero initial condition stated in the wavemaker problem (5).

Importantly, we address a uniqueness issue that arises for temporal frequencies below a certain threshold (the subcritical regime) and manifests itself through the presence of solely real singularities among which only one specific choice gives rise to the correct D-N map, despite the fact that, at first sight, all of them qualify as removable. In particular, by deriving the leading-order long-time asymptotics directly from the explicit solution formulae obtained via the Fokas method, we rigorously justify the empirical "rule of thumb" used previously in the literature, according to which, in the subcritical regime, the singularities that must be made removable are the ones along the boundary of the region 𝐷 -defined by ( 27) below.

General framework

Define the half-line Fourier transform of 𝑢 with respect to 𝑥 by

We emphasize that, unlike the traditional whole-line Fourier transform, the half-line transform (17) is valid not only for 𝑘 ∈ ℝ but also for Im(𝑘) < 0, due to the fact that 𝑥 ≥ 0. Upon taking the transform (17) and integrating by parts, Equation (5a) yields

where we have introduced the following notation for the various (known or unknown) boundary values of 𝑢 at 𝑥 = 0:

Consider now the case of an asymptotically periodic Dirichlet boundary datum, that is, suppose that 𝑔 0 is given and satisfies

Furthermore, suppose that the remaining (unknown) boundary values 𝑔 1 and 𝑔 2 , as well as the solution 𝑢, are also asymptotically time-periodic with the same periodicity as 𝑔 0foot_1 , that is,

Under these assumptions, the characterization of the D-N map for Equation (5a) on the half-line amounts to expressing the unknown coefficients 𝑏 𝑛 and 𝑐 𝑛 in terms of the coefficients 𝑎 𝑛 , which are known since they correspond to the prescribed boundary datum. In view of (21), we write 𝑢(𝑥, 𝑡) = ∑ 𝑛∈ℤ 𝑈 𝑛 (𝑥)𝑒 𝑖𝑛𝜔 0 𝑡 + 𝐸(𝑥, 𝑡) where 𝐸(𝑥, 𝑡) = 𝑜(1) as 𝑡 → ∞, that is, lim 𝑡→∞ 𝐸(𝑥, 𝑡) = 0 for each 𝑥 ∈ (0, ∞). Taking the half-line Fourier transform of this equation yields û(𝑘, 𝑡) = ∫ ∞ 0 𝑒 -𝑖𝑘𝑥 ∑ 𝑛∈ℤ 𝑈 𝑛 (𝑥)𝑒 𝑖𝑛𝜔 0 𝑡 𝑑𝑥 + Ê(𝑘, 𝑡). Thus, assuming that ∑ 𝑛∈ℤ 𝑈 𝑛 (𝑥)𝑒 𝑖𝑛𝜔 0 𝑡 is uniformly convergent in 𝑥, we can interchange the infinite sum with the integral to obtain û(𝑘, 𝑡) = ∑ 𝑛∈ℤ Û𝑛 (𝑘)𝑒 𝑖𝑛𝜔 0 𝑡 + Ê(𝑘, 𝑡). Assuming further that 𝐸(𝑥, 𝑡) ∈ 𝐿 1 𝑥 (0, ∞) for each 𝑡 ≥ 𝑡 0 > 0, the dominated convergence theorem implies lim 𝑡→∞ Ê(𝑘, 𝑡) = lim 𝑡→∞ ∫ ∞ 0 𝑒 -𝑖𝑘𝑥 𝐸(𝑥, 𝑡)𝑑𝑥 = ∫ ∞ 0 𝑒 -𝑖𝑘𝑥 lim 𝑡→∞ 𝐸(𝑥, 𝑡)𝑑𝑥 = 0, from which we deduce that the half-line Fourier transform û(𝑘, 𝑡) of 𝑢(𝑥, 𝑡) is also asymptotically time periodic with

Then, under the additional assumption that the time derivatives û𝑡 (𝑘, 𝑡), 𝑔 ′ 0 (𝑡) and 𝑔 ′ 1 (𝑡) are also asymptotically time-periodic with corresponding expressions given by the term-by-term differentiation in 𝑡 of the series for û(𝑘, 𝑡), 𝑔 0 (𝑡) and 𝑔 1 (𝑡), namely,

Equation (18) yields

where we have introduced the notation

Consequently, 𝜔(𝑘) = Ω(𝑘)∕(1 + 𝐴 -2 𝑘 2 ) is the dispersion relation (6). We deduce from ( 24)

Zeros of the denominator correspond to roots of 𝜔(𝑘) = -𝑛𝜔 0 . Note that the negative sign is due to the fact that the dispersion relation (6) was defined for waves in the form (1), that is, opposite that of the frequency defined in (22). As mentioned above, û(𝑘, 𝑡) is analytic for Im(𝑘) < 0 and well defined in the 𝐿 2 -sense for Im(𝑘) ≤ 0 (e.g., see Theorem 7.2.4 in [24]). Hence, it would be reasonable to expect that the expression (26) for Û𝑛 (𝑘) is well defined for Im(𝑘) ≤ 0. In turn, this would imply that any zeros of the quantity Ω + 𝑛𝜔 0 (1 + 𝐴 -2 𝑘 2 ) that arise in {Im(𝑘) ≤ 0} are actually removable singularities of (26). However, as we shall see later in the special cases of the linear KdV and linear BBM equations (Sections 2.2 and 2.3), this is true only for the zeros of Ω + 𝑛𝜔 0 (1 + 𝐴 -2 𝑘 2 ) that lie on the boundary 𝜕𝐷 -of the domain 𝐷 -, where

Observe that, since 𝑛𝜔 0 ∈ ℝ, the zeros of Ω + 𝑛𝜔 0 (1 + 𝐴 -2 𝑘 2 ) always satisfy the equation Re( 𝑖Ω 1+𝐴 -2 𝑘 2 ) = 0 and hence lie on 𝜕𝐷 + ∪ 𝜕𝐷 -, as ℝ ⊂ (𝜕𝐷 + ∪ 𝜕𝐷 -). To see this, note that Re

and so, Re( 𝑖Ω 1+𝐴 -2 𝑘 2 ) vanishes for 𝑘 ∈ ℝ and is odd with respect to Im(𝑘). Thus, ℝ ⊂ (𝜕𝐷 + ∪ 𝜕𝐷 -) since (𝜕𝐷 + ∪ 𝜕𝐷 -) ∩ ℝ ≠ ∅ and the conjugate of any point in 𝐷 + belongs to {Im(𝑘) < 0} ⧵ 𝐷 -and the conjugate of any point in 𝐷 -belongs to {Im(𝑘) > 0} ⧵ 𝐷 + .

Remark. The departure from the original intuition that the expression (26) for Û𝑛 (𝑘) should make sense for all Im(𝑘) ≤ 0 is due to the fact that, although related to each other, the quantities Û𝑛 (𝑘) and û(𝑘, 𝑡) are not quite the same in terms of analyticity; it seems that the manipulations between infinite series and integrals used for deriving (26) do not fully preserve analyticity in 𝑘 from û(𝑘, 𝑡) to Û𝑛 (𝑘).

Explicitly, Ω + 𝑛𝜔 0 (1 + 𝐴 -2 𝑘 2 ) = 0 corresponds to the cubic equation

The requirement that the zeros of (29) along 𝜕𝐷 -are removable singularities of the right-hand side of (26) implies that these zeros must also satisfy the equation

Then, for 𝐴 3 ≠ 0, assuming that precisely two roots 𝑘 - 1 = 𝑘 - 1 (𝑛) and 𝑘 - 2 = 𝑘 - 2 (𝑛) of the cubic equation (29) lie on 𝜕𝐷 -while the third root 𝑘 0 = 𝑘 0 (𝑛) lies on 𝜕𝐷 + , we have the following system for the coefficients 𝑏 𝑛 and 𝑐 𝑛 :

Solving this system while noting that, by Vieta's formulae, 𝑘 - 1 + 𝑘 - 2 + 𝑘 0 = -(𝐴 2 + 𝑛𝜔 0 𝐴 -2 )∕𝐴 3 and 𝑘 - 1 𝑘 - 2 + 𝑘 - 2 𝑘 0 + 𝑘 0 𝑘 - 1 = -𝐴 1 ∕𝐴 3 , we obtain the D-N map for the third-order model (5a) as

Similarly, for 𝐴 3 = 0, which corresponds to the special case where (5a) reduces to a second-order model, assuming that precisely, one root 𝑘 - 1 = 𝑘 - 1 (𝑛) of the now quadratic equation ( 29) lies on 𝜕𝐷 -while the second root 𝑘 0 = 𝑘 0 (𝑛) lies on 𝜕𝐷 + , we obtain the first of the relations (32), which amounts to the D-N map for this second-order reduction.

Remark (Series representation of the leading-order asymptotics via the D-N map). Under the assumption of sufficient smoothness, the leading-order long-time asymptotics of the solution of the general model (5a) on the half-line can be directly constructed via the D-N map (32) by means of the Maclaurin series 𝑢(𝑥, 𝑡) = ∑ ∞ 𝑗=0 𝜕 𝑗 𝑥 𝑢(0,𝑡) 𝑗! 𝑥 𝑗 . Indeed, for 𝐴 3 ≠ 0, (32) yields 𝑢 𝑥 (0, 𝑡) ∼ 𝑖 ∑ 𝑛∈ℤ 𝑘 0 (𝑛)𝑎 𝑛 𝑒 𝑖𝑛𝜔 0 𝑡 , 𝑢 𝑥𝑥 (0, 𝑡) ∼ -∑ 𝑛∈ℤ [𝑘 0 (𝑛)] 2 𝑎 𝑛 𝑒 𝑖𝑛𝜔 0 𝑡 , which can be combined with (5a) evaluated at the boundary 𝑥 = 0 and the fact that 𝑘 0 (𝑛) satisfies (29) to obtain 𝑢 𝑥𝑥𝑥 (0, 𝑡) ∼ -𝑖 ∑ 𝑛∈ℤ [𝑘 0 (𝑛)] 3 𝑎 𝑛 𝑒 𝑖𝑛𝜔 0 𝑡 .

Moreover, along the same lines and using induction while assuming that the 𝑥-derivatives of the solution are asymptotically time-periodic at 𝑥 = 0, we have the general formula

Similarly, for 𝐴 3 = 0, we can use the first of the relations (32) to establish (33) via induction. Therefore, in view of the Maclaurin series of 𝑢 with respect to 𝑥, and assuming that we can interchange the order of summation, we deduce

which shows that, under the asymptotic periodicity assumptions ( 21)-( 23) and also the aforementioned assumption of sufficient smoothness needed for the Maclaurin series, for asymptotically periodic boundary data of the form (20), the leading-order asymptotics of the solution of (5a) is a superposition of waves that either propagate (if 𝑘 0 (𝑛) ∈ ℝ) or spatially decay (if Im 𝑘 0 (𝑛) > 0).

To rigorously justify the D-N map (32), one must study the location of the roots of the cubic equation ( 29) with respect to the curves 𝜕𝐷 -and 𝜕𝐷 + . For general parameters 𝐴 𝑗 , this task involves the study of the quartic polynomial in Re(𝑘) and Im(𝑘) that appears in the expression (28) for Re( 𝑖Ω 1+𝐴 -2 𝑘 2 ). Hence, below we focus on two particular sets of parameters, namely, {𝐴 -2 = 0, 𝐴 3 < 0} and {𝐴 3 = 0, 𝐴 -2 > 0}, which still correspond to generalizations of the linear KdV and the linear BBM equations, respectively, that are more tractable due to the fact that the aforementioned quartic reduces to a quadratic.

2.2

A general local third-order model: 𝑨 -𝟐 = 𝟎 and 𝑨 𝟑 < 𝟎

For 𝐴 -2 = 0 and 𝐴 3 < 0, Equation (5a) no longer involves a mixed spatiotemporal derivative as it takes the form

and the denominator in (26) reduces to 𝑖(Ω + 𝑛𝜔 0 ) with Equation ( 29) simplifying to

Before discussing the roots of (36), it is useful to give the precise characterization of the domains 𝐷 -and 𝐷 + . Starting from (25), we compute Re

so that, according to the definition (27),

Thus, for real solutions of ( 36) that lie along the boundary 𝜕𝐷 -, we must have

Similarly, the region 𝐷 + is described by

and therefore, any real solution of ( 36) that lies along the boundary 𝜕𝐷 + must satisfy

We now turn our attention to the solution of Equation (36). The roots of this cubic are given by

where 𝜈 𝑗 are the solutions of the depressed cubic

The discriminant of this depressed cubic is Δ = -4𝑝 3 -27𝑞 2 . According to the classical result by del Ferro and Tartaglia, also known as Cardano's formula, there are three cases:

Case 1: Δ < 0. This is equivalent to

Note that this inequality is satisfied for all 𝑛 ∈ ℤ if 𝑝 > 0. On the other hand, if 𝑝 ≤ 0, then taking square roots and defining

we must have either 𝑛𝜔 0 > 𝜔 + 𝑐𝑟 or 𝑛𝜔 0 < 𝜔 - 𝑐𝑟 (note that 𝜔 - 𝑐𝑟 ≤ 𝜔 ≤ 𝜔 + 𝑐𝑟 since 𝐴 3 < 0 and 𝑝 ≤ 0). That is, Δ < 0 if either 𝑝 > 0, or 𝑝 ≤ 0 and 𝑛𝜔 0 > 𝜔 + 𝑐𝑟 or 𝑛𝜔 0 < 𝜔 - 𝑐𝑟 . Then, the depressed cubic ( 43) has the real solution 𝜈 1 = (-

𝑞 2 + √ 𝑞 2 4 + 𝑝 3 27 ) 1 3 + (-𝑞 2 -√ 𝑞 2 4 + 𝑝 3 27 ) 1 3 and the complex conjugate solutions 𝜈 2 = (-𝑞 2 + √ 𝑞 2 4 + 𝑝 3 27 ) 1 3 𝑒 𝑖 2𝜋 3 + (-𝑞 2 -√ 𝑞 2 4 + 𝑝 3 27 ) 1 3 𝑒 -𝑖 2𝜋 3 and 𝜈 3 = (-𝑞 2 + √ 𝑞 2 4 + 𝑝 3 27 ) 1 3 𝑒 -𝑖 2𝜋 3 + (-𝑞 2 -√ 𝑞 2 4 + 𝑝 3 27 ) 1 3 𝑒 𝑖 2𝜋 3

. In turn, the cubic (36) has one real root and two complex conjugate roots, namely,

Recalling that Re(𝑖Ω) = 0 at the roots 𝑘 𝑗 , it follows that precisely one of the complex conjugate roots 𝑘 2 and 𝑘 3 lies on 𝜕𝐷 -. Moreover, we claim that 𝑘 1 ∈ 𝜕𝐷 -. Indeed, we begin by observing that

. This is clearly true if 𝑝 ≥ 0. In addition, if 𝑝 < 0, then noting that for 𝑞 ≥ 0, we have ) 1∕3 for

as desired. In turn, using (42) and substituting for 𝑝 while recalling that 𝐴 3 < 0, we obtain precisely the condition (39) for 𝑘 1 and thus prove that 𝑘 1 ∈ 𝜕𝐷 -as claimed.

Case 2: Δ > 0. This amounts to 𝑞 2 4 + 𝑝 3

27

< 0, which yields 𝜔 - 𝑐𝑟 < 𝑛𝜔 0 < 𝜔 + 𝑐𝑟 with 𝜎 ± given by (44) (note that 𝜔 + 𝑐𝑟 > 𝜔 - 𝑐𝑟 because 𝐴 3 < 0 and for Δ > 0, we must have 𝑝 < 0).

Then, the depressed cubic ( 43

3 , which correspond to the following distinct real roots for the cubic (36):

) 14679590, 2024, 3, Downloaded from https://onlinelibrary.wiley.com/doi/10.1111/sapm.12665 by University Of Colorado Librari, Wiley Online Library on [28/04/2024]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License

We will now show that 𝑘 1 , 𝑘 2 ∈ 𝜕𝐷 -while 𝑘 3 ∈ 𝜕𝐷 + . Writing -

𝑞 2 + 𝑖 √ -𝑞 2 4 -𝑝 3 27 = 𝑟𝑒 𝑖𝜃 with 𝑟 = √ -𝑝 3 27 and tan 𝜃 = √ -𝑞 2 4 -𝑝 3 27 ∕(-𝑞 2

), we have

Since 𝜔 - 𝑐𝑟 < 𝜔 < 𝜔 + 𝑐𝑟 for 𝑝 < 0 and, in addition, 𝜔 = -𝐴 3 𝑞 + 𝑛𝜔 0 and 𝐴 3 < 0, we have two subcases:

]. Hence,

, and -

. Thus, since by (39), the portion of 𝜕𝐷 - along ℝ is given by (-∞, -

, ∞), we deduce that 𝑘 1 , 𝑘 2 ∈ 𝜕𝐷 - and 𝑘 3 ∈ 𝜕𝐷 + . In particular, inequalities ( 39) and ( 41) hold strictly at 𝑘 1 , 𝑘 2 and at 𝑘 3 , respectively.

(ii) If 𝜔 - 𝑐𝑟 < 𝑛𝜔 0 ≤ 𝜔, then 𝑞 ≥ 0 and so 𝜃 ∈ [ 𝜋 2

, 𝜋). Hence,

, and 0 ≤ cos

, which imply

, and -

. Thus, as in the previous subcase, we have 𝑘 1 , 𝑘 2 ∈ 𝜕𝐷 -and 𝑘 3 ∈ 𝜕𝐷 + with inequalities ( 39) and ( 41) holding strictly at 𝑘 1 , 𝑘 2 and at 𝑘 3 , respectively.

Case 3: Δ = 0. This is equivalent to 4𝑝 3 + 27𝑞 2 = 0, thus either 𝑛𝜔 0 = 𝜔 + 𝑐𝑟 or 𝑛𝜔 0 = 𝜔 - 𝑐𝑟 . In the former case, the cubic (36) admits the real solutions

When 𝑛𝜔 0 = 𝜔 - 𝑐𝑟 , the cubic (36) admits the real solutions

Hence, in view of the fact that 𝜕𝐷 -∩ ℝ = (-∞, -

, ∞), we deduce that if 𝑝 < 0, then 𝑘 1 ∈ 𝜕𝐷 -with strict inequality in (39) and 𝑘 2 = 𝑘 3 ∈ 𝜕𝐷 -∩ 𝜕𝐷 + with equality in (41), while if 𝑝 = 0, then 𝑘 1 = 𝑘 2 = 𝑘 3 ∈ 𝜕𝐷 -∩ 𝜕𝐷 + with equality in (39) and (41).

Overall, our analysis leads to the following conclusions:

• If 3𝐴 1 𝐴 3 + 𝐴 2 2 < 0, then the solutions of ( 36) are given by ( 45) for all 𝑛 ∈ ℤ. • If 3𝐴 1 𝐴 3 + 𝐴 2 2 > 0, then the solutions of ( 36) are given by ( 45) for 𝑛𝜔 0 > 𝜔 + 𝑐𝑟 or 𝑛𝜔 0 < 𝜔 - 𝑐𝑟 , by (46) for 𝜔 - 𝑐𝑟 < 𝑛𝜔 0 < 𝜔 + 𝑐𝑟 , and by (48) for 𝑛𝜔 0 = 𝜔 + 𝑐𝑟 or (49) for 𝑛𝜔 0 = 𝜔 - 𝑐𝑟 , where 𝜔 ± 𝑐𝑟 are defined by (44).

• If 3𝐴 1 𝐴 3 + 𝐴 2 2 = 0, then 𝜔 + 𝑐𝑟 = 𝜔 - 𝑐𝑟 = 𝜔 defined by (44) and the solutions of ( 36) are given by ( 45) for 𝑛𝜔 0 ≠ 𝜔 and by (48) for 𝑛𝜔 0 = 𝜔.

Importantly, in all three of the above cases, the contour 𝜕𝐷 -contains precisely two solutions of (36), whereas the contour 𝜕𝐷 + contains only one, which we label 𝑘 0 . That is, we have shown that for 𝐴 -2 = 0 and 𝐴 3 < 0 is satisfied, and hence, the D-N map for the general third-order model ( 35) is indeed given by (32).

This proves Theorem 1 when 𝐴 -2 = 0.

Remark (Connection to the radiation condition ( 4)). The group velocity associated with the thirdorder equation ( 35) is 𝑐 𝑔 = Ω ′ (𝑘) = 3𝐴 3 𝑘 2 + 2𝐴 2 𝑘 -𝐴 1 . Thus, for a real solution of the cubic (36), the condition that this solution lies on 𝜕𝐷 + (Equation ( 41)) is equivalent to the group velocity being nonnegative at that solution. In particular, recall that the scenario 𝑘 0 ∈ 𝜕𝐷 + ∩ ℝ is possible if either 3𝐴 1 𝐴 3 + 𝐴 2 2 > 0 and 𝜔 - 𝑐𝑟 < 𝑛𝜔 0 < 𝜔 + 𝑐𝑟 , or 3𝐴 1 𝐴 3 + 𝐴 2 2 > 0 and either 𝑛𝜔 0 = 𝜔 + 𝑐𝑟 or 𝑛𝜔 0 = 𝜔 - 𝑐𝑟 , or 3𝐴 1 𝐴 3 + 𝐴 2 2 = 0 and 𝑛𝜔 0 = 𝜔. Hence, we have a strictly positive group velocity in the first scenario (see "Case 2" above) and zero group velocity in the other three scenarios (see "Case 3" above), which deals with the boundary case 𝑐 𝑔 = 0. Furthermore, if 𝑘 0 ∈ 𝜕𝐷 + ⧵ ℝ, then it follows from the definition (40) of 𝐷 + that Im(𝑘) > 0.

This proves Corollary 1 when 𝐴 -2 = 0.

The analysis of the roots of the cubic Ω + 𝑛𝜔 0 (i.e., of Equation ( 36)) for the general third-order model ( 35) provided a complete justification of the assumption used in the derivation of the D-N map (32) that precisely two of those three roots lie on the contour 𝜕𝐷 -. However, there is one more assumption that was needed for deducing (32) and is yet to be justified, namely, that the right-hand side of (26) should have removable singularities only at those roots of Ω + 𝑛𝜔 0 that occur along 𝜕𝐷 -. Above a certain threshold of values for |𝑛𝜔 0 |, it is easy to see why this latter assumption must be imposed; below that threshold, however, the situation becomes more delicate.

More precisely, for those values of 𝑛𝜔 0 for which the roots of Ω + 𝑛𝜔 0 involve a complex conjugate pair (i.e., in Case 1: Δ < 0 above, 𝑛𝜔 0 ∈ [𝜔 - cr , 𝜔 + cr ] 𝑐 ), one of the complex roots lies in the upper half of the 𝑘-plane where û(𝑘, 𝑡) (and hence Û𝑛 (𝑘)) is not well defined (recall discussion below (17)). Thus, this root cannot be a removable singularity of (26). On the other hand, the remaining two roots do lie in the domain of definition {Im(𝑘) ≤ 0} of û(𝑘, 𝑡) and so, it is reasonable to expect that they must be removable singularities of (26), thereby leading to the system (31) and hence to the D-N map (32). The situation is analogous in the degenerate case of one distinct and two repeated real roots (i.e., in Case 3: Δ = 0 above, 𝑛𝜔 0 ∈ {𝜔 - 𝑐𝑟 , 𝜔 + 𝑐𝑟 }) because, counting without multiplicities, there are exactly two roots that belong to {Im(𝑘) ≤ 0}, and hence, it is reasonable to expect that must be removable singularities of (26).

Nevertheless, in the case of three distinct real roots (i.e., Case 2: Δ > 0 above, 𝑛𝜔 0 ∈ (𝜔 - 𝑐𝑟 , 𝜔 + 𝑐𝑟 )), we encounter the following issue: while all three roots, namely, the two lying on 𝜕𝐷 -and the one lying on 𝜕𝐷 + , belong to {Im(𝑘) ≤ 0}, and hence, according to our previous reasoning, should be removable singularities of (26), requiring that all three of them are removable results via (30) in

for the linear KdV equation as given by ( 51)-( 52). Inside these domains Re(𝑖Ω) < 0, while outside of them Re(𝑖Ω) > 0. In the case of the wavemaker data (53

3), then the three roots of Ω + 𝑛𝜔 0 are real and distinct with 𝑘 - 1 ∈ 𝜕𝐷 - 1 , 𝑘 - 2 ∈ 𝜕𝐷 - 2 , 𝑘 0 ∈ 𝜕𝐷 + as shown, and the contour of integration in the Fokas method solution (54) is deformed from the boundary 𝜕𝐷 + of 𝐷 + to the contour 𝜕 D+ that bypasses the real root 𝑘 0 ∈ 𝜕𝐷 + from below. The resulting formula (55) is used for deriving the D-N map of the wavemaker problem by computing the long-time asymptotics via Cauchy's residue theorem, exponential decay and, when necessary, the deformation of 𝜕 D+ to the half-lines  1 ,  2 given by (57), along which Re(𝑖Ω) > 0

where Re(𝑖Ω) = 0).

three equations for the two unknowns 𝑏 𝑛 and 𝑐 𝑛 , that is, an overdetermined system. Therefore, the question is: which two of the three distinct real roots of Ω + 𝑛𝜔 0 are actually removable singularities of ( 30)? Note that each of the three possible choices leads to a different D-N map since, via Vieta's formulae, the root that is not a removable singularity is the one that eventually appears in the D-N relations (32). From this point of view, the challenge we are facing is directly related to the uniqueness of the D-N map.

Our conjecture is that, even in the case of three distinct real roots, our earlier assumption that led to (32) is correct, namely, the roots of Ω + 𝑛𝜔 0 that must be removable singularities of (26) are precisely the two roots that lie along 𝜕𝐷 -(and not the one along 𝜕𝐷 + -see remark below (28) for a possible explanation of this fact in connection with the domains of definition of û(𝑘, 𝑡) and Û𝑛 (𝑘)). In what follows, we prove this conjecture for the special case of the linear KdV equation by exploiting the explicit solution formula produced by the Fokas method.

2.2.1

The linear KdV equation on the half-line

The linear KdV equation

is included in the four-parameter family of third-order models (35) (it corresponds to 𝐴 0 = 𝐴 2 = 0 and 𝐴 1 = 𝐴 3 = -1), and thus, the general framework developed earlier can readily be employed.

In particular, the dispersion relation is Ω(𝑘) = 𝑘 -𝑘 3 , and the regions 𝐷 + and 𝐷 -are given by

where (see Figure 1)

14679590, 2024, 3, Downloaded from https://onlinelibrary.wiley.com/doi/10.1111/sapm.12665 by University Of Colorado Librari, Wiley Online Library on [28/04/2024]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License

and the discriminant of the cubic Ω + 𝑛𝜔 0 is Δ = 4 -27(𝑛𝜔 0 ) 2 , which defines symmetric critical interval endpoints 𝜔 ± 𝑐𝑟 = ±𝜔 𝑐𝑟 = ±2∕(3 √ 3). Moreover, in all three scenarios Δ < 0, Δ > 0 and Δ = 0, each of the contours 𝜕𝐷 + , 𝜕𝐷 - 1 , and 𝜕𝐷 - 2 contains precisely one of the three roots 𝑘 0 , 𝑘 - 1 , and 𝑘 - 2 of the cubic Ω + 𝑛𝜔 0 , respectively. As explained above, the challenge is how to justify that, among the three distinct real roots that arise in the case Δ > 0, that is,

3), the two roots 𝑘 - 1 ∈ 𝜕𝐷 - 1 and 𝑘 - 2 ∈ 𝜕𝐷 - 2 that lie on 𝜕𝐷 -= 𝜕𝐷 - 1 ∪ 𝜕𝐷 - 2 should be removable singularities of ( 26) while the remaining root 𝑘 0 ∈ 𝜕𝐷 + should not, despite the fact that it belongs to the domain of definition {Im(𝑘) ≤ 0} of û(𝑘, 𝑡). Therefore, we consider Equation (50) in the context of the wavemaker problem associated with the data

and restrict our attention to the scenario

3) that, in the case of the data (53), corresponds to Ω + 𝑛𝜔 0 having three distinct real roots since 𝑢(0, 𝑡) = -∑ 𝑛=±1 sgn(𝑛)

2𝑖

𝑒 𝑖𝑛𝜔 0 𝑡 and, regardless of which two roots end up satisfying (30), the resulting D-N relations will imply that the Fourier series coefficients 𝑏 𝑛 , 𝑐 𝑛 of the boundary values 𝑢 𝑥 (0, 𝑡), 𝑢 𝑥𝑥 (0, 𝑡) vanish whenever the coefficient 𝑎 𝑛 of 𝑢(0, 𝑡) does.

In what follows, we will show that choosing 𝑘 - 1 and 𝑘 - 2 as the two removable singularities of (26) does indeed produce the correct D-N map, namely, the first relation given in (32). We will do so by independently deriving that D-N relation directly from the solution formula obtained for the linear KdV wavemaker problem (50)-(53) via the Fokas method, which is given by (see Ref. 10, Example 1.12)

Note that the three roots 𝑘 0 , 𝑘 - 1 , 𝑘 - 2 of Ω + 𝑛𝜔 0 are removable singularities of the above integrand, which is therefore analytic in 𝑘. Hence, in the particular case that 𝑘 0 , 𝑘 - 1 , 𝑘 - 2 are real and distinct, we invoke Cauchy's theorem to deform the contour of integration 𝜕𝐷 + in (54) to the modified contour 𝜕 D+ that bypasses the removable singularity 𝑘 0 ∈ 𝜕𝐷 + from below as shown in Figure 1, and thus, obtain

Since 𝑘 0 no longer lies on the contour of integration, we can split the fraction in the integrand to write

Furthermore, for 𝑘 ∈ {Im(𝑘) ≥ 0 ∶ 𝑘 ∉ 𝐷 + }, we have |𝑒 -𝑖Ω𝑡 | ≤ 1 so if, in addition, 𝑘 is large in modulus and stays away from 𝑘 - 1 and 𝑘 - 2 , then

. This decay bound combined with analyticity allows us to employ Cauchy's theorem and Jordan's lemma in the second integral of (56) in order to deform 𝜕 D+ to the union  1 ∪  2 of the half-lines (with orientation as shown in Figure 1)

and hence express (56) in the form

We emphasize the significance of the bound |

, which provides uniform decay in 𝑘 along the circular arcs connecting 𝜕 D+ with  1 and  2 and thus allows us to argue that the relevant integrals vanish by invoking Jordan's lemma (or, alternatively, by a direct estimation of those integrals along the lines of the integral 𝐼 𝑅 below).

We will now use the modified solution formula (58) to compute the leading-order long-time asymptotics of 𝑢 𝑥 (𝑥, 𝑡) for any finite 𝑥 > 0. Differentiating with respect to 𝑥 yields

Let 𝜕 D+ 𝑅 denote the finite portion of the contour 𝜕 D+ that lies below the circular arc

, 𝑅 > 1}. With this notation, the first integral on the right-hand side of (59) can be regarded as the limit

Therefore, by Cauchy's residue theorem,

where 𝐼 𝑅 (𝑥, 𝑡) ∶= ∫ 𝐶 𝑅 𝑒 𝑖𝑘𝑥 𝑘Ω ′ Ω+𝑛𝜔 0 𝑑𝑘. We claim that lim 𝑅→∞ 𝐼 𝑅 (𝑥, 𝑡) = 0 for any 𝑥 > 0, 𝑡 > 0, and hence, calculating the residue in (60),

Indeed, for 𝑘 ∈ 𝐶 𝑅 , we have

14679590, 2024, 3, Downloaded from https://onlinelibrary.wiley.com/doi/10.1111/sapm.12665 by University Of Colorado Librari, Wiley Online Library on [28/04/2024]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License and, using the symmetry of sin 𝜃 about 𝜋 2 as well as the well-known inequality sin 𝜃 ≥ 2𝜃 𝜋 , 0 ≤ 𝜃 ≤ 𝜋 2

, we find

) , which implies the claim for any 𝑥 > 0, 𝑡 > 0. In turn, combining (59) with (61), we obtain

Remark. We could have chosen the contour 𝜕 D+ to bypass 𝑘 0 from above. This would have made the first integral of (59) equal to zero and would have moved the residue contribution over to the second integral.

Next, we will prove that, for any 𝑥 > 0, the integral on the right-hand side of (62) vanishes as 𝑡 → ∞. We emphasize that the deformation from 𝜕 D+ to  1 ∪  2 in (58) was possible because Re(𝑖Ω) ≥ 0 in the upper half-plane and outside 𝐷 + . In fact, Re(𝑖Ω) ≥ 0 also in the lower halfplane and outside 𝐷 - 1 and 𝐷 - 2 and, more precisely, Re(𝑖Ω) > 0 for 𝑘 ∈  1 ∪  2 ⧵ {±1∕ √ 3} so that the exponential 𝑒 -𝑖Ω𝑡 is strictly decaying there. This fact will allow us to show that the integral on the right-hand side of (62) decays like 𝑂(1∕ √ 𝑡) as 𝑡 → ∞. We remark that our argument will not work without the presence of the exponential 𝑒 -𝑖Ω𝑡 , which decays as |𝑘| → ∞ along  1 ∪  2 . Thus, the deformation from 𝜕 D+ to  1 ∪  2 is only possible and helpful for the second integral in (56) and not for the first one, which does not involve 𝑒 -𝑖Ω𝑡 .

Since the estimation along the half-lines  1 and  2 is entirely analogous, we provide the details only for  2 . Then, 𝑘 = -𝑖 + 𝑟𝑒

𝑖 𝜋 3 = 𝑟 2 + 𝑖( √ 3𝑟 2 -1) and so Re(𝑖Ω) = 3 2 (𝑟 -2 √ 3

) 2 . Importantly, observe that due to the definition of  2 the expression for Re(𝑖Ω) is quadratic (as opposed to cubic) in 𝑟. Also, similar to the proof of (62), for 𝑟 sufficiently large, say 𝑟 ≥ 10, we have |Ω + 𝑛𝜔 0 | ≳ 𝑟 3 . Furthermore, |Ω + 𝑛𝜔 0 | is strictly positive by the definition of  2 and the fact that

3), that is, |Ω + 𝑛𝜔 0 | ≥ 𝐶(𝑛𝜔 0 ) > 0 for some positive constant depending on 𝑛𝜔 0 . Hence,

(63) A similar result holds for the integral along  1 and hence for (62).

Overall, we have shown that, in the scenario

3) of three distinct real roots of Ω + 𝑛𝜔 0 , for any 𝑥, 𝑡 > 0, the spatial derivative of the solution to the linear KdV wavemaker problem

In other words, for each 𝑥 ∈ [0, 𝐿], some 𝐿 > 0,

As we will see in Section 3, the bound (64) is not tight and can be refined to 𝑂(𝑡 -3∕2 ). Taking the limit, we find lim 𝑥→0 + 𝑢 𝑥 (𝑥, 𝑡) ∼ ∑ 𝑛=±1 𝑖𝑘 0 𝑎 𝑛 𝑒 𝑖𝑛𝜔 0 𝑡 as 𝑡 → ∞, in agreement with the first of the D-N relations (32). A similar analysis of the solution representation (58) can be used to prove that the second spatial derivative 𝑢 𝑥𝑥 satisfies the second of the D-N relations (32).

We have shown that the D-N relations are obtained from the solution formula (58). In fact, essentially the same analysis implies the asymptotic form of the solution

By linearity, this result extends to the more general boundary data 𝑔 0 (𝑡) ∼ ∑ 𝑛∈ℤ 𝑎 𝑛 𝑒 𝑖𝑛𝜔 0 𝑡 . In particular, we have also proven that the solution is asymptotically time periodic for each fixed 𝑥 ≥ 0. Thus, for the linear KdV wavemaker problem, we have proven Theorem 1 and the radiation condition (4) (Corollary 1) for any frequency 𝜔 0 ∈ ℝ without the assumption of asymptotic time periodicity of the solution.

The asymptotic solution (66) is only valid for fixed 𝑥 as 𝑡 → ∞. In Section 3, we extend the analysis of the linear KdV wavemaker problem to obtain an asymptotic expansion of the solution for large 𝑡 that is uniform in 𝑥 ≥ 0.

Remark. The above arguments and, in particular, estimate (63), crucially rely on the fact that Re(𝑖Ω) ≥ 0 along  1 and  2 . Attempting to deform the contour 𝜕𝐷 + differently in order to also enclose one or both of the roots 𝑘 1 , 𝑘 2 would force us to include contours in regions where Re(𝑖Ω) < 0 (see Figure 1) and hence would result in exponentially growing terms, thus an estimate like (63) would no longer hold. This is why the only root contributing to the D-N map (through the residue calculation (60)) is the pole along 𝜕𝐷 + , namely, 𝑘 0 .

2.3

A general nonlocal second-order model: 𝑨 -𝟐 > 𝟎 and 𝑨 𝟑 = 𝟎

For 𝐴 -2 > 0 and 𝐴 3 = 0, the third-order equation (5a) becomes of second order but still involves the mixed spatiotemporal derivative term, namely,

Moreover, the equation Ω + 𝑛𝜔 0 (1 + 𝐴 -2 𝑘 2 ) = 0, which, in general, corresponds to the cubic (29), now reduces to the quadratic

In addition, according to definition (27),

Thus, for real solutions of (68) that lie along the boundary 𝜕𝐷 -, we must have

Likewise, the region 𝐷 + is given by

and therefore, any real solution of (68) that lies along the boundary 𝜕𝐷 + must satisfy

Next, we focus on the solution of Equation (68). We begin with the scenario in which the parameters 𝐴 -2 , 𝐴 2 , 𝜔 0 are such that 𝐴 2 + 𝑛𝜔 0 𝐴 -2 ≠ 0 for all 𝑛 ∈ ℤ. Then, (68) is a quadratic equation with solutions

where the discriminant Δ is given by

Observe that Δ is itself a quadratic in the quantity 𝑛𝜔 0 with corresponding discriminant

] .

Note that Δ ≥ 0 since 𝐴 -2 > 0.

Case 1: Δ > 0. This case amounts to either 𝐴 1 ≠ 0 or 𝐴 2 + 𝐴 0 𝐴 -2 ≠ 0. Then, Δ has two distinct real roots

which determine a bifurcation in solution behavior. If 𝑛𝜔 0 > 𝜔 + 𝑐𝑟 or 𝑛𝜔 0 < 𝜔 - 𝑐𝑟 then Δ < 0 and the solutions 𝑘 + , 𝑘 -of equation (68) given by (73) form a complex conjugate pair; if 𝜔 - 𝑐𝑟 < 𝑛𝜔 0 < 𝜔 + 𝑐𝑟 , then Δ > 0 and 𝑘 + , 𝑘 -are real and distinct; and if 𝑛𝜔 0 = 𝜔 + 𝑐𝑟 or 𝑛𝜔 0 = 𝜔 - 𝑐𝑟 , then 𝑘 + , 𝑘 -reduce to a double real solution.

Importantly, for each 𝑛 ∈ ℤ, there is precisely one solution of (68) that lies on 𝜕𝐷 -. This is straightforward to see if 𝑛𝜔 0 > 𝜔 + 𝑐𝑟 or 𝑛𝜔 0 < 𝜔 - 𝑐𝑟 since, by the definition of the regions 𝐷 + and 𝐷 -, any nonreal solution of (68) must lie either on 𝜕𝐷 + or on 𝜕𝐷 -. On the other hand, if 𝑘 ± ∈ ℝ, then more details are needed:

(i) If 𝐴 1 = 0 and 𝐴 2 + 𝐴 0 𝐴 -2 > 0, then (70) requires 𝑘 ≤ 0. However, by Vieta's formulae, 𝑘 + + 𝑘 -= 0. Thus, precisely one (possibly double, if equal to zero) solution of (68) lies on 𝜕𝐷 - and, since ℝ ⊂ (𝜕𝐷 + ∩ 𝜕𝐷 -), the remaining solution lies on 𝜕𝐷 + . (ii) If 𝐴 1 = 0 and 𝐴 2 + 𝐴 0 𝐴 -2 < 0, then (70) requires 𝑘 ≥ 0 and, as in (i), precisely one (possibly double, if equal to zero) solution of (68) lies on 𝜕𝐷 -and the remaining solution lies on 𝜕𝐷 + . (iii) If 𝐴 1 > 0, then, since 𝐴 -2 > 0, (70) requires that 𝛽 -≤ 𝑘 ≤ 𝛽 + , where

First, we claim that 𝛽 -≤ 𝑘 -≤ 𝛽 + (and hence, via (70), that 𝑘 -∈ 𝜕𝐷 -). Indeed, note that 𝜔 ± 𝑐𝑟 =

𝐴 1 𝛽 ± 2 + 𝐴 0 and, since 𝑘 ± ∈ ℝ amounts to having

and 𝐴 1 + √ Δ > 0, we infer that

, which, in view of the inequalities 0 <

≤ 1, 𝛽 + > 0, 𝛽 -< 0, allows us to conclude that 𝛽 -≤ 𝑘 -≤ 𝛽 + as claimed. Note that the equalities are attained only when Δ = 0, that is, only when 𝑛𝜔 0 = 𝜔 + 𝑐𝑟 or 𝑛𝜔 0 = 𝜔 - 𝑐𝑟 , in which case 𝑘 -(𝜔 + 𝑐𝑟 ) = 𝑘 + (𝜔 + 𝑐𝑟 ) = 𝛽 + or 𝑘 -(𝜔 - 𝑐𝑟 ) = 𝑘 + (𝜔 - 𝑐𝑟 ) = 𝛽 -, respectively. In addition, we claim that 𝑘 + ≥ 𝛽 + or 𝑘 + ≤ 𝛽 -, with equalities attained only when Δ = 0 (and hence that 𝑘 + ∉ 𝜕𝐷 -except when 𝑘 + = 𝑘 -). Indeed, in view of (73), (75), and the fact that 𝐴 1 𝐴 -2 > 0, if 𝐴 2 + 𝑛𝜔 0 𝐴 -2 > 0, then the inequality

and, in turn, 𝑘 + ≤ 𝛽 -, establishing the claim. (iv) Finally, if 𝐴 1 < 0, then, since 𝐴 -2 > 0, we have 𝛽 + < 0 and 𝛽 -> 0, and (70) requires either 𝑘 ≥ 𝛽 -or 𝑘 ≤ 𝛽 + . Via entirely analogous steps as in (iii), using the fact that 𝜔 - 𝑐𝑟 ≤ 𝑛𝜔 0 ≤ 𝜔 + 𝑐𝑟 , we obtain 𝛽 + ≤ 𝑘 + ≤ 𝛽 -and 𝑘 -≥ 𝛽 -(for 𝐴 2 + 𝑛𝜔 0 𝐴 -2 < 0) or 𝑘 -≤ 𝛽 + (for 𝐴 2 + 𝑛𝜔 0 𝐴 -2 > 0), with equalities possible only when Δ = 0, that is, only when 𝑛𝜔 0 = 𝜔 + 𝑐𝑟 or 𝑛𝜔 0 = 𝜔 - 𝑐𝑟 , in which case 𝑘 -(𝜔 + 𝑐𝑟 ) = 𝑘 + (𝜔 + 𝑐𝑟 ) = 𝛽 + or 𝑘 -(𝜔 - 𝑐𝑟 ) = 𝑘 + (𝜔 - 𝑐𝑟 ) = 𝛽 -, respectively. Thus, once again, we see that there exists precisely one (possibly double) root along 𝜕𝐷 -, namely, 𝑘 -, while the other root, that is, 𝑘 + , lies on 𝜕𝐷 + . Therefore, in the case Δ > 0, we have shown that for each 𝑛 ∈ ℤ, there is precisely one (possibly double) solution of (68) that lies on 𝜕𝐷 -, with the remaining solution lying on 𝜕𝐷 + .

Case 2: Δ = 0. This case is much simpler because it requires 𝐴 1 = 𝐴 2 + 𝐴 0 𝐴 -2 = 0, and so, it reduces our second-order model (67) to the equation

which can be solved by elementary techniques. Indeed, the ordinary differential equation in 𝑥 along with the boundary conditions 𝑢(0, 𝑡) = 𝑔 0 (𝑡) and lim 𝑥→∞ 𝑢(𝑥, 𝑡) = 0 can be integrated to yield an expression for the quantity 𝑢 𝑡 -𝑖𝐴 0 𝑢, which can be integrated in 𝑡 to imply, upon imposing the initial condition 𝑢(𝑥, 0) = 𝑢 0 (𝑥), 𝑢(𝑥, 𝑡) = 𝑒 𝑖𝐴 0 𝑡 𝑢 0 (𝑥) +

(77)

The analysis of the scenario 𝐴 2 + 𝑛𝜔 0 𝐴 -2 ≠ 0 for all 𝑛 ∈ ℤ is complete and leads to the following conclusions:

• If either 𝐴 1 ≠ 0 or 𝐴 2 + 𝐴 0 𝐴 -2 ≠ 0, then one of the solutions 𝑘 + , 𝑘 -of the quadratic (68) given by (73) lies on the contour 𝜕𝐷 -, whereas the other solution does not lie on 𝜕𝐷 -(with the exception of the special case 𝑘 + = 𝑘 -). Denote the solution along 𝜕𝐷 -by 𝑘 - 1 and the solution off 𝜕𝐷 - (which, as explained above, lies on 𝜕𝐷 + ) by 𝑘 0 . Then, setting 𝑘 = 𝑘 - 1 into the D-N map condition (30), which for the second-order model (67) reads

and using the Vieta formula

for Equation (68), we obtain the D-N relation 𝑏 𝑛 = 𝑖𝑘 0 𝑎 𝑛 that was conjectured in (32). • If 𝐴 1 = 𝐴 2 + 𝐴 0 𝐴 -2 = 0, then the solution to the Dirichlet problem of Equation (67) on the half-line is given by the simple formula (77) and the D-N map can be computed by directly differentiating this formula.

Finally, we consider the scenario of 𝐴 2 + 𝑛 * 𝜔 0 𝐴 -2 = 0 for some 𝑛 * ∈ ℤ, that is, the possibility that the parameters 𝐴 -2 , 𝐴 2 , 𝜔 0 are such that 𝑛 * ∶= -𝐴 2 ∕(𝜔 0 𝐴 -2 ) ∈ ℤ. Then, Equation (30) does not yield information about 𝑏 𝑛 * , as either it becomes an identity (if 𝐴 1 = 0) or it forces 𝑎 𝑛 * = 0 (if 𝐴 1 ≠ 0).

If 𝐴 1 = 0, then Equation (68) evaluated at 𝑛 * implies 𝐴 0 𝐴 -2 + 𝐴 2 = 0 so we recover the simple model (76) with the explicit solution formula (77).

If 𝐴 1 ≠ 0, then we note that under the transformation 𝑢 ↦ 𝑢 + 𝑎 𝑛 * 𝑒 𝑖𝑛 * 𝜔 0 𝑡 , which in view of (20) eliminates the term corresponding to 𝑛 * from the series expansion of 𝑢(0, 𝑡), Equation ( 67) becomes ( 1 -𝐴 -2 𝜕 2 𝑥 )( 𝑢 𝑡 + 𝑖𝑛 * 𝜔 0 𝑎 𝑛 * 𝑒 𝑖𝑛 * 𝜔 0 𝑡 ) = 𝑖𝐴 0 ( 𝑢 + 𝑎 𝑛 * 𝑒 𝑖𝑛 * 𝜔 0 𝑡 ) + 𝐴 1 𝑢 𝑥 + 𝑖𝐴 2 𝑢 𝑥𝑥 .

If 𝐴 0 𝐴 -2 + 𝐴 2 = 0, then 𝑛 * 𝜔 0 = 𝐴 0 and Equation ( 79) is identical to Equation (67), so without loss of generality, we can take 𝑎 𝑛 * = 0. In particular, this is the case for the linear BBM equation considered below. On the other hand, if 𝐴 0 𝐴 -2 + 𝐴 2 ≠ 0, then 𝑛 * 𝜔 0 ≠ 𝐴 0 and (79) is different from (67), so the case 𝑛 * = -𝐴 2 ∕(𝜔 0 𝐴 -2 ) ∈ ℤ and 𝐴 1 (𝐴 0 𝐴 -2 + 𝐴 2 ) ≠ 0 is not covered by our analysis unless the associated Dirichlet datum actually satisfies 𝑎 𝑛 * = 0.

Concerning the value of 𝑏 𝑛 * , if the Dirichlet datum is actually periodic (and not just asymptotically periodic), then using the compatibility between the initial and boundary values at 𝑥 = 𝑡 = 0 and assuming that 𝑢(𝑥, 𝑡) is continuously differentiable up to the boundary, we find 𝑏 𝑛 * = 𝑢 ′ 0 (0) -∑ 𝑛∈ℤ,𝑛≠𝑛 * 𝑏 𝑛 , with the coefficients 𝑏 𝑛 for 𝑛 ≠ 𝑛 * given by the first of the D-N relations (32). Furthermore, if 𝐴 2 = 0 and hence 𝑛 * = 0, then 𝑏 𝑛 * = 𝑏 0 can be determined even if the solution is only asymptotically periodic. In particular, evaluating equation (67) at 𝑥 = 0 and then integrating over [𝑡, 𝑡 + 2𝜋 𝜔 0 ] with 𝑡 arbitrary, we exploit the periodicity of the relevant exponentials to deduce, by taking the limit 𝑡 → ∞,

Note that, in the particular case of the linear BBM equation, the above expression yields 𝑏 0 = 0. This proves Theorem 1 when 𝐴 3 = 0.

Remark (Connection to the radiation condition ( 4)). The group velocity associated with the second-order equation ( 67) is

. Thus, for a real solution of the quadratic (68), the condition that this solution lies on 𝜕𝐷 + (Equation ( 72)) amounts precisely to the group velocity being nonnegative at that solution. In particular, in the nontrivial scenario of either 𝐴 1 ≠ 0 or 𝐴 2 + 𝐴 0 𝐴 -2 ≠ 0, recall that 𝑘 0 ∈ 𝜕𝐷 + ∩ ℝ is possible if either 𝜔 - 𝑐𝑟 < 𝑛𝜔 0 < 𝜔 + 𝑐𝑟 or 𝑛𝜔 0 = 𝜔 ± 𝑐𝑟 . In the first case, we have a strictly positive group velocity, while in the second case, the group velocity is zero. Furthermore, if 𝑘 0 ∈ 𝜕𝐷 + ⧵ ℝ, then Im(𝑘 0 ) > 0.

This proves Corollary 1 when 𝐴 3 = 0.

The linear BBM on the half-line

Recall that the main idea behind deriving the D-N map (32) is that each root of Ω + 𝑛𝜔 0 (1 + 𝐴 -2 𝑘 2 ) that occurs along 𝜕𝐷 -should be a removable singularity of the right-hand side of ( 26) and hence a solution of (78). In this connection, the analysis provided above for the quadratic (68) associated with the second-order model (67) justifies the assumption made in the derivation of (32) that precisely one of the two roots of Ω + 𝑛𝜔 0 (1 + 𝐴 -2 𝑘 2 ) lies on the contour 𝜕𝐷 -. Nevertheless, the additional assumption that only this root must satisfy (78) has not yet been justified and, although this is an easy task for 𝑛𝜔 0 outside the interval [𝜔 - 𝑐𝑟 , 𝜔 + 𝑐𝑟 ]-since then the root off 𝜕𝐷 - is in the upper half-plane where ( 17) is generally not defined-it becomes more subtle for values within the interval. In what follows, we address this outstanding uniqueness issue in the case of the wavemaker problem (53) for the linear BBM equation, which reads 𝑢 𝑡 + 𝑢 𝑥 -𝑢 𝑥𝑥𝑡 = 0 (81) and corresponds to the choice of parameters 𝐴 0 = 𝐴 2 = 0, 𝐴 -2 = 1, 𝐴 1 = -1 in the general model (67). We remark that, despite the absence of the second spatial derivative term, the mixed derivative term-an essential element of the model (67) when compared to the model ( 35)is still present in Equation (81) and makes the analysis quite different from that of the linear KdV equation.

In the case of the linear BBM equation (81), the regions 𝐷 -and 𝐷 + of (69) and ( 71) are given by

that is, 𝐷 + is the interior of the upper half of the unit disk centered at the origin and 𝐷 -is the exterior of that disk in the lower half of the complex plane, as shown in Figure 2. The critical interval endpoints are symmetric 𝜔 ± 𝑐𝑟 = ±𝜔 𝑐𝑟 = ±1∕2. Furthermore,

when (cf. Equations ( 68) and ( 73))

(84)

14679590, 2024, 3, Downloaded from https://onlinelibrary.wiley.com/doi/10.1111/sapm.12665 by University Of Colorado Librari, Wiley Online Library on [28/04/2024]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License

F I G U R E 2

The regions 𝐷 + and 𝐷 -for the linear BBM equation as given by (82). In the case of the wavemaker data (53), for 0 < |2𝜔 0 | < 1, the two roots of 𝜔 + 𝑛𝜔 0 are real and distinct with 𝑘 - 1 ∈ 𝜕𝐷 -, 𝑘 0 ∈ 𝜕𝐷 + as shown. Then, to compute the long-time asymptotics of the solution, we deform the contour of integration in the Fokas formula (85) from  to the contour  along which Re(𝑖𝜔) > 0 (except for 𝑘 = ±1 where Re(𝑖𝜔) = 0). This deformation yields a residue contribution from 𝑘 0 ∈ 𝜕𝐷 + (as opposed to 𝑘 - 1 ∈ 𝜕𝐷 -), justifying the assumption that only 𝑘 - 1 (and not 𝑘 0 ) should satisfy (78).

If |2𝑛𝜔 0 | > 1, then (84) form a pair of complex conjugate roots that lie on the unit circle. Thus, due to the domain of (17), only the root lying in the lower half of the complex plane-and hence on 𝜕𝐷 --must satisfy (78). The same is true if |2𝑛𝜔 0 | = 1, since then the roots (84) reduce to a single double root equal to ∓1 ∈ 𝜕𝐷 -. On the other hand, if 0 < |2𝑛𝜔 0 | < 1, then the roots (84) are real and distinct, with one of them, denoted by 𝑘 0 , lying inside (-1, 1) ⊂ 𝜕𝐷 + and the other one, denoted by 𝑘 - 1 , lying in [-1, 1] 𝑐 ⊂ 𝜕𝐷 -(see Figure 2). In this case, as noted above, since both roots are in the domain of (17), it is not immediately clear why only 𝑘 - 1 , and not 𝑘 0 , should satisfy (78) as a removable singularity of (26). Next, we will show that this is indeed the case by taking advantage of the powerful solution representation provided by the Fokas method.

In particular, let us consider the linear BBM equation (81) in the context of the wavemaker problem associated with the data (53). In this case, observe that the scenario of two distinct real roots of (83) corresponds to the range 0 < |2𝜔 0 | < 1, since 𝑢(0, 𝑡) = -∑ 𝑛=±1 sgn(𝑛)

2𝑖

𝑒 𝑖𝑛𝜔 0 𝑡 and, regardless of which root ends up satisfying (78), the resulting D-N relation will imply that the Fourier series coefficient 𝑏 𝑛 of the Neumann value 𝑢 𝑥 (0, 𝑡) vanishes whenever the coefficient 𝑎 𝑛 of 𝑢(0, 𝑡) does. Hence, we work with 0 < |2𝜔 0 | < 1 and two distinct real roots as shown in Figure 2. The solution formula is given by 21

where the contour  is a positively oriented, closed contour in the upper half of the complex plane and enclosing 𝑘 = 𝑖, as shown in Figure 2 and 𝜔 = 𝜔(𝑘) ∶= 𝑘 1 + 𝑘 2 (86) with 𝜔 ′ denoting the derivative of 𝜔. Note that the integrand in (85) has essential singularities at 𝑘 = ±𝑖 due to the exponential 𝑒 -𝑖𝜔𝑡 . In addition, for 𝑘 ≠ ±𝑖, we have Re( 𝑖Ω 1+𝑘 2 ) = Re(𝑖𝜔), and so, inside the regions 𝐷 -and 𝐷 + of (82), we have Re(𝑖𝜔) < 0, whereas outside those regions, we have Re(𝑖𝜔) > 0.

Thanks to analyticity and Cauchy's theorem, the contour of integration  in (85) can be deformed to the contour  depicted in Figure 2. Thus, differentiating in 𝑥, we have

Since  does not cross any singularities, we can split the integral with respect to 𝑘 in two parts:

The contour  encloses the points 𝑘 = 𝑖 and 𝑘 = 𝑘 0 , which are simple poles of the integrand of the first of the above integrals. Hence, by Cauchy's residue theorem and the fact that 𝑘 0 satisfies (84),

In particular, first taking the limit 𝑥 → 0 + and then performing the above asymptotic analysis as 𝑡 → ∞, we obtain lim 𝑥→0 + 𝑢 𝑥 (𝑥, 𝑡) ∼ ∑ 𝑛=±1 𝑖𝑘 0 𝑎 𝑛 𝑒 𝑖𝑘 0 𝑥+𝑖𝑛𝜔 0 𝑡 as 𝑡 → ∞, in agreement with the D-N map (32).

A similar analysis of the solution formula (85) implies the asymptotic solution for fixed 𝑥 ≥ 0 𝑢(𝑥, 𝑡) ∼ ∑ 𝑛=±1 𝑎 𝑛 𝑒 𝑖𝑘 0 𝑥+𝑖𝑛𝜔 0 𝑡 , 𝑡 → ∞, 𝑥 = (1).

By linearity, this result extends to the more general boundary data 𝑔 0 (𝑡) ∼ ∑ 𝑛∈ℤ 𝑎 𝑛 𝑒 𝑖𝑛𝜔 0 𝑡 . We have proven the asymptotic time periodicity of the solution to the linear BBM wavemaker problem so that Theorem 1 and the radiation condition (4) (Corollary 1) hold for any frequency 𝜔 0 ∈ ℝ without the need to assume asymptotic time periodicity of the solution.

The asymptotic solution (92) is only valid for fixed 𝑥 as 𝑡 → ∞. In Section 4, we extend the analysis of the linear BBM wavemaker problem to obtain an asymptotic expansion of the solution for large 𝑡 that is uniform in 𝑥 ≥ 0. We then observe that this analysis quantitatively agrees with viscous core-annular experiments in Section 5.

Remark. To establish estimate (90), it is crucial for  to lie in regions where Re(𝑖𝜔) ≥ 0. This requirement implies that only 𝑘 0 (and not 𝑘 - 1 ) can be enclosed by , because attempting to also enclose 𝑘 - 1 would push C into regions where Re(𝑖𝜔) < 0. Consequently, only 𝑘 0 (and not 𝑘 - 1 ) contributes in the asymptotics (91) (via the residue calculation (89)).

LONG-TIME ASYMPTOTICS FOR THE LINEAR KdV EQUATION WITH A SINUSOIDAL BOUNDARY CONDITION

Section 2.2 derives the D-N map and the asymptotic solution to the linear KdV wavemaker problem ( 50)-(53) for 𝑡 → ∞ and a fixed 𝑥 > 0. Herein, we study the KdV wavemaker problem in a different asymptotic limit where 𝑡 → ∞ while 𝜉 = 𝑥∕𝑡 = (1) is fixed, that is, we consider the solution along rays in the 𝑥-𝑡 plane. This limit provides information about how waves propagate away from the origin in the wavemaker problem.

Recall the wavemaker initial-boundary value problem (53) in which 𝑢(𝑥, 0) = 0, 𝑥 ≥ 0 and 𝑢(0, 𝑡) = -sin(𝜔 0 𝑡), 𝑡 ≥ 0 where 𝜔 0 > 0 is the input frequency. To obtain the long-time asymptotics, we use the solution formula for the linear KdV wavemaker problem ( 50)-(53) given by (54), where the contour of integration 𝜕𝐷 + is shown in Figure 1. Note that the integrand in (54) has removable singularities at the six simple zeros 𝑘 1 , … , 𝑘 6 of the two cubic equations Ω + 𝑛𝜔 0 = 0, 𝑛 = ±1 where Ω is defined in (25). For 0 < 𝜔 0 < 2 3 √ 3

, we have

)

)

)

) -sin

)

14679590, 2024, 3, Downloaded from https://onlinelibrary.wiley.com/doi/10.1111/sapm.12665 by University Of Colorado Librari, Wiley Online Library on [28/04/2024]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License ). Because all the zeros are real as shown in Figure 3(a), we call this the subcritical regime. On the other hand, in the supercritical regime 𝜔 0 > 2 3 √ 3

, we have

where 𝑟 1 = 9𝜔 0 -√ 81𝜔 2 0 -12 and 𝑟 2 = 9𝜔 0 + √ 81𝜔 2 0 -12 are real numbers. In this case, the location of the zeros is shown in Figure 3(b) and described as follows: 𝑘 1 is complex and on the boundary of 𝐷 - 2 , 𝑘 2 ∈ (-∞, -2∕ √ 3) is real, 𝑘 3 is complex and on the right boundary of 𝐷 + , 𝑘 4 ∈ (2∕ √ 3, ∞) is real, 𝑘 5 is complex and on the left boundary of 𝐷 + , and 𝑘 6 is also complex and on the boundary of 𝐷 - 1 . The poles are symmetric with respect to the imaginary and real axes and lie on the ellipse In addition, 𝑘 1,3,5,6 are also on the hyperbola 3Re(𝑘) 2 -Im(𝑘) 2 = 1, 𝑗 = 1, 3, 5, 6, given as 𝜕𝐷 + . We define the critical frequency to be

Importantly, once the two 𝑡-exponentials in the solution formula (54) are separated, the removable singularities 𝑘 1 , … , 𝑘 6 turn into simple poles for the respective integrands. Therefore, in order to avoid crossing these poles, prior to splitting (54), we first indent 𝜕𝐷 + so that the poles are excluded from the interior of 𝜕𝐷 + , as shown in Figure 3. In what follows, 𝜕𝐷 + denotes the indented contour shown in Figure 3. Then, by substituting the linear dispersion relation Ω = 𝑘 -𝑘 3 and expressing the complex exponentials as trigonometric functions, we are able to rewrite the solution (54) as 𝑢(𝑥, 𝑡) = 𝑢 1 (𝑥, 𝑡) + 2 (𝑥, 𝑡), (96a)

Note that the dependence of 𝑓 on 𝑡 has been suppressed. We begin by observing that, since all of the poles 𝑘 1 , … , 𝑘 6 lie outside the contour 𝜕𝐷 + , by Cauchy's theorem inside the region enclosed by 𝜕𝐷 + and 𝐶 𝑅 (the semicircular arcs of radius 𝑅 in Figure 3), we have 𝑢 1 (𝑥, 𝑡) = 1 2𝜋 ∫ 𝑘∈𝐶 0 𝑓(𝑘)𝑒 𝑖𝑘𝑥 𝑑𝑘. To evaluate along the contour 𝐶 𝑅 , we first notice that the integrand 𝑓(𝑘) in (96b) is uniformly bounded in |𝑘| (𝑓(𝑘) = (𝑘 -1 ), 𝑘 → ∞). Since 𝑥 ≥ 0, we apply Jordan's lemma in the upper half plane by sending 𝑅 → ∞ and conclude that 𝑢 1 (𝑥, 𝑡) ≡ 0 → 𝑢(𝑥,𝑡) = 𝑢 2 (𝑥, 𝑡).

(97)

The remaining term 𝑢 2 (𝑥, 𝑡) defined by (96c) will be asymptotically evaluated in the large-time limit by the method of steepest descent. The difference is that 𝑒 𝑖𝑘𝑥 in 𝑢 1 (𝑥, 𝑡) decays in 𝐷 + so that one can apply Jordan's lemma along 𝐶 𝑅 , whereas 𝑢 2 (𝑥, 𝑡) involves 𝑒 𝑖𝑘𝑥-𝑖(𝑘-𝑘 3 )𝑡 that does not necessarily decay in 𝐷 + . We now asymptotically approximate 𝑢 2 (𝑥, 𝑡) in (96c) using the method of steepest descent. Define 𝜉 = 𝑥∕𝑡 and the phase 𝜙(𝑘) = 𝑖𝑘𝜉 -𝑖(𝑘 -𝑘 3 ), then 𝑢 2 = 𝜔 0 2𝜋 ∫ 𝜕𝐷 + 𝑔(𝑘)𝑒 𝜙(𝑘)𝑡 𝑑𝑘. The saddle points are found by requiring 𝜙 ′ (𝜌 𝑗 ) = 0 so that

The saddle points 𝑘 1,2 are real for 0 ≤ 𝜉 ≤ 1 and imaginary for 𝜉 > 1. Let 𝑘 = 𝑝 + 𝑖𝑞, 𝑝, 𝑞 ∈ ℝ, so that 𝜙(𝑘) = 𝑞 -3𝑝 2 𝑞 + 𝑞 3 -𝑞𝜉 + 𝑖(-𝑝 + 𝑝 3 -3𝑝𝑞 2 + 𝑝𝜉). We now deform the original contour 𝜕𝐷 + onto the steepest descent paths that are defined by Im{𝜙(𝑘)} = Im{𝜙(𝜌 𝑗 )} and Re{𝜙(𝑘) -𝜙(𝜌 𝑗 )} < 0, 𝑗 = 1, 2. 98) and the steepest descent contour (red) for the solution (96c) of the linear KdV wavemaker problem in the case of (A) 0 < 𝜉 < 1 and (B) 𝜉 > 1, against the original contour 𝜕𝐷 + (black). The contours 𝐶 1 ∪ 𝐶 2 ∪ 𝐶 3 and 𝐶 4 ∪ 𝐶 5 ∪ 𝐶 6 in (A) are given by Equations ( 99a) and (99b), respectively, and the contours 𝐶 1 ∪ 𝐶 2 ∪ 𝐶 3 ∪ 𝐶 4 in (B) are described by (108). The color map indicates levels of Im(𝜙(𝑘)) with Im(𝜙(𝑘)) constant along the steepest descent path. Recall that the contour 𝜕𝐷 + avoids the poles 𝑘 1 , … , 𝑘 6 , as shown in Figure 3.

Im(𝜙(𝑘)) = Im(𝜙(𝜌 2 )) ⇒ -𝑝 + 𝑝 3 -3𝑝𝑞 2 + 𝑝𝜉 = -2 9

respectively. Note that the exponential term satisfies |𝑒 𝜙(𝑘)𝑡 | = 1 on the steepest descent contour for 0 < 𝜉 < 1, so it is purely oscillatory. On the right steepest descent path 𝐶 1 ∪ 𝐶 2 ∪ 𝐶 3 in Figure 4(a), the dominant contribution comes from the neighborhood of 𝜌 1 . Parameterizing this neighborhood by 𝑘 = 𝜌 1 + (1 + 𝑖)𝑠 with small 𝑠 ∈ (-𝜖, 𝜖), we obtain

, 𝑡 → ∞.

(100)

Similarly, on the left steepest descent path 𝐶 4 ∪ 𝐶 5 ∪ 𝐶 6 in Figure 4(a), we parameterize 𝑘 = 𝜌 2 + (1 -𝑖)𝑠, 𝑠 ∈ (-𝜖, 𝜖). The dominant contribution is the complex conjugate of (100)

, 𝑡 → ∞.

(101)

14679590, 2024, 3, Downloaded from https://onlinelibrary.wiley.com/doi/10.1111/sapm.12665 by University Of Colorado Librari, Wiley Online Library on [28/04/2024]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License , 0], their relative ordering is important. The boundary between Regions I and II.a is defined by the coalescence 𝑘 3 = 𝜌 1 and 𝑘 6 = 𝜌 2 so that

We are now in a position to deform the original contour 𝜕𝐷 + onto the steepest descent paths. There are two cases to consider. First, the region 0 < 𝜉 < 𝜔 ′ (𝑘 3 ), which corresponds to the poles lying between the saddle points as shown in Region I of Figure 5.

Notice that the integrand 𝑔(𝑘)𝑒 𝜙(𝑘)𝑡 of (96c) decays for 𝑘 → ∞ in the region between 𝐷 1 and 𝐶 1 (also 𝐷 5 and 𝐶 6 ) with |𝑒 𝜙(𝑘)𝑡 | = 𝑒 -Im(𝑘)𝜉𝑡+Im(𝑘-𝑘 3 )𝑡 , 𝜉 > 0, 𝑡 > 0, where Im(𝑘) > 0 and Im(𝑘 -𝑘 3 ) < 0 in this region.

Thus, one can simply deform 𝐷 1 onto 𝐶 1 and 𝐷 5 onto 𝐶 6 by Cauchy's theorem and Jordan's lemma so that ∫ 𝐷 1 𝑔(𝑘)𝑒 𝜙(𝑘)𝑡 𝑑𝑘 = ∫ 𝐶 1 𝑔(𝑘)𝑒 𝜙(𝑘)𝑡 𝑑𝑘 and ∫ 𝐷 5 𝑔(𝑘)𝑒 𝜙(𝑘)𝑡 𝑑𝑘 = -∫ 𝐶 6 𝑔(𝑘)𝑒 𝜙(𝑘)𝑡 𝑑𝑘. The absence of poles between 𝐷 2 and 𝐶 2 (𝐷 4 and 𝐶 5 ) and Cauchy's theorem imply ∫ 𝐷 2 𝑔(𝑘)𝑒 𝜙(𝑘)𝑡 𝑑𝑘 = ∫ 𝐶 2 𝑔(𝑘)𝑒 𝜙(𝑘)𝑡 𝑑𝑘 and ∫ 𝐷 4 𝑔(𝑘)𝑒 𝜙(𝑘)𝑡 𝑑𝑘 = -∫ 𝐶 5 𝑔(𝑘)𝑒 𝜙(𝑘)𝑡 𝑑𝑘. Finally, since the poles 𝑘 3 and 𝑘 6 lie inside the region enclosed by 𝐷 3 , 𝐶 3 and 𝐶 4 , using Cauchy's residue theorem, we compute ∫ 𝐷 3 +𝐶 3 -𝐶 4 𝑔(𝑘)𝑒 𝜙(𝑘)𝑡 𝑑𝑘 = -2𝜋𝑖(Res(𝑘 3 ) + Res(𝑘 6 )), Res(𝑘 3 ) = 1 2𝜔 0 𝑒 𝑖(𝑘 3 𝜉𝑡-𝜔 0 𝑡) , Res(𝑘 6 ) = -Res(𝑘 3 ) * , where 𝑘 3 = -𝑘 6 > 0 satisfies the linear dispersion relation 𝑘 3 -𝑘 3 3 = 𝜔 0 . Combining these results, we obtain 𝑢 2 (𝜉; 𝑡) = sin[(𝑘 3 𝜉 -𝜔 0 )𝑡] + 𝜔 0

2𝜋

∫ 𝐶 1 +𝐶 2 -𝐶 3 +𝐶 4 -𝐶 5 -𝐶 6 𝑔(𝑘)𝑒 𝜙(𝑘)𝑡 𝑑𝑘. The remaining integral on the right-hand side corresponds to the saddle point contributions (100) and ( 101), which are dominated by the pole contributions. Therefore, the large 𝑡 solution to the wavemaker problem is

where 𝑘 0 ≡ 𝑘 3 in this case, which satisfies the radiation condition (4).

In contrast, the case 𝜔 ′ (𝑘 0 ) < 𝜉 < 1 occurs when the saddles are in between the poles depicted in Region II.a of Figure 5. Using Cauchy's theorem, we can deform 𝜕𝐷 + to the steepest descent contour ∪ 6 𝑗=1 𝐶 𝑗 with no pole contributions. Hence, we obtain 𝑢 2 (𝑥, 𝑡) = 𝜔 0

2𝜋

∫ 𝐶 1 +𝐶 2 -𝐶 3 +𝐶 4 -𝐶 5 -𝐶 6 𝑔(𝑘)𝑒 𝜙(𝑘)𝑡 𝑑𝑘, so the saddle point contributions (100) and ( 101) are asymptotically dominant

In summary, when 0 < 𝜔 0 < 2∕(3 √ 3), the leading order solution (103) is a traveling wave with real wave number 𝑘 3 and frequency 𝜔 0 satisfying the dispersion relation 𝜔 0 = 𝑘 3 -𝑘 3 3 as 𝑡 → ∞ when 0 < 𝜉 = 𝑥∕𝑡 < 𝜔 ′ (𝑘 3 ). This regime is dominated by the pole contributions from contour deformation. When 𝑥∕𝑡 > 𝜔 ′ (𝑘 3 ), that is, ahead of the boundary-generated, invading wavefront, the leading order solution (104) exhibits algebraic decay in time with 𝑢(𝑥, 𝑡) = (𝑡 -1∕2 ), 𝜉 = 𝑥∕𝑡 = (1), 𝑡 → ∞. Here, the solution is dominated by saddle point contributions. In contrast to the previous case (103) when 0 < 𝜉 < 𝜔 ′ (𝑘 0 ), the asymptotic solution (104) is not time periodic.

Note that the obtained asymptotic approximations are not uniform in 𝜉. For example, Equation ( 104) is singular when 𝜉 = 𝜔 ′ (𝑘 3 ). The cause of nonuniformity is the collisions of the poles 𝑘 3,6 with the saddle points 𝜌 1,2 . Uniform asymptotic approximations involving the error function with complex argument can be obtained. 26-28 (a.2) Supercritical regime: 𝜔 0 > 2∕(3

3). In this case, some poles are complex, including 𝑘 3,5 , which are avoided by semicircle indentations of the contour 𝜕𝐷 + (cf. Figure 3A). To determine the pole locations relative to the steepest descent contours 𝐶 𝑗 , 𝑗 = 1, … , 6 depicted in red in Figure 4(A) (case 0 < 𝜉 < 1), notice that 𝐶 1,2 and 𝐶 5,6 intersect with 𝜕𝐷 + at

Downloaded from https://onlinelibrary.wiley.com/doi/10.1111/sapm.12665 by University Of Colorado Librari, Wiley Online Library on [28/04/2024]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License (𝑝 2 , 𝑞 2 ) = ( -√ 2 + 𝜉 6 cos 𝜑 3 , √ -1 + 2 + 𝜉 2 cos 2 𝜑 3 ) , where 𝜑 = arctan[ 3 √ 𝜉(4-2𝜉+𝜉 2 ) 2 √ 2(𝜉-1) √ 1-𝜉 ]. Here, (𝑝 1 , 𝑞 1 ) is associated with 𝐶 1,2 ∩ 𝜕𝐷 + and (𝑝 2 , 𝑞 2 ) is associated with 𝐶 5,6 ∩ 𝜕𝐷 + . We now determine for what value of 𝜉 the poles 𝑘 3,5 coincide with (𝑝 1,2 , 𝑞 1,2

). Letting Re(𝑘 3 ) = 𝑝 1 , Re(𝑘 5 ) = 𝑝 2 and simplifying, results in

where 𝑟 1 is defined below (94). Equation ( 105) determines the curve in the 𝜔 0 -𝜉 plane that splits Regions II.b and III in Figure 5.

When 𝑙 1 (𝜉, 𝜔 0 ) < 0, the poles 𝑘 3,5 are avoided by the corner contours 𝐷 2,4 with Im(𝑘 3,5 ) > 0 (Region II.b in Figure 5). Using Cauchy's theorem and Jordan's lemma where appropriate, we deform each sector of 𝜕𝐷 + onto the steepest descent path with no pole contributions. This gives 𝑢 2 (𝜉; 𝑡) = 𝜔 0

2𝜋

∫ 𝐶 1 +𝐶 2 -𝐶 3 +𝐶 4 -𝐶 5 -𝐶 6 𝑔(𝑘)𝑒 𝜙(𝑘)𝑡 𝑑𝑘 so that the leading-order behavior of 𝑢(𝜉𝑡, 𝑡) is determined by the saddle contributions as in (104), with algebraic decay (𝑡 -1∕2 ).

When 𝑙 1 (𝜉, 𝜔 0 ) > 0 (see Region III in Figure 5), the poles 𝑘 3 and 𝑘 5 lie between the contours 𝐷 1 , 𝐶 1 and 𝐷 5 , 𝐶 6 , respectively. Thus, when invoking Cauchy's theorem, the residue contributions from these poles must be taken into account. Computing

because 𝑘 3 = -𝑘 * 5 . The pole 𝑘 3 corresponds to the singularity on 𝜕𝐷 + that determines the D-N map 𝑘 0 ∶= 𝑘 3 and satisfies the dispersion relation 𝑘 3 -𝑘 3 3 = 𝜔 0 . We find 𝑢 2 (𝜉𝑡, 𝑡) = 𝑒 -Im(𝑘 0 )𝜉𝑡 sin[(Re(𝑘 0 )𝜉 -𝜔 0 𝑡)] + 𝜔 0

2𝜋

∫ 𝐶 1 +𝐶 2 -𝐶 3 +𝐶 4 -𝐶 5 -𝐶 6 𝑔(𝑘)𝑒 𝜙(𝑘)𝑡 𝑑𝑘. In light of (100) and (101), we deduce

The second term in (107) is the dominant contribution, with (𝑡 -1∕2 ) for 𝜉 = (1). The first term is beyond all orders. In other words, the asymptotic behavior of the solution in Region III is the same as the asymptotic behavior (104) in Regions II.a and II.b. However, as we now analyze, both terms in (107) are required for obtaining an asymptotic expansion that is also valid for 𝑥 fixed and 𝑡 → ∞. Rewriting (107) in terms of 𝑥 = 𝜉𝑡 and 𝑡, we obtain 𝑢(𝑥, 𝑡) ∼𝑒 -Im(𝑘 0 )𝑥 sin(Re(𝑘 0 )𝑥 -𝜔 0 𝑡)

For 𝑥 ∈ [0, 𝐿], some 𝐿 > 0 and 𝑡 → ∞. This expression is asymptotically equivalent to 𝑢(𝑥, 𝑡) = 𝑒 -Im(𝑘 0 )𝑥 sin(Re(𝑘 0 )𝑥 -𝜔 0 𝑡) + (𝑡 -3∕2 ), 𝑡 → ∞.

The asymptotics of 𝑢(𝑥, 𝑡) match our earlier analysis using the D-N map for fixed 𝑥 and large 𝑡 given in Equation (34). Therefore, keeping both terms in (107) gives an asymptotic expansion that matches the limits 𝜉 = (1), 𝑡 → ∞ and 𝜉𝑡 = (1), 𝑡 → ∞.

is unbounded as 𝑡 → ∞. On the other hand, since 𝜙(𝜌 1 ) = -2 9 √ 3𝜉 -3(𝜉 -1) < 0, 𝑒 𝑡𝜙(𝜌 1 ) decays as 𝑡 → ∞. The deformed contour, which follows the steepest descent directions emanating from the saddle 𝜌 1 , is shown in red in Figure 4(B) and is described by setting 𝑘 = 𝑝 + 𝑖𝑞, 𝑝, 𝑞 ∈ ℝ and

Knowing that the main contribution along the contours 𝐶 1 , … , 𝐶 4 in Figure 4(B) comes from the neighborhood of the saddle 𝜌 1 , we parameterize 𝑘 = 𝜌 1 + 𝑠 with small 𝑠 ∈ (-𝜖, 𝜖) so that, using the Laplace method, we find

Now we deform 𝜕𝐷 + onto 𝐶 1 , … , 𝐶 4 taking the poles into consideration. (b.1) Subcritical regime:

3). This case corresponds to Region IV.a in Figure 5. The poles 𝑘 3,6 are real and lie below 𝜕𝐷 + . Consequently, they lie outside the region enclosed by 𝐶 2 , 𝐶 3 , 𝐷 2 , 𝐷 3 , 𝐷 4 , so that Cauchy's theorem implies 𝑢 2 (𝜉𝑡, 𝑡) =

𝑔(𝑘)𝑒 𝜙(𝑘)𝑡 𝑑𝑘. Using (109), the leading-order solution in this region decays exponentially

This case corresponds to Regions IV.b and IV.c in Figure 5. In particular, the poles 𝜅 3,5 are now near the right and left edge of the contour 𝜕𝐷 + , respectively, with Re(𝑘 3 ) = -Re(𝑘 5 ) and Im(𝑘 3 ) = Im(𝑘 5 ) > 0. Thus, we need to distinguish how they are located relative to the deformed contours 𝐶 1 , … , 𝐶 4 . Note from Figure 4(B) that these contours intersect with 𝜕𝐷 + when 𝑘 = 𝑝 + 𝑖𝑞, -𝑝 + 𝑝 3 -3𝑝𝑞 2 + 𝑝𝜉| 𝑞= √ 3𝑝 2 -1 = 0, that is, at the points

)

.

The poles 𝑘 3,5 coincide with these intersection points when 𝑝 3 = Re(𝑘 3 ) or, equivalently, 𝑝 4 = Re(𝑘 5 ), that is, when

where 𝑟 1 is defined below (94). Equation (111) determines the curve in the 𝜔 0 -𝜉 plane that splits Regions IV.b and IV.c in Figure 5. When 𝑙 2 (𝜉, 𝜔 0 ) < 0, then the poles 𝜅 3,5 lie outside the region enclosed by 𝐶 2 , 𝐶 3 , 𝐷 2 , 𝐷 3 , 𝐷 4 (see Region IV.b in Figure 5), thus 𝑢 2 (𝜉𝑡, 𝑡) = 𝜔 0

2𝜋

∫ 𝐶 1 +𝐶 2 -𝐶 3 -𝐶 4 𝑔(𝑘)𝑒 𝜙(𝑘)𝑡 𝑑𝑘 so that the leading-order behavior of 𝑢(𝜉𝑡, 𝑡) is given by (110). On the other hand, when 𝑙 2 (𝜉, 𝜔 0 ) > 0, then 𝑘 3 and 𝑘 5 are inside the regions enclosed by 𝐶 1 , 𝐷 1 and 𝐶 4 , 𝐷 5 , respectively (Region IV.c in Figure 5), so in view of (106), we have the additional residue contribution 𝑢 2 (𝜉𝑡, 𝑡) = 𝑒 -Im(𝑘 3 )𝜉𝑡 sin(Re(𝑘 3 )𝜉𝑡 -𝜔 0 𝑡) + 𝜔 0

2𝜋

∫ 𝐶 1 +𝐶 2 -𝐶 3 -𝐶 4 𝑔(𝑘)𝑒 𝜙(𝑘)𝑡 𝑑𝑘. Hence, using once again (109), the leading-order asymptotics of the solution in this region is described by 𝑢(𝜉𝑡, 𝑡) ∼𝑒 -Im(𝑘 3 )𝜉𝑡 sin(Re(𝑘 3 )𝜉𝑡 -𝜔 0 𝑡)

Observe that both terms on the right-hand side decay exponentially. In order to figure out which one of them is dominant, we notice that when 𝑙 2 > 0, that is,

√ 3𝜉 -3(𝜉 -1). Thus, (112) is dominated by the second term. Therefore, in the regime 𝜉 > 1, the leading-order behavior of 𝑢(𝜉𝑡, 𝑡) is always given by (110) in both the subcritical and the supercritical regime, that is, regardless of the range of values of 𝜔 0 . When 𝜉 passes through 1, the asymptotic expansion abruptly changes from oscillatory to exponential, for example, (107) to (112). This nonuniformity is due to the collision of the saddle points 𝜌 1,2 , which can be cast as a uniform asymptotic expansion by use of Airy-type integrals. 28 The large 𝑡 solution given in Equations ( 103), ( 104), (107), (110), and (112) proves Theorem 2. We now carry out a much simpler calculation of a subset of these results using modulation theory.

Remark (Linear modulation theory). Wave modulation theory is one of the primary techniques to describe slow modulations of periodic waves. Herein, we use linear modulation theory to obtain the wave's amplitude, wave number, and frequency for the linear KdV boundary value problem ( 50)-(53) and compare with the long-time asymptotic approximation just obtained. The modulation equations of wave conservation are 29 𝑘(𝜔) 𝑡 + 𝜔 𝑥 = 0, (113a)

where 𝜔 is the angular frequency and 𝑘(𝜔) is one of the roots of the linear dispersion relation 𝑘 -𝑘 3 = 𝜔. Translating the initial-boundary value problem ( 50)-(53) for 𝑢(𝑥, 𝑡) to the initial-boundary value problem for 𝑘(𝑥, 𝑡) and 𝑎(𝑥, 𝑡), we have

Equations (113) are a system of quasi-linear equations with the double characteristic speed (eigenvalue) 𝑐 𝑔 (𝜔 0 ) that has a one-dimensional eigenspace. Equation (113a) and the initialboundary conditions (114) are invariant under the hydrodynamic scaling symmetry 𝑥 → 𝑏𝑥, 𝑡 → 𝑏𝑡 for any 𝑏 > 0. Consequently, we seek a self-similar solution that preserves this symmetry: 𝜔 = 𝜔(𝜉), where 𝜉 = 𝑥∕𝑡. Since 𝑘(𝜔) 𝑡 = 𝜔 𝑡 ∕𝑐 𝑔 (𝜔), either 𝜔 ′ (𝜉) = 0 or 𝑐 𝑔 (𝜔) = 𝜉, where 𝑐 𝑔 (𝜔) = 1∕𝑘 ′ (𝜔) = 1 -3𝑘(𝜔) 2 is the group velocity. Combining these in a piecewise fashion, we obtain the solution

where 𝑘 0 = 𝑘(𝜔 0 ) satisfies the linear dispersion relation 𝑘 0 -𝑘 3 0 = 𝜔 0 . Here, we encounter the issue of selecting 𝑘(𝜔 0 ) from the three possible wave numbers that satisfy the cubic polynomial linear dispersion relation. In the context of modulation theory, there are two admissibility conditions for the solution (115) to be well defined: (i) 𝑐 𝑔 (𝜔 0 ) > 0 so that 𝜉 = 𝑥∕𝑡 varies across nonnegative values and (ii) Im 𝑘(𝜔 0 ) = 0 so that the solution is real-valued. These are equivalent to the radiation condition (4) when 0 < 𝜔 0 < 2∕(3

3). We proved in Section 2.2 the adherence of the linear KdV equation wavemaker problem to the radiation condition (4), so a unique 𝑘(𝜔 0 ) branch exists with these properties (𝑘 0 = 𝑘 3 from Equation ( 93)). With an input frequency 𝜔 0 ∈ (0, 2∕(3

3)), the solution consists of a plane wave with constant frequency 𝜔 0 and wave number 𝑘 0 for 0 < 𝜉 < 𝑐 𝑔 (𝜔 0 ). The group velocity demarcates the invading front of the plane wave into the quiescent region.

By integrating the generalized wave number 𝜃 𝑥 = 𝑘(𝑥∕𝑡) and generalized frequency 𝜃 𝑡 = -𝜔(𝑥∕𝑡) relations, as well as invoking the boundary data 𝑢(0, 𝑡) = -sin(𝜔 0 𝑡) (53), we can obtain the phase 𝜃(𝑥, 𝑡) up to an integration constant 𝜃 0 ∈ ℝ (more generally, 𝜃 0 can be slowly varying 30 )

𝑡(1 -𝑥) 3∕2 + 𝜃 0 , 𝑐 𝑔 (𝜔 0 )𝑡 < 𝑥 ≤ 𝑡.

(116)

14679590, 2024, 3, Downloaded from https://onlinelibrary.wiley.com/doi/10.1111/sapm.12665 by University Of Colorado Librari, Wiley Online Library on [28/04/2024]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License

In addition, solving (113b) yields the amplitude

The self-similar solution (117) for 𝑎 does not preserve the hydrodynamic scaling symmetry when 𝑐 𝑔 (𝜔 0 ) < 𝜉 ≤ 1 and 𝐹(𝜉) can be any smooth function. This is due to the degeneracy of the double characteristic speed 𝑐 𝑔 (𝜔 0 ). These results agree with the large 𝑡 asymptotic solutions in Equations ( 103) and (104). Note that neither the phase 𝜃 nor the amplitude 𝑎 are continuous at 𝜉 = 𝑐 𝑔 (𝜔 0 ). Here, 𝜃 0 and 𝐹(𝜉) cannot be determined through leading order modulation theory. However, we can appeal to our long-time asymptotic result (104) in order to identify

It is interesting to note that, for 𝑐 𝑔 (𝜔 0 )𝑡 < 𝑥 ≤ 𝑡, the phase 𝜃(𝑥, 𝑡) is independent of the boundary data. In this region, the phase is determined by a structural property of the linear KdV equation, namely, the saddle points 𝜌 1,2 , which are independent of 𝜔 0 . See the discussion leading up to (104). In summary, modulation theory is a simple method for determining some features of the longtime asymptotic solution. However, its leading order manifestation does not provide information on the integrated phase 𝜃 0 , the amplitude function 𝐹(𝜉), nor any information on the supercritical regime 𝜔 0 > 2∕(3

3) or when 𝜉 > 1. Yet, the properties that can be ascertained from modulation theory such as the frequency 𝜔, wave number 𝑘, and amplitude 𝑎 agree with the leading order, long-time asymptotic solutions (103) and (104) of Regions I and II.a, respectively in Figure 5.

LONG-TIME ASYMPTOTICS FOR THE LINEAR BBM EQUATION WITH A SINUSOIDAL BOUNDARY CONDITION

This section studies the linear BBM wavemaker problem in the asymptotic limit 𝑡 → ∞ and 𝜉 = 𝑥∕𝑡 = (1). The analysis is similar to that for the linear KdV wavemaker problem studied in the previous section. Consider the linear BBM initial-boundary value problem (81) posed on the semi-infinite line with initial data 𝑢(𝑥, 0) = 0, 𝑥 ≥ 0 and the time-periodic sinusoidal boundary condition 𝑢(0, 𝑡) = -sin(𝜔 0 𝑡), 𝑡 ≥ 0, where 𝜔 0 > 0. The solution of the problem can be derived via the Fokas method given by Equation (85), where the contour of integration  is a small, closed contour around 𝑘 = 𝑖 as shown in Figure 2. We rewrite the solution (85) in terms of trigonometric functions 𝑢(𝑥, 𝑡) = 𝑢 1 (𝑥, 𝑡) + 𝑢 2 (𝑥, 𝑡) + 𝑒 -𝑥 sin(-𝜔 0 𝑡), (119a)

The integrand in Equation (119b) has the simple pole 𝑘 = 𝑖 inside the contour , so computing the corresponding residue gives 𝑢 1 (𝑥, 𝑡) = sin(𝜔 0 𝑡)𝑒 -𝑥 . Therefore, the solution of the linear BBM wavemaker problem is reduced to evaluating (119c) whose integrand possesses an essential singularity at 𝑘 = 𝑖. Then, the solution is simply 𝑢(𝑥, 𝑡) = 𝑢 2 (𝑥, 𝑡).

We appeal to the method of steepest descent for studying the long-time asymptotic behavior of the integral (119c). The poles of ℎ(𝑘) occur when 𝑘 2 -(1 + 𝑘 2 ) 2 𝜔 2 0 = 0 and are given by 𝑘 = 𝑘 𝑗 where

When 0 < 𝜔 0 < 1∕2, all the poles are real and this corresponds to the subcritical regime. When 𝜔 0 > 1∕2, the poles are complex and corresponds to the supercritical regime. We therefore define

as the critical frequency. Defining the phase 𝜙(𝑘) = 𝑖𝑘𝜉 -𝑖 𝑘 1+𝑘 2 , where 𝜉 = 𝑥∕𝑡, then the saddle points satisfy 𝜙 ′ (𝑘) = 𝑖𝜉 -𝑖 1-𝑘 2 (1+𝑘 2 ) 2 = 0 and hence are located at

If 0 < 𝜉 < 1, we have 𝜌 1 , 𝜌 2 ∈ ℝ with 𝜌 1 ∈ (0, 1), 𝜌 2 ∈ (-1, 0), and 𝜌 3 , 𝜌 4 are purely imaginary with Im(𝜌 3 ) > √ 3 and Im(𝜌 4 ) < -√ 3. If 𝜉 > 1, all four saddle points are purely imaginary and symmetrically distributed along the imaginary axis with Im(𝜌 1 ) ∈ (0, 1), Im(𝜌 2 ) ∈ (-1, 0), Im(𝜌 3 ) ∈ (1, √ 3), and Im(𝜌 4 ) ∈ (-√ 3, -1). Let 𝑘 = 𝑝 + 𝑖𝑞 so that 𝜙(𝑘) = (-𝜉𝑞 +

. The steepest descent paths occur when Im{𝜙(𝑘)} = Im{𝜙(𝜌 𝑗 )} is a constant and Re{𝜙(𝑘) -𝜙(𝜌 𝑗 )} < 0, 𝑗 = 1, 2, 3, 4.

Now we select the steepest descent contours  ′ passing through the saddles that avoid the essential singularities at 𝑘 = ±𝑖. Then, we deform  onto  ′ and consider the large 𝑡 behavior. Our general strategy will be to first consider the limit 𝑡 → ∞ and then show that the remaining contributions from avoiding the essential singularities are negligible. Depictions of the contours in Figure 6 will be discussed after splitting the range of 𝜉 into different intervals.

(a) Regime 0 < 𝜉 < 1.

The deformed heart-shaped contour (reminiscent of the Catalpa tree leaf) is shown in Figure 6(A), which passes through the saddles 𝜌 1 and 𝜌 2 . It is not necessary to deform onto 𝜖 → 0. Therefore, we can refine the deformed contour  ′ as

so that the solution 𝑢 in Equation ( 119) can be represented as

We have just proven that the third integral on the right-hand side vanishes in the limit 𝜖 → 0 for any 𝑡 > 0. We now consider the second integral on the right-hand side and prove that its dominant contributions come from the saddles on  1 ∪  2 and  3 ∪  4 . The first integral may give pole contributions, which we consider afterward. Now, we require 0 < 𝜖 < 1. On the contour  1 ∪  2 of Figure 6(A), the neighborhood of 𝜌 1 dominates. In a neighborhood of 𝜌 1 , the path  1 ∪  ⨙ is approximately parameterized by 𝑘 = 𝜌 1 ± (1 + 𝑖)𝑠, 𝑠 ∈ (-𝜖, 𝜖), so that by the Laplace method, we obtain

Similarly, on the contour  3 ∪  4 of Figure 6(A), we obtain the dominant contribution by integrating over a neighborhood of 𝜌 2 . Since 𝜌 2 1 = 𝜌 2 2 , Im(ϕ(ρ 1 )) = -Im(ϕ(ρ 2 )) and |𝜙 ′′ (𝜌 1 )| = |𝜙 ′′ (𝜌 2 )|, we sum up the saddle point contributions as

which is (𝑡 -1∕2 ). To determine how the poles 𝑘 𝑗 , 𝑗 = 1, 2, 3, 4, in Equation (120) enter the contour integral, we need to consider different cases for 𝜔 0 . (a.1) Subcritical regime: 0 < 𝜔 0 < 1∕2.

In this case, all the poles lie on the real line with 𝑘 1 ∈ (1, ∞), 𝑘 2 ∈ (0, 1), 𝑘 3 ∈ (-1, 0) and 𝑘 4 ∈ (-∞, -1). While it is easy to see that 𝑘 1,4 are outside of the heart-shaped contour  ′ , 𝑘 2,3 are inside  ′ when 𝑘 2 < 𝜌 1 and 𝑘 3 > 𝜌 2 . Observe that

When 𝜉 < 𝜔 ′ (𝑘 2 ), 𝑘 2,3 are in  ′ (see Region I in Figure 7). Applying Cauchy's residue theorem and using we obtain the asymptotic solution to the linear BBM sinusoidal initial-boundary value problem

where the saddle contribution in Equation ( 126) is the (𝑡 -1∕2 ) correction. We notice

corresponds to the unique root that is given by the D-N map and satisfies the radiation condition (4).

When 𝜉 > 𝜔 ′ (𝑘 2 ), 𝑘 2,3 are outside  ′ (see Region II in Figure 7) so that there are no pole contributions and the saddle contributions dominate. Using Equation (126), we obtain the leading-order asymptotic solution as 𝑡 → ∞ 𝑢(𝜉𝑡, 𝑡) ∼ cos

where Ξ ∶= √ 8𝜉 + 1. The solution (129) is (𝑡 -1∕2 ), indicating the typical algebraic decay in time of dispersive waves propagating faster than the group velocity. (a.2) Supercritical regime: 𝜔 0 > 1∕2.

The poles 𝑘 𝑗 , 𝑗 = 1, 2, 3, 4, are all complex and lie on the unit circle Re(𝑘 𝑗 ) 2 + Im(𝑘 𝑗 ) 2 = 1. Furthermore, 𝑘 2,4 are in the lower half plane and are not enclosed by  ′ , whereas 𝑘 1,3 are inside  ′ when Im𝜙(𝑘 1 ) < Im𝜙(𝜌 1 ). Setting Im𝜙(𝑘 1 ) = Im𝜙(𝜌 1 ) defines

14679590, 2024, 3, Downloaded from https://onlinelibrary.wiley.com/doi/10.1111/sapm.12665 by University Of Colorado Librari, Wiley Online Library on [28/04/2024]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License the curve 𝑙 3 (𝜉, 𝜔 0 ) ∶= 𝜉 2𝜔 0 -𝜔 0 + 2 √ 2𝜉(-2𝜉 + Ξ(𝜉) -1) 3∕2

(Ξ(𝜉) -1) 2 = 0, 𝜔 0 > 1∕2 and 0 < 𝜉 < 1,

in the 𝜔 0 -𝜉 plane. When 𝑙 3 (𝜉, 𝜔 0 ) > 0 (Region III in Figure 7), 𝑘 1,3 are inside the heartshaped contour  ′ and give the pole contributions

The pole 𝑘 1 is consistent with the singularity 𝑘 0 that lies on 𝜕𝐷 + in the D-N map and satisfies the radiation condition (4). Combining the pole and saddle contributions, we obtain the asymptotic solution

where 𝑢 𝑠1 is the contribution from the saddles given in (129). Similar to the long-time asymptotic solution of the linear KdV wavemaker problem (107), both terms in (131) are required for an asymptotic expansion that matches the result in the limit 𝜉𝑡 = (1), 𝑡 → ∞.

On the other hand, when 𝑙 3 (𝜉, 𝜔 0 ) < 0 (Region II in Figure 7), that is, no poles are inside the deformed contour  ′ , the neighborhoods of the saddles 𝜌 1,2 dominate and the leading-order solution is given by (129). The transition from regions I or III to region II in 𝜉 corresponds to the collision of the poles 𝑘 2,3 with saddle points 𝜌 1,2 and can be made uniform in 𝜉. 28 (b) Regime 𝜉 > 1.

Under the condition 𝜉 > 1, all the saddle points 𝜌 𝑗 are imaginary and the second derivative 𝜙 ′′ (𝜌 𝑗 ) at each saddle 𝜌 𝑗 , 𝑗 = 1, 2, 3, 4, is real. As shown in Figure 6(b), the curve  ′ =  1 ∪  2 passing through 𝜌 1 and 𝜌 3 has constant Im(𝜙(𝑘)) defined by

where 𝑘 = 𝑝 + 𝑖𝑞. The steepest descent direction follows the circular contour  1 ∪  ⨙ , and therefore, the neighborhood of the saddle point 𝜌 1 gives the dominant contribution for the integral. With 𝜙(𝜌 1 ), 𝜙 ′′ (𝜌 1 ) ∈ ℝ and 𝜙(𝜌 1 ) < 0, 𝜙 ′′ (𝜌 1 ) < 0, we apply the Laplace method and obtain

Now we deform the original small closed contour  onto the steepest descent curve  ′ ; in doing so, we incorporate the pole contributions from 𝑘 𝑗 in (120) when necessary. (b.1) Subcritical regime: 0 < 𝜔 0 < 1∕2.

All the poles (120) are on the real line and so are not enclosed in  ′ (see Region IV.a in Figure 7). This indicates that ∫  ℎ(𝑘)𝑒 𝑡𝜙(𝑘) 𝑑𝑘 = ∫  ′ ℎ(𝑘)𝑒 𝑡𝜙(𝑘) 𝑑𝑘. We use the saddle contribution from (133) to obtain the leading-order solution

where

The solution in this region has rapid decay in time of order (𝑡 -1∕2 𝑒 -𝛼(𝜉)𝑡 ). (b.2) Supercritical regime: 𝜔 0 > 1∕2.

The poles 𝑘 1,3 are in the first and second quadrants on a unit circle centered at the origin, so are enclosed in  ′ when Im(𝜙(𝑘 1 )) < Im(𝜙(𝜌 1 )) = 0 and, equivalently, Im(𝜙(𝑘 3 )) > Im(𝜙(𝜌 1 )) = 0. This implies that 𝑘 1,3 are inside  ′ when 𝜔 0 > 1 √ 2

, 1 < 𝜉 < 2𝜔 2 0 (see Region IV.b in Figure 7) and outside  ′ when

, 𝜉 > 2𝜔 2 0 (Region IV.a in Figure 7). The poles 𝑘 2,4 are in the third and fourth quadrants so are outside of the deformed contour  ′ .

For 𝑘 1,3 outside the contour in Region IV.a, the asymptotic solution is given by (134). In Region IV.b, the contour integral over  can be expressed as

as shown by the deformation of  to  ′ in Figure 7. Evaluating the residues, we obtain

Then, the large 𝑡 asymptotics are

where 𝑢 𝑠2 is the saddle contribution given in (134). Because the exponential decay rate |𝜙(𝜌 1 )| of 𝑢 𝑠2 is monotone increasing, in Region IV.b where 1 < 𝜉 < 2𝜔 2 0 , we have

. Then, the solution (137) is dominated by the saddle contribution 𝑢 𝑠2 . Therefore, for large 𝑡, the leading-order contribution for all 𝜔 0 and 𝜉 > 1 is the same and is given by 𝑢 𝑠2 .

The large 𝑡 solution given in Equations ( 128), ( 129), ( 131), (134), and (137) proves Theorem 3.

Remark. At a qualitative level, the large 𝑡 solutions of the linear BBM and KdV wavemaker problems exhibit similar bifurcation diagrams in Figures 5 and 7. Consequently, the solution behavior is qualitatively similar, with traveling periodic waves in Region I, algebraically decaying waves in Region II, and exponential decay in Regions III and IV of the 𝜔 0 -𝜉 plane.

Remark (Linear modulation theory). The modulation equations (113) are applicable for general linear dispersive partial differential equations. Following a path similar to that of the previous section, we now solve (113) with the associated boundary values (114) where, in this case, 𝜔 = 𝑘∕(1 + 𝑘 2 ). For 𝜉 = 𝑥∕𝑡, the self-similar solution is

where 𝑘 0 = 𝑘(𝜔 0 ) > 0 satisfies the linear dispersion relation 𝜔 0 = 𝑘 0 ∕(1 + 𝑘 2 0 ). There are two choices for 𝑘 0 . Admissibility of the simple wave solution in the subcritical regime 0 < 𝜔 0 < 1∕2 implies the radiation condition 𝑐 𝑔 (𝜔 0 ) > 0. Then, 𝑘(𝜔 0 ) is uniquely determined (𝑘 0 = 𝑘 2 in Equation ( 120)) and

Suitable integration of 𝜃 𝑥 = 𝑘(𝑥∕𝑡), 𝜃 𝑡 = -𝜔(𝑥∕𝑡) determines the phase

where 𝜃 0 ∈ ℝ. The amplitude 𝑎 has the solution (117) with 𝑐 𝑔 (𝜔 0 ) given in (139) and 𝐹(𝜉) an arbitrary, smooth function. As before, these results agree with the asymptotic solutions (128) and (129) if we take

where Ξ = √ 8𝜉 + 1.

COMPARISON WITH CONDUIT EXPERIMENTS

In the long-wavelength, small-amplitude regime, the linear BBM equation (equivalent to the linearized conduit equation 22 ) describes the interfacial wave dynamics of a particular type of viscous two-fluid core-annular flow 31,32 that models magma migration through the upper mantle 33 and channelized water flow through glaciers. 34 This flow system comprises a pressure-driven, buoyant, Stokes core fluid within an annulus of another, miscible, heavier Stokes fluid, with a small core-to-annular fluid viscosity ratio. At the injection site, a piston pump acts as a wavemaker by injecting a steady stream of core fluid with an added, small-amplitude, time-periodic injection to generate interfacial waves (see Ref. 22). The system can be realized in a simple laboratory setting and provides a quantitative experimental platform for the study of linear and nonlinear wave propagation. 22,31,35-42 The linear BBM equation was found to more accurately capture the dispersive wave properties of this viscous core-annular flow across a wider range of frequencies than the asymptotically equivalent long wave, linear KdV equation. 22 In this section, we investigate the linear time-periodic wavemaker problem in this viscous core-annular system and compare our experimental observations with the long-time asymptotic solutions of the linear BBM equation, derived in Section 4.

Consider the linearized BBM equation initial-boundary value problem

𝑢(𝑥, 0) = 0, 𝑥 > 0,

𝑢(0, 𝑡) = sin(-𝜔 0 𝑡), 𝑡 > 0.

The dependent variable 𝑢 represents small-amplitude fluctuations of the two-fluid interfacial cross-sectional area scaled by the steady, constant cross-sectional area of the established, vertical, straight conduit; 𝑥 is proportional to the vertical distance from the boundary, and 𝑡 is proportional to time. See Ref. 22 for the nondimensionalization. Here, we will use the dimensional length 𝐿 and time 𝑇 scalings as fitting parameters, which amounts to fitting the fluid viscosities and density difference, in order to compare with experiment. Past work has demonstrated quantitative agreement between the nonlinear extension of the linearized BBM equation (the conduit equation) and experiment with no fitting parameters. 37,39,41 A linear time-periodic sinusoidal time series of unit amplitude and frequency 𝜔 0 is generated at the boundary and propagates into the constant background. The experimental configuration is identical to that employed in Ref. 22. The experiments are conducted within a vertical square acrylic column that has the cross-section 5 cm × 5 cm and a height of 180 cm. The core fluid is composed of dyed, diluted glycerin that exhibits lower density and significantly less viscosity than the reservoir fluid into which it is injected. The core fluid rises buoyantly in an external pure glycerin reservoir to establish a circularly symmetric free interface between the two fluids. A computer-programmed flow rate pump is used to regulate the flow rate Q (cm 3 /s) of the core fluid, with the cross-sectional area à (cm 2 ) at the injection site satisfying the Hagen-Poiseulle pipe flow law

with a constant of proportionality involving the ratio of the fluid density difference to the interior viscosity. The imposition of the time-periodic boundary data in (142c) is realized by the injection flow rate

where Q0 (cm 3 /s) is the background flow rate, 𝑎 is the (nondimensional) amplitude, ω0 (rad/s) denotes the dimensional angular input frequency, and t (s) is time. A high-resolution camera with a 3 Hz sampling rate is used for data acquisition above the injection boundary to capture the evolving wavetrains. The conduit edges are extracted from camera images using edge detection and smoothing algorithms. At each vertical image pixel, the conduit diameter is calculated as the number of horizontal pixels between the two conduit edges. The diameter is then converted to the dimensionless conduit cross-sectional area 𝐴 by dividing by the number of horizontal pixels in the undisturbed conduit prior to commencing time-periodic injection and squaring the result. For experimental comparison with theoretical predictions, we appeal to the nondimensionalization parameters 𝐿 (cm) as the vertical length scale, 𝑈 (cm/s), as the vertical velocity scale and 𝑇 = 𝐿∕𝑈 (s) as the time scale. For each experiment presented below, approximately 15 sets of linear periodic traveling waves were first generated for calibration, and the dimensional angular frequency and wave number ( ω, k) were measured. The scales 𝐿 and 𝑈 were determined by fitting the dimensional quantities ( ω, k) to the BBM linear dispersion relation

The scales associated with each displayed data set are reported. We scale the dimensional space and time as 𝑥 = x∕𝐿, 𝑡 = t∕𝑇, respectively. Upon linearization of the cross-sectional area 𝐴 on the uniform background

and using the boundary condition (144) and the Hagen-Poiseulle law (143) as

where the dimensionless input frequency is 𝜔 0 = ω0 𝑇, we obtain the initial, boundary value problem (142a), which we compare with experiment. Previous viscous core-annular flow experiments examined the accuracy of the linear dispersion relation (145), the existence of the critical frequency 𝜔 𝑐𝑟 separating sub-and supercritical regimes, and the sinusoidal spatial profiles of linear waves, in addition to their nonlinear counterparts. 22 The BBM dispersion relation (145) was shown to be an accurate model of observed linear waves, with discrepancies due to a weak recirculating flow in the outer, annular fluid. In what follows, we compare new experiments with asymptotic predictions of (i) the spatial structure of waves along rays 𝑥∕𝑡 = constant in the sub-and supercritical regimes and (ii) the subcritical group velocity transition between Regions I and II in Figure 7.

The space-time contour plot in Figure 8(A) depicts the dimensionless cross-sectional area 𝐴 for the wavemaker problem in the subcritical regime 0 < 𝜔 0 < 𝜔 𝑐𝑟 = 1∕2. The color scale denotes wave amplitude with the red-blue region representing waves of substantial amplitude, whereas the lighter-colored, blue region corresponds to the small amplitude, noisy background. The method for measuring the parameters of the periodic traveling wave is described in Ref. 22. Measurement of this experimental data set yields the periodic wave amplitude, angular frequency, and wave number (𝑎, 𝜔, 𝑘) = (0.11 ± 0.02, 0.41 ± 0.02, 0.52 ± 0.02), respectively. Figure 8(B) displays example time-lapse images captured by a high-resolution camera during the experiment. At 𝑡 = 150 and 𝑡 = 300, the interfacial wave is observed as it enters the uniform conduit, resulting in decay at the leading edge. Subsequently, as the wave propagates, for example, at 𝑡 = 600, it attains local periodicity with small amplitude relative to the average cross-sectional diameter or area. In Figure 8(C), we compare the experimentally extracted traveling wave data from Figure 8(A) along rays 𝑥∕𝑡 = constant with the long-time asymptotic solutions of the linear BBM wavemaker problem. The comparison utilizes the transformation (146) and the solutions (128, 129) in Regions I and II of Figure 7, with 𝑎 = 0.11, 𝜔 0 = 0.41, and the previously fitted scaling parameters 𝐿 and 𝑇. However, the experiment identifies 𝑥 = 0 with the starting point of the cropped image and 𝑡 = 0 with the initiation of imaging, which differ slightly from the location of the injection site and the initiation of the periodic time series for the injection rate. Consequently, a shift is applied to 𝑥 and 𝑡 everywhere they appear in the solution formulas. For example,

𝜉 =

𝑥 -𝑥 0 𝑡 -𝑡 0 .

The experiment in Figure 8(C) is subject to the shifts 𝑥 0 = 6 (≈ 2 cm) and 𝑡 0 = 6 (≈ 2 s), which are selected by comparing the experimental data with the theoretical solutions at various 𝜉. The shift in 𝑥 is consistent with the experimental setup, in which case the injection site is, while not shown, on the left of the experimental images in Figure 8(B). The shift 𝑡 0 is the best fit time shift.

In Figure 8(C), good agreement is achieved for the comparison between the experiment and the long-time asymptotic solutions in Regions I (0 < 𝜉 < 𝑐 𝑔 ) and II (𝑐 𝑔 < 𝜉 < 1). Region I contains the solutions asymptoting to periodic traveling waves, and the solutions in Region II asymptote to algebraically time-decaying waves. According to the linear BBM theoretical prediction, the corresponding threshold between Regions I and II, that is, the group velocity, at 𝜔 0 = 0.41 is predicted to be 𝑐 𝑔 = 𝜔 ′ (𝑘 2 ) = 0.45 (cf. Equation ( 127)). For small 𝑡, the asymptotic solution deviates somewhat from experiment. However, this is expected because the formula represents the large 𝑡 asymptotics of the solution. In particular, the observed oscillation amplitude, period, and phase for sufficiently large 𝑡 are very close to experiment. Note that the small, noisy oscillations correspond to the resolution of our imaging system. In particular, we observe diameter fluctuations on the order of four pixels, which, for a 30 pixel diameter established conduit, corresponds to (4∕30) 2 ≈ 0.02 fluctuations in the nondimensional cross-sectional area 𝐴. The reason for this resolution is the 200:1 vertical to horizontal aspect ratio of the imaging system used to capture the data in Figure 8(A). As shown in Figure 8, wave propagation is periodic when 𝜉 is less than a threshold 𝜉 𝑔 and decays when 𝜉 > 𝜉 𝑔 . To obtain a value for this threshold from experiment, we employ the following method. Time slices are extracted from the data set for 𝑡 ∈ [30,150]. In each time slice, the rightmost location 𝑥 𝑔 (𝑡) is defined to be the last spatial peak with a value greater than 80% of the averaged periodic wave amplitude at the leading edge where the wave starts to decay. The fitted slope of 𝑥 𝑔 (𝑡), 𝜉 𝑔 , is a reliable approximation of the threshold between Regions I and II. This measurement method is applied to thirteen distinct data sets of linear waves generated from varied input frequencies in the subcritical regime. Figure 9 compares the measurement of 𝜉 𝑔 with the group velocity 𝑐 𝑔 (𝜔) (139) at different input angular frequencies 𝜔. The angular frequencies are measured and averaged over the periodic wave region. The error bars presented in Figure 9 represent the uncertainties resulting from the linear fit of 𝑥 𝑔 (𝑡) and the variation in the measured wave frequency. Figure 9 verifies the threshold 𝜉 𝑔 as the linear BBM group velocity prediction experimentally.

Figure 10 depicts an experiment in the supercritical regime with an input frequency greater than the critical value, that is, 𝜔 0 > 1∕2 (see Ref. 22 for an experimental study of the critical value). The wave appears to be time-periodic near the boundary and rapidly decays as it propagates. Near

CONCLUSIONS

This work presents the long-time, asymptotic analysis of initial-boundary value problems for linear, scalar evolution equations on the half-line with time-periodic boundary conditions. We obtain the D-N map for a general third-order model with Dirichlet boundary data that are sufficiently smooth and asymptotically periodic in time. The uniqueness of the D-N map is established (Theorem 1) if and only if the well-known radiation condition in physics (Corollary 1) holds. By employing the D-N map, we obtain the leading-order asymptotic solution for fixed 𝑥 and 𝑡 → ∞. These results hinge upon the assumption of asymptotic time periodicity of the solution.

This assumption is proven in the specific case of the wavemaker problems for the linear KdV and linear BBM equations with zero initial condition and sinusoidal boundary condition with frequency 𝜔 0 (Theorems 2 and 3). Using the Fokas method and steepest descent analysis, we obtain the leading-order, asymptotic solution as 𝑡 → ∞ that is uniform in 𝑥 ≥ 0. Crossing the interval of frequencies (-𝜔 𝑐𝑟 , 𝜔 𝑐𝑟 ) identifies a bifurcation in solution behavior. When |𝜔 0 | < 𝜔 𝑐𝑟 , the solution is asymptotically the propagating wave 𝑢(𝑥, 𝑡) ∼ sin(𝑘 0 𝑥 -𝜔 0 𝑡), 𝑡 → ∞ for 𝑥 ∈ [0, 𝑐 𝑔 (𝜔 0 )𝑡) where 𝑐 𝑔 (𝜔 0 ) = 𝜔 ′ (𝑘 0 ) is the group velocity and 𝑘 0 = 𝑘 0 (𝜔 0 ) is the unique branch of the dispersion relation 𝜔 0 = 𝜔(𝑘 0 ) such that 𝑐 𝑔 (𝜔 0 ) > 0 and Im 𝑘 0 = 0. This choice of 𝑘 0 corresponds to the radiation condition. When |𝜔 0 | > 𝜔 𝑐𝑟 , the asymptotic solution is the time-periodic, spatially decaying, causal solution 𝑢(𝑥, 𝑡) ∼ exp(-Im(𝑘 0 )𝑥) sin(Re(𝑘 0 )𝑥 -𝜔 0 𝑡) for 𝑥 = (1), 𝑡 → ∞ uniformly in 𝑥 like 𝑂( 1 𝑘 4 )), we infer lim 𝑥→0 + 𝑢(𝑥, 𝑡) = 𝜔 0 2𝜋 ∫ 𝜕 D+ 𝑒 -𝑖(𝑘-𝑘 3 )𝑡 1 -3𝑘 2 (𝑘 -𝑘 3 ) 2

[ 𝜔 0 g0 (𝑘 -𝑘 3 , 𝑇) + 1 ] 𝑑𝑘.

But note that, via integration by parts, 𝜔 0 g0 (𝑘 -𝑘 3 , 𝑇) + 1 = 𝑒 𝑖(𝑘-𝑘 3 )𝑇 cos(𝜔 0 𝑇) -𝑖(𝑘 -𝑘 3 ) 𝜔 0 𝑒 𝑖(𝑘-𝑘 3 )𝑇 sin(𝜔 0 𝑇) + (𝑘 -𝑘 3 ) 2 𝜔 0 g0 (𝑘 -𝑘 3 , 𝑇).

Therefore, employing (A.6) (which is valid for 𝑥 = 0 and 𝑛 = 2), we find (after deforming 𝜕 D+ back to 𝜕𝐷 + once the singular terms are gone) lim 𝑥→0 + 𝑢(𝑥, 𝑡) = 1 2𝜋 ∫ 𝜕𝐷 + 𝑒 -𝑖(𝑘-𝑘 3 )𝑡 (1 -3𝑘 2 ) g0 (𝑘 -𝑘 3 , 𝑇)𝑑𝑘 = 1 2𝜋 ∫ 𝜕𝐷 + 𝑒 -𝑖(𝑘-𝑘 3 )𝑡 (1 -3𝑘 2 ) {𝐺 0 }(-(𝑘 -𝑘 3 ))𝑑𝑘,

where  {𝐺 0 }(𝑘) is the (temporal) Fourier transform on the whole line of the piecewise-defined function 𝐺 0 (𝑡) equal to 𝑔 0 (𝑡) on (0, 𝑇) and to 0 on [0, 𝑇] 𝑐 . Making the change of variable 𝜆 = -(𝑘 -𝑘 3 ), which maps 𝜕𝐷 + to -ℝ, and then using the Fourier inversion theorem, we obtain lim 𝑥→0 + 𝑢(𝑥, 𝑡) = 1 2𝜋 ∫ ℝ 𝑒 𝑖𝜆𝑡  {𝐺 0 }(𝜆)𝑑𝜆 = 𝐺 0 (𝑡) = 𝑔 0 (𝑡), 0 < 𝑡 < 𝑇, which is the prescribed boundary condition at 𝑥 = 0. Examination of the formula (A.4) for either Im𝑘 > 0 and 𝑘 ∈ 𝜕𝐷 + or 𝑘 ∈ ℝ and 𝑘 ∈ 𝜕𝐷 + shows that the limit as 𝑥 → ∞ of 𝑢(𝑥, 𝑡) is zero for each fixed 𝑡 > 0 by either exponential decay or the Riemann-Lebesgue lemma, respectively. We can pass differentiation by 𝑥 through the integral in (A.4) owing to the dominated convergence theorem and thereby obtain lim 𝑥→∞ 𝑢 𝑥 (𝑥, 𝑡) = 0 by the same argument. This step concludes our proof that the Fokas method formula (A.1) satisfies the linear KdV wavemaker problem (A.2).

Having shown that formula (A.1) satisfies problem (A.2), it remains to show that problem (A.2) admits at most one solution, so that (A.1) is the unique solution to (A.2). Uniqueness follows via a standard energy estimate argument. Indeed, suppose that 𝑢 and 𝑣 both satisfy problem (A.2). Then, their difference 𝑤 ∶= 𝑢 -𝑣 satisfies the homogeneous problem 𝑤 𝑡 + 𝑤 𝑥 + 𝑤 𝑥𝑥𝑥 = 0, 0 < 𝑥 < ∞, 𝑡 > 0, 𝑤(𝑥, 0) = 𝑤(0, 𝑡) = 0.

(A.9) Multiplying the linear KdV equation in (A.9) by 𝑤 and integrating with respect to 𝑥, we have

Thus, in view of the boundary condition satisfied by 𝑤 as well as of the assumption that the solution and its derivatives vanish at infinity, we infer 𝑑 𝑑𝑡 ∫ ∞ 0 𝑤(𝑥, 𝑡) 2 𝑑𝑥 = -𝑤 𝑥 (0, 𝑡) 2 .

14679590, 2024, 3, Downloaded from https://onlinelibrary.wiley.com/doi/10.1111/sapm.12665 by University Of Colorado Librari, Wiley Online Library on [28/04/2024]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License

In the context of real solutions, which is the one relevant to the KdV equation, it follows that the 𝐿 2 𝑥 (0, ∞)-norm of 𝑤 is nonincreasing in 𝑡. Hence, in particular, using also the initial condition for 𝑤, we obtain ∫ ∞ 0 𝑤(𝑥, 𝑡) 2 𝑑𝑥 ≤ ∫ ∞ 0 𝑤(𝑥, 0) 2 𝑑𝑥 = 0 which implies that 𝑤 ≡ 0. Equivalently, 𝑢 ≡ 𝑣 and problem (A.2) admits at most one solution. Since we have already shown that formula (A.1) satisfies (A.2), we conclude that (A.1) provides the unique solution to (A.2).