# Achieving Optimal Locomotion using Self-Generated Waves

Daire O'donovan, Miguel Bustamante, Olivier Devauchelle, Graham P Benham

► To cite this version:

Daire O'donovan, Miguel Bustamante, Olivier Devauchelle, Graham P Benham. Achieving Optimal Locomotion using Self-Generated Waves. Journal of Fluid Mechanics, 2026, 1029, pp.A4. 10.1017/jfm.2026.11133. hal-05502420

HAL Id: hal-05502420

https://hal.science/hal-05502420v1

Submitted on 10 Feb 2026

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# Achieving Optimal Locomotion using Self-Generated Waves

D O'Donovan¹†, MD Bustamante¹, O Devauchelle² and GP Benham¹

¹School of Mathematics and Statistics, University College Dublin, Belfield, Dublin 4, Ireland

²Université de Paris, Institut de Physique du Globe de Paris, CNRS, F-75005 Paris, France

(Received xx; revised xx; accepted xx)

An oscillating body floating at the water surface produces a wave-field of self-generated waves. When the oscillation induces a difference in fore-aft wave amplitude squared, these self-generated waves can be used as a mechanism to propel the body horizontally across the surface (Longuet-Higgins & Stewart 1964). The optimisation of this wave-driven propulsion (WDP) is the interest of this work. To study the conditions necessary to produce optimal thrust we will utilise a shallow water set-up where a periodically oscillating pressure source acts as the body. In this framework, an expression for the thrust is derived by relation to the aforementioned difference in fore-aft amplitude squared. The conditions on the source for maximal thrust are explored both analytically and numerically in two optimal control problems. The first case is where a bound is imposed on the norm of the control function to regularise it. Secondly, a more physically motivated case is studied where the power injected by the source is bounded. The body is permitted to have a drift velocity *U*. When scaled with the wave speed *c*, the dimensionless velocity *v* = *U*/*c* divides the study into subcritical, critical and supercritical regimes and the optimal conditions are presented for each. The result in the bounded power case is then used to demonstrate how the modulation of power injected can slowly change the cruising velocity from rest to supercritical velocities.

## 1. Introduction

Horizontal propulsion along the water surface is seen in the natural world with examples including various insects (e.g water striders and whirligig beetles) (Bush & Hu 2006) as well as ducks (Yuan et al. 2021). Motivated by a curiosity to understand the underpinning fluid mechanics describing the actions of these creatures, various studies have been conducted such as high speed imaging of insects like the water strider (Hu et al. 2003; Steinmann et al. 2021) and the stricken honeybee (Roh & Gharib 2019). Inspired by these creatures, robots that mimic their movements have been developed, with *Robostrider* (Hu et al. 2003) in the case of the water strider and *SurferBot* in the case of the honeybee (Rhee et al. 2022). In a similar vein to such movements we can achieve horizontal locomotion by riding waves produced by gunwale bobbing on a canoe (Benham et al. 2022).

The water strider was first thought to use Wave-Driven Propulsion (WDP), but the observation of subsurface vortex generation (Hu et al. 2003) shifted our understanding to the use of WDP in combination with vortex-driven propulsion (VDP) (Bühler 2007; Steinmann et al. 2018; Gao & Feng 2011; Steinmann et al. 2021). We note that although multiple inertial propulsion mechanisms can be used including vortices (VDP), waves (WDP) or a combination of both (as in the water strider case), in this study we will focus on optimising thrust using WDP only.

In the case of the honeybee trapped on the water surface, it beats its wings to create an asymmetry in the fore-aft surface wave-field which is indicative of the use of radiation stresses described by Longuet-Higgins & Stewart (1964) where the former later applied these ideas

<footnote>† Email address for correspondence: daire.odonovan@ucdconnect.ie</footnote>
to propel a small craft (Longuet-Higgins 1977). Radiation stress is defined as proportional to the difference in the square of the fore-aft amplitudes generating momentum and hence thrust. Mathematically in the shallow limit the thrust force is (Longuet-Higgins & Stewart 1964)

$$F_T = \frac{3}{4} \rho g [h^2]_+^+ \quad (1.1)$$

where h is the wave amplitude, ρ is density, g is acceleration due to gravity and +,− denotes the fore and aft waves respectively. Radiation stress can be used to find the thrust generated by waves to propel bodies on a vibrating surface (Ho et al. 2023; Pucci et al. 2015; Barotta et al. 2023) as well as a body producing its own waves in WDP (Longuet-Higgins 1977; Benham et al. 2024). Motivated by the latter, we will characterise the thrust due to WDP by starting from the wave equation and recovering the difference in fore-aft amplitude squared. To reduce complexity, we will work in the shallow water approximation and focus only on surface waves ignoring any subsurface effects. We also need to characterise the body which we will be studying. Rather than a rigid body on the surface (Benham et al. 2024), we will keep things general and consider a pressure region influencing the surface similar to the characterisation of an air-cushion vehicle (Doctors & Sharma 1972).

Once we have built the mathematical problem, we can look at posing an optimal control problem to maximise the thrust. Optimisation studies are motivated not only by understanding the natural world with examples including ducklings wave-riding in their mother's wake (Yuan et al. 2021) but also for industrial applications. Optimising WDP can lead to improved fuel efficiency of boats by extracting propulsive energy from the waves (Bøckmann et al. 2018). Meanwhile in sports, optimisation of WDP can be used to gain an advantage by using opponents' waves to improve propulsion (Tuck & Lazauskas 1998; Dode et al. 2022). This paper will consider an optimal control problem from a different perspective compared to the examples in Yuan et al. (2021) and Doctors (1997). Rather than minimising wave resistance, we are interested in maximising the thrust due to WDP. The only constraints on the problem are the shallow water setup and appropriate bounds on the pressure source.

The optimal control problem is divided into multiple cases according to the value of the drift velocity U. The case where the drift velocity U is zero will be studied in Sections 2 and 3. This will simplify the expressions to allow us to focus on the concepts such as deriving the thrust and the power. Following the completion of the start-up case, Section 4 will be concerned with the case where the drift velocity is less than the wave speed (subcritical). This section will introduce the application of a Galilean transformation to operate in the rest frame of the body. The case where the body moves at the same velocity as the waves emitted (critical) will follow in Section 5. The intuition will be useful when considering the case where the body is travelling faster than the wave speed (supercritical) in Section 6 requiring us to think carefully about the left and right wave since they are both left-travelling after the transformation. This will complete the template of velocity regimes from which we can study the injection of power necessary to accelerate from rest to supercritical velocities in Section 7, given that the acceleration is small enough to maintain the periodic assumptions that are made. Finally, Section 8 discusses the scope and possible improvements of this work.

## 2. 1D shallow water waves due to a source

Given a single body of length L propelling itself forward using WDP, we wish to study the conditions to produce maximal forward thrust. We will consider a surface region where there is a pressure disturbance induced by the motion of the body. To begin, we will consider zero horizontal velocity U of the body. This will be referred to as the start-up phase where a body is beginning to propel itself forward from rest. Following this, the horizontal velocity will be
<figure><img src="image_2.png" /><figcaption>Figure 1: Diagram of the set-up where there is a pressure source $P(x, t)$, over a region of length $L$, acting on the surface of the shallow water, with depth $H$. This induces waves on the fluid surface, denoted by $h(x, t)$, resulting in locomotion with a drift velocity $U$.</figcaption></figure>

introduced in three different regimes: subcritical velocities ($0 < v < 1$), critical velocities ($v = 1$) and finally supercritical velocities ($v > 1$), where $v = U/c$ is the drift velocity divided by the wave speed $c$.

To add a simplification which allows us to focus on the surface of the water, we will operate within the shallow water approximation and a one-dimensional set-up as shown in figure 1. A pressure source $P(x, t)$ is applied at the surface of an otherwise undisturbed body of shallow water where the fluid is in the region $0 \le z \le H$ in the steady state. The source induces a perturbation so the surface of the water is at $z = H + h(x, t)$, where $h(x, t) \ll H$. The shape of the surface perturbed by the pressure is encapsulated in $h(x, t)$. Unlike the honeybee and the water strider which produce gravity-capillary waves, here we focus on gravity waves and ignore the effect of surface tension.

As described in Appendix A.1, the perturbed wave field $z = h(x, t)$ satisfies the one-dimensional wave equation,

$$ \frac{\partial^2 h}{\partial x^2} - \frac{1}{c^2} \frac{\partial^2 h}{\partial t^2} = -Q(x, t), \quad -\infty < x < \infty, \qquad (2.1) $$

where $c = \sqrt{gH}$ is the wave speed and $Q(x, t) = P_{xx}(x, t)/\rho g$ is the pressure source term. The source is defined to be non-zero in the region $x \in [-\frac{L}{2}, \frac{L}{2}]$ and zero outside this. Since this investigation is concerned with a single body undergoing periodic oscillations, with a frequency $\omega$, the variables will be rescaled as follows:

$$ h, x \sim \frac{c}{\omega}, \qquad t \sim \frac{1}{\omega}, \qquad Q \sim \frac{\omega}{c}. \qquad (2.2) $$

After considering these scalings, we arrive at the dimensionless wave equation with a source,

$$ \frac{\partial^2 h}{\partial x^2} - \frac{\partial^2 h}{\partial t^2} = -Q(x, t), \quad -\infty < x < \infty. \qquad (2.3) $$

In dimensionless form, the source $Q$ is defined to be non-zero on the interval $[-\frac{l}{2}, \frac{l}{2}]$, where $l$ is the dimensionless body length $l = L\omega/c$. A periodically oscillating source and wave field are assumed, since there is a single body and the surface is otherwise undisturbed. Therefore, the source and wave field can be decomposed as follows:

$$ h(x, t) = \mathrm{Re} \left[ \hat{h}(x)e^{it} \right], \qquad Q(x, t) = \mathrm{Re} \left[ \hat{Q}(x)e^{it} \right]. \qquad (2.4) $$

The above expressions in (2.4) are used to reduce the PDE for $h(x, t)$ in (2.3) to an ODE for the complex function $\hat{h}(x)$,

$$ \hat{h}''(x) + \hat{h}(x) = -\hat{Q}(x). \qquad (2.5) $$
The ODE has a general solution that is calculated using Green's functions,

$$ \hat{h}(x) = Ae^{ix} + Be^{-ix} + \frac{i}{2} \int_{-\frac{l}{2}}^{x} \hat{Q}(X) \left( e^{i(x-X)} - e^{i(X-x)} \right) dX, \quad (2.6) $$

where $A, B \in \mathbb{C}$ are constants to be found. We note that although there is a dependence on $x$ in the integral, for $x > \frac{l}{2}$ the upper bound will be $\frac{l}{2}$ since $\hat{Q}(X > \frac{l}{2}) = 0$. Likewise for $x < -\frac{l}{2}$ the integral is zero since $\hat{Q}(X < -\frac{l}{2}) = 0$. The following Sommerfeld boundary conditions will be used, where $k_{\pm}$ is the wavenumber:

$$ \hat{h}' \left( \pm \frac{l}{2} \right) = \mp ik_{\pm} \hat{h} \left( \pm \frac{l}{2} \right), \quad (2.7) $$

and $k_{\pm} = 1$ for the start-up case. This is equivalent to imposing that waves only radiate away from the source. Through application of the boundary conditions, $A$ and $B$ are found to be

$$ A = -\frac{i}{2} \int_{-\frac{l}{2}}^{\frac{l}{2}} \hat{Q}(X) e^{-iX} dX, \quad B = 0, \quad (2.8) $$

and the resulting solution for the wave-field is

$$ \hat{h}(x) = -\frac{i}{2} \int_{x}^{\frac{l}{2}} \hat{Q}(X) e^{i(x-X)} dX - \frac{i}{2} \int_{-\frac{l}{2}}^{x} \hat{Q}(X) e^{i(X-x)} dX. \quad (2.9) $$

This solution is therefore composed only of a left-travelling wave for $x < -\frac{l}{2}$ and a right travelling wave for $x > \frac{l}{2}$.

The problem we aim to solve is to find a $\hat{Q}$ such that it gives maximal thrust. Therefore, an expression for the thrust is required. The thrust force is not immediately obvious to define but can be derived through manipulation of the wave equation (2.3) to recover the “radiation stress” as described by Longuet-Higgins & Stewart (1964), i.e. an asymmetrical wave-field resulting in a difference in fore-aft amplitude that generates thrust. Multiplying (2.3) by $h_x$ and taking an integral with respect to $x$ (see Appendix A.2), the radiation stress result can be recovered where the difference in fore-aft amplitude squared ($|\hat{h}|^2$) produces a thrust. As such we derive the following equivalence:

$$ -\frac{1}{2} [k |\hat{h}|^2]_{-}^{+} = \langle \hat{Q}, \hat{h}' \rangle, \quad (2.10) $$

where the + and - symbols in (2.10) represent values for $x$ that are on the right and left hand side of the $[-\frac{l}{2}, \frac{l}{2}]$ region respectively and $\langle \cdot, \cdot \rangle$ is the time-averaged inner product defined as

$$ \langle f, g \rangle := \frac{1}{T} \int_{0}^{T} \int_{-\frac{l}{2}}^{\frac{l}{2}} \mathrm{Re} [f(x) e^{it}] \mathrm{Re} [g(x) e^{it}] \, dx \, dt = \frac{1}{4} \int_{-\frac{l}{2}}^{\frac{l}{2}} f g^{*} + f^{*} g \, dx. \quad (2.11) $$

The left hand side of (2.10) has included differing fore-aft wavenumbers $k_{\pm}$, but in the start-up case $k_{+} = k_{-} = 1$. This will be more important when velocity is introduced to apply a correction due to the Doppler shift introduced through operating in the moving frame. Relating back to the radiation stress, we interpret (2.10) as an equivalence relation between the force injected and the resulting asymmetrical wave-field that induces thrust. It is noted that there is a different prefactor compared to (1.1) which can easily be fixed through scaling the equations by 3/2, but this will not be done since it has no effect on the work to follow. Hence, motivated by (2.10) the time-averaged thrust is taken to be

$$ \bar{F}_T(\hat{Q}(x)) := \int_{-\frac{l}{2}}^{\frac{l}{2}} \overline{Qh_x} dx = \langle \hat{Q}, \hat{h}' \rangle. \quad (2.12) $$
Some constraints will be required in the optimal control problem. One case that will be investigated is when the time-averaged power is bounded. To derive an expression for the power, we begin with the energy of the wave as the sum of the dimensionless kinetic and potential energy integrated over the whole domain

$$E(t) = \int_{-\frac{L}{2}}^{\frac{L}{2}} \left( \frac{1}{2} h_t^2 + \frac{1}{2} h_x^2 \right) dx. \qquad (2.13)$$

Power is given as Pow = $\frac{dE}{dt}$. Since we are working with periodic oscillations, the time-averaged power is $\overline{Pow} = E(T) - E(0) = 0$, which means the energy is conserved over a period. The right hand side of (2.13) can be differentiated with respect to time and (2.3) is substituted and then time-averaged to generate another equivalence relationship (see Appendix A.3):

$$\int_{-\frac{L}{2}}^{\frac{L}{2}} \overline{Qh_t} dx = - \left[ \overline{h_t h_x} \right]_{-}^{+}. \qquad (2.14)$$

The average power injected by the pressure source on the left-hand side of (2.14) is related by equality to the average power radiated at the boundaries of the pressure source on the right hand side of (2.14). Since this is an idealised system where no power inputted is lost due to external effects, this is expected. Motivated by (2.14) we take the time-averaged power injected to be

$$\overline{Pow} := \int_{-\frac{L}{2}}^{\frac{L}{2}} \overline{Qh_t} dx = \langle \hat{Q}, i\hat{h} \rangle. \qquad (2.15)$$

Furthermore, we can substitute (2.4) into the right hand side of (2.14) to show

$$\frac{1}{2} [k_+ |\hat{h}_+|^2 + k_- |\hat{h}_-|^2] = \langle \hat{Q}, i\hat{h} \rangle. \qquad (2.16)$$

As such, while the thrust is the difference in fore-aft amplitude squared, the power is the sum.

Now, we can proceed with posing the optimal control problem. In the work so far, the body is taken to be at rest. The thrust would introduce movement itself but, for now, we will restrict the problem to that of “start-up,” with the horizontal drift velocity close to zero, $U \approx 0$.

### 3. Optimal locomotion on start-up

We wish to pose an optimal control problem to maximise the thrust as it is defined in (2.12). The control function $\hat{Q}(x)$ is in the expression for $\hat{h}(x)$ in (2.9) which implicitly enforces that the wave ODE (2.5) and its boundary conditions (2.7) must be satisfied. Some sort of regularisation to the problem is required. In this paper we will use two different constraints. Firstly we will bound the time-averaged norm of the control,

$$\langle \hat{Q}, \hat{Q} \rangle \le \varepsilon, \qquad (3.1)$$

for some dimensionless $\varepsilon > 0$†. Secondly, we will bound the time-averaged power, as it is seen in (2.15),

$$\langle \hat{Q}, i\hat{h} \rangle \le \delta, \qquad (3.2)$$

for some dimensionless $\delta > 0$. It is worth noting that, neither constraint requires that $\hat{Q}$ be continuous at the boundary and jumps may be present at $x = \pm \frac{L}{2}$.

Along with an analytical approach, a numerical optimisation approach is taken to validate the results using the Ipopt solver that is part of the JuMP package in Julia (Lubin et al. 2023). To pose the problem numerically, the inputs are divided into variables, objective and constraints.

<footnote>† We define $\varepsilon = \hat{\varepsilon}c/\omega$, $\delta = \delta\omega/\rho g H^2 c$ where $\hat{\varepsilon}, \delta$ are dimensional bounds.</footnote>
There are four arrays of variables to represent the real and imaginary parts of $\hat{Q}$ and $\hat{h}$ which are discretised in space. The thrust (2.12), discretised using the trapezium rule, is the objective function we wish to maximise. As for the constraints, these are the boundary conditions (2.7), the discretised ODE for $\hat{h}(x)$ (2.5), and either the bounded norm of $\hat{Q}$, or the bounded power. In the interest of comparing the analytical and numerical work, we take $\varepsilon, \delta = 1$. These will be dealt separately in the subsections that follow.

### 3.1. Case 1 - Bounded Norm

To be physically plausible, the source $Q$ can be regularised by bounding the $L_2$ norm which will provide us with well behaved solutions for $\hat{Q}$. In (3.1) the value $\varepsilon = 1$ is taken to compare the analytical and numerical approaches. The variational calculus approach in Appendix B.1 results in the following condition to be an optimal solution,

$$ \mathcal{L}_0 \hat{Q} = 0, \qquad (3.3) $$

where the (1D Helmholtz) differential operator $\mathcal{L}_0 = \partial_{xx} + \mathbb{I}$, applied to $\hat{Q}$ implies wave-like solutions,

$$ \hat{Q}(x) = D_L e^{ix} + D_R e^{-ix}, \qquad (3.4) $$

where $D_L, D_R \in \mathbb{C}$. In Appendix B.1, we go a step further and find the constants $D_L$ and $D_R$ in terms of $\lambda$, the Lagrange multiplier, which is also found. The constraint on the modulus squared of the constants is

$$ |D_L|^2 = 4 \left( l + \frac{l \sin^2(l)}{(2\lambda - l)^2} + \frac{2 \sin^2(l)}{2\lambda - l} \right)^{-1}, \qquad |D_R|^2 = |D_L|^2 \frac{\sin^2(l)}{(2\lambda - l)^2}, \qquad (3.5) $$

where $\lambda$ is

$$ \lambda_{\pm} = \pm \frac{\sqrt{l^2 - \sin^2(l)}}{2}. \qquad (3.6) $$

The $(\pm)$ in the expression for $\lambda$ corresponds to the choice of thrust direction. To move in the positive direction (right) corresponds to taking the negative part in $\lambda$, thus maximising positive thrust. Since $\lambda$ is the Lagrange multiplier, $\lambda \in \mathbb{R}$ is required. We have the condition

$$ l^2 - \sin^2(l) > 0, \qquad (3.7) $$

which is true for all real values of $l$, so $\lambda$ is always real. Hence, in (B 9) we note $D_L$ and $D_R$ are related via multiplication by a real number meaning the two waves are in phase.

A periodic source $\hat{Q}$ that is composed of a left travelling and/or right travelling waves of this form will produce optimum thrust. There are infinitely many such solutions which differ by an arbitrary choice of phase Arg($D$). An example of one such $\hat{Q}(x)$, found numerically, is displayed in figure 2 along with the resulting wave-field $\hat{h}(x)$. In figure 2b we see a larger wave on the left hand side which corresponds to positive thrust and hence start-up locomotion in the right-travelling direction is achieved.

The expressions for $|D_L|^2$ and $|D_R|^2$ in (3.5) are validated numerically for a given $l$. Using (B 15) allows us to use the right and left waves to relate $|D_L|^2, |D_R|^2$ to the information we have on $\hat{Q}$ from the numerical result which is one of infinite solutions that satisfy

$$ \left| \frac{1}{2} \int_{-\frac{l}{2}}^{\frac{l}{2}} \hat{Q}(X) e^{-iX} dX \right|^2 = \lambda^2 |D_L|^2, \qquad \left| \frac{1}{2} \int_{-\frac{l}{2}}^{\frac{l}{2}} \hat{Q}(X) e^{iX} dX \right|^2 = \lambda^2 |D_R|^2. \qquad (3.8) $$

We note that $\lambda^2|D_L|^2$ and $\lambda^2|D_R|^2$ are related to the amplitude squared of the left and right travelling waves ($|\hat{h}_L|^2, |\hat{h}_R|^2$) respectively. For $l = 3\pi/2$, the relative error between the analytical
<figure><img src="image_3.png" /><figcaption>Figure 2: (a) An example of a source, $\hat{Q}(x)$, that results in the optimal thrust under the bounded norm constraint where $v = 0$ and $l = 3\pi/2$. (b) The corresponding wave-field, $\hat{h}(x)$ under the same conditions.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>Figure 3: Plot of the time-averaged thrust $\bar{F}_T$ vs the dimensionless length scale $l$ resulting from the bounded norm optimisation on start-up $v \approx 0$. It can also be seen that as $l$ increases, the dimensional thrust tends to scale with $l^2$ which is the same behaviour found in (3.9) for $l \gg 1$</figcaption></figure>

and numerical expressions was found to be 0.003% for $\lambda|D_L|^2$ and 0.014% for $\lambda^2|D_R|^2$. Given the aforementioned arbitrary choice of phase, the numerical scheme will converge to one of an infinite set of solutions. If we use (B 15), we find the numerical values for $D_L$ and $D_R$ and use them in conjunction with (3.4) to demonstrate correspondence between the numerical and analytical results in figure 2.

The dependence of the expressions on the length scale $l$, which is related to $\omega$, motivates that we investigate how the resulting thrust depends on $l$. The thrust can be written (see Appendix B.1) in terms of the constants that have been found,

$$ \bar{F}_T = \frac{1}{4} (2\lambda^2|D_L|^2 - 2\lambda^2|D_R|^2). \qquad (3.9) $$

Plotting the force that results from the numerical optimisation, we see that it increases with $l$ (figure 3), which is a property we will return to in Section 4.

## 3.2. Case 2 - Bounded power

Rather than bounding the norm of $Q$, let us instead bound the power (2.15). A body on the water propelling itself through vertical oscillations has control over the power it inputs to generate thrust. Mathematically, the power is bounded by $\overline{Pow} \le \delta$, where $\delta = 1$ is chosen for the following work. The variational calculus approach in Appendix B.2 is taken and the optimal condition is a
first order homogeneous ODE,

$$ (\partial_x + i\lambda) \left( -\frac{i}{2} \int_{-\frac{1}{2}}^{\frac{1}{2}} \hat{Q}(X) e^{i(X-x)} \, dX - \frac{i}{2} \int_{-\frac{1}{2}}^{\frac{1}{2}} \hat{Q}(X) e^{i(x-X)} \, dX \right) = 0, \quad (3.10) $$

which results in two solutions depending on the eigenvalue, $\lambda \in \mathbb{R}$, which again corresponds to the Lagrange multiplier in the variational calculus. The two solutions for $\lambda$ will be denoted as $\lambda_1$ and $\lambda_2$ and the corresponding solutions to the ODE (3.10) are given as follows:

$$ C_L e^{ix} + C_R e^{-ix} = C e^{ix} \quad \text{for } \lambda_1 = -1, \quad (3.11) $$

$$ C_L e^{ix} + C_R e^{-ix} = C e^{-ix} \quad \text{for } \lambda_2 = 1, \quad (3.12) $$

where $C_L, C_R \in \mathbb{C}$. This gives us the option of a purely left or purely right travelling wave and to illuminate this fact, the left and right wave amplitudes have been labelled $C_L$ and $C_R$ respectively. The constants are found by equating coefficients to be

$$ C_L = -\frac{i}{2} \int_{-\frac{1}{2}}^{\frac{1}{2}} \hat{Q}(X) e^{-iX} \, dX, \qquad C_R = -\frac{i}{2} \int_{-\frac{1}{2}}^{\frac{1}{2}} \hat{Q}(X) e^{iX} \, dX. \quad (3.13) $$

Similar to the bounded norm case, we can note that $C_L$ and $C_R$ correspond the fore-aft wave heights respectively. Before making a choice of $\lambda$ and discussing that solution, we recall that the bounded power (2.15), when written in terms of the above constants, gives a condition on $|C_L|^2$ and $|C_R|^2$,

$$ \frac{1}{4} (2 |C_L|^2 + 2 |C_R|^2) = 1. \quad (3.14) $$

To achieve a positive thrust, a purely left travelling wave is chosen ($\lambda = -1$), meaning $C_L = C$ and $C_R = 0$ (see Appendix B.2). Using (3.14) the condition for optimal thrust is

$$ \left| -\frac{i}{2} \int_{-\frac{1}{2}}^{\frac{1}{2}} \hat{Q}(X) e^{-iX} \, dX \right|^2 = 2, \quad \text{and} \quad \left| -\frac{i}{2} \int_{-\frac{1}{2}}^{\frac{1}{2}} \hat{Q}(X) e^{iX} \, dX \right|^2 = 0. \quad (3.15) $$

These equations are two linear constraints on the function $\hat{Q}$, constituting thus an underdetermined system. Therefore there are infinitely many possible solutions that satisfy this. To verify that there is indeed a positive thrust, the conditions (3.15) are substituted into the thrust formula (2.12) to give

$$ \bar{F}_T = 1. \quad (3.16) $$

This makes physical sense, since to get a motion in a right travelling direction, a left travelling wave would be required. The absence of a wave in the rightward direction agrees with this intuition since a wave on the right hand side would reduce the difference in fore-aft amplitude which would reduce the radiation stress and hence reduce the thrust by wasting some power injected.

To validate the numerical result with the analytical expression, we want to do something similar to what was done in Section 3.1. This time, we have less information because the constraints do not give us information on $\hat{Q}$ like (3.8). Therefore, as an ansatz we will prescribe $\hat{Q}$ to be a step function of the form

$$ \hat{Q} = \begin{cases} \alpha, & \text{if } x \in [-\frac{l}{2}, 0), \\ \beta, & \text{if } x \in [0, \frac{l}{2}], \\ 0, & \text{otherwise,} \end{cases} \quad (3.17) $$

and use the numerical optimisation code to find $\alpha, \beta \in \mathbb{C}$ which can then be used to validate $\hat{h}$. This will show that there is no wave on the right hand side as is the case in (3.15). Given the step function in figure 4a, the wave-field is found analytically and numerically and their comparison
<figure><img src="image_5.png" /><figcaption>Figure 4: (a) Resulting plots for the real and imaginary parts of the source $\hat{Q}$ for the bounded power case where $\frac{U}{c} = v = 0$ and $\frac{L\omega}{c} = l = 3\pi/2$. (b) A similar plot for the wave-field $\hat{h}$ for the same case. There is no wave on the right implying this is the optimum result to start and travel rightward.</figcaption></figure>

is demonstrated in figure 4b. To fully validate the result, the condition for the left-travelling wave in (3.15) and the thrust (3.16) are compared using relative error, which were within 0.008% and 0.0069% respectively. To validate the right hand wave, the relative error is undefined since the value we are comparing to is zero. We will instead present the absolute difference in the numerical value, which is $5.9 \times 10^{-8}$. Given the prescribed step function can be substituted into (3.15), we can get constraints on $\alpha$ and $\beta$. However, we leave this calculation for Section 4.2 in the case of general subcritical $v$.

## 4. Subcritical motion $0 < v < 1$

Up to now, all that has been considered is the start-up case. This section will introduce a velocity $U$ as in figure 1. In dimensionless form, the source is moving at a velocity $v = \frac{U}{c}$. The body velocity is scaled by the velocity of the waves emitted $c$ such that the Froude number can be interpreted like a Mach number, $Fr = \frac{U}{\sqrt{gH}} = \frac{U}{c} = Ma$. With this scaling, the dimensionless velocity is split into three different regimes: subcritical velocities ($0 < v < 1$) where the source is travelling slower than the waves it emits, supercritical velocities ($v > 1$) where the source moves quicker than the waves it produces, and the critical velocity ($v = 1$). We will discuss the latter two in Sections 6 and 5 respectively. To begin, we will rewrite the optimisation problem to include a subcritical velocity and discuss the two choices of bound similar to the previous section.

As was the case in Section 3, the assumption of a periodic source is applied. Due to a moving body in a static frame, expressions for the wave-field and the source are

$$h(x, t) = \operatorname{Re} [\hat{h}(x - vt)e^{it}], \qquad Q(x, t) = \operatorname{Re} [\hat{Q}(x - vt)e^{it}]. \qquad (4.1)$$

A Galilean transformation is applied to work in the rest frame of the source. The Galilean transformation will define new spatial and temporal coordinates $\tilde{x}$ and $\tilde{t}$,

$$\tilde{x} = x - vt \qquad \tilde{t} = t. \qquad (4.2)$$

As such, the PDE in (2.3) is reduced to a modified ODE,

$$(1 - v^2)\hat{h}''(\tilde{x}) + 2iv\hat{h}'(\tilde{x}) + \hat{h}(\tilde{x}) = -\hat{Q}(\tilde{x}). \qquad (4.3)$$
This is solved with the same Green's function approach as before:

$$ \hat{h}(\tilde{x}) = A e^{\frac{i \tilde{x}}{1+\nu}} + B e^{-\frac{i \tilde{x}}{1-\nu}} + \frac{i}{2} \int_{-\frac{1}{2}}^{\tilde{x}} \hat{Q}(X) \left( e^{\frac{i(\tilde{x}-X)}{1+\nu}} - e^{\frac{i(X-\tilde{x})}{1-\nu}} \right) dX, \quad (4.4) $$

where $A, B \in \mathbb{C}$ are constants to be found. Sommerfeld boundary conditions will be used again, but they now take a Doppler shifted form,

$$ \hat{h}' \left( \pm \frac{l}{2} \right) = \mp ik_{\pm} \hat{h} \left( \pm \frac{l}{2} \right), \quad (4.5) $$

where $k_+ = \frac{1}{1-\nu}$, $k_- = \frac{1}{1+\nu}$ are the right- and left-travelling Doppler shifted wavenumbers respectively. The resulting wave field is given by

$$ \hat{h}(\tilde{x}) = -\frac{i}{2} \int_{\tilde{x}}^{\frac{l}{2}} \hat{Q}(X) e^{\frac{i(\tilde{x}-X)}{1+\nu}} dX - \frac{i}{2} \int_{-\frac{l}{2}}^{\tilde{x}} \hat{Q}(X) e^{\frac{i(X-\tilde{x})}{1-\nu}} dX. \quad (4.6) $$

An equivalent expression for the thrust can be derived (see Appendix A.2) using the same method as was done to get (2.10) and $\bar{F}_T$ will be as in (2.12). In this case $k_\pm$ are the same as in (4.5) and are interpreted as a correction due to the Doppler shift in the moving frame. Similarly, the time-averaged power in the moving case can be found. Following the steps taken in Appendix A.3, the resulting expression for the time-averaged injected power is

$$ \overline{\text{Pow}} = \langle \hat{Q}, i\hat{h} \rangle - v \langle \hat{Q}, \hat{h}' \rangle, \quad (4.7) $$

where the second term is proportional to the thrust as in (2.12).

Given that we have the power and thrust, the efficiency $\eta$ can be discussed next. The efficiency will tell us how much of the injected power is used as thrust. In dimensionless form, all that is required is (2.12) and (4.7) to provide this information

$$ \eta = \frac{\bar{F}_T}{\overline{\text{Pow}}} = \frac{\langle \hat{Q}, \hat{h}' \rangle}{\langle \hat{Q}, i\hat{h} \rangle - v \langle \hat{Q}, \hat{h}' \rangle}. \quad (4.8) $$

Note that in dimensional form, we see that the efficiency is given as $\eta = \bar{F}_T \cdot c/\overline{\text{Pow}}$, because the energy is being radiated at speed $c$. In the start-up case, where the power was bounded to be $\overline{\text{Pow}} = 1$ the resulting force was $\bar{F}_T = 1$. Hence, the efficiency is $\eta = 1$ when $\hat{Q}$ is optimised. If we had a purely symmetric source, the difference in fore-aft amplitude squared would be zero resulting in $\bar{F}_T = 0$ and hence $\eta = 0$. We will revisit the efficiency in the discussion of the results for the bounded norm and power in the subcritical case.

### 4.1. Case 1 - Bounded norm

Similar to the start-up case (Section 3.1), an analytical approach shows us that to be an optimal solution the source $\hat{Q}(\tilde{x})$ must satisfy

$$ \mathcal{L}_v \hat{Q} = 0, \quad (4.9) $$

where $\mathcal{L}_v = (1-v^2)\partial_{\tilde{x}}\tilde{x} + 2iv\partial_{\tilde{x}} + I$. This means that the source must be composed of a rightward and leftward Doppler-shifted wave

$$ \hat{Q}(\tilde{x}) = D_L e^{\frac{i \tilde{x}}{1+\nu}} + D_R e^{-\frac{i \tilde{x}}{1-\nu}}, \quad (4.10) $$

where conditions on $D_L, D_R \in \mathbb{C}$ are found in Appendix B.1. These are summarised as follows:

$$ |D_L|^2 = 4 \left( l + \frac{l(1-v^2)^2}{[2\lambda(1-v)-l]^2} \sin^2\left(\frac{l}{1-v^2}\right) + \frac{2(1-v^2)^2}{2\lambda(1-v)-l} \sin^2\left(\frac{l}{1-v^2}\right) \right)^{-1}, \quad (4.11) $$
<figure><img src="image_6.png" /><figcaption>Figure 5: (a) Comparison between the analytically and numerically calculated left and (b) right wave amplitudes squared for $l = 3\pi/2$ over both subcritical and supercritical velocities. For the right wave, the numerical scheme becomes unstable close to $v = 1$ and hence it has been hatched.</figcaption></figure>

$$|D_R|^2 = \frac{|D_L|^2(1-v^2)^2}{(2(1-v)\lambda - l)^2} \sin^2\left(\frac{l}{1-v^2}\right). \quad (4.12)$$

Again, $\lambda \in \mathbb{R}$ is the Lagrange multiplier associated with the variational calculus, with eigenvalues

$$\lambda_{\pm} = \frac{vl \pm \sqrt{l^2 - (1-v^2)^3 \sin^2\left(\frac{l}{1-v^2}\right)}}{2(1-v^2)}. \quad (4.13)$$

It is worth noting than given $\lambda \in \mathbb{R}$, the constants $D_L$ and $D_R$ are in phase as they are related via multiplication by a real scalar (see (B 9)). As was the case on start-up, the negative sign is taken to maximise the thrust for a body moving in the positive (right) direction. We validate that $\lambda \in \mathbb{R}$, since

$$l^2 - (1 - v^2)^3 \sin^2\left(\frac{l}{1 - v^2}\right) > 0, \quad (4.14)$$

is satisfied for all real $l, v > 0$.

A study over the dimensionless velocity $v$ is demonstrated in figures 5a and 5b for $l = 3\pi/2$. The forward wave appears to have zero amplitude at certain resonant values of $v$ and $l$, which is discussed later in Section 6. The expressions involving $D_L$ and $D_R$ can be compared with numerics. Using (B 15) to recover $\lambda^2|D_L|^2$ and $\lambda^2|D_R|^2$ in a similar way to (3.8), we show correspondence within 0.221% and 0.0005 respectively. The latter error is quoted as an absolute value since the right wave has zero amplitude at resonance points.

Since there is no forward wave in some conditions, we ask, is there an optimal $l$ for a given $v$? The contour plots in figure 6 best illuminate that the same trend as in figure 3 continues for subcritical values of $v$. For a given $v$, maximal thrust is achieved at maximal $l$. However, there is an interesting behavioural difference between smaller and larger length scales. For a fixed $l$ that is large, the force is maximal close to $v = 0$, but as we shrink the length scale, the optimum begins to move towards $v = 1$, which is demonstrated in figure 6b. For a given short $l$, the time-averaged thrust may be non-monotone with velocity, which is demonstrated in the subcritical part of figure 7b. The efficiency is also plotted showing that $\eta = 1$ is achieved for a small set of velocities close to $v = 1$ for small $l$. This will be shown in greater detail in Section 6.

## 4.2. Case 2 - Bounded power

The general approach taken while moving is the same as it is in the start-up case. The variational calculus approach (Appendix B.2) produces two results depending on the choice of the solution
<figure><img src="image_7.png" /><figcaption>Figure 6: (a) Contour plot where colour corresponds to the magnitude of the force. For a given $v$, a larger time-averaged thrust can be achieved with a larger $l$. (b) The same heat map over a range of smaller $l$ denoted by the red box in (a) to convey the trend for smaller length scales, demonstrating maximal thrust at dimensionless velocities closer to 1.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Figure 7: (a) Time-averaged thrust over velocities from $v = 0$ to $v = 3$ for $l = 3\pi/2$. The efficiency is also plotted on the right axis showing the reduction after the last resonance point. (b) A similar plot for $l = \pi/6$ demonstrates that the region of maximal efficiency reduces due to the last and first resonance points being closer to $v = 1$ for smaller values of $l$.</figcaption></figure>

for the Lagrange multiplier, $\lambda$:

$$C_L e^{\frac{i\bar{x}}{1+v}} + C_R e^{-\frac{i\bar{x}}{1-v}} = C e^{\frac{i\bar{x}}{1+v}} \quad \text{for } \lambda_1 = -1, \qquad (4.15)$$

$$C_L e^{\frac{i\bar{x}}{1+v}} + C_R e^{-\frac{i\bar{x}}{1-v}} = C e^{-\frac{i\bar{x}}{1-v}} \quad \text{for } \lambda_2 = 1, \qquad (4.16)$$

where

$$C_L = -\frac{i}{2} \int_{-\frac{l}{2}}^{\frac{l}{2}} \hat{Q}(X) e^{-\frac{iX}{1+v}} dX, \qquad C_R = -\frac{i}{2} \int_{-\frac{l}{2}}^{\frac{l}{2}} \hat{Q}(X) e^{\frac{iX}{1-v}} dX. \qquad (4.17)$$

The choice of $\lambda$ governs whether we are going to maximise the thrust in the positive or negative direction. Given that the subcritical velocity is positive ($0 < v < 1$), to maximise thrust in the positive direction, $\lambda_1$ is chosen. Using (4.15), $C_L = C$ and $C_R = 0$. The power constraint (4.7) is then used to give

$$\left| -\frac{i}{2} \int_{-\frac{l}{2}}^{\frac{l}{2}} \hat{Q}(X) e^{-\frac{iX}{1+v}} dX \right|^2 = 2(1 + v), \qquad \left| -\frac{i}{2} \int_{-\frac{l}{2}}^{\frac{l}{2}} \hat{Q}(X) e^{\frac{iX}{1-v}} dX \right|^2 = 0 \qquad (4.18)$$
and this can be used to show that for all subcritical $v$,

$$ \bar{F}_T = \frac{1}{4} \left( \frac{2}{1+v} |C_L|^2 - \frac{2}{1-v} |C_R|^2 \right) = 1. \quad (4.19) $$

The above results can be verified numerically over the interval $0 \le v < 1$. The condition for the left travelling wave (4.18) was tested and the infinity norm of the relative error at each point sampled was taken to show that all the points were within 0.0010%. Through a similar test the force (4.19) was shown to be within 0.0011%. As discussed in Section 3.2, the right hand wave condition of zero (4.18) does not permit the calculation of the relative error so we will just provide the infinity norm of the absolute difference between the analytical and numerical solution, which is $1.11 \times 10^{-5}$.

This is a very general solution since there are infinite solutions for $\hat{Q}$ that satisfy (4.18). It is also interesting that we are free to choose any $l$ and still maintain an optimal solution that has $\bar{F}_T = 1$. This corresponds to an efficiency of $\eta = 1$ since the power and force are found to be equal.

To shed some light onto what an example solution would look like, we will return to the example used in Section 3.2 where $\hat{Q}(\tilde{x})$ is defined as a step function in the non-zero domain (3.17). The conditions (4.18) become

$$ \left| -\frac{i}{2} \left( \alpha \int_{-\frac{l}{2}}^{0} e^{-\frac{iX}{1+v}} dX + \beta \int_{0}^{\frac{l}{2}} e^{-\frac{iX}{1+v}} dX \right) \right|^2 = 2(1+v), \quad (4.20) $$

$$ \left| -\frac{i}{2} \left( \alpha \int_{-\frac{l}{2}}^{0} e^{\frac{iX}{1-v}} dX + \beta \int_{0}^{\frac{l}{2}} e^{\frac{iX}{1-v}} dX \right) \right|^2 = 0. \quad (4.21) $$

Equation (4.21) gives us

$$ \alpha = \beta \left( \frac{1 - e^{\frac{il}{2(1-v)}}}{1 - e^{-\frac{il}{2(1-v)}}} \right) = \beta \zeta, \quad (4.22) $$

where $\zeta \in \mathbb{C}$. This can be substituted back into (4.20) to show that

$$ |\beta|^2 = 4 \left( \left| \zeta \left( 1 - e^{\frac{il}{2(1+v)}} \right) + e^{-\frac{il}{2(1+v)}} - 1 \right|^2 \right)^{-1} . \quad (4.23) $$

If we wish to prescribe that $\hat{Q}$ be a step function of this form, the constants $\alpha$ and $\beta$ are constrained to satisfy the condition on their modulus squared, and are only defined up to an arbitrary phase, similar to the bounded norm case. In contrast to the bounded norm case, there is a phase difference between $\alpha$ and $\beta$ since $\zeta \in \mathbb{C}$.

## 5. The critical case ($v = 1$)

In order to study the whole range of velocities, the critical velocity case ($v = 1$) must be studied separately. This is a simple case because (4.3) reduces to a first order ODE,

$$ 2i\hat{h}'(\tilde{x}) + \hat{h}(\tilde{x}) = -\hat{Q}(\tilde{x}). \quad (5.1) $$

The sole boundary condition required is that the fluid to the right of the body is unperturbed, since the body is now moving as quickly as the waves that it emits. The boundary condition takes the form

$$ \hat{h} = 0: \quad \tilde{x} > \frac{l}{2}. \quad (5.2) $$
The resulting particular solution for the wave-field $\hat{h}(\tilde{x})$ is

$$ \hat{h}(\tilde{x}) = -\frac{i}{2} \int_{\tilde{x}}^{\frac{l}{2}} \hat{Q}(X) e^{\frac{i(\tilde{x}-X)}{1+v}} dX, \qquad (5.3) $$

where $v = 1$. Given this expression for $\hat{h}$, we can recover the thrust through radiation stresses

$$ \frac{1}{2} k_{-} |\hat{h}_{-}|^{2} = \langle \hat{Q}, \hat{h}' \rangle, \qquad (5.4) $$

where $k_{-} = \frac{1}{1+v}$ is the wavenumber for the left travelling wave with amplitude $|\hat{h}_{-}|$. This is the same form as (2.10) but the fore wave is zero.

The same analytical and numerical methods are used as the subcritical case and we will briefly outline the results. In the bounded norm case there will only be a wave in the aft direction which from the variational calculus is generated by $\hat{Q}$ such that,

$$ \hat{Q}(\tilde{x}) = D_L e^{\frac{i\tilde{x}}{2}}, \qquad (5.5) $$

where $D_L$ is found to satisfy

$$ |D_L|^2 = \frac{4}{l}. \qquad (5.6) $$

Since there is no forward wave, efficiency of $\eta = 1$ is achieved in this case. The Lagrange multiplier is found to be $\lambda = -\frac{l}{4}$ and the condition in (5.6) can be tested using expressions found in a similar way to (B 15) and $\lambda^2|D_L|^2$ is found to have a relative error of 0.0002%.

For the bounded power case, the condition is the first equation in (4.18) when $v = 1$ (the second condition is automatically satisfied due to the impossibility of a wave on the right hand side). Hence our condition is

$$ \left| -\frac{i}{2} \int_{-\frac{l}{2}}^{\frac{l}{2}} \hat{Q}(X) e^{-\frac{iX}{2}} dX \right|^2 = 4. \qquad (5.7) $$

The resulting force is again $\bar{F}_T = 1$. These two results are validated with the numerical approach with relative errors of 0.037% and $6.07 \times 10^{-12}$% respectively.

## 6. Supercritical motion $v > 1$

The final part of our study is the supercritical velocity regime. If we recall, the dimensionless velocity, $v = U/c$, is the ratio of the body velocity and the emitted wave velocity. The $v > 1$ case, corresponds to the body moving faster than the waves that it emits. Hence, any wave emitted will immediately be overtaken by the body. Similar to the subcritical case, a Galilean transformation into the moving frame is done using the transformation in (4.2).

The radiative boundary conditions applied in Sections 3 and 4 are no longer suitable because no wave can possibly propagate on the right hand side and the left hand side will be composed of two waves of differing wavenumbers. We can use the fact that there is no wave on the right to our advantage by applying the following boundary conditions

$$ \hat{h} = \hat{h}' = 0: \quad \tilde{x} > \frac{l}{2}. \qquad (6.1) $$

The Green's function approach to solving (4.3) is applied again to calculate the solution for the wave-field, $\hat{h}(\tilde{x})$:

$$ \hat{h}(\tilde{x}) = -\frac{i}{2} \int_{\tilde{x}}^{\frac{l}{2}} \hat{Q}(X) e^{\frac{i(\tilde{x}-X)}{1+v}} dX + \frac{i}{2} \int_{\tilde{x}}^{\frac{l}{2}} \hat{Q}(X) e^{\frac{i(\tilde{x}-X)}{v-1}} dX, \qquad (6.2) $$
where we have swapped $-(1 - v)$ with $v - 1$ to emphasise that this solution is now composed of two left travelling waves in the rest frame of the body.

A similar approach of recovering the radiation stress is applied to get the expression for the thrust as in Appendix A.2. The following equivalence relation is recovered:

$$ -\frac{1}{2} [k|\hat{h}|^2]_{\leftarrow} = \langle \hat{Q}, \hat{h}' \rangle, \qquad (6.3) $$

where $\rightarrow$ and $\leftarrow$ denote the right and left travelling waves in the rest frame respectively. This choice of notation represents a subtlety of the supercritical case where both waves are left-travelling waves in the moving frame. But once we translate back into the rest frame, there is a forward and backwards wave. As such, $k$ represents the Doppler shifted wavenumber in the moving frame where $k_{\leftarrow} = \frac{1}{1+v}$ and $k_{\rightarrow} = \frac{1}{v-1}$.

Variational calculus is undertaken in a similar way to the subcritical case. The bounded norm and power cases will be studied again, but since there are only minor changes, rather than being split into sections, we simply summarise the results.

In the bounded norm case the optimal condition again is that $\hat{Q}(\tilde{x})$ is composed of wave-like parts as was the case in (4.10). This time the corresponding $|D_L|^2$, $|D_R|^2$ and $\lambda$ are

$$ |D_L|^2 = 4 \left( l + \frac{l(v^2 - 1)^2}{[2\lambda(v-1)-l]^2} \sin^2\left(\frac{l}{v^2-1}\right) + \frac{2(v^2-1)^2}{2\lambda(v-1)-l} \sin^2\left(\frac{l}{v^2-1}\right) \right)^{-1}, \quad (6.4) $$

$$ |D_R|^2 = \frac{|D_L|^2 (v^2 - 1)^2}{(2(v-1)\lambda - l)^2} \sin^2\left(\frac{l}{v^2-1}\right), \qquad (6.5) $$

$$ \lambda_{\pm} = \frac{l \pm \sqrt{v^2 l^2 - (v^2 - 1)^3 \sin^2\left(\frac{l}{v^2-1}\right)}}{2(v^2 - 1)}. \qquad (6.6) $$

Similar to before, to maximise the force in the positive direction, the negative sign is taken in the expression for $\lambda \in \mathbb{R}$. The analytical and numerical results are verified as they were in the subcritical case. Looking again at the supercritical part of figure 5b, the forward wave demonstrates a resonance behaviour where at certain velocities it has zero amplitude. Turning to the analytical expressions above, we can see that the $\sin^2(l/(v^2 - 1))$ term in (6.5) is important. To verify this, we look at the cases where $|D_R|^2 = 0$. To also link to the subcritical case we will use the fact that $\sin^2(l/(v^2 - 1)) = \sin^2(l/(1 - v^2))$ and find

$$ \sin^2\left(\frac{l}{1-v^2}\right) = 0 \implies v = \left(1 - \frac{l}{n\pi}\right)^{\frac{1}{2}}, \qquad (6.7) $$

where $n \in \mathbb{Z}\{0\}$. The positive square root is taken on the right hand side since $v \ge 0$. In this current form, we see that the subcritical resonances will be seen for all possible $0 < n < \frac{1}{\pi}$ while supercritical resonances for all $n < 0$. Equation (6.7) is used to find the velocities where the right hand wave is zero and is validated in figure 5b. For the supercritical resonances, the largest resonant velocity will be found when $n = -1$.

We have shown that for a given $l$, we can find values of $v$ for which the front wave will be zero. Inspired by the bounded power results, we can use this to our advantage in the interest of efficiency. For supercritical velocities this will be

$$ \eta = \frac{\bar{F}_T}{\text{Pow}} = \frac{-2(v-1)\lambda^2|D_R|^2 + 2(1+v)\lambda^2|D_L|^2}{2(v-1)\lambda^2|D_R|^2 + 2(1+v)\lambda^2|D_L|^2}. \qquad (6.8) $$

It is clear from (6.8) that $|D_R| = 0$ will produce maximal efficiency (see figure 7a at the values of $v$ where $\eta = 1$) which is achieved by wasting no power on the wave that inhibits thrust.
We note that $\eta$ drops for larger $v$ because the body is too short to generate a substantial wave-field to generate thrust. This is demonstrated in figure 7b where the range of values that efficiency is maximal appears to narrow down for small $l$. This manifests itself in (6.7). The last velocity value before efficiency diminishes at the $n = -1$ resonance will reduce as $l$ reduces.

In contrast, the bounded power case will continue to follow the same result in (4.18) where the only change is a $\frac{l}{2}$ rather than $-\frac{l}{2}$ preceding the second integral. Since the absolute value is taken, this has no effect. After inserting these optimal conditions into the expressions for the force, we find that $\bar{F}_T = 1$ still holds. Hence the bounded power case results in $\eta = 1$ for all $v > 0$. The ansatz (3.17) can be used again here but we do not include the results since they are similar.

## 7. Quasi-periodic theory for a changing $v$

We have studied the optimal solution for a given velocity and found that in the bounded power case the time-averaged thrust is a constant ($\bar{F}_T = 1$). All this work has been done in the absence of acceleration. We can now ask if we could create a model that introduces a small acceleration to move between different velocities while maintaining the optimal solution and the assumption that the waves are periodic? A dimensional acceleration such that $\dot{U} \ll \omega c$, permits the use of the previous results where acceleration can be interpreted as changing from one moving frame to the next.

In order to make sense of moving between different velocities, the cruising velocity $v_\star$ must be known. For a given set of conditions, the balance between thrust $F_T$ and drag $F_D$ will return the equilibrium velocity. Given the expression for the time-averaged dimensionless drag force, the result from bounding the power will be used to create the aforementioned force balance and find the intercept which will be referred to as the cruising velocity. The form of $v_\star$ can then be used to understand how slowly varying the parameters of the problem can be used to move from $v = 0$ to $v > 1$.

To find $v_\star$, we need to introduce a drag force. It is important to note that the true drag may contain contributions from wave drag, form drag due to shape, and shallow water resonance, but we skip these details here since they aren't the focus of the study. Instead, we use a simple drag law to demonstrate how the force changes with speed. Therefore the dimensional drag will be considered to be,

$$F_D = \frac{1}{2} C_D \rho U^2 L, \qquad (7.1)$$

where $C_D$ is a general drag coefficient. Similar to Benham et al. (2024), we will neglect oscillatory drag and focus solely on horizontal drag. To get the drag in dimensionless form, the dimensional drag is set equal to the dimensional thrust,

$$\bar{F}_T = \rho g H^2 \langle \hat{Q}, \hat{h}' \rangle. \qquad (7.2)$$

Hence, by dividing by $\rho g H^2$, the corresponding time-averaged dimensionless drag is found to be

$$\bar{F}_D = \frac{1}{2} C_D v^2 \frac{L}{H}. \qquad (7.3)$$

In the previous sections we chose $\delta = 1$ for the bounded power case to compare to the numerical work. In this section we will return to the general result $\overline{Pow} \le \delta$. The resulting optimal time-averaged thrust is $\bar{F}_T = \delta$ which gives us the cruising velocity $v_\star$,

$$\delta = \frac{1}{2} C_D v^2 \frac{L}{H} \implies v_\star(\delta) = \sqrt{\frac{2H\delta}{LC_D}}. \qquad (7.4)$$

The drag coefficient $C_D$, the body length $L$ and the fluid height $H$ are fixed which leaves the power injected $\delta$ as the only parameter to control the velocity. The velocity and power injected
<figure><img src="image_9.png" /><figcaption>Figure 8: (a) Demonstration of the one period of oscillation where $v_{\star} = 0.2$ and $l = 3\pi/2$. The body is defined as the region between the vertical lines ($\tilde{x} \in [-l/2, l/2]$). The period is divided into 10 points where the first half period is in blue and the latter in red. (b) A similar plot showing the time evolution at $v_{\star} = 0.3$ demonstrating that the amplitude must increase to maintain the optimal solution at the larger velocity.</figcaption></figure>

are related by $v_{\star} \sim \delta^{1/2}$. Using the thrust result $\bar{F}_T = \delta$, the optimal conditions on the source $\hat{Q}$ can be shown to also increase with $\delta$ (see (7.5) below). In terms of $\delta$ and $v_{\star}$, the variational calculus produces the following conditions,

$$ \left| -\frac{i}{2} \int_{-\frac{l}{2}}^{\frac{l}{2}} \hat{Q}(X) e^{-\frac{iX}{1+v_{\star}}} dX \right|^2 = 2\delta(1 + v_{\star}(\delta)), \quad \left| -\frac{i}{2} \int_{-\frac{l}{2}}^{\frac{l}{2}} \hat{Q}(X) e^{\frac{iX}{1-v_{\star}}} dX \right|^2 = 0. \quad (7.5) $$

Hence, as the velocity increases, we need a larger wave to maintain the optimal thrust.

Given a small acceleration we can vary the cruising velocity $v_{\star}$ with the time-averaged power injected $\delta$. From here we can make two observations. To move quicker, one has to provide more power and to remain moving efficiently, the aft wave gets larger while the fore wave remains zero. This makes physical sense because the faster velocity requires that we generate more thrust which requires a bigger difference in fore-aft amplitude as is the case in (2.10). Modifications can be made to this to delve further into the conditions for accelerating from start-up. Motivated by previous sections, we could prescribe $\hat{Q}$. For example, if we used a step function like (3.17), we could study how the conditions for $\alpha$ and $\beta$ evolve in terms of the changing speed $v_{\star}$. We could also study adding in another constraint on top of the bounded power. If the norm is bounded simultaneously, there will be an upper bound on the magnitude of the source and hence there will come a point where the thrust will diminish from achieving $\bar{F}_T = \delta$ and efficiency will suffer. Since this section is a demonstration of a simple theory, these modifications will not be investigated here.

To conclude this section, we will take this opportunity to look at the time-evolution of the solutions. We will reintroduce the temporal part and study the wave-field $h(x, t)$ as it is in (4.1) to produce a graph showing one whole period of oscillation. The numerically optimised bounded power solution is plotted over one period in figure 8a (note, this is just one optimal $\hat{Q}$ out of infinitely many that satisfy (4.18)). To visualise the results of (7.5), a comparison between a faster cruising velocity is shown in figure 8b. Here we can see that the amplitude of oscillation is required to be greater to maintain optimal efficiency at the larger velocity.
# 8. Discussion

In this paper, we have used gravity waves to generate a thrust to optimally push a body along the shallow water surface. Given one of two constraints where a bound is placed on either the power or the norm of the source, the conditions for optimal thrust are found. The analytical work sheds light on the choice of solution to optimise thrust while the numerical work is a validation step through solving the problem in a different way. A dimensionless velocity based on the Froude number ($v = U/c = U/\sqrt{gH}$) is used to divide the work into subcritical, critical and supercritical velocities. Studying each of these velocities using a Galilean transformation to operate in the rest frame of the body allowed us to build a template for the optimal solution at each velocity. In the case of the bounded power, we use this template to study how the power can be modulated to vary the cruising velocity from start-up to supercritical velocities. A small acceleration $\dot{U} \ll \omega c$ is required to maintain the periodic assumption. In the bounded norm case, the intuition from the bounded power result is applied to find resonance points in the bounded norm case to minimise the forward wave and thus increase efficiency.

One point to note is that this work is not restricted to a fluid dynamics setting. The present work could also be used as a template to study the optimal conditions associated with variations of the wave equation and similar PDEs such as the heat, Schrödinger and KdV equations. It would be interesting to replace the Galilean transformation with a relativistic Lorentz transformation to allow for time dilation. One major difference will be that the critical and supercritical velocities will be disallowed due to the constraint that $U/c < 1$.

These results are of course limited by the assumptions made in the shallow water limit, but they do provide some general guidance. In order to achieve maximal efficiency, the goal based off of the intuition built here is to minimise the wave that is being produced ahead of the body. To refine this intuition, we wish to extend to the case where the surface of the shallow water is two dimensional. The expectation is to see a symmetry on either side of the centreline that is parallel to the direction of motion, much like what is observed when a honeybee propels itself along the water surface (Roh & Gharib 2019). Working in the same set of velocity regimes, a similar template for optimal solutions in the two dimensional case can be used to generate an updated description of the quasiperiodic movement between different velocities under the restrictions of $\dot{U} \ll \omega c$.

We can also study the problems outside of the shallow limit to build a model for the subsurface dynamics and dispersive waves. Some work with potential flow as used by Benham et al. (2024) could be applied with the addition of a moving reference frame to operate at different velocities. We could increase the complexity further with the full Navier-Stokes equations. This will introduce the potential for vortices to be generated by the movement of the body. Similar to Gao & Feng (2011) we could then estimate the role that the momentum transferred to the vortices plays in propulsion as well as the viscous and/or turbulent dissipation which are properties not present in the current framework. We could then study the difference in propulsion of the body from waves and vortices (WDP and VDP). Similar to Steinmann et al. (2021) we could estimate the values of the different forces at play to understand what dominates in our set-up at different scales.

Unlike what is done by Benham et al. (2024), the work presented in this paper only considers the raft as a pressure region. An interesting study would be to couple this work to a rigid or flexible beam and investigate the optimal solution for WDP in this case both theoretically and experimentally. The example of the step function sheds some light onto the type of optimal pressure distribution that may be required. Similarly, we could also look at optimal solutions for pressure distributions used in the literature like in the study of air-cushion vehicles (Doctors & Sharma 1972).

On our way to considering more physical setups, we can broaden our horizon past studying a periodic source. This will allow the implementation of multiple bodies on the fluid surface as
well as the study of different incident waves on a single body. Removal of the assumption of a periodically oscillating body will permit the body to interact with external waves as well as allowing us to generate a more complete picture of accelerating between different velocities. We can then study an optimal control problem that could apply the wave-riding and wave-passing that ducklings demonstrate as described by Yuan et al. (2021). The interactions between the body and the waves will have to be studied in more detail since the body is always riding the waves produced by another. We could also use this version of the problem to investigate ways of optimising the distance between bodies to keep them stationary.

This paper has helped us mathematically validate the key elements such as reducing the forward wave to maximise efficiency as well as build an understanding on the representation of WDP in the different velocity regimes. What is clear is there is plenty of more work to do, but this is a beginning that builds physical intuition in a simplified case that presents the opportunity for analytical as well as numerical study.

## Funding

This work was partially funded by the European Research Council (ERC) under the EU Horizon 2020 research and innovation program (Grant Agreement 833125-HIGHWAVE)

A github repository of example code used in this paper can be found at https://github.com/daireodonovan/Optimal-locomotion-using-self-generated-waves.git

# Appendix A. Derivations of the shallow water wave equation, thrust and power

## A.1. Shallow water wave equation

Beginning with the Euler equations, we will work in a one dimensional shallow water set up where the fluid is of height $H$ with characteristic length scale $L_w$. Through the use of the kinematic, dynamic and impermeable base boundary conditions we have a set of equations that can be solved. A pressure source perturbs the waves such that the dynamic boundary condition is,

$$p = P(x,t) : z = H + h(x,t), \qquad (A1)$$

where $h(x,t) \ll H$ is the height of the resulting wave from the perturbation. We can begin by solving the rest case where $p = p_a$ (atmospheric pressure) followed by a perturbation with a pressure $p = \epsilon P(x,t)$, where $\epsilon = H/L_w \ll 1$ is the aspect ratio. The system of equations is taken to the leading order and the resulting horizontal momentum equation can be substituted into the vertically integrated conservation of mass equation to arrive at,

$$\frac{\partial^2 h}{\partial x^2} - \frac{1}{c^2} \frac{\partial^2 h}{\partial t^2} = - \frac{1}{\rho g} Q(x,t), \qquad -\infty < x < \infty. \qquad (A2)$$

## A.2. Time-averaged thrust

We begin by manipulating the wave equation (in the subcritical moving frame) through multiplication by $\hat{h}_{\bar{x}}$ and integrating, such that,

$$\int_{-}^{+} (h_{\bar{t}\bar{t}} - 2v h_{\bar{x}\bar{t}} - (1 - v^2)h_{\bar{x}\bar{x}}) h_{\bar{x}} d\bar{x} = \int_{-}^{+} Q h_{\bar{x}} d\bar{x}, \qquad (A3)$$

where we assume a periodic source such that the wave-field is also periodic,

$$h = \mathrm{Re} \left( \hat{h}(\bar{x}) e^{i \bar{t}} \right). \qquad (A4)$$
Substituting this in and taking the time average, we note that the second term on the left of (A 3) does not contribute to the thrust after time-averaging. This results in,

$$ \frac{1}{4}(v^2 - 1)[\hat{h}'\hat{h}'']_{-}^{+} - \frac{1}{4}[\hat{h}\hat{h}^{*}]_{-}^{+} = \frac{1}{4}\int_{-\frac{l}{2}}^{\frac{l}{2}} (\hat{Q}^{*}\hat{h}' + \hat{Q}\hat{h}'') d\tilde{x}. \quad (A5) $$

Applying the boundary conditions on the right and left hand side of the body (+ and – respectively)

$$ \hat{h}'(+) = -\frac{i}{1-v}\hat{h}(+), \qquad \hat{h}'(-) = \frac{i}{1+v}\hat{h}(-), \quad (A6) $$

results in

$$ \frac{1}{4} \left( \frac{v^2 - 1}{(1 - v)^2} |\hat{h}|^2 \Big|_+ - \frac{v^2 - 1}{(1 + v)^2} |\hat{h}|^2 \Big|_- - |\hat{h}|^2 \Big|_+ + |\hat{h}|^2 \Big|_- \right) = \langle \hat{Q}, \hat{h}' \rangle, \quad (A7) $$

$$ \Rightarrow -\frac{1}{2} \left( \frac{1}{1-v} |\hat{h}|^2 \Big|_+ - \frac{1}{1+v} |\hat{h}|^2 \Big|_- \right)^+ = -\frac{1}{2} [k |\hat{h}|^2]_-^+ = \langle \hat{Q}, \hat{h}' \rangle, \quad (A8) $$

where we have used the notation defined in (2.11). This can be defined as the thrust because the left hand side is equivalent to the radiation stress, with a difference in fore-aft amplitude. The k acts as a correction to the wavenumber due to the waves getting Doppler shifted in the moving frame. Hence, we define the force as,

$$ \bar{F}_T := \langle \hat{Q}, \hat{h}' \rangle. \quad (A9) $$

### A.3. Time-averaged power

Two routes can be taken to derive the time-averaged power. We could follow steps analogous to Appendix A.2 and instead of multiplying by $\hat{h}_{\tilde{x}}$ we instead multiply by $\hat{h}_{\tilde{t}}$. Here we will derive the power by taking the time derivative of the energy as it is commonly done.

Beginning with the dimensionless energy,

$$ E = \int_{(-)+v\tilde{t}}^{(+)+v\tilde{t}} \left( \frac{1}{2} h_{\tilde{t}}^2 + \frac{1}{2} h_{\tilde{x}}^2 \right) dx, \quad (A10) $$

we can transform the energy equation into the moving frame using the Galilean transformation in (4.2),

$$ E = \int_{-}^{+} \left( \frac{1}{2} h_{\tilde{t}}^2 + \frac{1}{2} (1 + v^2) h_{\tilde{x}}^2 - v h_{\tilde{x}} h_{\tilde{t}} \right) d\tilde{x}. \quad (A11) $$

Taking the time derivative will give us power,

$$ \frac{dE}{d\tilde{t}} = \int_{-}^{+} \left( (h_{\tilde{t}} - v h_{\tilde{x}}) h_{\tilde{t}\tilde{t}} + (1 + v^2) h_{\tilde{x}\tilde{t}} h_{\tilde{x}} - v h_{\tilde{t}} h_{\tilde{x}\tilde{t}} \right) d\tilde{x}. \quad (A12) $$

We can use the wave equation in the moving frame to substitute in for $h_{\tilde{t}\tilde{t}}$,

$$ \frac{dE}{d\tilde{t}} = \int_{-}^{+} \left[ Q h_{\tilde{t}} - v Q h_{\tilde{x}} + (1 - v^2)(h_{\tilde{x}} h_{\tilde{x}\tilde{t}} + h_{\tilde{t}} h_{\tilde{x}\tilde{x}}) - v(1 - v^2)h_{\tilde{x}} h_{\tilde{x}\tilde{x}} + v h_{\tilde{t}} h_{\tilde{x}\tilde{t}} \right] d\tilde{x}. \quad (A13) $$

The first two elements of the integral will be left as they are because they include the control function $Q$, but the rest of the terms can be written in terms of derivatives with respect to $\tilde{x}$ and result in boundary terms

$$ \frac{dE}{d\tilde{t}} = \int_{-}^{+} \left( Q h_{\tilde{t}} - v Q h_{\tilde{x}} \right) d\tilde{x} + (1 - v^2) \left[ h_{\tilde{x}} h_{\tilde{t}} \right]_{-}^{+} - \frac{1}{2} v (1 - v^2) \left[ h_{\tilde{x}}^2 \right]_{-}^{+} + \frac{1}{2} v \left[ h_{\tilde{t}}^2 \right]_{-}^{+}. \quad (A14) $$

We are working with periodic waves and if we take the time-average over a period, the result
is $E(T) - E(0) = 0$. Thus we make the equivalence relation between the power radiated by the waves and the power injected by the source,

$$ -(1 - v^2) [\overline{\hat{h}_{\tilde{x}} \hat{h}_{\tilde{t}}}]_-^+ + \frac{1}{2} v (1 - v^2) [\overline{\hat{h}_{\tilde{x}}^2}]_-^+ - \frac{1}{2} v [\overline{\hat{h}_{\tilde{t}}^2}]_-^+ = \int_-^+ \overline{(Q \hat{h}_{\tilde{t}} - v Q \hat{h}_{\tilde{x}})} d\tilde{x}. \quad (A 15) $$

We can see that the second term on the right hand side is precisely the time-averaged thrust as derived in A.2. Using the periodic assumption to decompose into spatial and temporal parts, we write the time-averaged power injected as,

$$ \overline{\text{Pow}} = \langle \hat{Q}, i\hat{h} \rangle - v \langle \hat{Q}, \hat{h}' \rangle. \quad (A 16) $$

Motivated by the thrust being written as the difference in fore-aft amplitude squared, we can expand out the left hand side of (A 15) in a similar way to what is done in A.2 to show,

$$ \frac{1}{2} \left( \frac{1}{1-v} |\hat{h}_+|^2 + \frac{1}{1+v} |\hat{h}_-|^2 \right) = \langle \hat{Q}, i\hat{h} \rangle - v \langle \hat{Q}, \hat{h}' \rangle. \quad (A 17) $$

Hence, the power is the sum of the fore-aft amplitude squared.

# Appendix B. Variational Calculus Approach

## B.1. Case 1 - Bounded norm variational calculus

The time-averaged thrust force from (2.12) is used with $\hat{h}(\tilde{x})$ (4.6) substituted in for the subcritical case. This will simplify the following steps in the analytical optimisation approach where the thrust is given by,

$$ \bar{F}_T = \frac{1}{4} \left( - \frac{1}{2(1-v)} \int_{-\frac{1}{2}}^{\frac{1}{2}} \hat{Q}(\tilde{x}) \int_{-\frac{1}{2}}^{\frac{1}{2}} \hat{Q}^*(X) e^{\frac{i(\tilde{x}-X)}{1-v}} dX d\tilde{x} \right. \\ \left. + \frac{1}{2(1+v)} \int_{-\frac{1}{2}}^{\frac{1}{2}} \hat{Q}(\tilde{x}) \int_{-\frac{1}{2}}^{\frac{1}{2}} \hat{Q}^*(X) e^{\frac{i(X-\tilde{x})}{1+v}} dX d\tilde{x} \right). \quad (B 1) $$

The action $f(\hat{Q})$ can be written with the thrust as the objective and the norm of the control function $\hat{Q}(\tilde{x})$ as the constraint

$$ f(\hat{Q}) = \frac{1}{4} \left( - \frac{1}{2(1-v)} \int_{-\frac{1}{2}}^{\frac{1}{2}} \hat{Q}(\tilde{x}) \int_{-\frac{1}{2}}^{\frac{1}{2}} \hat{Q}^*(X) e^{\frac{i(\tilde{x}-X)}{1-v}} dX d\tilde{x} \right. \\ \left. + \frac{1}{2(1+v)} \int_{-\frac{1}{2}}^{\frac{1}{2}} \hat{Q}(\tilde{x}) \int_{-\frac{1}{2}}^{\frac{1}{2}} \hat{Q}^*(X) e^{\frac{i(X-\tilde{x})}{1+v}} dX d\tilde{x} \right) + \frac{\lambda}{4} \left( \int_{-\frac{1}{2}}^{\frac{1}{2}} \hat{Q}^* \hat{Q} d\tilde{x} \right), \quad (B 2) $$

where $\lambda$ is the Lagrange multiplier. Consider a perturbation $\epsilon\alpha(\tilde{x})$ to the function $\hat{Q}(\tilde{x})$. Variation of the action is defined as

$$ \delta f = \lim_{\epsilon \to 0} \frac{f(\hat{Q}(\tilde{x}) + \epsilon\alpha(\tilde{x})) - f(\hat{Q}(\tilde{x}))}{\epsilon}. \quad (B 3) $$

The perturbed action is

$$ \delta f(\hat{Q}) = \frac{1}{4} \left( - \frac{1}{2(1-v)} \int_{-\frac{1}{2}}^{\frac{1}{2}} \alpha(\tilde{x}) \int_{-\frac{1}{2}}^{\frac{1}{2}} \hat{Q}^*(X) e^{\frac{i(\tilde{x}-X)}{1-v}} dX d\tilde{x} \right. \\ \left. + \frac{1}{2(1+v)} \int_{-\frac{1}{2}}^{\frac{1}{2}} \alpha(\tilde{x}) \int_{-\frac{1}{2}}^{\frac{1}{2}} \hat{Q}^*(X) e^{\frac{i(X-\tilde{x})}{1+v}} dX d\tilde{x} \right) + \frac{\lambda}{4} \left( \int_{-\frac{1}{2}}^{\frac{1}{2}} \alpha(\tilde{x}) \hat{Q}^*(\tilde{x}) d\tilde{x} \right) + c.c., \quad (B 4) $$
where c.c. denotes the complex conjugate.

When the variation is set to zero, the integrand of the $\tilde{x}$ integral can be interpreted as equalling to zero. The expression is further reduced by only taking the coefficients of $\alpha(\tilde{x})$ since the coefficients of $\alpha^*(\tilde{x})$ will also be zero since it is just the complex conjugate. Taking the complex conjugate of the results gives,

$$ -\frac{1}{2(1-\nu)} \int_{-\frac{l}{2}}^{\frac{l}{2}} \hat{Q}(X) e^{\frac{i(X-\tilde{x})}{1-\nu}} dX + \frac{1}{2(1+\nu)} \int_{-\frac{l}{2}}^{\frac{l}{2}} \hat{Q}(X) e^{\frac{i(\tilde{x}-X)}{1+\nu}} dX + \lambda \hat{Q}(\tilde{x}) = 0. \quad (B 5) $$

The optimal condition can be found through operating on (B 5) with $\mathcal{L}_v = (1 - \nu^2)\partial_{\tilde{x}}\tilde{x} + 2iv\partial_{\tilde{x}} + \mathbb{I}$. This has a null space proportional to $e^{\frac{i\tilde{x}}{1+\nu}}$ and $e^{-\frac{i\tilde{x}}{1-\nu}}$ as shown previously. Therefore, application of $\mathcal{L}_v$ on (B 5) results in

$$ \mathcal{L}_v \hat{Q} = 0. \quad (B 6) $$

This equation is easily solved and it tells us that $\hat{Q}(\tilde{x})$ is composed of the sum of two wave-like parts

$$ \hat{Q}(\tilde{x}) = D_L e^{\frac{i\tilde{x}}{1+\nu}} + D_R e^{-\frac{i\tilde{x}}{1-\nu}}, \quad (B 7) $$

where the constants $D_L, D_R \in \mathbb{C}$. This information can be substituted back into (B 5), where the coefficients to $e^{\pm \frac{i\tilde{x}}{1+\nu}}$ are matched to generate two equations. Another equation is needed since there are three unknowns, $D_L, D_R$ and $\lambda$. The norm constraint in (3.1) is therefore bounded to an arbitrary value (we will choose 1) to get another equation to make the following system solvable,

$$ \frac{1}{2(1+\nu)} \left( D_L l + D_R (1-\nu^2) \sin\left(\frac{l}{1-\nu^2}\right) \right) = -\lambda D_L, \quad (B 8a) $$

$$ -\frac{1}{2(1-\nu)} \left( D_L (1-\nu^2) \sin\left(\frac{l}{1-\nu^2}\right) + D_R l \right) = -\lambda D_R, \quad (B 8b) $$

$$ |D_L|^2 l + |D_R|^2 l + D_L^* D_R (1-\nu^2) \sin\left(\frac{l}{1-\nu^2}\right) + D_L^* D_R (1-\nu^2) \sin\left(\frac{l}{1-\nu^2}\right) = 4. \quad (B 8c) $$

Equation (B 8b) demonstrates the relation between the constants $D_L$ and $D_R$,

$$ D_R = D_L \frac{(1-\nu^2)}{2(1-\nu)\lambda - l} \sin\left(\frac{l}{1-\nu^2}\right). \quad (B 9) $$

Since the Lagrange multiplier $\lambda$ is real by definition, $D_R$ is multiplied by a real number. This means that $D_L$ and $D_R$ are in phase. Equation (B 8a) is used along with (B 9) to show that

$$ \lambda_{\pm} = \frac{\nu l \pm \sqrt{l^2 - (1-\nu^2)^3 s^2}}{2(1-\nu^2)}, \quad (B 10) $$

where $s = \sin(\frac{l}{1-\nu^2})$. As mentioned already, $\lambda \in \mathbb{R}$, which means that $l^2 - (1-\nu^2)^3 s^2 \ge 0$ must hold. It can be shown that this is satisfied for all real $l, \nu > 0$. Information on $D_L$ and therefore $D_R$ can be found using (B 8c),

$$ |D_L|^2 = 4 \left( l + \frac{l(1-\nu^2)^2}{[2\lambda(1-\nu)-l]^2} \sin^2\left(\frac{l}{1-\nu^2}\right) + \frac{2(1-\nu^2)^2}{2\lambda(1-\nu)-l} \sin^2\left(\frac{l}{1-\nu^2}\right) \right)^{-1}. \quad (B 11) $$

Finally (B 9) is used to get information on $D_R$,

$$ |D_R|^2 = \frac{|D_L|^2 (1-\nu^2)^2}{(2(1-\nu)\lambda - l)^2} \sin^2\left(\frac{l}{1-\nu^2}\right). \quad (B 12) $$
All that is recovered is information on the modulus squared of $D_L, D_R \in \mathbb{C}$. Thus, there are infinite possible solutions for $\hat{Q}$ defined up to an arbitrary phase.

In the interest of studying efficiency, we can write the force and power in terms of the constants $D_L$ and $D_R$. The force is given by

$$ \bar{F}_T = \frac{1}{4} \left( -\frac{1}{2(1-\nu)} \left| \int_{-\frac{1}{2}}^{\frac{1}{2}} \hat{Q}(X) e^{\frac{iX}{1-\nu}} dX \right|^2 + \frac{1}{2(1+\nu)} \left| \int_{-\frac{1}{2}}^{\frac{1}{2}} \hat{Q}(X) e^{-\frac{iX}{1+\nu}} dX \right|^2 \right), \quad (B13) $$

$$ \Rightarrow \bar{F}_T = \frac{1}{4} \left( -2(1-\nu)\lambda^2 |D_R|^2 + 2(1+\nu)\lambda^2 |D_L|^2 \right). \quad (B14) $$

We have used the result in (B 7) and substituted in for $\lambda \hat{Q}$ in (B 5) to compare coefficients,

$$ \frac{1}{2(1+\nu)} \int_{-\frac{1}{2}}^{\frac{1}{2}} \hat{Q}(X) e^{\frac{i(\bar{x}-X)}{1+\nu}} dX = -\lambda D_L, \quad \frac{1}{2(1-\nu)} \int_{-\frac{1}{2}}^{\frac{1}{2}} \hat{Q}(X) e^{\frac{i(X-\bar{x})}{1-\nu}} dX = \lambda D_R, \quad (B15) $$

and the modulus squared of (B 15) is taken to write the thrust in terms of $D_L, D_R$. Similarly the power can be found to be,

$$ \overline{\text{Pow}} = \frac{1}{4} \left( 2(1-\nu)\lambda^2 |D_R|^2 + 2(1+\nu)\lambda^2 |D_L|^2 \right). \quad (B16) $$

## B.2. Case 2 - Bounded power variational calculus

A similar approach will be taken to the bounded norm case. Beginning with the expression for the thrust in (B 1), the constraint in this case will be the power. Similar to what was done with the thrust, the power is simplified by substituting in for $\hat{h}(\bar{x})$

$$ \overline{\text{Pow}} = \frac{1}{4} \left( \frac{1}{2} \int_{-\frac{1}{2}}^{\frac{1}{2}} \hat{Q}(\bar{x}) \int_{-\frac{1}{2}}^{\frac{1}{2}} \hat{Q}^*(X) e^{\frac{i(X-\bar{x})}{1+\nu}} dX d\bar{x} + \frac{1}{2} \int_{-\frac{1}{2}}^{\frac{1}{2}} \hat{Q}(\bar{x}) \int_{-\frac{1}{2}}^{\frac{1}{2}} \hat{Q}^*(X) e^{\frac{i(\bar{x}-X)}{1-\nu}} dX d\bar{x} \right) \\ - \nu \bar{F}_T. \quad (B17) $$

The resulting action is

$$ f(\hat{Q}) = (1 - \lambda\nu)\bar{F}_T + \frac{\lambda}{4} \left( \frac{1}{2} \int_{-\frac{1}{2}}^{\frac{1}{2}} \hat{Q}(\bar{x}) \int_{-\frac{1}{2}}^{\frac{1}{2}} \hat{Q}^*(X) e^{\frac{i(X-\bar{x})}{1+\nu}} dX d\bar{x} + \frac{1}{2} \int_{-\frac{1}{2}}^{\frac{1}{2}} \hat{Q}(\bar{x}) \int_{-\frac{1}{2}}^{\frac{1}{2}} \hat{Q}^*(X) e^{\frac{i(\bar{x}-X)}{1-\nu}} dX d\bar{x} \right). \quad (B18) $$

Perturbing the action as it was defined in (B 3) and setting the variation to be zero results in

$$ \left( -\frac{1}{2(1-\nu)} \int_{-\frac{1}{2}}^{\frac{1}{2}} \hat{Q}^*(X) e^{\frac{i(\bar{x}-X)}{1-\nu}} dX + \frac{1}{2(1+\nu)} \int_{-\frac{1}{2}}^{\frac{1}{2}} \hat{Q}^*(X) e^{\frac{i(X-\bar{x})}{1+\nu}} dX \right) + \frac{\lambda}{1-\lambda\nu} \left( \frac{1}{2} \int_{-\frac{1}{2}}^{\frac{1}{2}} \hat{Q}^*(X) e^{\frac{i(\bar{x}-X)}{1-\nu}} dX + \frac{1}{2} \int_{-\frac{1}{2}}^{\frac{1}{2}} \hat{Q}^*(X) e^{\frac{i(X-\bar{x})}{1+\nu}} dX \right) = 0. \quad (B19) $$

The complex conjugate of (B 19) is then taken to work in terms of $\hat{Q}$,

$$ \left( -\frac{1}{2(1-\nu)} \int_{-\frac{1}{2}}^{\frac{1}{2}} \hat{Q}(X) e^{\frac{i(X-\bar{x})}{1-\nu}} dX + \frac{1}{2(1+\nu)} \int_{-\frac{1}{2}}^{\frac{1}{2}} \hat{Q}(X) e^{-\frac{i(X-\bar{x})}{1+\nu}} dX \right) + \frac{\lambda}{1-\lambda\nu} \left( \frac{1}{2} \int_{-\frac{1}{2}}^{\frac{1}{2}} \hat{Q}(X) e^{\frac{i(X-\bar{x})}{1-\nu}} dX + \frac{1}{2} \int_{-\frac{1}{2}}^{\frac{1}{2}} \hat{Q}(X) e^{-\frac{i(X-\bar{x})}{1+\nu}} dX \right) = 0. \quad (B20) $$
The first bracket of equation (B 20) is related to the second bracket by differentiation, leaving an eigenvalue problem to solve

$$ \left( \partial_{\tilde{x}} + \frac{i\lambda}{1-\lambda v} \right) \left( -\frac{i}{2} \int_{-\frac{1}{2}}^{\frac{1}{2}} \hat{Q}(X) e^{\frac{i(X-\tilde{x})}{1-v}} dX - \frac{i}{2} \int_{-\frac{1}{2}}^{\frac{1}{2}} \hat{Q}(X) e^{\frac{i(\tilde{x}-X)}{1+v}} dX \right) = 0. \quad (B 21) $$

Equation (B 21) has solutions of the form $Ce^{-\frac{i\lambda}{1-\lambda v}\tilde{x}}$, such that

$$ -\frac{i}{2} \int_{-\frac{1}{2}}^{\frac{1}{2}} \hat{Q}(X) e^{\frac{i(X-\tilde{x})}{1-v}} dX - \frac{i}{2} \int_{-\frac{1}{2}}^{\frac{1}{2}} \hat{Q}(X) e^{\frac{i(\tilde{x}-X)}{1+v}} dX = Ce^{-\frac{i\lambda}{1-\lambda v}\tilde{x}}. \quad (B 22) $$

To find $\lambda$, the operator $\mathcal{L}_v = (1 - v^2)\partial_{\tilde{x}\tilde{x}} + 2iv\partial_{\tilde{x}} + \mathbb{I}$ is applied to (B 22). All terms proportional to $e^{\frac{i\tilde{x}}{1+v}}$ and $e^{-\frac{i\tilde{x}}{1-v}}$ are part of the null space of $\mathcal{L}_v\hat{h}(\tilde{x})$ and hence vanish. As a result we have,

$$ \mathcal{L}_v e^{-\frac{i\lambda}{1-\lambda v}\tilde{x}} = 0, \quad (B 23) $$

and two solutions for $\lambda$ are recovered,

$$ \lambda_1 = -1 \qquad \lambda_2 = 1. \quad (B 24) $$

There are two possible solutions to the eigenvalue problem,

$$ C_L e^{\frac{i\tilde{x}}{1+v}} + C_R e^{-\frac{i\tilde{x}}{1-v}} = C e^{\frac{i\tilde{x}}{1+v}} \quad \text{for } \lambda_1, \quad (B 25) $$

$$ C_L e^{\frac{i\tilde{x}}{1+v}} + C_R e^{-\frac{i\tilde{x}}{1-v}} = C e^{-\frac{i\tilde{x}}{1-v}} \quad \text{for } \lambda_2. \quad (B 26) $$

Each solution has one constant set to zero and the other non-zero. This corresponds to a purely left or right travelling wave. Expressions for $C_L$ and $C_R$ are found by matching wave coefficients in (B 22),

$$ C_L = -\frac{i}{2} \int_{-\frac{1}{2}}^{\frac{1}{2}} \hat{Q}(X) e^{-\frac{iX}{1+v}} dX, \qquad C_R = -\frac{i}{2} \int_{-\frac{1}{2}}^{\frac{1}{2}} \hat{Q}(X) e^{\frac{iX}{1-v}} dX. \quad (B 27) $$

To find conditions for $C_L$ and $C_R$ to satisfy, the power (A 17) is used,

$$ \frac{2}{(1+v)} |C_L|^2 + \frac{2}{(1-v)} |C_R|^2 = 4, \quad (B 28) $$

where (B 27) is used to substitute in $C_L$ and $C_R$. To maximise thrust in the positive direction, we will choose a purely left travelling wave ($\lambda_1$). This choice sets $C_L = C$ and $C_R = 0$,

$$ \left| -\frac{i}{2} \int_{-\frac{1}{2}}^{\frac{1}{2}} \hat{Q}(X) e^{-\frac{iX}{1+v}} dX \right|^2 = 2(1+v), \qquad \left| -\frac{i}{2} \int_{-\frac{1}{2}}^{\frac{1}{2}} \hat{Q}(X) e^{\frac{iX}{1-v}} dX \right|^2 = 0. \quad (B 29) $$

This is the condition that has to be met for the optimal solution comprising of a left travelling wave only. Similarly, the thrust can be written in terms of $C_L$ and $C_R$ and when the conditions in (B 27) are substituted into (B 13), the optimum thrust for a velocity $v$ is recovered

$$ \bar{F}_T = \frac{1}{4} \left( \frac{2}{1+v} |C_L|^2 - \frac{2}{1-v} |C_R|^2 \right) = 1. \quad (B 30) $$

## REFERENCES

* BAROTTA, J.W, THOMSON, S.J., ALVENTOSA, L.F.L., LEWIS, M. & HARRIS, D.M. 2023 Bidirectional wave-propelled capillary spinners. *Commun. Physics* **6** (1), 87.

* BENHAM, G.P., DEVAUCHELLE, O., MORRIS, S.W. & NEUFELD, J.A. 2022 Gunwale bobbing. *Phys. Rev. Fluids* **7** (7), 074804.
BENHAM, G.P., DEVAUCHELLE, O. & THOMSON, S.J. 2024 On wave-driven propulsion. *J. Fluid Mech.* **987**, A44.

BUSH, J.W.M. & Hu, D.L. 2006 Walking on water: Biolocomotion at the interface. *Annu. Rev. Fluid Mech.* **38**, 339–369.

BÖCKMANN, E., YRKE, A. & STEEN, S. 2018 Fuel savings for a general cargo ship employing retractable bow foils. *App. Ocean Research* **76**, 1–10.

BÜHLER, O. 2007 Impulsive fluid forcing and water strider locomotion. *J. Fluid Mech.* **573**, 211–236.

DOCTORS, L.J 1997 Optimal pressure distributions for river-based air-cushion vehicles. *Ship Technology Research* **44**, 32–36.

DOCTORS, L.J. & SHARMA, S.D. 1972 The wave resistance of an air-cushion vehicle in steady and accelerated motion. *J. Ship Research* **16** (04), 248–260.

DODE, A., CARMIGNIANI, R., COHEN, C., CLANET, C. & BOCQUET, L. 2022 Wave drag during an unsteady motion. *J. Fluid Mech.* **951**, A15.

GAO, P. & FENG, J.J. 2011 A numerical investigation of the propulsion of water walkers. *J. Fluid Mech.* **668**, 363–383.

HO, I., PUCCI, G., OZA, A.U. & HARRIS, D.M. 2023 Capillary surfers: Wave-driven particles at a vibrating fluid interface. *Phys. Rev. Fluids* **8**, L112001.

HU, D.L., CHAN, B. & BUSH, J.W.M. 2003 The hydrodynamics of water strider locomotion. *Nature* **424** (6949), 663–666.

LONGUET-HIGGINS, M.S. & STEWART, R.W. 1964 Radiation stresses in water waves; a physical discussion, with applications. *Deep Sea Research and Oceanographic Abstracts* **11** (4), 529–562.

LONGUET-HIGGINS, M. S. 1977 The Mean Forces Exerted by Waves on Floating or Submerged Bodies with Applications to Sand Bars and Wave Power Machines. *Proc. R. Soc. Lon. A* **352** (1671), 463–480.

LUBIN, M., DOWSON, O., DIAS GARCIA, F., HUCHETTE, J., LEGAT, B. & VIELMA, J.P. 2023 JuMP 1.0: Recent improvements to a modeling language for mathematical optimization. *Math. Prog. Comp.* **15**, 581–589.

PUCCI, G., BEN AMAR, M. & COUDER, Y. 2015 Faraday instability in floating drops. *Phys. of Fluids* **27** (9), 091107.

RHEE, E., HUNT, R., THOMSON, S. & HARRIS, D.M. 2022 Surferbot: a wave-propelled aquatic vibrobot. *Bioinspir. and Biomim.* **17** (5), 055001.

ROH, C. & GHARIB, M. 2019 Honeybees use their wings for water surface locomotion. *Proc. Natl Acad. Sci. USA* **116** (49), 24446–24451.

STEINMANN, T., ARUTKIN, M., COCHARD, P., RAPHAËL, E., CASAS, J. & BENZAQUEN, M. 2018 Unsteady wave pattern generation by water striders. *J. Fluid Mech.* **848**, 370–387.

STEINMANN, T., CRIBELLIER, A. & CASAS, J. 2021 Singularity of the water strider propulsion mechanisms. *J. Fluid Mech.* **915**, A118.

TUCK, E. O. & LAZAUSKAS, L. 1998 Optimum hull spacing of a family of multihulls. *Ship Technol. Res.-Schiffstechnik* **45** (4), 180.

YUAN, Z.M, CHEN, M., JIA, L., JI, C. & INCECIK, A. 2021 Wave-riding and wave-passing by ducklings in formation swimming. *J. Fluid Mech.* **928**, R2.