# PMC-LLaMA: Further Finetuning LLaMA
on Medical Papers

Chaoyi Wu1,2, Xiaoman Zhang1,2, Ya Zhang1,2, Yanfeng Wang1,2, Weidi Xie1,2,🖂
{wtzxxxwcy02, xm99sjtu, ya_zhang, wangyanfeng<EMAIL_ADDRESS>
1Cooperative Medianet Innovation Center, Shanghai Jiao Tong University
2Shanghai AI Laboratory

###### Abstract

Large Language Models (LLMs) have showcased remarkable capabilities in natural
language understanding in various domains. These models can usually behave
well on daily dialog, or question answering scenarios, however, in areas that
value precision, for example, in medical applications, they often exhibit
unsatisfactory performance due to a lack of domain-specific knowledge. In this
report, we introduce PMC-LLaMA, an open-source language model that is acquired
by fine-tuning an open-source language model on a total of 4.8 million
biomedical academic papers for further injecting medical knowledge, enhancing
its capability in medical domain. Our preliminary evaluations are conducted on
three biomedical QA datasets, including PubMedQA, MedMCQA, and USMLE, showing
that the our model after finetuning, i.e., PMC-LLaMA, demonstrates better
understanding of biomedical domain-specific concepts, thus achieving high
performance on QA benchmarks. The model and codes, along with an online demo,
are publicly available111Huggingface Page: https://huggingface.co/chaoyi-
wu/PMC_LLAMA_7B,222Github Page: https://github.com/chaoyi-wu/PMC-LLaMA.

##  1 Introduction

The rapid advancement of large language models (LLMs), for example, OpenAI’s
ChatGPT (cha, 2023) and GPT-4 (OpenAI, 2023) has truly revolutionized
artificial intelligence in various domain, for example, natural language
processing, computer vision, and biomedical domain (Moor et al., 2023; Nori et
al., 2023; Singhal et al., 2022), unfortunately, the training details and
model architectures for ChatGPT, and its variants, still remain unclear. In
contrast, as an open-source foundational language model, LLaMA (Touvron et
al., 2023) often performs poorly on applications that require heavy domain
knowledge, which we conjecture is due to inadequate domain-specific data at
the model pre-training stage.

In the recent literature, there has been a growing interest in leveraging
open-source LLMs and adapting them towards specific applications or domains.
For instance, Alpaca (Taori et al., 2023) and Vicuna (Chiang et al., 2023)
have concentrated on enhancing the model’s interactive capabilities using
machine-generated instruction-following samples. In contrast to their goal on
daily dialogue, in this report, our focus is to steer the foundational
language model towards medical-specific corpus, by further injecting domain
knowledge into one pre-trained LLaMA.

Specifically, we introduce PMC-LLaMA, an open-source language model trained by
fine-tuning LLaMA-7B on 4.8 million biomedical academic papers. The whole
procedure can be seen in Figure 1. With preliminary evaluation, PMC-LLaMA has
demonstrated superior performance on various medical QA datasets, including
PubMedQA, MedMCQA, and USMLE, for end-to-end full fine-tuning, parameter-
efficient fine-tuning and data-efficient fine-tuning. We believe that such a
medical-specific foundational language model would be more suitable for
specialization in various healthcare sub-tasks, such as medical dialogue or
consultation.

Figure 1: PMC-LLaMA Training Pipeline.

##  2 Experiment

In this section, we begin by outlining the fine-tuning procedure (Sec. 2.1),
including the dataset used and specific training details; Subsequently, in
Sec. 2.2, we provide a detailed description of the evaluation benchmark,
composing three QA datasets, namely PubMedQA, MedMCQA, and UMLSE. Lastly, in
Sec. 2.3, we present three evaluation scenarios: full fine-tuning, parameter-
efficient fine-tuning, and data-efficient fine-tuning.

###  2.1 Fine-tuning Procedure

Dataset.  We start with the S2ORC (Lo et al., 2020) Datasets with $81.1$M
English-language academic papers, and filter them with PubMed Central
(PMC)-id. As a result, there are around $4.9$M papers left, that are highly
related to medical knowledge totaling over $75$B tokens.

Fine-tuning detail.  We fine-tune the LLaMA-7B model on these open-accessed
PMC papers, by optimising an autoregressive generation objective, as
introduced in GPT2 (Radford et al., 2019). Specifically, during our
finetuning, the max context length is set as 512, with a batch size to be 128,
and the model is trained with AdamW optimizer (Loshchilov and Hutter, 2017)
with learning rate 2e-5. To accelerate the training, we adopt the Fully
Sharded Data Parallel (FSDP) acceleration strategy and bf16 (Brain Floating
Point) data format. The model is trained for 5 epochs with 8 A100 GPUs in
around 7 days. Note that, in each epoch, we randomly sample 512 continuous
tokens per paper for training.

###  2.2 Benchmark Descriptions

In this section, we detail the three QA benchmarks for evaluation in the
following.

PubMedQA (Jin et al., 2019) contains questions on biomedical research, the
model is provided with paper abstracts from PubMed, and is required to
complete multiple-choice questions. It is split into three subsets: labeled
(PQA-L), unlabeled (PQA-U), and artificially generated (PQA-A). We use PQA-A
for training, which contains 211,269 Question-Answer pairs, and PQA-L for
testing, which contains 1000 pairs.

MedMCQA (Pal et al., 2022) is a dataset of multiple choice questions sourced
from mock exams and past exams of two Indian medical school entrance exams
called AIIMS and NEET-PG (Pal et al., 2022). The train set contains 182,822
questions, and the test set contains 4183 questions. Each question has 4
choices.

USMLE (Jin et al., 2021) is a dataset of multiple choice questions (4 choices
per question), based on the United States Medical License Exams. The dataset
is collected from the professional medical board exams, covering three
languages: English, simplified Chinese, and traditional Chinese, and contains
12,724, 34,251, and 14,123 questions for these three languages, respectively.
We only use the English parts and split it into 10,178 questions for training,
1273 for validation, and 1273 for testing, following the official data spits.

###  2.3 Evaluation Scenario

In this section, we evaluate PMC-LLaMA by fine-tuning on the above-mentioned
corresponding medical QA datasets, under three scenarios, namely full fine-
tuning, parameter-efficient fine-tuning, and data-efficient fine-tuning as
detailed in the following.

Full fine-tuning. We fine-tune PMC-LLAMA on the combination of the training
set from PubMedQA and MedMCQA, for evaluation, we treat the test set from
PubMedQA and MedMCQA as in-domain (ID) evaluation, while treating the USMLE
test set as out-of-domain (OOD) evaluation, similarly to LMFlow (Diao et al.,
2023).

Parameter-efficient fine-tuning (Mangrulkar et al., 2022) enables efficient
adaptation of pre-trained LLMs to various downstream applications without
fine-tuning all parameters, greatly reducing the time and computation cost.
Here we use the most widely used method, called PEFT Low-Rank Adaptation
(LoRA) (Hu et al., 2021). The hyper-parameter for LoRA is set by default
provided in the python package PEFT (Mangrulkar et al., 2022). The data
settings for training and testing remain the same as full fine-tuning.

Data-efficient fine-tuning. In addition to measuring the number of training
parameters, we also conduct experiments on the data side, i.e., data-efficient
fine-tuning. Specifically, we conduct experiments on USMLE dataset, _i.e._ ,
simply train and test on its own splits, resembling an USMLE ID evaluation.
However, as the training set in USMLE is much smaller compared with that of
combining PubMedQA and MedMCQA ($10k$ vs. $400k$), less medical knowledge can
be learned from the downstream training data, it turns out to be rather more
challenging than the previous OOD test, as will be shown by our experimental
results.

Implementation details. For each downstream task, we fine-tune the network
using the same setting as training on PMC papers, specifically, we set the max
context length as 512 and fine-tune the models with AdamW optimizer with
learning rate 2e-5, the batch size is set to be 128, and the model is trained
for 3 epochs.

##  3 Results

In Fig.2 we show the fine-tuning training curve under different settings and
compare the convergence speed. The final evaluation results are shown in
Tab.1, generally speaking, in almost all medical QA benchmarks, our proposed
PMC-LLaMA converges faster and exhibits better performance than the original
LLaMA (Touvron et al., 2023).

Model convergence. In Fig.2, we name each training curve using the convention
{Model Name}-{Training Setting}-{Dataset}”, for example, {LLaMA-7B}-{Full-
Finetune}-{PM&MedMC}” refers to the full fine-tuning of the LLaMA-7B model on
the PubMedQA and MedMCQA training set. As illustrated in the figure,
generally, PMC-LLaMA converges more quickly and achieves lower loss compared
to the original LLaMA, especially when the training dataset is larger (with
PM&MedMC) and has more trainable parameters (Full-Finetune). This suggests
that PMC-LLaMA offers a better initialization for medical tasks.

Figure 2: The training curve compared with LLaMA under two settings with
different training data. The curve is smoothed with $0.6$ Exponential Moving
Average. Depending on the amount of training data, the curves have different
ending steps, thus showing seesaw shapes with different periods.

[b]  Method Setting USMLE(OOD/ID) MedMCQA(ID) PubMedQA(ID) Human (pass)
Manual* 50.0 – 60.0 Human (expert) 87.0 90.0 78.0 InstructGPT-175B Zero-shot*
46.0 44.0 73.2 ChatGPT 57.0 44.7 63.9 LLaMA-7B 27.1 24.3 5.2 LLaMA-33B 43.4
30.3 1.8 LLaMA-7B${}_{\text{Full}}$ Full fine-tuning 44.55/35.66 48.15 73.41
PMC-LLaMA-7B${}_{\text{Full}}$ 44.70/40.61 50.54 69.53
LLaMA-7B${}_{\text{PEFT}}$ PEFT 29.38/27.34 32.37 65.81 PMC-
LLaMA-7B${}_{\text{PEFT}}$ 30.64/28.52 34.33 68.23

Table 1: Comparison between LLaMA-7B and PMC-LLaMA-7B under different
settings. ACC in percentages is reported in the table. Note that, the manual
and zero-shot results with * are referred from LMFlow (Diao et al., 2023).

Full fine-tuning results. As shown in Tab. 1, PMC-LLaMA-7B${}_{\text{Full}}$
exceeds the LLaMA-7B${}_{\text{Full}}$ in two of three test sets, improving
the performance from $44.55\%$ to $44.70\%$ and $35.66\%$ to $40.61\%$ on
USMLE for OOD and ID settings respectively, and $48.15\%$ to $50.74\%$ on
MedMCQA.

Parameter-efficient fine-tuning results. As shown in Tab. 1, PMC-
LLaMA-7B${}_{\text{PEFT}}$ demonstrates superior performance than
LLaMA-7B${}_{\text{PEFT}}$, particularly on the in-domain datasets, $1.22\%$
improvement on USMLE, $1.96\%$ improvement on MedMCQA and $2.42\%$ on
PubMedQA. These results demonstrate that the original LLaMA only provides
suboptimal embedding spaces for Medical QA and further fine-tuning on the
biomedical corpus is beneficial for model domain adaptation.

Data-efficient fine-tuning results. We refer to USMLE ID test results as the
Data-efficient fine-tuning comparison. As shown in Tab. 1, whatever the
training setting is “Full fine-tuning” or “PEFT”, PMC-LLaMA-7B can both
achieve better results on the USMLE ID test. Specifically, we improve the ACC
scores from $35.66\%$ to $40.61\%$ under full fine-tuning and from $27.34\%$
to $28.52\%$ under PEFT.

Compare to ChatGPT and its variants. As is evidently shown in Tab.1,
InstructGPT and chatGPT achieve fantastic results even under zero-shot
evaluation, however, we argue that, as they are not open-sourced, and do not
reveal any training details, there is no guarantee on whether the evaluation
data has been exploited for training their model, hence their results can only
be used as a reference here, and should not be taken for a fair comparison.

##  4 Zero-shot Case Study

In this section, we show zero-shot generation cases in Fig.3 with the
judgement from GPT-4, and present the differences between original LLaMA and
PMC-LLaMA. The GPT-4 Judgement is inspired by (Chiang et al., 2023), and uses
the following prompt:

Considering the two outputs, which one do you think is more related to the
input? Input: …PMC-LLaMA Output: …Original LLaMA Output: …

Figure 3: Zero-shot evalution on medical-related sentence completion, we
compare the results between PMC-LLaMA and original LLaMA-7B with the judgement
provided by GPT-4.

Generally speaking, we find PMC-LLaMA has more expertise in medical knowledge
and responds better to some cutting-edge professional medical concepts.

Case 1: COPD. According to World Health Organization (WHO111WHO
https://www.who.int), “Chronic obstructive pulmonary disease (COPD) is a
common and treatable disease characterized by progressive airflow limitation
and tissue destruction and is the third leading cause of death worldwide”,
which is consistent with the output of PMC-LLaMA, in contrast, although LLaMA
has pointed out the fact that COPD can be identified by spirometry results,
quickly lose the attention on COPD and generate irrelevant sentences. In GPT-4
Judgement, it also points out that the original LLaMA “goes on to discuss
unrelated topics”.

Case 2: Robotic Cardiac Surgery. According to Johns Hopkins Medicine222Johns
Hopkins Medicine: https://www.hopkinsmedicine.org/, Robotic cardiac surgery
involves performing heart surgery through very small incisions in the chest.
Utilizing miniature instruments and robot-assisted tools, surgeons can conduct
heart surgery in a manner that is significantly less invasive than traditional
open-heart surgery”. The assertion made by PMC-LLaMA closely aligns with this
statement, while the output from the original LLaMA appears to be unclear. Its
phrase through the placement of surgical robotic tools” is inaccurate and
should instead refer to “robot-assisted tools”. Although GPT-4 considered the
original LLaMA output to be satisfactory, it failed to detect the inaccuracies
in LLaMA’s descriptions, focusing merely on their relevance.

Case 3 & 4: Diabetes and Pneumonia. In cases 3 and 4, PMC-LLaMA consistently
provides more context about the input, suggesting that PMC-LLaMA possesses a
deeper understanding of medical background knowledge. In contrast, the
original LLaMA tends to introduce irrelevant content in its output. As GPT-4
notes, the rest of the text is less clear and informative compared to the PMC-
LLaMA output ”and the response goes on to discuss PCT concentrations, which
are not directly related to the input statement”.

##  5  Conclusion and Future Work

In conclusion, in this report, we investigate the existing open-source
language model, namely, LLaMA, on medical applications, showing that the model
performs unsatisfactorily on QA tasks. To inject domain knowledge into the
pre-trained model, We conduct a preliminary investigation by fine-tuning LLaMA
with 4.8 million medical papers, resulting in a new foundation model for the
medical domain PMC-LLaMA. Experimental evaluation demonstrates that PMC-LLaMA
is more suitable for medical tasks compared with LLaMA and after instruction
tuning it can get better performance compared to other models instruction
tuned from LLaMA in the medical domain.

One limitation of PMC-LLaMA is that we have only trained 5 epochs in the
current version and have not seen every token in the 4.8 million papers. In
future work, we will continuously train PMC-LLaMA and update our base model on
the hugging face page and gradually train PMC-LLaMA models with more
parameters.

## References

  * cha [2023] Openai. introducing chatgpt.  https://openai.com/blog/chatgpt/, 2023.
  * Chiang et al. [2023] Wei-Lin Chiang, Zhuohan Li, Zi Lin, Ying Sheng, Zhanghao Wu, Hao Zhang, Lianmin Zheng, Siyuan Zhuang, Yonghao Zhuang, Joseph E. Gonzalez, Ion Stoica, and Eric P. Xing.  Vicuna: An open-source chatbot impressing gpt-4 with 90%* chatgpt quality, March 2023.  URL https://vicuna.lmsys.org.
  * Diao et al. [2023] Shizhe Diao, Rui Pan, Hanze Dong, KaShun Shum, Jipeng Zhang, Wei Xiong, and Tong Zhang.  Lmflow: An extensible toolkit for finetuning and inference of large foundation models.  https://optimalscale.github.io/LMFlow/, 2023.
  * Hu et al. [2021] Edward J Hu, Yelong Shen, Phillip Wallis, Zeyuan Allen-Zhu, Yuanzhi Li, Shean Wang, Lu Wang, and Weizhu Chen.  Lora: Low-rank adaptation of large language models.  _arXiv preprint arXiv:2106.09685_ , 2021.
  * Jin et al. [2021] Di Jin, Eileen Pan, Nassim Oufattole, Wei-Hung Weng, Hanyi Fang, and Peter Szolovits.  What disease does this patient have? a large-scale open domain question answering dataset from medical exams.  _Applied Sciences_ , 11(14):6421, 2021.
  * Jin et al. [2019] Qiao Jin, Bhuwan Dhingra, Zhengping Liu, William W Cohen, and Xinghua Lu.  Pubmedqa: A dataset for biomedical research question answering.  _arXiv preprint arXiv:1909.06146_ , 2019.
  * Lo et al. [2020] Kyle Lo, Lucy Lu Wang, Mark Neumann, Rodney Kinney, and Daniel Weld.  S2ORC: The semantic scholar open research corpus.  In _Proceedings of the 58th Annual Meeting of the Association for Computational Linguistics_ , pages 4969–4983, Online, July 2020. Association for Computational Linguistics.  doi: 10.18653/v1/2020.acl-main.447.  URL https://www.aclweb.org/anthology/2020.acl-main.447.
  * Loshchilov and Hutter [2017] Ilya Loshchilov and Frank Hutter.  Decoupled weight decay regularization.  _arXiv preprint arXiv:1711.05101_ , 2017.
  * Mangrulkar et al. [2022] Sourab Mangrulkar, Sylvain Gugger, Lysandre Debut, Younes Belkada, and Sayak Paul.  Peft: State-of-the-art parameter-efficient fine-tuning methods.  https://github.com/huggingface/peft, 2022.
  * Moor et al. [2023] Michael Moor, Oishi Banerjee, Zahra Shakeri Hossein Abad, Harlan M Krumholz, Jure Leskovec, Eric J Topol, and Pranav Rajpurkar.  Foundation models for generalist medical artificial intelligence.  _Nature_ , 616(7956):259–265, 2023.
  * Nori et al. [2023] Harsha Nori, Nicholas King, Scott Mayer McKinney, Dean Carignan, and Eric Horvitz.  Capabilities of gpt-4 on medical challenge problems.  _arXiv preprint arXiv:2303.13375_ , 2023.
  * OpenAI [2023] OpenAI.  Gpt-4 technical report, 2023.
  * Pal et al. [2022] Ankit Pal, Logesh Kumar Umapathi, and Malaikannan Sankarasubbu.  Medmcqa: A large-scale multi-subject multi-choice dataset for medical domain question answering.  In _Conference on Health, Inference, and Learning_ , pages 248–260. PMLR, 2022.
  * Radford et al. [2019] Alec Radford, Jeffrey Wu, Rewon Child, David Luan, Dario Amodei, Ilya Sutskever, et al.  Language models are unsupervised multitask learners.  _OpenAI blog_ , 1(8):9, 2019.
  * Singhal et al. [2022] Karan Singhal, Shekoofeh Azizi, Tao Tu, S Sara Mahdavi, Jason Wei, Hyung Won Chung, Nathan Scales, Ajay Tanwani, Heather Cole-Lewis, Stephen Pfohl, et al.  Large language models encode clinical knowledge.  _arXiv preprint arXiv:2212.13138_ , 2022.
  * Taori et al. [2023] Rohan Taori, Ishaan Gulrajani, Tianyi Zhang, Yann Dubois, Xuechen Li, Carlos Guestrin, Percy Liang, and Tatsunori B. Hashimoto.  Stanford alpaca: An instruction-following llama model.  https://github.com/tatsu-lab/stanford_alpaca, 2023.
  * Touvron et al. [2023] Hugo Touvron, Thibaut Lavril, Gautier Izacard, Xavier Martinet, Marie-Anne Lachaux, Timothée Lacroix, Baptiste Rozière, Naman Goyal, Eric Hambro, Faisal Azhar, et al.  Llama: Open and efficient foundation language models.  _arXiv preprint arXiv:2302.13971_ , 2023.


# Interfacial waves in the presence of wind and current: an asymptotic study

A. F. Bonfils  Nordita, KTH Royal Institute of Technology and Stockholm
University, Stockholm 10691, Sweden  Dhrubaditya Mitra  Nordita, KTH Royal
Institute of Technology and Stockholm University, Stockholm 10691, Sweden  W.
Moon  Department of Environmental Atmospheric Sciences, Pukyong National
University, 48513 Pusan, South Korea  J. S. Wettlaufer  Nordita, KTH Royal
Institute of Technology and Stockholm University, Stockholm 10691, Sweden
Yale University, New Haven, Connecticut 06520-8109, USA

###### Abstract

We use asymptotic methods to study the evolution of short wavelength
interfacial waves driven by the combined action of wind and current. We solve
the Rayleigh equation for the stability of the shear flow, and construct a
uniformly valid approximation for the perturbed streamfunction, or
eigenfunction. We then expand the real part of the eigenvalue, the phase
speed, in a power series of the inverse wavenumber and show that the imaginary
part is exponentially small. We give expressions for the growth rates of the
Miles (1957) and rippling (e.g., Young and Wolfe, 2013) instabilities that are
valid for an arbitrary shear flow. The accuracy of the results is demonstrated
by a comparison with the exact solution of the eigenvalue problem in the case
when both the wind and the current have an exponential profile.

##  I Introduction

Figure 1: Schematic of wind driving a current in the water. Both are modeled
as a parallel shear flow $U=U(z)$ and perturbed by a wave with wavenumber $k$.

Waves at the interface between two fluids with different densities are
ubiquitous in nature. A natural question concerns how a flow in either fluid
affects these waves? Here, we consider fluid layers of infinite extent. A
canonical example is the wind flowing over the ocean, both of which can be
modeled as parallel flows of the form $\boldsymbol{U}=U(z)\
\boldsymbol{\hat{x}}$, with $z$ the vertical coordinate and
$\boldsymbol{\hat{x}}$ a horizontal unit vector (Figure 1). The stability of
$\boldsymbol{U}$ under the influence of small two-dimensional perturbations
has been studied extensively over the last seventy years. In the absence of a
current in the water, Miles (1957) found an instability of the wind leading to
the growth of water waves. The theory of Miles is inviscid and assumes that
the wind and the waves are weakly coupled because the air/water density ratio,
$r$, is small. Hence, the wind transfers energy to the waves within a critical
layer in the air, where the wind speed equals the phase speed of free surface
waves. The growth rate can be calculated analytically in only a few cases
(Bonfils  _et al._ , 2022).

When there is a current in the water, the position of the critical layer
responsible for the Miles instability is unknown because the phase speed is
itself unknown. This is associated with the fact that the surface waves are no
longer free in the sense that the current modifies their dispersion relation
in a non-trivial manner. Water currents are often wind-induced and decay with
depth. In such cases, Young and Wolfe (2013) showed that sheared surface waves
propagating slower than the water surface, which is dragged by the wind,
undergo the so-called rippling instability. In fact, there is another critical
layer in the water, where the unknown phase speed of the waves matches the
speed of the current. Analytical progress on the rippling instability in deep
water has only been made for piece-wise linear or exponential currents; see
Young and Wolfe (2013) for a review. Finally, Kadam  _et al._ (2022) gave an
exact analytical treatment of the stability of an exponential wind profile
over a finite-depth water layer that is either quiescent or has linear or
quadratic current profiles. For the quadratic current, they introduced
spheroidal wave functions to assess the stability of a shear flow. Here,
however, we obtain general results on both the Miles and rippling
instabilities in the framework of asymptotic analysis.

In the absence of a current in the water and for $r$ less than unity, and yet
not small, the Miles instability is difficult to study. Indeed, the waves are
then strongly coupled to the wind so that the shear substantially changes the
dispersion relation, yielding yet another situation in which the position of
the critical layer is unknown. Not only is such a strong wind–wave coupling a
central process on Earth, but it may play an important role in a number of
astrophysical settings, including white dwarfs, one of the end products of
stellar evolution, neutron stars and black holes (Shapiro and Teukolsky,
2008). For example, although most recent high-resolution three-dimensional
simulations (Casanova  _et al._ , 2011) suggest that Kelvin-Helmholtz
instabilities driven by buoyant fingering may be able to explain the
composition of white dwarfs, an alternative proposal by Rosner  _et al._
(2001) and Alexakis  _et al._ (2004a), is that the Miles instability is
responsible.

Because most ocean waves have wavelengths much larger than the characteristic
length scale of the wind profile, the Miles instability can be treated using
asymptotic methods in the long-wave-limit (Bonfils  _et al._ , 2022). Here, we
focus on short waves with the goal of capturing the combined influences of an
underlying current and a moderate density ratio. Whereas short waves may not
be central in terrestrial oceanography, Alexakis  _et al._ (2004b) argued that
they are at the core of the astrophysical setting described in the previous
paragraph. Thus, we treat a sheared two-fluid interface where the upper fluid
is less dense than the lower fluid. For clarity of discussion, we refer to the
upper and lower fluids as air and water, respectively, and refer to wind as
the shear flow in the air, and current as the shear flow in the water.

In §II, we describe the linear stability analysis of a sheared two-fluid
interface. The eigenfunction satisfies the Rayleigh equation and the
eigenvalue is a complex intrinsic phase speed. We focus on the case of a wind-
induced current in §III, and obtain the real part of the eigenvalue as a
series in powers of the dimensionless inverse wavenumber. Moreover, we show
that the imaginary part of the eigenvalue is small and obtain general formulae
for the growth rates of the Miles and the rippling instabilities. In §IV, we
treat the case of the wind and the current governed by an exponential profile
and compare our asymptotic results to the exact eigenvalue. In §V, we treat
another type of current wherein a mode can have two critical layers, one in
the air and one in the water. We then derive in §VI the asymptotic solution of
the Rayleigh equation that was used in §III for the calculation of the
eigenvalue. In particular, we show that an internal boundary layer emerges
from the singularity at the critical level. Finally, before concluding, in
§VII we demonstrate the limitations of the Wentzel-Kramers-Brillouin (WKB)
approach to solving the Rayleigh equation.

##  II Linear stability of a sheared two-fluid interface

Young and Wolfe (2013) derived the eigenvalue problem associated with the
linear stability of an inviscid parallel shear flow across a two-fluid
interface as sketched in Fig. 1. First we outline the governing Rayleigh
equation and boundary conditions, after which we describe the non-
dimensionalization of the problem.

###  II.1 Eigenvalue problem

We consider a parallel shear flow, $U=U(z)$, monotonic in both air and water.
The air to water density ratio is $r\equiv\rho_{\rm{a}}/\rho_{\rm{w}}<1$, and
the gravitational acceleration and surface tension are $g$ and $\sigma$
respectively. Incompressibility ensures that a perturbation of the flow is
entirely determined by the streamfunction $\psi=\psi(x,z,t)$, where $t$ is the
time. We use normal modes in the form

$\psi(x,z,t)=\Re\big{\\{}\hat{\psi}(z)\ e^{ik(x-ct)}\big{\\}},$ (1)

where $k$ is a real wavenumber and $c$ a complex phase speed to be determined
(the eigenvalue), conservation of vorticity yields the Rayleigh equation as

$\big{[}U(z)-c\big{]}\big{[}\hat{\psi}^{\prime\prime}(z)-k^{2}\hat{\psi}(z)\big{]}-U^{\prime\prime}(z)\hat{\psi}(z)=0.$
(2)

We require the function $U(z)$ to be continuous at $z=0$, which excludes a
Kelvin-Helmholtz type of instability and ensures the continuity of
$\hat{\psi}(z)$ (Drazin and Reid, 1981). Note that the derivative
$U^{\prime}(z)$ may have a finite jump at $z=0$ and the perturbation must
decay in the far field. Provided that

$\lim\limits_{z\to\pm\infty}\frac{U^{\prime\prime}(z)}{U(z)-c}=0,$ (3)

we can impose an exponential decay, viz.,

$\lim\limits_{z\to\pm\infty}\psi^{\prime}(z)\pm k\psi(z)=0.$ (4)

Finally, we impose continuity of the normal stress at $z=0$ and require the
air-water interface to be a streamline, which yields the boundary condition

$\frac{\sigma}{\rho_{\rm{w}}}\
k^{2}+(1-r)g+r\Big{[}(U_{\rm{S}}-c)^{2}\frac{\hat{\psi}^{\prime}(0^{+})}{\hat{\psi}(0)}-(U_{\rm{S}}-c)U^{\prime}(0^{+})\Big{]}-(U_{\rm{S}}-c)^{2}\frac{\hat{\psi}^{\prime}(0^{-})}{\hat{\psi}(0)}+(U_{\rm{S}}-c)U^{\prime}(0^{-})=0,$
(5)

where $U_{\rm{S}}=U(z=0)$ is the surface drift.

###  II.2 Non-dimensionalization

The length and velocity scales, $L$ and $V$, differ in the air and water,
which we denote with the subscripts $\rm{a}$ and $\rm{w}$ respectively.
Whereas when there is no current in the water, we use the scales of the wind,
$L_{\rm{a}}$ and $V_{\rm{a}}$, to non-dimensionalize the problem, because
water is the primary medium of surface waves we use $L_{\rm{w}}$ and
$V_{\rm{w}}$ when a current is present. Thus, the ratios

$\mathcal{R}_{1}\equiv\frac{V_{\rm{a}}}{V_{\rm{w}}}\qquad\text{and}\qquad\mathcal{R}_{2}\equiv\frac{L_{\rm{a}}}{L_{\rm{w}}}$
$None$

act as additional control parameters. In a frame moving at speed $U_{\rm{S}}$,
we use the dimensionless profile

$\mathcal{U}\equiv\frac{U-U_{\rm{S}}}{V_{\rm{w}}}\
,\qquad\text{where}\qquad\mathcal{U}(0)=0,$ (7)

and the dimensionless intrinsic phase speed is

$\mathcal{C}\equiv\frac{c-U_{\rm{S}}}{V_{\rm{w}}}\ .$ (8)

The dimensionless wavenumber and vertical coordinate are $\mathcal{k}\equiv
kL_{\rm{w}}$ and $\mathcal{z}\equiv z/L_{\rm{w}}$, respectively, and we define
the dimensionless gravity and surface tension as

$\mathcal{G}\equiv\frac{gL_{\rm{w}}}{V_{\rm{w}}^{2}}\qquad\text{and}\qquad\mathcal{S}\equiv\frac{\sigma}{\rho_{\rm{w}}V_{\rm{w}}^{2}L_{\rm{w}}}\
.$ (9)

From equation (4), the far field behavior of the streamfunction has the form
$e^{\pm\mathcal{k}\mathcal{z}}$. Short waves are characterized by
$\mathcal{k}\gg 1$, so the exponential decay is captured as follows,

$\frac{\hat{\psi}(z)}{\hat{\psi}(0)}\equiv\begin{cases}e^{-\mathcal{k}\mathcal{z}}f(\mathcal{z})\qquad\text{if}\qquad\mathcal{z}>0,\\\
e^{\mathcal{k}\mathcal{z}}\
h(\mathcal{z})\qquad\text{if}\qquad\mathcal{z}<0.\end{cases}$ (10)

We introduce the small parameter $\epsilon\equiv 1/\mathcal{k}$ and use (10)
in the Rayleigh equation (2), which gives

$\displaystyle\epsilon
f^{\prime\prime}(\mathcal{z})-2f^{\prime}(\mathcal{z})-\epsilon\
\frac{\mathcal{U}^{\prime\prime}(\mathcal{z})}{\mathcal{U}(\mathcal{z})-\mathcal{C}}\
f(\mathcal{z})$ $\displaystyle=0,\quad f(0)=1,$ (11)
$\displaystyle\text{and}\qquad\epsilon
h^{\prime\prime}(\mathcal{z})+2h^{\prime}(\mathcal{z})-\epsilon\
\frac{\mathcal{U}^{\prime\prime}(\mathcal{z})}{\mathcal{U}(\mathcal{z})-\mathcal{C}}\
h(\mathcal{z})$ $\displaystyle=0,\quad h(0)=1.$ (12)

We assume that

$f(\mathcal{z})=O(1)\quad\text{as}\quad\mathcal{z}\to+\infty\quad\text{and}\quad
h(\mathcal{z})=O(1)\quad\text{as}\quad\mathcal{z}\to-\infty,$ (13)

the veracity of which we check a posteriori. The dimensionless form of
equation (5) is

$\frac{\mathcal{S}}{\epsilon}+\mathcal{G}(1-r)\epsilon+\big{[}r\mathcal{U}^{\prime}(0^{+})-\mathcal{U}^{\prime}(0^{-})\big{]}\epsilon\mathcal{C}-(1+r)\mathcal{C}^{2}+\big{[}rf^{\prime}(0^{+})-h^{\prime}(0^{-})\big{]}\epsilon\mathcal{C}^{2}=0.$
(14)

For a given profile $\mathcal{U}(\mathcal{z})$, our main task is to solve
equations (11) and (12) subject to the boundary conditions (13) and (14).

We consider the canonical situation where the wind blows over the water in
which it induces a current (fig. 1). Hence,
$\mathcal{U}^{\prime}\big{(}\mathcal{z}\lessgtr 0\big{)}>0$,
$\mathcal{U}(\mathcal{z}>0)>0$, and $\mathcal{U}(\mathcal{z}<0)<0$. We explore
another situation in §V.

###  II.3 Examples of profiles

A typical example of a wind and a wind-induced current is the double
exponential profile,

$U(z)=\begin{cases}U_{\infty}+(U_{\rm{S}}-U_{\infty})\
e^{-z/d_{\rm{a}}}\qquad\text{if}\qquad z>0,\\\ U_{\rm{S}}\
e^{z/d_{\rm{w}}}\qquad\text{if}\qquad z<0,\end{cases}$ (15)

where $U_{\infty}$ is the free-stream air velocity, and $d_{\rm{a}}$ and
$d_{\rm{w}}$ denote the thicknesses of the air and water shear boundary layers
respectively. In the frame of the water surface, the dimensionless form of Eq.
(15) is

$\mathcal{U}(\mathcal{z})=\begin{cases}(\mathcal{R}_{1}-1)(1-\
e^{-\mathcal{z}/\mathcal{R}_{2}})\qquad\text{if}\qquad\mathcal{z}>0,\\\
e^{\mathcal{z}}-1\qquad\text{if}\qquad\mathcal{z}<0,\end{cases}$ (16)

with $\mathcal{R}_{1}=U_{\infty}/U_{\rm{S}}$ and
$\mathcal{R}_{2}=d_{\rm{a}}/d_{\rm{w}}$. Young and Wolfe (2013) showed that
for such a profile the eigenvalue problem can be solved exactly in terms of
hypergeometric functions.

In the context of physical oceanography, the double log profile,

$U(z)=\begin{cases}U_{\rm{S}}+u_{\star\rm{a}}\ln(1+z/z_{0\rm{a}})/\kappa\qquad\text{if}\qquad
z>0,\\\
U_{\rm{S}}-u_{\star\rm{w}}\ln(1-z/z_{0\rm{w}})/\kappa\qquad\text{if}\qquad
z<0,\end{cases}$ (17)

may be more realistic (Wu, 1975), where $u_{\star\rm{a}}$ and
$u_{\star\rm{w}}$ are the friction velocities of air and water, respectively,
$\kappa=0.4$ is the von Kármán constant, and $z_{0\rm{a}}$ and $z_{0\rm{w}}$
are air and water roughness lengths, accounting for the presence of waves.
Note that the velocity of the logarithmic current is negative for
$z<z_{\rm{min}}$, where

$z_{\rm{min}}\equiv-z_{0\rm{w}}\big{(}e^{\kappa
U_{\rm{S}}/u_{\star\rm{w}}}-1\big{)}.$ (18)

Such a change of sign is unphysical so we take $U(z)=0$ for $z<z_{\rm{min}}$.
The dimensionless form of (17) in the frame of the water surface is

$\mathcal{U}(\mathcal{z})=\begin{cases}\mathcal{R}_{1}\ln(1+\mathcal{z}/\mathcal{R}_{2})/\kappa\qquad\text{if}\qquad\mathcal{z}>0,\\\
-\ln(1-\mathcal{z})/\kappa\qquad\text{if}\qquad\mathcal{z}<0,\end{cases}$ (19)

with $\mathcal{R}_{1}=u_{\star\rm{a}}/u_{\star\rm{w}}$ and
$\mathcal{R}_{2}=z_{0\rm{a}}/z_{0\rm{w}}$. In this case, exact analytical
solutions are unknown.

##  III Short wavelength expansions and exponential asymptotics

We treat the eigenvalue problem described in §II.1 perturbatively, where the
small parameter is the dimensionless inverse wavenumber, $\epsilon\ll 1$. We
draw intuition for the approach from Miles (1957), who considered a simpler
version of our problem, with the small parameter $r\ll 1$ and in the absence
of a current. Hence, the leading order eigenvalue corresponding to $r=0$ is
real and equal to $c_{s}$, the phase speed of free surface waves. Moreover,
due to the critical layer at height $z_{c}$, such that $U(z_{c})=c_{s}$, the
eigenvalue has an imaginary part of order $r$ as well as real corrections,
also of order $r$.

In contrast to the treatment of Miles (1957), the leading order eigenvalue is
unknown. In fact, even the nature of the lowest order behavior in $\epsilon$
is murky and hence we only assume that

$\mathcal{C}=\mathcal{C}_{r}(\epsilon)+\rm{i}\
\mathcal{C}_{i}(\epsilon),\qquad\text{with}\qquad\mathcal{C}_{i}(\epsilon)\ll\mathcal{C}_{r}(\epsilon),\qquad\epsilon\to
0.$ (20)

We use the notation of Bender and Orszag (1999) for
“$\mathcal{C}_{i}(\epsilon)$ is much smaller than $\mathcal{C}_{r}(\epsilon)$
as $\epsilon$ tends to $0$”, namely that

$\lim\limits_{\epsilon\to
0}\frac{\mathcal{C}_{i}(\epsilon)}{\mathcal{C}_{r}(\epsilon)}=0.$ (21)

The purpose of such an assumption is to replace $\mathcal{C}\in\mathbb{C}$ by
$\mathcal{C}_{r}\in\mathbb{R}$ in equations (11) and (12). Recalling that
$\mathcal{U}(\mathcal{z}>0)>0$ and $\mathcal{U}(\mathcal{z}<0)<0$, if
$\mathcal{C}_{r}$ is positive (negative) then there is a critical layer in the
air (water), associated with a singular point in Eq. 11 (Eq. 12). However, at
this stage in the development, we do not know the possible values for
$\mathcal{C}_{r}$. Indeed, the term $rf^{\prime}(0^{+})-h^{\prime}(0^{-})$
depends on $\mathcal{C}$ in a non-trivial manner so that Eq. (14) is not
necessarily quadratic in $\mathcal{C}$. Physically, $\mathcal{C}_{r}$ is the
phase speed of sheared interfacial waves in the reference frame of the water
surface. In that frame, a wave propagating in the direction of (against) the
current has a positive (negative) $\mathcal{C}_{r}$ and Young and Wolfe (2013)
refer to these as prograde (retrograde) modes. In the case of a constant (zero
shear) current $U$, there is one prograde and one retrograde mode, which are
simply the Doppler-shifted forward and backward interfacial waves. We
generalize this result to the case of arbitrary shear, which we check a
posteriori. Thus we assume that there are two solutions, $\mathcal{C}_{r+}>0$
and $\mathcal{C}_{r-}<0$ that correspond to two critical levels,
$\mathcal{z}_{\rm{c}+}>0$ and $\mathcal{z}_{\rm{c}-}<0$, such that

$\mathcal{U}(\mathcal{z}_{\rm{c}\pm})=\mathcal{C}_{r\pm}.$ (22)

Therefore, the prograde (retrograde) mode undergo the Miles (rippling)
instability, and hence $\mathcal{C}_{i}\neq 0$ for both modes. We stress that,
provided that the values of $\mathcal{C}_{r\pm}$ are within the bounds of the
function $\mathcal{U}(\mathcal{z})$, critical layers actually exist. When
$\mathcal{C}_{r\pm}$ are equal to these bounds, the system is marginally
stable. In Appendix C, we derive a general asymptotic formula for the large
wavenumber cut-off of the rippling instability.

Hence, although $\mathcal{C}_{r\pm}$ are unknown, we replace $\mathcal{C}$ by
$\mathcal{C}_{r+}$ in Eq. (11) and by $\mathcal{C}_{r-}$ in Eq. (12). We then
solve these equations using boundary layer theory in §VI, where we construct
uniformly valid composite solutions. However, here we need only substitute the
derivatives at $\mathcal{z}=0^{\pm}$ into Eq. (14). In Appendix B, we show
that for the prograde mode,

$\displaystyle f^{\prime}(0^{+})$ $\displaystyle=-i\pi\
\frac{\mathcal{U}^{\prime\prime}(\mathcal{z}_{\rm{c}+})}{|\mathcal{U}^{\prime}(\mathcal{z}_{\rm{c}+})|}\
e^{-2\mathcal{z}_{\rm{c}+}/\epsilon}+\epsilon\
\frac{\mathcal{U}^{\prime\prime}(0^{+})}{2\mathcal{C}_{r+}(\epsilon)}-\frac{\mathcal{U}^{\prime\prime}(\mathcal{z}_{\rm{c}+})}{\mathcal{U}^{\prime}(\mathcal{z}_{\rm{c}+})}\
\frac{\epsilon}{2\mathcal{z}_{\rm{c}+}}\ ,$ (23)
$\displaystyle\text{and}\qquad h^{\prime}(0^{-})$ $\displaystyle=-\epsilon\
\frac{\mathcal{U}^{\prime\prime}(0^{-})}{2\mathcal{C}_{r+}(\epsilon)}\ ,$ (24)

whereas for the retrograde mode,

$\displaystyle f^{\prime}(0^{+})$ $\displaystyle=\epsilon\
\frac{\mathcal{U}^{\prime\prime}(0^{+})}{2\mathcal{C}_{r-}(\epsilon)}\ ,$ (25)
$\displaystyle\text{and}\qquad h^{\prime}(0^{-})$ $\displaystyle=i\pi\
\frac{\mathcal{U}^{\prime\prime}(\mathcal{z}_{\rm{c}-})}{|\mathcal{U}^{\prime}(\mathcal{z}_{\rm{c}-})|}\
e^{2\mathcal{z}_{\rm{c}-}/\epsilon}-\epsilon\
\frac{\mathcal{U}^{\prime\prime}(r^{-})}{2\mathcal{C}_{r-}(\epsilon)}+\frac{\mathcal{U}^{\prime\prime}(\mathcal{z}_{\rm{c}-})}{\mathcal{U}^{\prime}(\mathcal{z}_{\rm{c}-})}\
\frac{\epsilon}{2\mathcal{z}_{\rm{c}-}}\ .$ (26)

We can now solve Eq. (14). With the advent of results (23) and (26), our key
assumption, Eq. (20), is made more precise by letting

$\mathcal{C}_{\pm}=\mathcal{C}_{r\pm}(\epsilon)+\rm{i}\ A_{\pm}(\epsilon)\
\epsilon\ e^{-2|\mathcal{z}_{\rm{c}\pm}|/\epsilon}.$ (27)

Therefore, in order to have solutions with a non-zero imaginary part, we must
include exponentially small terms. When first discarding those terms, we
obtain

$\frac{\mathcal{S}}{\epsilon}+\mathcal{G}(1-r)\epsilon+\big{[}r\mathcal{U}^{\prime}(0^{+})-\mathcal{U}^{\prime}(0^{-})\big{]}\epsilon\mathcal{C}_{r}-(1+r+\rm{h.o.t.})\mathcal{C}_{r}^{2}=0,\qquad\epsilon\ll
1,$ (28)

where h.o.t denotes higher order terms, such as terms of order
$\epsilon/\mathcal{C}_{r}$ and $\epsilon/z_{\rm{c}}$.

We seek solutions of Eq. (28) as a series in powers of $\epsilon$. The nature
of the leading order term depends on the value of the dimensionless surface
tension $\mathcal{S}$. When gravity is the only restoring force, we find

$\mathcal{C}_{r\pm}^{\rm{grav}}=\pm\sqrt{\frac{1-r}{1+r}\
\mathcal{G}\epsilon}\
+\frac{r\mathcal{U}^{\prime}(0^{+})-\mathcal{U}^{\prime}(0^{-})}{2(1+r)}\
\epsilon\
\pm\frac{[r\mathcal{U}^{\prime}(0^{+})-\mathcal{U}^{\prime}(0^{-})]^{2}}{8\sqrt{\mathcal{G}(1-r)}}\bigg{[}\frac{\epsilon}{1+r}\bigg{]}^{\frac{3}{2}}+\rm{h.o.t}.$
(29)

When capillary forces are included, we instead obtain

$\mathcal{C}_{r\pm}^{\rm{cap-
grav}}=\pm\sqrt{\frac{\mathcal{S}}{(1+r)\epsilon}}+\frac{r\mathcal{U}^{\prime}(0^{+})-\mathcal{U}^{\prime}(0^{-})}{2(1+r)}\
\epsilon\ \pm\frac{\mathcal{G}(1-r)}{2\sqrt{\mathcal{S}(1+r)}}\
\epsilon^{\frac{3}{2}}+\rm{h.o.t}.$ (30)

Hence, apart from the effect of shear, for short waves gravity acts as a high
order correction to the effect of surface tension. Note that the third term in
Eq. (30) can be derived by expanding the dispersion relation of interfacial
capillary-gravity waves,

$\mathcal{C}(\epsilon)=\pm\sqrt{\frac{\mathcal{S}}{(1+r)\epsilon}+\frac{1-r}{1+r}\
\mathcal{G}\epsilon}\ .$ (31)

Armed with $\mathcal{C}_{r\pm}$, we can use assumption (27), and Eqs. (23) and
(26), in Eq. (14). We collect terms of order $\epsilon\
e^{-2|\mathcal{z}_{\rm{c}\pm}|/\epsilon}$ to find the amplitudes $A_{\pm}$,
and infer the growth rates of the prograde (+) and retrograde (-) modes as

$\displaystyle\Im\\{\mathcal{k}\mathcal{C}_{+}\\}$
$\displaystyle=-\mathcal{C}_{r+}\ \frac{\pi}{2}\ \frac{r}{1+r}\
\frac{\mathcal{U}^{\prime\prime}(\mathcal{z}_{\rm{c}+})}{|\mathcal{U}^{\prime}(\mathcal{z}_{\rm{c}+})|}\
e^{-2\mathcal{k}\mathcal{z}_{\rm{c}+}},$ (32)
$\displaystyle\text{and}\qquad\Im\\{\mathcal{k}\mathcal{C}_{-}\\}$
$\displaystyle=-\mathcal{C}_{r-}\ \frac{\pi}{2}\ \frac{1}{1+r}\
\frac{\mathcal{U}^{\prime\prime}(\mathcal{z}_{\rm{c}-})}{|\mathcal{U}^{\prime}(\mathcal{z}_{\rm{c}-})|}\
e^{2\mathcal{k}\mathcal{z}_{\rm{c}-}},$ (33)

where $\mathcal{k}=1/\epsilon$. We emphasize two aspects of the present
results. Firstly, equation (32) for the prograde mode is a generalization of
the growth rate obtained by Miles in an Appendix of Morland and Saffman
(1992). That result was obtained from short wavelength asymptotic analysis of
the exact solution of the Rayleigh equation for an exponential wind profile.
Secondly, equation (33) for the retrograde mode is a generalization of the
short wavelength limit of the growth rate of the rippling instability found by
Young and Wolfe (2013) for an exponential wind-induced current.

(a) (b) (c) (d) (e) (f)

Figure 2: Comparison of the short wavelength asymptotic results (dashed lines)
for the prograde mode with the exact solution (solid lines) in the case of an
exponential wind profile (no current), $U(z)=U_{\infty}\
(1-e^{-z/d_{\rm{a}}})$, a density ratio $r=0.001$, and different values (grey
and black) of the dimensionless gravity,
$\mathcal{G}=gd_{\rm{a}}/U_{\infty}^{2}$, and surface tension,
$\mathcal{S}=\sigma/(\rho_{\rm{w}}U_{\infty}^{2}d_{\rm{a}})$. The phase speed
as a function of the wavenumber is shown for gravity (a), capillary (b), and
capillary-gravity waves (c); the exact and the asymptotic solutions are
indistinguishable. The growth rate as a function of the wavenumber for gravity
(d), capillary (e), and capillary-gravity waves (f).

(a) (b) (c) (d) (e) (f)

Figure 3: Same as Fig. 2 but with a density ratio of $r=0.9$. Whereas in Figs.
2(a)-(c) the asymptotic results (dashed lines) and the exact solution (solid
lines) are indistinguishable, such is not the case here, where we also see the
phase speed of free surface waves (dotted lines).

##  IV Interpretation

(a) (b) (c) (d) (e) (f)

Figure 4: Comparison of the short wavelength asymptotic results (dashed lines)
for the retrograde mode with the exact solution (solid lines) in the case of
an exponential current (no air), $U(z)=U_{\rm{S}}\ e^{z/d_{\rm{w}}}$, and
different values (grey and black) of the dimensionless gravity,
$\mathcal{G}=gd_{\rm{w}}/U_{\rm{S}}^{2}$, and surface tension,
$\mathcal{S}=\sigma/(\rho_{\rm{w}}U_{\rm{S}}^{2}d_{\rm{w}})$. The dotted lines
depict the corresponding phase speed of Doppler-shifted surface waves. The
growth rate as a function of the wavenumber for gravity (d), capillary (e),
and capillary-gravity waves (f).

The behavior of short waves in presence of a wind and a wind-induced current
depends on the profile of the latter solely through the derivatives at
$\mathcal{z}=0$ and at the critical levels,
$\mathcal{z}=\mathcal{z}_{\rm{c}\pm}$. In consequence, the results are
qualitatively the same for the two profiles introduced in §II.3. We note that,
as shown by Young and Wolfe (2013), the eigenvalue problem can be solved
exactly in terms of hypergeometric functions in the case of double exponential
profiles, whereas exact analytical solutions for double log profiles are
unknown.

In figures 2 and 3 we show the phase speed and the growth rate of the prograde
mode undergoing the Miles instability. We compare the results of our short
wavelength asymptotic analysis with the exact solution for an exponential wind
profile, in the case of gravity (a and d), capillary (b and e), and capillary-
gravity waves (c and f), for different values of the density ratio $r$ and the
control parameters $\mathcal{G}$ and $\mathcal{S}$. Figure 4 shows a similar
comparison for the retrograde mode undergoing the rippling instability in the
absence of air, that is $r=0$. The asymptotic results are in accord with the
exact results.

The phase speed of Doppler-shifted interfacial waves in presence of a
constant, zero shear, current is

$c(k)=U_{\rm{S}}\pm\sqrt{\frac{1-r}{1+r}\
\frac{g}{k}+\frac{\sigma}{\rho_{\rm{w}}(1+r)}\ k}\ ,$ (34)

which is the dimensional form of Eq. (31), where the plus (minus) sign
corresponds to the prograde (retrograde) mode. Equations (29) and (30) predict
that the non-uniformity of the current increases the phase speed by
$rU^{\prime}(0^{+})/2k$ and reduces it by $U^{\prime}(0^{-})/2k$ respectively.
Indeed, Figures 3(a)-(c) show that when $r$ is close to $1$, and $\mathcal{k}$
approaches unity, the phase speed of sheared waves is significantly larger
than the phase speed of free surface waves. Moreover, Figures 4(a)-(c) show
the predicted reduction of the retrograde mode in the absence of air when
$\mathcal{k}$ is of order unity.

Generally, for small values of $r$, the effect of the wind on dispersion is
insignificant, as shown by the phase speed of free surface in Figures
2(a)-(c). However, for large values of $r$, the wind has a major influence.
For instance, as shown in Figure 3(b), the wind is responsible for a minimum
phase speed of capillary waves. Moreover, as $r$ increases so does the growth
rate of the prograde mode. Thus, we interpret the density ratio as a wind-wave
coupling constant.

Finally, for both modes the growth rates increase when $\mathcal{G}$ and
$\mathcal{S}$ are small, which is the case for large velocity scales (c.f.,
Figs 2 and 4). Thus, consistent with physical intuition, a small perturbation
grows faster in the presence of strong winds and/or currents.

##  V Prograde instability due to critical layers in both air and water

Figure 5: Schematic of a wind blowing in the air and a current in the water.

We now explore the case of a wind blowing in the air when there is also a
current in the water, as shown in Figure 5. Thus we have
$\mathcal{U}\big{(}\mathcal{z}\lessgtr 0\big{)}>0$,
$\mathcal{U}^{\prime}(\mathcal{z}>0)>0$, and
$\mathcal{U}^{\prime}(\mathcal{z}<0)<0$, and again assume an infinite depth
and proceed as in §III. The key difference is that the prograde mode has a
critical layer in both the air and the water, and the retrograde mode has no
critical layer. Therefore, there are levels $\mathcal{z}_{\rm{c}\pm}$ such
that

$\mathcal{U}(\mathcal{z}_{\rm{c}+})=\mathcal{C}_{r+}\qquad\text{and}\qquad\mathcal{U}(\mathcal{z}_{\rm{c}-})=\mathcal{C}_{r+},$
(35)

and there is no value of $\mathcal{z}$ such that $\mathcal{C}_{r-}<0$ is equal
to $\mathcal{U}(\mathcal{z})$. This implies that

$\displaystyle f^{\prime}(0^{+})$ $\displaystyle=-i\pi\
\frac{\mathcal{U}^{\prime\prime}(\mathcal{z}_{\rm{c}+})}{|\mathcal{U}^{\prime}(\mathcal{z}_{\rm{c}+})|}\
e^{-2\mathcal{z}_{\rm{c}+}/\epsilon}+\epsilon\
\frac{\mathcal{U}^{\prime\prime}(0^{+})}{2\mathcal{C}_{r+}(\epsilon)}-\frac{\mathcal{U}^{\prime\prime}(\mathcal{z}_{\rm{c}+})}{\mathcal{U}^{\prime}(\mathcal{z}_{\rm{c}+})}\
\frac{\epsilon}{2\mathcal{z}_{\rm{c}+}}\ ,$ (36)
$\displaystyle\text{and}\qquad h^{\prime}(0^{-})$ $\displaystyle=i\pi\
\frac{\mathcal{U}^{\prime\prime}(\mathcal{z}_{\rm{c}-})}{|\mathcal{U}^{\prime}(\mathcal{z}_{\rm{c}-})|}\
e^{2\mathcal{z}_{\rm{c}-}/\epsilon}-\epsilon\
\frac{\mathcal{U}^{\prime\prime}(0^{-})}{2\mathcal{C}_{r+}(\epsilon)}+\frac{\mathcal{U}^{\prime\prime}(\mathcal{z}_{\rm{c}-})}{\mathcal{U}^{\prime}(\mathcal{z}_{\rm{c}-})}\
\frac{\epsilon}{2\mathcal{z}_{\rm{c}-}}\ ,$ (37)

for the prograde mode, while

$\displaystyle f^{\prime}(0^{+})$ $\displaystyle=\epsilon\
\frac{\mathcal{U}^{\prime\prime}(0^{+})}{2\mathcal{C}_{r-}(\epsilon)}\ ,$ (38)
$\displaystyle\text{and}\qquad h^{\prime}(0^{-})$ $\displaystyle=-\epsilon\
\frac{\mathcal{U}^{\prime\prime}(0^{-})}{2\mathcal{C}_{r-}(\epsilon)}\ ,$ (39)

for the retrograde mode. Hence, the retrograde mode is neutral but the
prograde mode can undergo both the Miles and rippling instabilities, and has a
growth rate of

$\Im\\{\mathcal{k}\mathcal{C}_{+}\\}=-\frac{\pi}{1+r}\
\frac{\mathcal{C}_{r+}}{2}\bigg{[}r\
\frac{\mathcal{U}^{\prime\prime}(\mathcal{z}_{\rm{c}+})}{|\mathcal{U}^{\prime}(\mathcal{z}_{\rm{c}+})|}\
e^{-2\mathcal{k}\mathcal{z}_{\rm{c}+}}+\frac{\mathcal{U}^{\prime\prime}(\mathcal{z}_{\rm{c}-})}{|\mathcal{U}^{\prime}(\mathcal{z}_{\rm{c}-})|}\
e^{2\mathcal{k}\mathcal{z}_{\rm{c}-}}\bigg{]}.$ (40)

The real parts $\mathcal{C}_{r\pm}$ are still given by equations (29) and
(30). Such an enhancement of the Miles instability by the rippling instability
may be difficult to observe in geophysical flows but could appear in other
flows.

##  VI Asymptotic solution of the Rayleigh equation for short waves

Figure 6: Schematic of the domain of analysis of the Rayleigh equation for
short waves; $\mathcal{k}\gg 1$. A boundary layer of thickness
$\epsilon=1/\mathcal{k}$ emerges from the singularity at the critical level,
$\mathcal{z}=\mathcal{z}_{\rm{c}}$.

Here we solve equation (11) for $\mathcal{C}=\mathcal{C}_{r+}$ and note that a
similar procedure is applicable to equation (12) when
$\mathcal{C}=\mathcal{C}_{r-}$. For simplicity, we rewrite equation (11)

$\epsilon f^{\prime\prime}(\mathcal{z})-2f^{\prime}(\mathcal{z})-\epsilon\
\frac{\mathcal{U}^{\prime\prime}(\mathcal{z})}{\mathcal{U}(\mathcal{z})-\mathcal{C}_{r}}\
f(\mathcal{z})=0,\qquad f(0)=1,$

and we drop the subscript $+$ for the rest of this section. There is a regular
singularity at $\mathcal{z}=\mathcal{z}_{\rm{c}}$ such that

$\mathcal{U}(\mathcal{z}_{\rm{c}})=\mathcal{C}_{r}.$

Because the small parameter $\epsilon$ multiplies the highest order derivative
in equation (11), we expect the solution to have a boundary layer somewhere in
the domain, but do not know its location a priori. However, guidance is
provided by the presence of the singularity. The Frobenius exponents are $0$
and $1$, so that the solution of equation (11) is finite at
$\mathcal{z}=\mathcal{z}_{\rm{c}}$ whereas its derivative has a logarithmic
divergence (Drazin and Reid, 1981). Therefore, we assume that an internal
boundary layer emerges from the singularity. We assume that the point
$\mathcal{z}=\mathcal{z}_{\rm{c}}(\epsilon)$ is well separated from the lower
boundary, $\mathcal{z}=0$, as $\epsilon$ goes to zero. We check this a
posteriori once the dependence of $\mathcal{C}_{r}$ on $\epsilon$ is known. In
consequence, $\mathcal{C}_{r}$ can be treated as a constant in the following
analysis.

The boundary layer is an inner region where the solution of equation (11)
changes rapidly. We define two outer regions, where the solution changes
slowly. One spans $\mathcal{z}=0$ to the inner region; the other spans the
inner region to the far field. We call them ‘lower outer region’ and ‘upper
outer region’, respectively (see Figure 6).

###  VI.1 Outer solutions

Following Miles (1962), we seek outer solutions as a power series in
$\epsilon$ as

$f_{\rm{out}}(\mathcal{z})=f_{0}(\mathcal{z})+\epsilon\
f_{1}(\mathcal{z})+O\big{(}\epsilon^{2}\big{)},\qquad\epsilon\to 0^{+},$ (41)

and find that

$f_{0}=\text{constant}\quad\text{and}\quad
f_{1}(\mathcal{z})=-\frac{f_{0}}{2}\int^{\mathcal{z}}_{\mathcal{z}_{1}}d\tilde{z}\
\frac{\mathcal{U}^{\prime\prime}(\tilde{z})}{\mathcal{U}(\tilde{z})-\mathcal{C}_{r}}\
.$ $None$

The constant $f_{0}$ and the lower limit of the integral $\mathcal{z}_{1}$ are
not a priori the same for the two outer solutions. In particular, they cannot
be determined by the boundary conditions at $z=0$ and infinity, requiring us
to find an inner solution.

###  VI.2 Inner solution

Within the boundary layer, we introduce the stretched coordinate
$Z\equiv(\mathcal{z}-\mathcal{z}_{\rm{c}})/\delta$, where $0<\delta\ll 1$, and
seek an inner solution, $F_{\rm{in}}(Z)$, to equation (11). We approximate the
coefficients by their Taylor series expansions about
$\mathcal{z}=\mathcal{z}_{\rm{c}}$, so that equation (11) becomes

$\frac{\epsilon}{\delta^{2}}\ F_{\rm{in}}^{\prime\prime}(Z)-\frac{2}{\delta}\
F_{\rm{in}}^{\prime}(Z)-\frac{\epsilon}{\delta}\frac{\mathcal{U}^{\prime\prime}_{\rm{c}}}{\mathcal{U}^{\prime}_{\rm{c}}Z}\
F_{\rm{in}}(Z)=0\ ;\qquad Z=O(1),$ (43)

where the subscript ‘c’ denotes evaluation at the critical level,
$\mathcal{z}=\mathcal{z}_{\rm{c}}$. By balancing the two first terms, we
obtain the distinguished limit $\delta=\epsilon$ (Bender and Orszag, 1999),
and hence must solve

$F_{\rm{in}}^{\prime\prime}(Z)-2F_{\rm{in}}^{\prime}(Z)=\epsilon\
\frac{\mathcal{U}^{\prime\prime}_{\rm{c}}}{\mathcal{U}^{\prime}_{\rm{c}}Z}\
F_{\rm{in}}(Z),\qquad\epsilon\to 0^{+}.$ (44)

Because the boundary layer has an $O(\epsilon)$ thickness, the appropriate
inner expansion is

$F_{\rm{in}}(Z)=F_{0}(Z)+\epsilon\
F_{1}(Z)+O\big{(}\epsilon^{2}\big{)},\qquad\epsilon\to 0^{+}.$ (45)

The general solutions are

$\displaystyle F_{0}(Z)$ $\displaystyle=A\ e^{2Z}+B,\qquad A,B\in\mathbb{C},$
(46) $\displaystyle\text{and}\qquad F_{1}(Z)$
$\displaystyle=\frac{\mathcal{U}^{\prime\prime}_{\rm{c}}}{\mathcal{U}^{\prime}_{\rm{c}}}\
\int^{Z}_{b}dx\ e^{2x}\int^{x}_{a}dt\ \frac{e^{-2t}}{t}\ F_{0}(t).$ (47)

Next we determine the integration constants, $A$ and $B$, and the bounds of
integration, $a$ and $b$, by asymptotic matching.

###  VI.3 Asymptotic matching and uniformly valid composite solutions

We must match the inner solution and the two outer solutions, as $\epsilon\to
0^{+}$. We use the superscripts ‘$\ell$’ and ‘u’ for lower and upper,
respectively.

The lower outer solution, $f_{\rm{out}}^{\ell}(\mathcal{z})$, is defined for
$0~{}\leq~{}\mathcal{z}~{}\ll~{}\mathcal{z}_{\rm{c}}$, and must satisfy the
lower boundary condition;

$f_{0}^{\ell}(0)=1,\qquad\text{and}\qquad f_{1}^{\ell}(0)=0.$ $None$

The upper outer solution, $f_{\rm{out}}^{\rm{u}}(\mathcal{z})$, is defined for
$\mathcal{z}\gg\mathcal{z}_{\rm{c}}$ and hence must satisfy condition (4). We
use the following matching conditions,

$\displaystyle\lim\limits_{\mathcal{z}\to\mathcal{z}_{\rm{c}}^{-}}f_{\rm{out}}^{\ell}(\mathcal{z})$
$\displaystyle=\lim\limits_{Z\to-\infty}F_{\rm{in}}(Z),$ (49)
$\displaystyle\text{and}\quad\lim\limits_{\mathcal{z}\to\mathcal{z}_{\rm{c}}^{+}}f_{\rm{out}}^{\rm{u}}(z)$
$\displaystyle=\lim\limits_{Z\to+\infty}F_{\rm{in}}(Z),$ (50)

and then apply the Van Dyke additive rule (Bender and Orszag, 1999) to
construct uniformly valid composite solutions;

$\text{uniform approx}=\text{inner}+\text{outer}-\text{common part}.$ (51)

We stress that this is possible because the lower and upper outer solutions
happen to have a common analytical expression.

####  VI.3.1 Leading order

To leading order, the outer solutions are constant (see Eq. $None$a), and
because of the boundary condition ($None$a) we have

$f_{0}^{\ell}=1.$ (52)

The leading order inner solution is given by equation (46). To preempt
divergence as $Z\to+\infty$, we impose $A=0$, and the matching condition (49)
implies that $B=1$. Therefore,

$F_{0}=1,$ (53)

which, upon imposition of the matching condition (50), yields

$f_{0}^{\rm{u}}=1.$ (54)

We combine the solutions (52), (53) and (54) using the additive rule (51), and
thus obtain a uniformly valid composite solution at leading order;

$f_{\rm{unif},0}=1+O(\epsilon).$ (55)

The effect of the boundary layer appears only at the next order.

####  VI.3.2 Order $\epsilon$

Following equation ($None$b), the lower and upper outer solutions at order
$\epsilon$ are

$\displaystyle f_{1}^{\ell}(\mathcal{z})$
$\displaystyle=-\frac{1}{2}\int^{\mathcal{z}}_{\mathcal{z}_{1}^{\ell}}d\tilde{z}\
\frac{\mathcal{U}^{\prime\prime}(\tilde{z})}{\mathcal{U}(\tilde{z})-\mathcal{C}_{r}}\
,$ (56) $\displaystyle\text{and}\quad f_{1}^{\rm{u}}(\mathcal{z})$
$\displaystyle=-\frac{1}{2}\int^{\mathcal{z}}_{\mathcal{z}_{1}^{\rm{u}}}d\tilde{z}\
\frac{\mathcal{U}^{\prime\prime}(\tilde{z})}{\mathcal{U}(\tilde{z})-\mathcal{C}_{r}}\
,$ (57)

respectively. From the Laurent series expansion

$\frac{\mathcal{U}^{\prime\prime}(\mathcal{z})}{\mathcal{U}(\mathcal{z})-\mathcal{C}_{r}}=\frac{\mathcal{U}^{\prime\prime}_{\rm{c}}}{\mathcal{U}^{\prime}_{\rm{c}}(\mathcal{z}-\mathcal{z}_{\rm{c}})}+\frac{1}{2}\bigg{[}\frac{\mathcal{U}^{\prime\prime}_{\rm{c}}}{\mathcal{U}^{\prime}_{\rm{c}}}\bigg{]}^{2}+O(\mathcal{z}-\mathcal{z}_{\rm{c}}),$
(58)

we deduce that

$\displaystyle f_{1}^{\ell}(\mathcal{z})$
$\displaystyle\sim-\frac{\mathcal{U}^{\prime\prime}_{\rm{c}}}{2\mathcal{U}^{\prime}_{\rm{c}}}\operatorname{Log}\bigg{(}\frac{\mathcal{z}-\mathcal{z}_{\rm{c}}}{\mathcal{z}_{1}^{\ell}-\mathcal{z}_{\rm{c}}}\bigg{)},\quad\mathcal{z}\to\mathcal{z}_{\rm{c}}^{-},$
(59) $\displaystyle\text{and}\quad f_{1}^{\rm{u}}(\mathcal{z})$
$\displaystyle\sim-\frac{\mathcal{U}^{\prime\prime}_{\rm{c}}}{2\mathcal{U}^{\prime}_{\rm{c}}}\operatorname{Log}\bigg{(}\frac{\mathcal{z}-\mathcal{z}_{\rm{c}}}{\mathcal{z}_{1}^{\rm{u}}-\mathcal{z}_{\rm{c}}}\bigg{)},\quad\mathcal{z}\to\mathcal{z}_{\rm{c}}^{+},$
(60)

where $\operatorname{Log}$ denotes a continuation of the natural logarithm to
the negative real numbers,

$\operatorname{Log}(\mathcal{z}-\mathcal{z}_{\rm{c}})=\begin{cases}\ln|\mathcal{z}-\mathcal{z}_{\rm{c}}|-i\pi\quad\text{if
}\mathcal{U}^{\prime}_{\rm{c}}>0,\\\
\ln|\mathcal{z}-\mathcal{z}_{\rm{c}}|+i\pi\quad\text{if
}\mathcal{U}^{\prime}_{\rm{c}}<0,\end{cases}\quad\text{for
}\mathcal{z}<\mathcal{z}_{\rm{c}}.$ (61)

The choice of the branch cut – just above the negative real axis if
$\mathcal{U}^{\prime}_{\rm{c}}>0$, just below otherwise – follows from Lin
(1955).

From equation (47), the inner solution at order $\epsilon$ is

$F_{1}(Z)=\frac{\mathcal{U}^{\prime\prime}_{\rm{c}}}{\mathcal{U}^{\prime}_{\rm{c}}}\
\int^{Z}_{b}dx\ e^{2x}\int^{x}_{a}dt\ \frac{e^{-2t}}{t}.$ (62)

We split the integral over $t$ as

$\int^{x}_{a}=\int^{x}_{+\infty}+\int_{a}^{+\infty}.$ (63)

Noting that the second integral is a constant, $C\in\mathbb{C}$, we obtain

$F_{1}(Z)=-\frac{\mathcal{U}^{\prime\prime}_{\rm{c}}}{\mathcal{U}^{\prime}_{\rm{c}}}\bigg{\\{}\int^{Z}_{b}dx\
e^{2x}E_{1}(2x)+\frac{C}{2}\Big{(}e^{2Z}-e^{2b}\Big{)}\bigg{\\}}.$ (64)

where

$E_{1}(x)\equiv\int_{x}^{+\infty}dt\ \frac{e^{-t}}{t},\qquad|\arg(x)|<\pi,$
(65)

is the exponential integral. For large values of $|x|$, the divergent series
(Bender and Orszag, 1999)

$E_{1}(x)=\frac{e^{-x}}{x}\ \sum_{n=0}^{N-1}\frac{n!}{(-x)^{n}}\ ,\quad
x\to\infty\ ;\quad|\arg(x)|<\frac{3\pi}{2},$ (66)

yields

$e^{2x}E_{1}(2x)\sim\frac{1}{2x}\ ,\qquad x\to\pm\infty.$ (67)

After integration, we find

$F_{1}(Z)\sim-\frac{\mathcal{U}^{\prime\prime}_{\rm{c}}}{2\mathcal{U}^{\prime}_{\rm{c}}}\bigg{\\{}\operatorname{Log}\bigg{[}\frac{Z}{b}\bigg{]}+C\Big{(}e^{2Z}-e^{2b}\Big{)}\bigg{\\}},\quad
Z\to\pm\infty.$ (68)

Recalling that $Z=(\mathcal{z}-\mathcal{z}_{\rm{c}})/\epsilon$, the matching
of the inner limits, given by equations (59) and (60), and the outer limits,
given by equation (68), yields

$C=0\qquad\text{and}\qquad\mathcal{z}_{1}^{\ell}-\mathcal{z}_{\rm{c}}=\epsilon
b=\mathcal{z}_{1}^{\rm{u}}-\mathcal{z}_{\rm{c}}.$ $None$

Because the boundary condition ($None$b) requires that

$\mathcal{z}_{1}^{\ell}=0,$ (70)

the outer solution at order $\epsilon$ has the same expression in the lower
and upper outer regions;

$f_{1}(\mathcal{z})=-\frac{1}{2}\int_{0}^{\mathcal{z}}d\tilde{z}\
\frac{\mathcal{U}^{\prime\prime}(\tilde{z})}{\mathcal{U}(\tilde{z})-\mathcal{C}_{r}}\
,\quad\mathcal{z}\ll\mathcal{z}_{\rm{c}}\ \text{ and }\
\mathcal{z}\gg\mathcal{z}_{\rm{c}}.$ (71)

The lower outer solution is real, because the path of integration along the
real axis does not reach the branch point $\mathcal{z}=\mathcal{z}_{\rm{c}}$.
However, for the upper outer solution we have to make a detour in the complex
plane. We deform the path of integration in the upper part if
$\mathcal{U}^{\prime}_{\rm{c}}>0$, and in the lower part otherwise. In
Appendix B, we show that

$\Im\\{f_{1}^{\rm{u}}\\}=-\frac{\pi}{2}\
\frac{\mathcal{U}^{\prime\prime}_{\rm{c}}}{|\mathcal{U}^{\prime}_{\rm{c}}|}\
.$ (72)

Using the matching conditions ($None$) together with equation (70), the inner
solution (64) becomes

$F_{1}(Z)=-\frac{\mathcal{U}^{\prime\prime}_{\rm{c}}}{\mathcal{U}^{\prime}_{\rm{c}}}\int^{Z}_{-\mathcal{z}_{\rm{c}}/\epsilon}dx\
e^{2x}E_{1}(2x).$ (73)

Finally, we use the Van Dyke rule (51) to construct a uniformly valid
composite solution at order $\epsilon$;

$f_{\rm{unif},1}(\mathcal{z})=1+\epsilon\Bigg{\\{}-\frac{\mathcal{U}^{\prime\prime}_{\rm{c}}}{\mathcal{U}^{\prime}_{\rm{c}}}\
\int^{Z}_{-\mathcal{z}_{\rm{c}}/\epsilon}dx\
e^{2x}E_{1}(2x)-\frac{1}{2}\int_{0}^{\mathcal{z}}d\tilde{z}\
\frac{\mathcal{U}^{\prime\prime}(\tilde{z})}{\mathcal{U}(\tilde{z})-\mathcal{C}_{r}}+\frac{\mathcal{U}^{\prime\prime}_{\rm{c}}}{2\mathcal{U}^{\prime}_{\rm{c}}}\
\operatorname{Log}\bigg{(}1-\frac{\mathcal{z}}{\mathcal{z}_{\rm{c}}}\bigg{)}\Bigg{\\}}+O\big{(}\epsilon^{2}\big{)}.$
(74)

In Appendix B, we verify that
$\chi(\mathcal{z})=e^{-\mathcal{k}\mathcal{z}}f_{\rm{unif},1}(\mathcal{z})$
satisfies the global property;

$\Im\big{\\{}\chi^{\prime}(0^{+})\big{\\}}=-\pi\
\frac{\mathcal{U}^{\prime\prime}_{\rm{c}}}{|\mathcal{U}^{\prime}_{\rm{c}}|}\
|\chi_{\rm{c}}|^{2},$ (75)

at leading order. In Appendix A we show that

$f(\mathcal{z}_{\rm{c}})=1+\frac{\epsilon}{2}\
\frac{\mathcal{U}^{\prime\prime}_{\rm{c}}}{\mathcal{U}^{\prime}_{\rm{c}}}\Bigg{(}\gamma_{\rm{E}}+\ln\bigg{[}\frac{2\mathcal{z}_{\rm{c}}}{\epsilon}\bigg{]}-i\pi\Bigg{)}+O\big{(}\epsilon^{2}\big{)}.$
(76)

Equation (76) generalizes a result that Miles derived in an Appendix to
Morland and Saffman (1992) using an exact solution of the Rayleigh equation
for an exponential wind profile.

###  VI.4 Comparison with the numerical solution

(a) (b) (c) (d)

Figure 7: Comparison of the uniformly valid composite solution (74) with the
numerical solution of the Rayleigh equation for the exponential (a,b) and
logarithmic (c,d) wind profiles, with $\mathcal{C}_{r}=0.1$ and different
values of $\mathcal{k}$. The dots and the stars denote the real and imaginary
parts of the numerical solution, respectively. The solid line shows the real
part of (74) and the dashed line the imaginary part.

We numerically solve equation (11) for the two standard wind profiles

$\mathcal{U}(\mathcal{z})=1-e^{-\mathcal{z}}\qquad\text{and}\qquad\mathcal{U}(\mathcal{z})=\ln(1+\mathcal{z})/\kappa,$
$None$

where, as before, $\kappa=0.4$ is the von Kármán constant. We compare the
numerical solution with our uniformly valid composite solution (74) for fixed
values of $\mathcal{k}$ and $\mathcal{C}_{r}$ in Figure 7. (Note that any
dispersion relation can be retrieved with a proper choice of the control
parameters.) For both profiles the composite solution and the numerical
solution agree very well. We distinguish the lower (upper) outer solutions
with their imaginary part being equal to zero (a positive constant).
Consistent with the Frobenius solution (Drazin and Reid, 1981), the solution
within the inner layer depends on the wind profile and the dispersion relation
solely through the scale factor
$\mathcal{U}^{\prime\prime}_{\rm{c}}/\mathcal{U}^{\prime}_{\rm{c}}$ and the
bound of integration $\mathcal{z}_{\rm{c}}/\epsilon$. Since the phase of the
solution of the Rayleigh equation changes only within the boundary layer, we
conclude that the interaction of short waves with the wind principally occurs
therein. In contrast, we showed that for $\mathcal{k}\ll 1$ the phase varies
from $\mathcal{z}=0$ to $\mathcal{z}=\mathcal{z}_{\rm{c}}$, so that long waves
interact with the wind all the way from the mean water surface to the critical
level (Bonfils  _et al._ , 2022).

###  VI.5 Similarity solution

We have non-dimensionalized the variables using external parameters
characterizing the shear in the air. However, the general short wave solution
(74) can also be written in terms of $\xi\equiv kz$ and $\mathfrak{U}\equiv
U/\tilde{c}_{r}$, where $\tilde{c}_{r}$ is the dimensional version of
$\mathcal{C}_{r}$. The wind-generated ripples have a continuous wavenumber
spectrum so that $k$, and subsequently $\tilde{c}_{r}$, are variables rather
than parameters. In that sense, $\xi$ is a similarity variable introduced by
Miles (1957), and the Rayleigh equation has a self-similar solution,
$\phi=\phi(\xi)$, which for short waves is

$\phi(\xi)=e^{-\xi}\Bigg{\\{}1-\frac{\mathfrak{U}^{\prime\prime}(\xi_{\rm{c}})}{\mathfrak{U}^{\prime}(\xi_{\rm{c}})}\
\int^{\xi-\xi_{\rm{c}}}_{-\xi_{\rm{c}}}dx\
e^{2x}E_{1}(2x)-\frac{1}{2}\int_{0}^{\xi}d\tilde{z}\
\frac{\mathfrak{U}^{\prime\prime}(\tilde{z})}{\mathfrak{U}(\tilde{z})-1}+\frac{\mathfrak{U}^{\prime\prime}(\xi_{\rm{c}})}{2\mathfrak{U}^{\prime}(\xi_{\rm{c}})}\
\operatorname{Log}\bigg{(}1-\frac{\xi}{\xi_{\rm{c}}}\bigg{)}\Bigg{\\}}.$ (78)

Figure 8: Comparison of the growth rate of gravity waves obtained by Alexakis
_et al._ (2004b) with the predictions of our short wavelength asymptotic
treatment (cf. Eq. 32) and the growth rate calculated from the exact
(hypergeometric) solution for a wind profile $U(z)=U_{\infty}\
(1-e^{-z/d_{\rm{a}}})$ and a density ratio $r=0.5$.

##  VII The WKB method in the shortwavelength limit

Alexakis  _et al._ (2004b) proposed a solution of the Rayleigh equation for
short waves using the WKB method. However, because their solution does not
satisfy the global property (75), their growth rate has an extra factor
$\pi/\mathcal{k}$; see their equation (A31). In Figure 8, we show the growth
rate of gravity waves for an exponential wind profile and compare the results
of Alexakis  _et al._ (2004b) with our short wave asymptotic solution and the
exact solution. Clearly the WKB method deviates from the others for all
$kd_{a}\gtrsim 25$, and in what follows we explain the origin of the
deviation.

Setting $\epsilon=1/\mathcal{k}$, we write the Rayleigh equation in the form

$\epsilon^{2}\chi^{\prime\prime}(\mathcal{z})=Q(\mathcal{z},\epsilon)\chi(\mathcal{z}),\quad\text{with}\quad
Q(\mathcal{z},\epsilon)\equiv 1+\epsilon^{2}\
\frac{\mathcal{U}^{\prime\prime}(\mathcal{z})}{\mathcal{U}(\mathcal{z})-\mathcal{C}_{r}}\
.$ (79)

A WKB expansion is (Bender and Orszag, 1999)

$\chi(\mathcal{z})\sim\exp\bigg{[}\frac{1}{\epsilon}\
\sum_{n=0}^{+\infty}\epsilon^{n}S_{n}(\mathcal{z})\bigg{]},\qquad\epsilon\to
0.$ (80)

In the standard case, $Q(\mathcal{z},\epsilon)$ does not depend on the small
parameter $\epsilon$. Then, we can take only two terms in the series (80):
this, physical optics approximation, yields

$\chi(\mathcal{z})\sim
C_{1}\big{[}Q(\mathcal{z})\big{]}^{-\frac{1}{4}}\exp\bigg{[}\frac{1}{\epsilon}\int_{a}^{\mathcal{z}}d\tilde{z}\sqrt{Q(\tilde{z})}\bigg{]}+C_{2}\big{[}Q(\mathcal{z})\big{]}^{-\frac{1}{4}}\exp\bigg{[}-\frac{1}{\epsilon}\int_{a}^{\mathcal{z}}d\tilde{z}\sqrt{Q(\tilde{z})}\bigg{]}$
(81)

for $Q(\mathcal{z})>0$, and

$\chi(\mathcal{z})\sim
C_{3}\big{[}-Q(\mathcal{z})\big{]}^{-\frac{1}{4}}\exp\bigg{[}\frac{i}{\epsilon}\int_{b}^{\mathcal{z}}d\tilde{z}\sqrt{-Q(\tilde{z})}\bigg{]}+C_{4}\big{[}-Q(\mathcal{z})\big{]}^{-\frac{1}{4}}\exp\bigg{[}-\frac{i}{\epsilon}\int_{b}^{\mathcal{z}}d\tilde{z}\sqrt{-Q(\tilde{z})}\bigg{]},$
(82)

for $Q(\mathcal{z})<0$, where the bounds of integration, $a$ and $b$, are
arbitrary. For $\mathcal{U}^{\prime\prime}(\mathcal{z})<0$ and
$Q(\mathcal{z},\epsilon)$ introduced in equation (79), Alexakis  _et al._
(2004b) used the physical optics approximation on three intervals, defined in
Figure 9, and obtained three solutions of the form of (81) or (82), in which
they neglected the fact that $\epsilon\ll 1$ within $Q(\mathcal{z},\epsilon)$.
Then, they had to match those solutions in order to determine the integration
constants $C_{1}$, $C_{2}$, etc. The common matching procedure takes place at
the simple turning point $\mathcal{z}_{\star}$ (e.g. Bender and Orszag, 1999).
However, they needed inner solutions near the critical level,
$\mathcal{z}_{\rm{c}}$, and it is at this juncture that the appeal of the WKB
method is understood. Indeed, they found that the solution of the Rayleigh
equation near the singularity can be represented in terms of Bessel functions
for $\mathcal{z}>\mathcal{z}_{\rm{c}}$, and in terms of modified Bessel
functions for $\mathcal{z}<\mathcal{z}_{\rm{c}}$. Moreover, the outer limits
of those inner solutions formally match the inner limits of the physical
optics approximations. Unfortunately, however, the distance between the
critical point $\mathcal{z}_{\rm{c}}$ and the turning point
$\mathcal{z}_{\star}$ actually shrinks as $\epsilon$ tends to zero. For
instance, in the case of the exponential profile
$\mathcal{U}(\mathcal{z})=1-e^{-\mathcal{z}}$, we have

$\mathcal{z}_{\star}-\mathcal{z}_{\rm{c}}=\ln(1+\epsilon).$ (83)

Figure 9: The three intervals in the WKB approach of Alexakis  _et al._
(2004b) where $\mathcal{z}_{\star}$ is a simple turning point.

Hence, the interval where $Q<0$ cannot be taken as an outer region. In fact,
as seen in Figure 7, the numerical solution of the Rayleigh equation does not
exhibit oscillatory behavior.

In the vicinity of the critical level, $\mathcal{z}_{\rm{c}}$, we introduce
the inner variable

$\mathzapf{Z}=-\frac{\mathcal{U}^{\prime\prime}_{\rm{c}}}{\mathcal{U}^{\prime}_{\rm{c}}}\
(\mathcal{z}-\mathcal{z}_{\rm{c}}),$ (84)

in terms of which the inner solutions are

$\displaystyle\chi_{\text{in+}}(\mathzapf{Z}>0)$
$\displaystyle=\sqrt{\mathzapf{Z}}\Big{[}A\
J_{1}\big{(}2\sqrt{\mathzapf{Z}}\big{)}+B\
Y_{1}\big{(}2\sqrt{\mathzapf{Z}}\big{)}\Big{]},$ (85)
$\displaystyle\chi_{\text{in-}}(\mathzapf{Z}<0)$
$\displaystyle=\sqrt{-\mathzapf{Z}}\Big{[}C\
I_{1}\big{(}2\sqrt{-\mathzapf{Z}}\big{)}+D\
K_{1}\big{(}2\sqrt{-\mathzapf{Z}}\big{)}\Big{]},$ (86)

where $A$, $B$, $C$, and $D$ are complex constants, and $J_{1}$ and $Y_{1}$ (
$I_{1}$ and $K_{1}$) are the Bessel functions (modified Bessel functions) of
order $1$. For the solution to be continuous and its derivative to have the
correct behavior at $\mathcal{z}=\mathcal{z}_{\rm{c}}$, we must establish some
relationships between these constants. Although this was not done explicitly
by Alexakis  _et al._ (2004b), we could use their final matching formulae to
determine such relations to find them to be incompatible with the following;

$I_{1}(\mathcal{z})=-i\ J_{1}(i\mathcal{z})\ \text{ and }\
K_{1}(\mathcal{z})=-\frac{\pi}{2}\Big{[}J_{1}(-i\mathcal{z})-i\
Y_{1}(-i\mathcal{z})\Big{]},$ (87)

for $-\pi<\arg(\mathcal{z})\leq\pi/2$. We find that

$C=-A+i\ B\qquad\text{and}\qquad D=-\frac{2}{\pi}\ B.$ (88)

To check the consistency of these results, we use equation (88) for the
numerical integration of the Rayleigh equation and retrieve the growth rates
calculated by Beji and Nadaoka (2004).

##  VIII Conclusion

We have studied the effect of a shear flow on interfacial capillary-gravity
waves when their wavelength is much smaller than the characteristic length
scale of the flows in the fluids bounding the interface. Using the
dimensionless inverse wavenumber as a small parameter,
$\epsilon=1/\mathcal{k}$, we asymptotically solved the eigenvalue problem for
the stability of an arbitrary parallel flow $\boldsymbol{U}=U(z)\
\boldsymbol{\hat{x}}$ through a two-fluid interface with a density ratio that
is not necessarily small. We constructed uniformly valid composite solutions
of the governing Rayleigh equation, where the real part of the eigenvalue is a
power series in $\epsilon$. We showed that including the effect of surface
tension changes the nature of the leading order solution, and that
exponentially small terms must be considered in order to have a non-zero
imaginary part of the solution. From a physical view-point, there is a
prograde mode whose phase speed is greater than the speed of the water
surface, $U_{\rm{S}}$, and a retrograde mode whose phase speed is smaller than
$U_{\rm{S}}$. When the velocity of the shear flow equals the phase speed of
one these modes, there is a critical layer where the flow transfers energy to
the wave. If the critical layer is in the air (water), this is called the
Miles (rippling) instability. The only case considered in the literature thus
far is when the prograde mode undergoes the Miles instability and the
retrograde mode undergoes the rippling instability. In § V we studied the
situation where the prograde mode can undergo both instabilities and the
retrograde mode is neutral. In the short wave limit, we found that (i) the
effect of the shear on the phase speed of the two modes depends only on the
derivatives of $U$ at $z=0^{\pm}$; and (ii) the Miles and rippling
instabilities have a growth rate of the same form. Indeed, the interaction
between the shear flow and the waves is mostly reduced to a narrow region
around the critical level, an internal boundary layer of thickness $\epsilon$
where the solution of the Rayleigh equation has a self-similar structure.
Heuristically speaking, the waves are barely influenced by the flow outside of
this region. Nonetheless, we showed that there are signifiant effects on
dispersion and growth rates when the characteristic velocity of the shear flow
is large and the density ratio is close to $1$. Finally, we showed how the WKB
approach of solving the Rayleigh equation breaks down.

## Acknowledgements

We thank Cristobal Arratia for suggesting the possibility of a prograde mode
with two critical layers, developed in §V. A. F. Bonfils thanks Neil Balmforth
and Robert Rosner for discussions of the WKB approach and the 2022 Geophysical
Fluid Dynamics Summer Study Program at the Woods Hole Oceanographic
Institution, which is supported by the National Science Foundation and the
Office of Naval Research under No. OCE-1332750. All authors acknowledge
Swedish Research Council grant no. 638-2013-9243. Nordita is supported in part
by NordForsk.

##  Appendix A Solution of equation (11) at the singular point
$\mathcal{z}=\mathcal{z}_{\rm{c}}$

The regular singular point $\mathcal{z}=\mathcal{z}_{\rm{c}}$ is located in
the inner region where the solution at order $\epsilon$ is given by equation
(73), which we reproduce here;

$F_{1}(Z)=-\frac{\mathcal{U}^{\prime\prime}_{\rm{c}}}{\mathcal{U}^{\prime}_{\rm{c}}}\int^{Z}_{-\mathcal{z}_{\rm{c}}/\epsilon}dx\
e^{2x}E_{1}(2x),\ \text{ where }\
Z=\frac{\mathcal{z}-\mathcal{z}_{\rm{c}}}{\epsilon}$ (89)

is the inner variable.

###  A.1 Imaginary part

From the series representation of the exponential integral

$E_{1}(x)=-\gamma_{\rm{E}}-\operatorname{Log}(x)-\sum_{n=1}^{+\infty}\frac{(-x)^{n}}{n\
n!},\quad|\arg(x)|<\pi,$ (90)

where $\gamma_{\rm{E}}=0.577$ is the Euler constant, we see that $x=0$ is a
logarithmic branch point for the integral in equation (89). The only
contribution to the imaginary part of $F_{1}$ arises from the integration of
$\operatorname{Log}(x)$ along the branch cut. Namely, if $\Theta$ is the
Heaviside step function then

$\Im\bigg{\\{}\int^{Z}_{-\mathcal{z}_{\rm{c}}/\epsilon}dx\
e^{2x}E_{1}(2x)\bigg{\\}}=\pm\Theta(-Z)\
\pi\int^{Z}_{-\mathcal{z}_{\rm{c}}/\epsilon}dx\ e^{2x}\pm\Theta(Z)\
\pi\int^{0}_{-\mathcal{z}_{\rm{c}}/\epsilon}dx\ e^{2x},$ (91)

where we select the plus sign when $\mathcal{U}^{\prime}_{\rm{c}}>0$, and the
branch cut is just above the negative real axis, and the minus sign otherwise.
From equation (91), we find

$\Im\big{\\{}F_{1}(Z)\big{\\}}=\begin{cases}-\frac{\pi}{2}\
\frac{\mathcal{U}^{\prime\prime}_{\rm{c}}}{|\mathcal{U}^{\prime}_{\rm{c}}|}\big{(}e^{2Z}-e^{-2\mathcal{z}_{\rm{c}}/\epsilon}\big{)}\quad\text{if
}Z<0,\\\ \\\ -\frac{\pi}{2}\
\frac{\mathcal{U}^{\prime\prime}_{\rm{c}}}{|\mathcal{U}^{\prime}_{\rm{c}}|}\big{(}1-e^{-2\mathcal{z}_{\rm{c}}/\epsilon}\big{)}\quad\text{if
}Z\geq 0.\end{cases}$ (92)

Then, upon taking the limit $Z\to 0$ we obtain

$\Im\big{\\{}f(\mathcal{z}_{\rm{c}})\big{\\}}=-\epsilon\ \frac{\pi}{2}\
\frac{\mathcal{U}^{\prime\prime}_{\rm{c}}}{|\mathcal{U}^{\prime}_{\rm{c}}|}\Big{(}1-e^{-2\mathcal{z}_{\rm{c}}/\epsilon}\Big{)}+O\big{(}\epsilon^{2}\big{)}.$
(93)

It is interesting to take the outer limits $Z\to\pm\infty$ of the expression
(92). In particular, if the matching was done properly, these should be equal
to the inner limits of the imaginary parts of the lower and upper outer
solutions, viz.,

$\displaystyle\lim\limits_{\mathcal{z}\to\mathcal{z}_{\rm{c}}^{-}}\Im\big{\\{}f_{1}^{\ell}(\mathcal{z})\big{\\}}$
$\displaystyle=\frac{\pi}{2}\
\frac{\mathcal{U}^{\prime\prime}_{\rm{c}}}{|\mathcal{U}^{\prime}_{\rm{c}}|}\
e^{-2\mathcal{z}_{\rm{c}}/\epsilon},$ (94)
$\displaystyle\text{and}\quad\lim\limits_{\mathcal{z}\to\mathcal{z}_{\rm{c}}^{+}}\Im\big{\\{}f_{1}^{\rm{u}}(\mathcal{z})\big{\\}}$
$\displaystyle=-\frac{\pi}{2}\
\frac{\mathcal{U}^{\prime\prime}_{\rm{c}}}{|\mathcal{U}^{\prime}_{\rm{c}}|}\Big{(}1-e^{-2\mathcal{z}_{\rm{c}}/\epsilon}\Big{)}.$
(95)

In each case, we shall check that

$\lim\limits_{\epsilon\to
0^{+}}\frac{\mathcal{z}_{\rm{c}}(\epsilon)}{\epsilon}=+\infty.,$ (96)

which is the condition of separation of the boundary layer of thickness
$\epsilon$ from the lower boundary. Thus, we consistently retrieve a zero
imaginary part in the lower outer region. We recall that the upper outer
solution is given by (Eq. 71)

$f_{1}^{\rm{u}}(\mathcal{z})=-\frac{1}{2}\int_{0}^{\mathcal{z}}d\tilde{z}\
\frac{\mathcal{U}^{\prime\prime}(\tilde{z})}{\mathcal{U}(\tilde{z})-\mathcal{C}_{r}}\
,\qquad\mathcal{z}\gg\mathcal{z}_{\rm{c}}.$ (97)

According to the Sokhotski-Plemelj theorem,

$\lim\limits_{\mathcal{C}_{i}\to
0^{+}}-\frac{1}{2}\int_{0}^{\mathcal{z}}d\tilde{z}\
\frac{\mathcal{U}^{\prime\prime}(\tilde{z})}{\mathcal{U}(\tilde{z})-(\mathcal{C}_{r}+i\mathcal{C}_{i})}=-\frac{1}{2}\
\mathzapf{P}\int_{0}^{\mathcal{z}}d\tilde{z}\
\frac{\mathcal{U}^{\prime\prime}(\tilde{z})}{\mathcal{U}(\tilde{z})-\mathcal{C}_{r}}-i\pi\
\frac{\mathcal{U}^{\prime\prime}_{\rm{c}}}{2|\mathcal{U}^{\prime}_{\rm{c}}|}\
,\qquad\forall\mathcal{z}>\mathcal{z}_{\rm{c}},$ (98)

where $\mathzapf{P}$ denotes the Cauchy principal value. Hence, the imaginary
part in the upper outer region is constant, equal to the upper outer limit of
$\Im\\{F_{1}\\}$ up to exponentially small terms.

###  A.2 Real part

Figure 10: The ‘Remainder’ $R(\zeta)$ given by equation (109).

At the leading order, we have

$f(\mathcal{z}_{\rm{c}})=1+O(\epsilon).$ (99)

Here we calculate the next order contribution to the real part. For
convenience, we define

$J(Z)\equiv\Re\bigg{\\{}\int^{Z}_{-\mathcal{z}_{\rm{c}}/\epsilon}dx\
e^{2x}E_{1}(2x)\bigg{\\}},$ (100)

and we express the exponential integral using the Ramanujan series as

$E_{1}(x)=-\gamma_{\rm{E}}-\operatorname{Log}(x)+e^{-\frac{x}{2}}\sum_{n=1}^{+\infty}\frac{s_{n}}{n!}\bigg{[}\frac{x}{2}\bigg{]}^{n},\quad
s_{n}=\sum_{j=0}^{\left\lfloor\frac{n-1}{2}\right\rfloor}\frac{1}{j+\frac{1}{2}}\
,$ (101)

where $\left\lfloor x\right\rfloor$ denotes the floor function. This series
converges faster than (90). We perform a direct integration as

$J(Z)=-\int^{Z}_{-\mathcal{z}_{\rm{c}}/\epsilon}dx\
e^{2x}\big{(}\gamma_{\rm{E}}+\ln(2)+\ln|x|\big{)}+\sum_{n=1}^{+\infty}\frac{s_{n}}{n!}\int^{Z}_{-\mathcal{z}_{\rm{c}}/\epsilon}dx\
e^{x}x^{n},$ (102)

and after some algebra we obtain

$J(Z)=\frac{1}{2}\Bigg{\\{}\Big{[}\gamma_{\rm{E}}+\ln(2)\Big{]}\Big{[}e^{-2\mathcal{z}_{\rm{c}}/\epsilon}-e^{2Z}\Big{]}+e^{-2\mathcal{z}_{\rm{c}}/\epsilon}\ln\bigg{[}\frac{\mathcal{z}_{\rm{c}}}{\epsilon}\bigg{]}-e^{2Z}\ln|Z|+E_{1}\bigg{[}\frac{2\mathcal{z}_{\rm{c}}}{\epsilon}\bigg{]}-E_{1}\big{(}2|Z|\big{)}\Bigg{\\}}+S(Z),$
(103)

with

$S(Z)\equiv\sum_{n=1}^{+\infty}(-1)^{n}\frac{\Gamma(n+1,-Z)-\Gamma(n+1,\mathcal{z}_{\rm{c}}/\epsilon)}{n!}\
s_{n},$ (104)

in which

$\Gamma(a,x)\equiv\int_{x}^{+\infty}dt\ t^{a-1}e^{-t},$ (105)

is the upper incomplete gamma function. Note that

$\Gamma(n+1,0)=n!\quad\text{and}\quad s_{2p-1}=s_{2p},\quad\forall
p\in\mathbb{N}^{\star},$ $None$

and that

$\lim\limits_{Z\to
0}E_{1}(2\mathcal{z}_{\rm{c}}/\epsilon)-E_{1}\big{(}2|Z|\big{)}-e^{2Z}\ln|Z|=-\ln\bigg{[}\frac{\mathcal{z}_{\rm{c}}}{\epsilon}\bigg{]}+e^{-\mathcal{z}_{\rm{c}}/\epsilon}\sum_{n=1}^{+\infty}\bigg{[}\frac{\mathcal{z}_{\rm{c}}}{\epsilon}\bigg{]}^{n}\frac{s_{n}}{n!}\
.$ (107)

Therefore,

$\lim\limits_{Z\to
0}J(Z)=-\frac{1}{2}\Bigg{(}\gamma_{\rm{E}}+\ln\bigg{[}\frac{2\mathcal{z}_{\rm{c}}}{\epsilon}\bigg{]}\Bigg{)}\Big{(}1-e^{-2\mathcal{z}_{\rm{c}}/\epsilon}\Big{)}+R\bigg{[}\frac{\mathcal{z}_{\rm{c}}}{\epsilon}\bigg{]}.$
(108)

We define the ‘remainder’ as

$R(\zeta)\equiv\frac{e^{-\zeta}}{2}\sum_{n=1}^{+\infty}\frac{\zeta^{n}}{n!}\
s_{n}-\sum_{n=1}^{+\infty}(-1)^{n}\frac{\Gamma(n+1,\zeta)}{n!}\ s_{n},$ (109)

which is a negligible quantity for $\zeta\gg 1$, as shown in Figure 10.
Discarding the exponentially small terms, we arrive at

$\Re\big{\\{}f(\mathcal{z}_{\rm{c}})\big{\\}}=1+\frac{\epsilon}{2}\
\frac{\mathcal{U}^{\prime\prime}_{\rm{c}}}{\mathcal{U}^{\prime}_{\rm{c}}}\Bigg{(}\gamma_{\rm{E}}+\ln\bigg{[}\frac{2\mathcal{z}_{\rm{c}}}{\epsilon}\bigg{]}\Bigg{)}+O\big{(}\epsilon^{2}\big{)}.$
(110)

##  Appendix B Global property of the solution of the Rayleigh equation

The dimensionless Rayleigh equation is

$\chi^{\prime\prime}(\mathcal{z})+\bigg{[}\mathcal{k}^{2}+\frac{\mathcal{U}^{\prime\prime}(\mathcal{z})}{\mathcal{U}(\mathcal{z})-\mathcal{C}_{r}}\bigg{]}\chi(\mathcal{z})=0,$
(111)

where $\mathcal{C}_{r}=\mathcal{C}_{r}(\mathcal{k})$ is a known function; the
boundary conditions are

$\chi(0)=1\quad\text{and}\quad\lim\limits_{\mathcal{z}\to+\infty}\chi(\mathcal{z})=0.$
(112)

Assuming that there is a unique level $\mathcal{z}_{\rm{c}}>0$ such that
$\mathcal{U}(\mathcal{z}_{\rm{c}})=\mathcal{C}_{r}$, Miles (1957) showed that

$\Im\big{\\{}\chi^{\prime}(0^{+})\big{\\}}=-\pi\
\frac{\mathcal{U}^{\prime\prime}_{\rm{c}}}{|\mathcal{U}^{\prime}_{\rm{c}}|}\
|\chi_{\rm{c}}|^{2},$ (113)

where the subscript ‘c’ denotes evaluation at
$\mathcal{z}=\mathcal{z}_{\rm{c}}$. The derivation of the global property
(113) is analogous to the canonical derivation of Rayleigh’s inflection point
theorem as follows. Multiply the Rayleigh equation (111) by the complex
conjugate of $\chi(\mathcal{z})$, integrate by parts, use the boundary
conditions (112) to evaluate the integrated term, and take the imaginary part;
the result follows from the Sokhotski-Plemelj theorem.

Here, we show that our asymptotic solution for short waves satisfies equation
(113). For $\mathcal{z}>0$, we let

$\chi(\mathcal{z})=e^{-\mathcal{z}/\epsilon}\
f(\mathcal{z}),\qquad\text{with}\qquad\epsilon=1/\mathcal{k}\ .$ (114)

We take the derivative of the composite solution,
$f_{\rm{unif},1}(\mathcal{z})$, obtained in equation (74) and obtain

$f_{\rm{unif},1}^{\prime}(0^{+})=-\frac{\mathcal{U}^{\prime\prime}_{\rm{c}}}{\mathcal{U}^{\prime}_{\rm{c}}}\bigg{[}E_{1}\bigg{(}\frac{2\mathcal{z}_{\rm{c}}}{\epsilon}\bigg{)}\pm
i\pi\bigg{]}e^{-2\mathcal{z}_{\rm{c}}/\epsilon}+\epsilon\
\frac{\mathcal{U}^{\prime\prime}(0^{+})}{2\mathcal{C}_{r}(\epsilon)}-\frac{\mathcal{U}^{\prime\prime}_{\rm{c}}}{\mathcal{U}^{\prime}_{\rm{c}}}\
\frac{\epsilon}{2\mathcal{z}_{\rm{c}}}\ ,$ (115)

where we select $+i\pi$ if $\mathcal{U}^{\prime}_{\rm{c}}>0$, and $-i\pi$
otherwise. The contribution to the real part of the exponential integral
multiplied by the decaying exponential is negligible, and hence is discarded
in equations (23) and (25).

We use equation (115) to calculate the left-hand side of the identity (113);

$\Im\big{\\{}\chi^{\prime}(0^{+})\big{\\}}=\Im\big{\\{}f_{\rm{unif},1}^{\prime}(0^{+})\big{\\}}+\text{h.o.t.}=-\pi\
\frac{\mathcal{U}^{\prime\prime}_{\rm{c}}}{|\mathcal{U}^{\prime}_{\rm{c}}|}\
e^{-2\mathcal{z}_{\rm{c}}/\epsilon}.$ (116)

Using the results of Appendix A, we have

$\chi_{\rm{c}}=e^{-\mathcal{z}_{\rm{c}}/\epsilon}\big{\\{}1+O(\epsilon)\big{\\}}.$
(117)

Using equation (117) in the right-hand side of equation (113), we recover the
left-hand side of equation (116).

##  Appendix C Short wavelength cut-off of the rippling instability

As shown in Figure 4, the rippling instability has a short wavelength or high
wavenumber cut-off when the effect of surface tension is taken into account.
Following Young and Wolfe (2013), we denote this cut-off $k^{+}_{\rm{neut}}$
where the subscript ‘neut’ denotes neutral. In other words, the wave with
wavenumber $k^{+}_{\rm{neut}}$ is marginally stable because its phase speed
equals the lower bound of the velocity profile in the water, namely

$c_{-}(k^{+}_{\rm{neut}})=\lim\limits_{z\to-\infty}U(z)\equiv U_{-\infty},$
(118)

where all variables are dimensional. Hence, using equation (30) for the phase
speed of short waves,
$k^{+}_{\rm{neut}}~{}L~{}\equiv~{}\mathcal{k}^{+}_{\rm{neut}}$ is solution of

$\frac{U_{\rm{S}}}{V}-\sqrt{\frac{\mathcal{S}\mathcal{k}}{1+r}}+\frac{r\mathcal{U}^{\prime}(0^{+})-\mathcal{U}^{\prime}(0^{-})}{2(1+r)\mathcal{k}}\
-\frac{\mathcal{G}(1-r)}{2\sqrt{\mathcal{S}(1+r)}}\
\mathcal{k}^{-\frac{3}{2}}=\frac{U_{-\infty}}{V}\ .$ (119)

For convenience, we introduce the following notation:

$\displaystyle\tilde{\mathcal{S}}\equiv\frac{\mathcal{S}}{1+r},\qquad$
$\displaystyle\tilde{\mathcal{G}}\equiv\frac{1-r}{1+r}\ \mathcal{G},$ (120)
$\displaystyle\mathcal{A}\equiv\frac{\mathcal{U}^{\prime}(0^{-})-r\mathcal{U}^{\prime}(0^{+})}{1+r}\quad\text{and}\quad$
$\displaystyle\mathcal{B}\equiv\frac{U_{\rm{S}}-U_{-\infty}}{V}\ .$ (121)

We assume a wind-induced current in the water so that $\mathcal{B}>0$.
Following Young and Wolfe (2013), we let the inverse Weber number be a small
parameter, $\tilde{\mathcal{S}}\ll 1$, and solve equation (119) in the form

$\sqrt{\tilde{\mathcal{S}}\mathcal{k}}\bigg{(}\mathcal{B}-\frac{\mathcal{A}}{2\mathcal{k}}\bigg{)}=\tilde{\mathcal{S}}\mathcal{k}+\frac{\tilde{\mathcal{G}}}{2\mathcal{k}},\qquad\tilde{\mathcal{S}}\ll
1.$ (122)

Taking the square of this result and letting
$X\equiv\tilde{\mathcal{S}}\mathcal{k}$ we obtain

$X\bigg{[}\mathcal{B}-\frac{\tilde{\mathcal{S}}\mathcal{A}}{2X}\bigg{]}^{2}=\bigg{[}X+\frac{\tilde{\mathcal{S}}\tilde{\mathcal{G}}}{2X}\bigg{]}^{2},\qquad\tilde{\mathcal{S}}\ll
1.$ (123)

Equation (123) can be readily solved using regular perturbation theory. We
seek solutions in the form

$X=X_{0}+\tilde{\mathcal{S}}\ X_{1}+O\big{(}\tilde{\mathcal{S}}^{2}\big{)}.$
(124)

We find

$X_{0}=\mathcal{B}^{2}\qquad\text{and}\qquad
X_{1}=-\frac{\mathcal{A}\mathcal{B}+\tilde{\mathcal{G}}}{2X_{0}-\mathcal{B}^{2}}\
,$ (125)

and infer

$\mathcal{k}^{+}_{\rm{neut}}=\frac{\mathcal{B}^{2}(1+r)}{\mathcal{S}}-\frac{\mathcal{U}^{\prime}(0^{-})-r\mathcal{U}^{\prime}(0^{+})}{(1+r)\mathcal{B}}-\frac{1-r}{1+r}\
\frac{\mathcal{G}}{\mathcal{B}^{2}}+O(\mathcal{S}).$ (126)

For the double exponential profile introduced in §II.3 we find

$\mathcal{k}^{+}_{\rm{neut}}=\frac{1+r}{\mathcal{S}}-\frac{1-r(\mathcal{R}_{1}-1)/\mathcal{R}_{2}}{1+r}-\frac{1-r}{1+r}\
\mathcal{G}+O(\mathcal{S}),$ (127)

which is a generalization of equation (5.10) in Young and Wolfe (2013).

## References

  * Miles (1957) J. W. Miles, J. Fluid Mech. 3, 185 (1957).
  * Young and Wolfe (2013) W. R. Young and C. L. Wolfe, J. Fluid Mech. 739, 276 (2013).
  * Bonfils  _et al._ (2022) A. F. Bonfils, D. Mitra, W. Moon, and J. S. Wettlaufer, J. Fluid Mech. 944, A8 (2022).
  * Kadam  _et al._ (2022) Y. Kadam, R. Patibandla, and A. Roy, J. Fluid Mech. (2022).
  * Shapiro and Teukolsky (2008) S. L. Shapiro and S. A. Teukolsky,  _Black holes, white dwarfs, and neutron stars: The physics of compact objects_ (John Wiley & Sons, 2008).
  * Casanova  _et al._ (2011) J. Casanova, J. José, E. Garcia-Berro, S. N. Shore, and A. C. Calder, Nature 478, 490 (2011).
  * Rosner  _et al._ (2001) R. Rosner, A. Alexakis, Y. Young, J. W. Truran, and W. Hillebrandt, Astrophys. J. Lett. 562, L177 (2001).
  * Alexakis  _et al._ (2004a) A. Alexakis, A. C. Calder, A. Heger, E. F. Brown, L. J. Dursi, J. W. Truran, R. Rosner, _et al._ , Astrophys. J. 602, 931 (2004a).
  * Alexakis  _et al._ (2004b) A. Alexakis, Y. Young, and R. Rosner, J. Fluid Mech. 503, 171 (2004b).
  * Drazin and Reid (1981) P. G. Drazin and W. H. Reid,  _Hydrodynamic stability_ (Cambridge University Press, 1981).
  * Wu (1975) J. Wu, J. Fluid Mech. 68, 49 (1975).
  * Bender and Orszag (1999) C. M. Bender and S. A. Orszag,  _Advanced Mathematical Methods for Scientists and Engineers_ (Springer, 1999).
  * Morland and Saffman (1992) L. C. Morland and P. G. Saffman, J. Fluid Mech. 252, 383 (1992).
  * Miles (1962) J. W. Miles, J. Fluid Mech. 13, 427 (1962).
  * Lin (1955) C. C. Lin,  _The Theory of Hydrodynamic Stability_ (Cambridge University Press, 1955).
  * Beji and Nadaoka (2004) S. Beji and K. Nadaoka, J. Fluid Mech. 500, 65 (2004).