# Periodic and solitary wave solutions of the long wave-short wave Yajima-
Oikawa-Newell model

Marcos Caso-Huerta
Department of Mathematics, Physics and Electrical Engineering
Northumbria University, Newcastle upon Tyne, UK
Email<EMAIL_ADDRESS>

Antonio Degasperis
Dipartimento di Fisica
“Sapienza” Università di Roma, Rome, Italy
Email<EMAIL_ADDRESS>

Priscila Leal da Silva
Department of Mathematical Sciences, School of Science
Loughborough University, Loughborough, UK
Email<EMAIL_ADDRESS>

Sara Lombardo
Department of Mathematical Sciences, School of Science
Loughborough University, Loughborough, UK
Email<EMAIL_ADDRESS>

Matteo Sommacal
Department of Mathematics, Physics and Electrical Engineering
Northumbria University, Newcastle upon Tyne, UK
Email<EMAIL_ADDRESS>
###### Abstract

Models describing long wave-short wave resonant interactions have many
physical applications from fluid dynamics to plasma physics. We consider here
the Yajima-Oikawa-Newell (YON) model, which has been recently introduced
combining the interaction terms of two long wave-short wave, integrable
models, one proposed by Yajima-Oikawa, and the other one by Newell.

The new YON model contains two arbitrary coupling constants and it is still
integrable – in the sense of possessing a Lax pair – for any values of these
coupling constants. It reduces to the Yajima-Oikawa or the Newell systems for
special choices of these two parameters. We construct families of periodic and
solitary wave solutions, which display the generation of very long waves. We
also compute the explicit expression of a number of conservation laws.

##  1 Introduction

The study of wave propagation poses quite a number of challenging different
mathematical and computational problems [1].

Wave motion of continuous media is generally represented by solutions of one
or more partial differential equations (PDEs). Typically, homogeneous wave
propagation equations have a linear part which is characterised by a
dispersion law, and a nonlinear one which is responsible for self and/or cross
interaction. While the linear part can be treated by decomposition into
Fourier harmonics, the nonlinear part – even if its dependence on the fields
is analytic – is generically treatable only by numerical methods. These two
parts, in terms of solvability, recombine in cooperation with each other for a
small set of very special wave equations, the so-called _integrable_ ones, by
allowing for a sort of nonlinear Fourier-like analysis [2, 3].

We have in mind (and consider below) only waves propagating in a one-
dimensional space, even if integrable wave equations in higher dimension are
known too. Since their first discovery more than half a century ago, the
firmament of integrable wave equations has been continuously growing. From the
very beginning, water waves have played a pivotal role in this research, with
the discovery of the integrability of the Boussinesq equation, the Korteweg-de
Vries (KdV) equation and the nonlinear Schrödinger (NLS) equation, which may
serve as approximate models of waves travelling on a water free surface
subject to various physical conditions.

The existence of solitons is the most celebrated and effective prediction from
the study of water wave equations of integrable type. Subsequent research has
pointed out that discovering integrable wave equations is not only due to a
lucky strike in the process of approximating complicated PDEs of physical
significance. In fact, integrable models can be obtained by means of a
perturbation approach, the so-called multiscale method (see [4, 5, 6]), when
applied to a given known nonlinear wave equation. Indeed, given a physically
relevant PDE, by appropriately rescaling both wave amplitudes and space-time
coordinates via the introduction of a small parameter $\epsilon$, and by
expanding in powers of $\epsilon$, one ends up with a different PDE which
models the amplitude modulation in the rescaled space-time coordinates. The
point is that this process preserves the integrability property of the
original wave equation and thus it yields a (possibly) _new_ integrable wave
equation. The rationale behind this method is based on physical arguments: one
considers the nonlinear terms as a perturbation of known Fourier-like
solutions (harmonics) of the linear dispersive part. Consequently, approximate
solutions of the original PDE look like superpositions of exponentials as, for
instance, $\sum_{j}\psi_{j}e^{i\,j\,(kx-\omega t)}$, whose amplitudes
$\psi_{j}$ depend on rescaled coordinates only. The parameters $k$ and
$\omega$ are the wave number and frequency on the dispersion curve
$\omega=\omega(k)$. If the superposition of harmonics contains more wave-
numbers, this multiscale method in general yields a system of coupled PDEs
which models wave-wave interactions.

A necessary condition for this to happen is the (weak) resonance relation:

$\omega(k_{1}+k_{2}+\cdots+k_{j})=\omega(k_{1})+\omega(k_{2})+\cdots+\omega(k_{j}),\quad
j>1\,.$

Several models of resonant interaction have been derived in this way (e.g.,
see [5, 9, 6]) to investigate their main features by means of known
integrability techniques, provided such models were selected as sufficiently
“close” to a wave equation of specific physical interest.

This perturbative approach is also appropriate to investigate the resonant
coupling of two quasi-monochromatic waves, one with very long wave-length, say
with wave number $k_{L}\approx 0$, and the second one with a much shorter
wave-length, say with wave number $k_{S}\gg k_{L}$. As originally pointed out
in [7], this interaction can be understood as a resonant triad $k_{1}$,
$k_{2}$, $k_{3}$, namely

$k_{1}+k_{2}=k_{3}\,,\quad\omega(k_{1})+\omega(k_{2})=\omega(k_{3})\,,$

with

$k_{1}=k_{L}\,,\qquad k_{2}=k_{S}-(1/2)k_{L}\,,\qquad
k_{3}=k_{S}+(1/2)k_{L}\,.$

Indeed, if the long wave is sufficiently long, say $k_{L}\rightarrow 0$, and
the dispersion function $\omega(k)$ is analytic at $k=0$, this condition is
equivalent to the stronger condition that the long wave and the short wave
have the same group velocity, $v_{L}=v_{S}$.

The search for integrable PDEs which reasonably model phenomena due to long
wave-short wave (L-S) interaction started in the early years of the soliton
era in fluid dynamics, plasma physics and optics. The Yajima-Oikawa (YO)
system [8]

$iS_{t}+S_{xx}-LS=0\,,\quad L_{t}=2(|S|^{2})_{x}$ (1)

was first derived in the one-way wave approximation in plasma. Here and in the
following $S$ and $L$ are the complex, and, respectively, real amplitudes of
the short and long waves. This system shows up also via multiscale technique
[9] in multiple ways. In fact, it proves to be a multiscale reduction of an
integrable equation of interest in water waves, namely the Boussinesq equation
[5]. A second and alternative integrable L-S wave system has been proposed by
Newell [10]. This one reads

$iS_{t}+S_{xx}+(iL_{x}+L^{2}-2\sigma|S|^{2})S=0\,,\quad
L_{t}=2\sigma(|S|^{2})_{x}\,,\quad\sigma^{2}=1\,,$ (2)

where, in addition to a long wave-short wave coupling, the short wave has the
same self-interaction as the NLS equation, which may be both defocusing
$(\sigma=1)$ or focusing $(\sigma=-1)$ according to the sign $\sigma$. We will
refer to the system (2) as to the N equation. However, the way to obtain this
system as a multiscale reduction of an integrable equation does not seem to be
known.

In [11], it has been shown that the YO and N equations, (1) and (2), need not
to be separately treated to investigate the short wave-long wave interaction.
Indeed, these two different model equations, remarkably enough, can be
combined in just one system, which we refer to as YON model, that is itself
integrable for _any_ real value of the two arbitrary parameters $\alpha$ and
$\beta$, namely

$iS_{t}+S_{xx}+\left(i\alpha L_{x}+\alpha^{2}L^{2}-\beta
L-2\alpha|S|^{2}\right)S=0\,,\quad L_{t}=2(|S|^{2})_{x}\,.$ (3)

This system coincides with the YO equation (1) for $\alpha=0$, $\beta=1$,
while it reads as the N equation (2) by setting $\alpha=\sigma$, $\beta=0$ and
by substituting the field $L$ with $\sigma\,L$. This unifying result provides
a greater flexibility in modelling the resonant interaction of long and short
waves, and allows to construct and analyse in one go special solutions of the
two models, (1) and (2). On the mathematical side, we note that the YON model
(3) turns out to be a reduction of a larger system of four coupled PDEs [12].

The preliminary step to inquire on the application of the YON system to
specific physical wave phenomena is the compatibility of the L-S resonance
condition (see the triad resonance above) with the dispersion law (for an
instance of such analysis in optics, see [15]). Having in mind a fluid
dynamical context, we recall a few elementary facts to show that water waves
on the free surface of a two-dimensional rectangular container require that
gravity be contrasted by surface tension [16]. Indeed, this effect allows the
resonance, which, however – at least in geophysical applications – is far from
being of experimental relevance, see [17]. As it happens, if the surface
tension is neglected, the dispersion law

$\omega^{2}=g\,k\,\tanh(h\,k)\,,$

where $g$ is the gravity acceleration and $h$ is the depth of the flat bottom,
does not allow for the L-S resonance since the group velocity
$\mathrm{d}\omega(k)/\mathrm{d}k$ is monotonically decreasing in the entire
range $0<k<+\infty$. Surface tension, if the wave length is sufficiently
short, may contrast gravity and change the dispersion law into

$\omega^{2}=g\,k\,\tanh(h\,k)\,\left[1+\frac{n\,k^{2}}{g\,\rho}\right]\,,$

where $n=0.074\,\mathrm{N}/\mathrm{m}$ is the surface tension constant and
$\rho$ is the water density. Although this change of the frequency dispersion
formally allows for the S-L resonance, the value of $k_{S}$ which satisfies
the strong resonance condition $v_{L}=v_{S}$ strongly depends on the depth
$h$. For instance, a short wave of about $1\,\mathrm{cm}$ length requires a
flat bottom of approximately $1\,\mathrm{cm}$ depth. For stratified fluids,
see for instance [18, 19].

Despite such resonant conditions lead to small effects in capillarity-gravity
wave propagation, we deem it of interest to investigate the YON model because
of its potential applicability while being integrable. Indeed it has been
shown [11] that there exists a Lax pair of equations for an auxiliary function
$\Psi(x,t,\lambda)$ (where $\lambda$ is a complex number called spectral
parameter or spectral variable)

$\Psi_{x}=X\,\Psi\,,\quad\Psi_{t}=T\,\Psi\,\,,$ (4)

whose compatibility condition

$X_{t}-T_{x}+[X\,,\,T]=0$ (5)

is equivalent to the YON system (3). Here $X(x,t,\lambda)$ and
$T(x,t,\lambda)$ are complex matrix-valued functions depending on the field
variables $S$ and $L$, and also polynomially on the spectral variable
$\lambda$,

$X=i\lambda\,X_{1}+X_{0}\,,\quad
T=(i\lambda)^{2}T_{2}+i\lambda\,T_{1}+T_{0}\,,$ (6)

where $X_{1}$ and $T_{2}$ denote constant, traceless, diagonal matrices

$X_{1}=\left(\begin{array}[]{ccc}1&0&0\\\ 0&0&0\\\
0&0&-1\end{array}\right)\,,\quad
T_{2}=\frac{i}{3}\left(\begin{array}[]{ccc}1&0&0\\\ 0&-2&0\\\
0&0&1\end{array}\right)\,,$ (7a) and the matrices $X_{0}$, $T_{1}$ and $T_{0}$
have the form $X_{0}=\left(\begin{array}[]{ccc}0&S&iL\\\ \alpha
S^{\ast}&0&S^{\ast}\\\ i\alpha^{2}L-i\beta&\alpha
S&0\end{array}\right)\,,\quad T_{1}=\left(\begin{array}[]{ccc}0&iS&0\\\
i\alpha S^{\ast}&0&-iS^{\ast}\\\ 0&-i\alpha S&0\end{array}\right)\,,$ (7b)
$T_{0}=\left(\begin{array}[]{ccc}-i\alpha|S|^{2}&-\alpha LS+iS_{x}&i|S|^{2}\\\
-\alpha^{2}LS^{\ast}+\beta S^{\ast}-i\alpha
S^{\ast}_{x}&2i\alpha|S|^{2}&-\alpha LS^{\ast}-iS^{\ast}_{x}\\\
i\alpha^{2}|S|^{2}&-\alpha^{2}LS+\beta S+i\alpha
S_{x}&-i\alpha|S|^{2}\end{array}\right)\,.$ (7c)

However, in the following, we will not make use of the Lax pair (4), and we
refer the interested reader to [11], where further details on the integrable
character of the system can be found.

In the next section we discuss several families of periodic and traveling wave
solutions of the YON system (3), including dark and bright solitons, as well
as rational solitons. Finally, in Section 3, using a set of multipliers, we
construct and exhibit a family of conservation laws. To this purpose, it is
instrumental the use of symmetry transformations of our system (3). Thus we
end this section by observing that the YON system (3) is trivially invariant
with respect to translations in space and time, as well as to rotations around
the origin in the $S$-plane. Moreover, the YON system is invariant under the
general transformation

$\left(x,\,t,\,S,\,L\right)\to\Big{(}\varepsilon^{-1}x,\,\varepsilon^{-2}t,\,\varepsilon\,\exp{\Big{[}i\frac{(\varepsilon^{2}-1)\beta^{2}t}{4\alpha^{2}}\Big{]}}\,\,S,\,\varepsilon
L-\frac{(\varepsilon-1)\beta}{2\alpha^{2}}\Big{)}$

where $\varepsilon\neq 0$ is an arbitrary parameter. Note that the limit
$\alpha\to 0$ of the above transformation can be obtained by mapping
$\varepsilon\to\exp{(-2\alpha^{2}\varepsilon)}$ first. If $\alpha=0$, we have
that (3), namely the YO system (1), is invariant under the additional scaling:

$\left(x,\,t,\,S,\,L\right)\to\left(\varepsilon^{2}\,x,\,\varepsilon^{4}\,t,\,\varepsilon^{-3}\,S,\,\varepsilon^{-4}\,L\right)\,,$

while, if $\beta=0$, say if (3) is the N system (2), the invariant scaling
transformation is

$\left(x,\,t,\,S,\,L\right)\to\left(\varepsilon\,x,\,\varepsilon^{2}\,t,\,\varepsilon^{-1}\,S,\,\varepsilon^{-1}\,L\right)\,,$

where $\varepsilon\neq 0$ is an arbitrary parameter.

##  2 Solutions of the YON model

Although system (3) is integrable and hence it allows for solutions to be
found by a variety of elegant and powerful solution techniques rooted into
integrability theory, for the purpose of this paper we use instead an Ansatz
to derive periodic and travelling wave solutions, without resorting to heavier
mathematical machineries. Our Ansatz naturally follows the form of the
periodic and travelling wave solutions of the nonlinear Schrödinger equation
and it reads:

$\displaystyle S(t,x)=s(z)e^{i(\phi(z)-\omega t)},\quad L(t,x)=\ell(z),$ (8)

where $z=x-Vt$ for $V\in\mathbb{R}$, $V\neq 0$, $\omega\in\mathbb{R}$ and $s$,
$\phi$, $\ell$ are real valued functions. After substitution of the Ansatz
into the second equation of (3), the resulting equation can be integrated in
order to obtain $\ell$ in terms of $s$:

$\displaystyle\ell(z)=-\frac{2}{V}s(z)^{2}+c_{1},$ (9)

where $c_{1}$ is an arbitrary integration constant. We now make use of this
new expression for $L$, substitute in the first equation of (3) and separate
it into real and imaginary parts in order to obtain the system

$\displaystyle
s(z)\phi^{\prime\prime}(z)+2s^{\prime}(z)\phi^{\prime}(z)-\left[V+\frac{4\alpha}{V}s(z)^{2}\right]s^{\prime}(z)=0,$
(10a) $\displaystyle
s^{\prime\prime}(z)-s(z)\phi^{\prime}(z)^{2}+V\,s(z)\phi^{\prime}(z)+(\alpha^{2}c_{1}^{2}-\beta
c_{1}+\omega)\,s(z)+\left(\frac{2\beta}{V}-\frac{4\alpha^{2}c_{1}}{V}-2\alpha\right)s(z)^{3}+\frac{4\alpha^{2}}{V^{2}}s(z)^{5}=0.$
(10b)

for the functions $s$ and $\phi$. Equation (10a) can be integrated with
respect to $z$ after multiplication by $s$, yielding an expression for the
first derivative of $\phi$ in terms of $s$:

$\displaystyle\phi^{\prime}(z)=\frac{\alpha}{V}\,s(z)^{2}+\frac{V}{2}+c_{2}\,s(z)^{-2},$
(11)

where $c_{2}$ is an integration constant. We now plug this expression for
$\phi^{\prime}$ into equation (10b), obtaining

$\displaystyle
s^{\prime\prime}(z)+\frac{1}{4V}\left[V^{3}+4V\left(\alpha^{2}c_{1}^{2}-\beta
c_{1}+\omega\right)-8\alpha c_{2}\right]\,s(z)+\frac{2}{V}(\beta-\alpha
V-2\alpha^{2}c_{1})\,s(z)^{3}+\frac{3\alpha^{2}}{V^{2}}\,s(z)^{5}-c_{2}^{2}\,s(z)^{-3}=0\,.$
(12)

Equation (12) can be integrated with respect to $z$ after multiplication by
$s^{\prime}$, resulting in a differential equation for $s^{\prime}(z)^{2}$:

$\displaystyle
s^{\prime}(z)^{2}=\displaystyle{-\frac{1}{4V}\left[V^{3}+4V\left(\alpha^{2}c_{1}^{2}-\beta
c_{1}+\omega\right)-8\alpha c_{2}\right]s(z)^{2}}+\frac{1}{V}(-\beta+\alpha
V+2\alpha^{2}c_{1})s(z)^{4}-\frac{\alpha^{2}}{V^{2}}s(z)^{6}+2c_{3}-c_{2}^{2}\,s(z)^{-2},$
(13)

where $c_{3}$ is an integration constant. Observe that for $\alpha=0$ the
coefficient of $s(z)^{6}$ becomes zero and the equation simplifies leading to
the Weierstraß elliptic function. The case $\alpha=0$ corresponds to the YO
model, well studied in the literature, e.g. see [14, 13] for some recent
results about periodic and rational solutions, and literature therein.
Here and thereafter we will assume $\alpha\neq 0$. Introducing the following
change of variable

$\dfrac{V}{2\alpha}u(z)=s(z)^{2}-u_{0}\,,\quad
u_{0}=\dfrac{V(\alpha\,V+2K_{1}\alpha^{2}-\beta)}{4\alpha^{2}}\,,\quad\alpha\neq
0$ (14)

equation (13) becomes, without any loss of generality,

$u^{\prime}(z)^{2}=-u(z)^{4}+\mu_{2}u(z)^{2}+\mu_{1}u(z)+\mu_{0}$ (15)

where

$\begin{split}\mu_{0}&=-\dfrac{V}{2\alpha^{3}}c^{3}+\dfrac{4\alpha^{2}\omega-8\alpha^{2}\,V^{2}-\beta^{2}}{4\alpha^{4}}c^{2}+\dfrac{\alpha^{2}\mu_{1}+2\alpha\,V^{2}(\mu_{2}+6\omega)-\alpha^{2}\,V^{3}-3\beta^{2}\,V}{2\alpha^{3}}c+\\\
&+\frac{1}{4}\left\\{2\mu_{1}\,V+V^{2}\left(4\omega-\dfrac{\beta^{2}}{\alpha^{2}}\right)-\dfrac{\left[\beta^{2}-\alpha^{2}(\mu_{2}+4\omega)\right]^{2}}{\alpha^{4}}\right\\}\end{split}$
(16)

whereas $c$, $\mu_{1}$, $\mu_{2}$ are arbitrary constants. Note that the
number of arbitrary constants is preserved. The integration constants $c_{1}$,
$c_{2}$ and $c_{3}$ can be rewritten in terms of $c$, $\mu_{1}$ and $\mu_{2}$
as follows

$c_{1}=\dfrac{c+\beta}{2\alpha^{2}},$ (17a)
$c_{2}=-\dfrac{V\left\\{c^{2}+6c\alpha\,V+2\beta^{2}+\alpha^{2}\left[V^{2}-2(\mu_{2}+4\omega)\right]\right\\}}{16\alpha^{3}}\,.$
(17b)
$c_{3}=\dfrac{V^{3}\left[-2\alpha^{3}\mu_{1}+2\alpha^{2}\mu_{2}(c+\alpha\,V)-(c+\alpha\,V)^{3}\right]}{32\alpha^{4}}\,.$
(17c)

As for $\phi(z)$, from (11) we obtain the quadrature

$\phi^{\prime}(z)=\dfrac{\alpha\left[\alpha\,V^{2}-c\,V+\alpha(\mu_{2}+4\omega)\right]-\beta^{2}+2\alpha
u(z)[c+2\alpha\,V+\alpha u(z)]}{2\alpha[c+\alpha\,V+2\alpha u(z)]}$ (18)

For $\mu_{1}=0$ equation (15) admits periodic solutions in terms of Jacobi
elliptic functions, which we are going to discuss below.

###  2.1 Jacobi elliptic sine solution

Let us assume that $u(z)$ has the form

$u(z)=\gamma_{0}+\gamma_{1}\operatorname{sn}\Big{(}a(z-z_{0}),m\Big{)}\,,$
(19)

where $\operatorname{sn}(z)$ denotes the Jacobi elliptic sine of $z$, where
$\gamma_{0}$, $\gamma_{1}$, $a$, $z_{0}$ and $m$ are real parameters, and
$0\leq m\leq 1$. Inserting (19) into (15) with $\mu_{1}=0$ and playing with
the properties of the Jacobi elliptic functions, one obtains a polynomial of
degree four in $\operatorname{sn}(a(z-z_{0}),m)$ equated to zero. Setting the
coefficients of each power of $\operatorname{sn}(z)$ to zero, one obtains a
set of algebraic equations for the parameters $\omega$, $m$, $a$,
$\gamma_{0}$, $\gamma_{1}$, $\mu_{2}$ and $c$. In particular, a relation for
$\omega$ in terms of the other parameters can be found by setting the constant
term of the polynomial to zero. This latter relation returns a non-real
$\omega$ for any choice of the other parameters, and therefore this excludes
the existence of a $\operatorname{sn}$-solution of the form (19) starting from
Ansatz (8).

###  2.2 Jacobi elliptic cosine solution

Proceeding as above but with the Jacobi elliptic cosine $\operatorname{cn}(z)$
replacing $\operatorname{sn}(z)$ in (19), we obtain the following solution to
(15) in:

$u(z)=ma\operatorname{cn}\Big{(}a(z-z_{0}),m\Big{)}\,,$ (20)

with

$\mu_{2}=(2m^{2}-1)a^{2}\,,\quad c=b-\alpha\,V\,,$ (21a) and

$\omega=\dfrac{2\beta^{2}+2\alpha^{2}\left[a^{2}(1-2m^{2})-2V^{2}\right]+4\alpha\,Vb+b^{2}\pm\sqrt{(b-2am\alpha)\,(b+2am\alpha)\,\left[b^{2}+4a^{2}(1-m^{2})\alpha^{2}\right]}}{8\alpha^{2}}\,,$
(21b)

where $m$, $z_{0}$, $a\neq 0$ and $b$ are real parameters, with $0\leq m\leq
1$. Replacing the expression of $u(z)$ into the original Ansatz (8), we obtain
the following solution of the YON system:

$\displaystyle S(x,t)$ $\displaystyle=\frac{1}{2}e^{i(\phi(z-\omega
t)}\sqrt{\dfrac{V\left[b+2m\alpha
a\operatorname{cn}\Big{(}a(z-z_{0}),m\Big{)}\right]}{\alpha^{2}}}\,,$ (22a)
$\displaystyle L(x,t)$ $\displaystyle=\dfrac{\beta-\alpha\,V-2m\alpha
a\operatorname{cn}\Big{(}a(z-z_{0}),m\Big{)}}{2\alpha^{2}}\,,\quad z=x-Vt\,,$
(22b) where $\phi(z)$ satisfies the quadrature

$\phi^{\prime}(z)=\dfrac{1}{4\alpha}\left\\{b+2\alpha\left[V+ma\operatorname{cn}\Big{(}a(z-z_{0}),m\Big{)}\right]\pm\dfrac{\sqrt{(b-2am\alpha)\,(b+2am\alpha)\,\left[b^{2}+4a^{2}(1-m^{2})\alpha^{2}\right]}}{2m\alpha
a\operatorname{cn}\Big{(}a(z-z_{0}),m\Big{)}+b}\right\\}\,,$ (22c)

where the sign in front of the square root is the same sign chosen for
$\omega$. Observe that, in addition to the coupling parameters $\alpha$ and
$\beta$ in the YON model (3), the solutions (22) features five real
parameters, namely $a$, $b$, $m$, $V$ and $z_{0}$, with a sixth real,
arbitrary parameter coming from the integration of (22c).

For this to work, we need that all $s$, $\phi$, $\omega$ be real. In order to
assure that, we need to check the sign inside all the square roots involved,
which gives us constraints on the parameters, namely

$\left|\dfrac{b}{\alpha a}\right|\geq 2m\,,\quad Vb\geq 0\,.$ (23)

Furthermore, in the special case $m=1$, the value $b=0$ is allowed as long as
$\alpha\,Va>0$, see Section 2.4.

The short wave $|S|$ oscillates between the values

$\dfrac{1}{2}\sqrt{\dfrac{V(b+2m\alpha
a)}{\alpha^{2}}}\qquad\mbox{and}\qquad\dfrac{1}{2}\sqrt{\dfrac{V(b-2m\alpha
a)}{\alpha^{2}}}\,,$

that is to say, $|S|^{2}$ oscillates with amplitude $\left|Vma/\alpha\right|$,
while the long wave $L$ oscillates between the values

$\dfrac{\beta-\alpha\,V-2m\alpha
a}{2\alpha^{2}}\qquad\mbox{and}\qquad\dfrac{\beta-\alpha\,V+2m\alpha
a}{2\alpha^{2}}\,,$

that is, with amplitude $2\left|ma/\alpha\right|$. The cn-solution is periodic
in $x$ with period $\left|4K(m)/a\right|$ and in $t$ with period
$\left|4K(m)/(aV)\right|$, where $K(m)$ is the complete elliptic integral of
the first kind of $m$,

$K(m)=\int_{0}^{\frac{\pi}{2}}\dfrac{d\theta}{\sqrt{1-m\sin^{2}\theta}}=\int_{0}^{1}\dfrac{dt}{\sqrt{(1-t^{2})(1-mt^{2})}}\,.$
(24)

Moreover, we observe that there are two special solutions corresponding to the
particular choices $m=0$ and $m=1$. If $m=0$, the elliptic cosine reduces to
the trigonometric cosine and the solution becomes a plane wave. If $m=1$ the
elliptic cosine reduces to the hyperbolic secant, leading to a localised
solution, treated in Section 2.4.

(a) Short wave $|S|$.

(b) Long wave $L$.

Figure 1: Elliptic cosine solution with $\alpha=1$, $\beta=2$, $V=1.2$, $b=5$,
$a=1.3$, $z_{0}=0$, $m=0.5$.

###  2.3 Jacobi delta amplitude solution

Proceeding as above with the Jacobi delta amplitude $\operatorname{dn}(z)$
replacing $\operatorname{sn}(z)$ in (19), we obtain the following solution to
(15)

$u(z)=\frac{a}{m}\operatorname{dn}\Big{(}a(z-z_{0}),m\Big{)}\,,$ (25)

with

$\mu_{2}=\left(\frac{2}{m^{2}}-1\right)a^{2}\,,\quad c=b-\alpha\,V\,,$ (26a)
and

$\omega=\dfrac{1}{8m^{2}\alpha^{2}}\Bigg{\\{}2\alpha^{2}\big{[}-2V^{2}m^{2}+(m^{2}-2)a^{2}\big{]}+4Vm^{2}\alpha
b+m^{2}(2\beta^{2}+b^{2})\pm\sqrt{(mb-2a\alpha)\,(mb+2a\alpha)\left[m^{2}b^{2}-4a^{2}(1-m^{2})\alpha^{2}\right]}\Bigg{\\}}\,,$
(26b)

where $m$, $z_{0}$, $a\neq 0$ and $b$, are real parameters, with $0\leq m\leq
1$.

The solution in this case has the form

$\displaystyle S(x,t)=\frac{1}{2}e^{i(\phi(z)-\omega
t)}\sqrt{\dfrac{V\left[mb+2\alpha
a\operatorname{dn}\Big{(}a(z-z_{0}),m\Big{)}\right]}{m\alpha^{2}}}\,,$ (27a)
$\displaystyle L(x,t)=\dfrac{m(\beta-\alpha\,V)-2\alpha
a\operatorname{dn}\Big{(}a(z-z_{0}),m\Big{)}}{2m\alpha^{2}}\,,\quad z=x-Vt$
(27b) where $\phi$ satisfies the quadrature

$\phi^{\prime}(z)=\dfrac{1}{4m\alpha}\left[2Vm\alpha+mb+2\alpha
a\operatorname{dn}\Big{(}a(z-z_{0}),m\Big{)}\pm\frac{\sqrt{(mb-2a\alpha)\,(mb+2a\alpha)\left[m^{2}b^{2}-4a^{2}(1-m^{2})\alpha^{2}\right]}}{2\alpha
a\operatorname{dn}\Big{(}a(z-z_{0}),m\Big{)}+mb}\right]\,,$ (27c)

where the sign in front of the square root is the same sign chosen for
$\omega$. Similarly to the Jacobi elliptic solutions (22), observe that the
solutions (27) features five real parameters, namely $a$, $b$, $m$, $V$ and
$z_{0}$, with a sixth real, arbitrary parameter coming from the integration of
(27c).

Again, checking the square roots that appear for having real solutions, we get
the following constraints on the parameters:

$2\sqrt{\dfrac{1-m}{m^{2}}}\leq\dfrac{b}{\alpha a}\leq
2\sqrt{\dfrac{1-m^{2}}{m^{2}}}\,,\quad\alpha\,Va>0\,,\quad 0<m\leq 1\,.$ (28)

It also allows the special values $\frac{b}{\alpha
a}=-2\sqrt{\frac{1-m}{m^{2}}}$ and $\frac{b}{\alpha a}=\frac{2}{m}$, for
$0<m\leq 1$. We discuss the special case $m=1$ in Section 2.4.

The dn-solution has periodicity for $L$ and $|S|$ with period
$\left|2K(m)/a\right|$, where $K(m)$ is the complete elliptic integral of the
first kind of $m$, while the phase of $S$ has a period $\left|4K(m)/a\right|$.

The short wave $|S|$ oscillates between the values

$\dfrac{1}{2}\sqrt{\dfrac{V[2\alpha
a(1-m)^{\nicefrac{{1}}{{2}}}+mb]}{m\alpha^{2}}}\qquad\mbox{and}\qquad\dfrac{1}{2}\sqrt{\dfrac{V(2\alpha
a+mb)}{m\alpha^{2}}}\,,$

that is to say, the oscillations in $|S|^{2}$ have an amplitude
$\left|Va(1-\sqrt{1-m})/(2\alpha)\right|$, while the long wave $L$ oscillates
between the values

$\dfrac{m(\beta-\alpha\,V)-2\alpha
a}{2m\alpha^{2}}\qquad\mbox{and}\qquad\dfrac{m(\beta-\alpha\,V)-2\alpha
a\sqrt{1-m}}{2m\alpha^{2}}\,,$

that is, with an amplitude $\left|a(1-\sqrt{1-m})/\alpha\right|$.

(a) Short wave $|S|$.

(b) Long wave $L$.

Figure 2: Delta amplitude solution with $\alpha=1$, $\beta=1$, $V=2$, $b=9$,
$a=3$, $z_{0}=-2$, $m=0.5$.

###  2.4 Travelling waves: solitons

The choice $m=1$ in (20) makes the period of the elliptic cosine diverge, so
the solution becomes localised. The corresponding solutions are solitons, both
of dark and bright type.

The solution for $m=1$, which for a generic choice of parameters corresponds
to a dark soliton, has the form

$\displaystyle S(x,t)=\frac{1}{2}e^{i(\phi(z)-\omega
t)}\sqrt{\dfrac{V\left[b+2\alpha
a\operatorname{sech}\Big{(}a(z-z_{0})\Big{)}\right]}{\alpha^{2}}}\,,$ (29a)
$\displaystyle L(x,t)=\dfrac{\beta-\alpha\,V-2\alpha
a\operatorname{sech}\Big{(}a(z-z_{0})\Big{)}}{2\alpha^{2}}\,,\quad z=x-Vt\,,$
(29b)
$\displaystyle\omega=\dfrac{1}{8\alpha^{2}}\Big{[}2\beta^{2}-2\alpha^{2}(2V^{2}+a^{2})+4V\alpha
b+b^{2}\pm\sqrt{b^{4}-4\alpha^{2}a^{2}b^{2}}\Big{]}\,,$ (29c)

with

$\displaystyle\phi(z)=$
$\displaystyle\dfrac{z-z_{0}}{4\alpha}\left(2\alpha\,V+b\pm\mathrm{sgn}(b)\,\sqrt{b^{2}-4\alpha^{2}a^{2}}\right)+\arctan\Big{(}\tanh\Big{(}\frac{a}{2}(z-z_{0})\Big{)}\Big{)}$
$\displaystyle\mp\mathrm{sgn}(b)\,\arctan\Bigg{(}\dfrac{(b-2\alpha
a)\tanh\Big{(}\frac{a}{2}(z-z_{0})\Big{)}}{\sqrt{b^{2}-4\alpha^{2}a^{2}}}\Bigg{)}+\phi_{0}\,,$
(29d)

where $\phi_{0}$ is an arbitrary phase, and the sign function $\mathrm{sgn}$
satisfies $\mathrm{sgn}(0)=0$. As a consequence of (23), and as it can be
observed from the formulae above, the general condition on the parameters for
the validity of the soliton solution, when $b\neq 0$, is

$\left|\dfrac{b}{\alpha a}\right|\geq 2\,.$ (30)

The special choice $b=0$ (see below) is also allowed by the system, though the
resulting solution has the phase:

$\phi(z)=\dfrac{V\,(z-z_{0})}{2}+\arctan\Big{(}\tanh\Big{(}\frac{a}{2}(z-z_{0})\Big{)}\Big{)}+\phi_{0}\,.$

The square of the short wave, $|S|^{2}$, has an amplitude
$\dfrac{1}{2}\left|\dfrac{Va}{\alpha}\right|$ over the background
$\left|\dfrac{Vb}{4\alpha}\right|$, while the long wave $L$ has an amplitude
$-\dfrac{a}{\alpha}$ over the background
$\dfrac{\beta-\alpha\,V}{2\alpha^{2}}$. Note that both amplitudes and the
background of $S$ do not depend on $\beta$ at all, while they all depend
inversely on $\alpha$.

By construction, both $S(x,0)$ and $L(x,0)$ are centred at $x=z_{0}$, while
$S(x,t_{0})$ and $L(x,t_{0})$ for a given $t_{0}$ are both centred at
$x=z_{0}+Vt_{0}$.

Whenever $b=0$, $S$ has zero background, and whenever $V=\beta/\alpha$, $L$
has zero background too. Both equalities being true means having a bright
soliton solution, both being false leads to a dark soliton solution, while the
cases where one is true and the other is not lead to mixed bright-dark
solutions.

The resulting formula for the bright case is

$\displaystyle S(x,t)=\frac{\sqrt{2}}{2}e^{i(\phi(z)-\omega
t)}\sqrt{\dfrac{\beta
a\operatorname{sech}\Big{(}a(z-z_{0})\Big{)}}{\alpha^{2}}}\,,$ (31a)
$\displaystyle
L(x,t)=-\dfrac{a\operatorname{sech}\Big{(}a(z-z_{0})\Big{)}}{\alpha}\,,\quad
z=x-\dfrac{\beta}{\alpha}t\,,$ (31b)
$\displaystyle\omega=-\dfrac{\beta^{2}+\alpha^{2}a^{2}}{4\alpha^{2}}\,,$ (31c)
$\displaystyle\phi(z)=\dfrac{\beta(z-z_{0})}{2\alpha}+\arctan\Big{(}\tanh\Big{(}\frac{a}{2}(z-z_{0})\Big{)}\Big{)}+\phi_{0}\,.$
(31d)

The same procedure can be carried out by taking $m=1$ in the dnoidal solution
(25). The solitons obtained in this way have the exact same formula as the
ones obtained from the cnoidal case.

(a) Short wave $|S|$.

(b) Long wave $L$.

Figure 3: Bright soliton solution, with $\alpha=0.5$, $\beta=2$, $a=0.25$,
$b=0$, $V=4$, $z_{0}=-20$.

###  2.5 Travelling waves: rational solution

By letting $c_{1}=c_{2}=c_{3}=0$ in (13) it reads

$s^{\prime}(z)^{2}=-\frac{\alpha^{2}}{V^{2}}s(z)^{2}\left[s(z)^{4}-\frac{V}{\alpha^{2}}(\alpha
V-\beta)+\frac{V^{2}}{3\alpha^{4}}(\alpha V-\beta)^{2}\right],$ (32)

where $\omega$ is given below. An integration process now yields

$\displaystyle\int_{s}^{0}\frac{3\sqrt{3}\,\alpha^{2}|\alpha|}{(V(\alpha
V-\beta)-3\alpha^{2}\zeta^{2})^{3/2}}d\zeta=\pm\sqrt{\frac{\alpha^{2}}{V^{2}}}z\,.$

Computing the integral, solving with respect to $s(z)$ and substituting into
the expressions for $\ell(z)$ and $\phi(z)$ leads to the solution

$\displaystyle S(x,t)=\frac{z}{\sqrt{3}\,\alpha}\,e^{i(\phi(z)-\omega
t)}\sqrt{\frac{V\,(\alpha V-\beta)^{3}}{9\alpha^{2}+(\alpha
V-\beta)^{2}z^{2}}},$ (33a) $\displaystyle L(x,t)=-\frac{2(\alpha
V-\beta)^{3}z^{2}}{3\alpha^{2}[9\alpha^{2}+(\alpha V-\beta)^{2}z^{2}]},\quad
z=x-Vt,$ (33b) $\displaystyle\omega=\frac{\alpha^{2}V^{2}-8\alpha\beta
V+4\beta^{2}}{12\alpha^{2}}\,,$ (33c)
$\displaystyle\phi(z)=\arctan\left(\frac{3\alpha}{(\alpha
V-\beta)z}\right)+\left(\frac{V}{2}+\frac{\alpha
V-\beta}{3\alpha}\right)z+\phi_{0},$ (33d)

where $\phi_{0}$ is an arbitrary phase, and the constraint $V(\alpha
V-\beta)>0$ is assumed.

The short wave $|S|$ is a dark rational solution with an amplitude depression
of $\sqrt{V(\alpha V-\beta)/(3\alpha^{2})}$ propagating on the non-vanishing
background $\sqrt{V(\alpha V-\beta)/(3\alpha^{2})}$, whereas the long wave $L$
has an amplitude $2(\alpha V-\beta)/(3\alpha^{2})$ over the asymptotic
background $-2(\alpha V-\beta)/(3\alpha^{2})$.

(a) Short wave $|S|$.

(b) Long wave $L$.

Figure 4: Rational solution with $\alpha=2$, $\beta=1$, $V=1$.

In the case where $\beta=0$ and we are back to the Newell system (2), the
solution is simplified to

$\displaystyle S(x,t)=\frac{z}{\sqrt{3}\,\alpha}\,e^{i(\phi(z)-\omega
t)}\sqrt{\frac{\alpha\,V^{4}}{9+V^{2}z^{2}}},$ (34a) $\displaystyle
L(x,t)=-\frac{2V^{3}z^{2}}{3\alpha[9+V^{2}z^{2}]},\quad z=x-Vt,$ (34b)
$\displaystyle\omega=\frac{V^{2}}{12}\,,$ (34c)
$\displaystyle\phi(z)=\arctan\left(\frac{3}{Vz}\right)+\left(\frac{V}{2}+\frac{\alpha
V}{3\alpha}\right)z+\phi_{0},$ (34d)

where $\alpha>0$. To the best of our knowledge, this is the first time such a
solution is derived for (2).

Another solution can be obtained for $V(\alpha V-\beta)>0$ but now for
$\omega=-c^{2}/4$, still with $c_{1}=c_{2}=c_{3}=0$. The quadrature (13) is
now rewritten as

$s^{\prime}(z)^{2}=-\frac{\alpha^{2}}{V^{2}}s(z)^{4}\left[s(z)^{2}-\frac{V}{\alpha^{2}}(\alpha
V-\beta)\right],$

and integration leads to the solutions

$\displaystyle S(x,t)=e^{i(\phi(z)-\omega t)}\sqrt{\frac{V(\alpha
V-\beta)}{\alpha^{2}+(\alpha V-\beta)^{2}z^{2}}},$ (35a) $\displaystyle
L(t,x)=-\frac{2(\alpha V-\beta)}{\alpha^{2}+(\alpha V-\beta)^{2}z^{2}},\quad
z=x-Vt,$ (35b) $\displaystyle\omega=-c^{2}/4,$ (35c)
$\displaystyle\phi(z)=\frac{V}{2}z+\arctan\left(\frac{(\alpha
V-\beta)z}{\alpha}\right)+\phi_{0}\,.$ (35d)

The short wave $|S|$ and long wave $L$ are bright rational solutions with
amplitudes of $\sqrt{V(\alpha V-\beta)/(\alpha^{2})}$ and $-2(\alpha
V-\beta)/(\alpha^{2})$, respectively, on a zero background. To the best of our
knowledge, this is also a novel solution of (2).

##  3 Conservation laws

The YON model (3) is integrable and therefore it allows infinitely many
conservation laws. In this section, we are interested in deriving a few
explicit ones, by finding convenient multipliers. A conservation law of (3)
corresponds to an expression

$\rho_{t}+f_{x}=0,$ (36)

where $\rho\equiv\rho(S,L,S_{x},L_{x},\dots)$ is the density and $f\equiv
f(S,L,S_{x},L_{x},\dots)$ the corresponding flux, respectively. For instance,
a trivial conservation law of the YON system (3) is

$\displaystyle\rho_{0}=L,$ (37a) $\displaystyle f_{0}=-2|S|^{2}$ (37b)

which coincides with the second equation in (3).

In [20, 21, 22], it was established that every conservation law corresponds to
a symmetry of the equation, but very often this symmetry is non-classical or
even non-local. For cases where Noether’s theorem is applicable, this
correspondence is explicit as Noether symmetries lead to (possibly trivial)
conserved vectors; however, we ruled out this approach as we have not been
able to find a variational formulation for system (3). Instead of using
Noether’s approach, we will consider the direct method illustrated in [20, 21,
22, 28], consisting in finding a vector $g=(g_{1},g_{2},g_{3})$, called
multiplier, that depends on $S$, $L$ and their derivatives up to a certain
fixed but arbitrary order, such that

$\frac{\delta g_{1}\mathcal{F}_{1}}{\delta S}=0,\quad\frac{\delta
g_{2}\mathcal{F}_{2}}{\delta S^{\ast}}=0,\quad\frac{\delta
g_{3}\mathcal{F}_{3}}{\delta L}=0,$ (38)

where

$\mathcal{F}_{1}=iS_{t}+S_{xx}+\left(i\alpha L_{x}+\alpha^{2}L^{2}-\beta
L-2\alpha|S|^{2}\right)S\,,\quad\mathcal{F}_{2}=\mathcal{F}^{\ast}_{1}\,,\quad\mathcal{F}_{3}=L_{t}-2(|S|^{2})_{x}\,,$
(39)

and where $\delta/\delta u$ is the variational derivative with respect to the
variable $u$. If such a vector $g$ can be found, then it ensures the existence
of a $\rho$ and $f$ such that

$\rho_{t}+f_{x}=g_{1}\mathcal{F}_{1}+g_{2}\mathcal{F}_{2}+g_{3}\mathcal{F}_{3}=0\,.$
(40)

By making use of GeM [23, 24, 25, 26, 27], for any choices of
$\alpha,\beta\in\mathbb{R}$, we can find some pairs of conserved densities and
fluxes depending on derivatives up to the second order, the simples ones of
which read

$\displaystyle\rho_{1}=\frac{\alpha}{2}L^{2}-|S|^{2},$ (41a) $\displaystyle
f_{1}=-2\alpha L|S|^{2}-2\,Im(S^{\ast}S_{x}),$ (41b) $\displaystyle\rho_{2}$
$\displaystyle=2\alpha^{2}L|S|^{2}-\beta|S|^{2}+2\alpha\,Im(S^{\ast}S_{x}),$
(42a) $\displaystyle f_{2}$
$\displaystyle=2\alpha|S_{x}|^{2}-2\alpha\,Re(S^{\ast}S_{xx})+4\alpha^{2}L\,Im(S^{\ast}S_{x})-2\beta\,Im(S^{\ast}S_{x}),$
(42b)

while the remaining ones can be found in Appendix A.

##  4 Conclusions

Models describing long wave-short wave resonant interactions arise in a
variety of physical contexts, from fluid dynamics to plasma physics. In this
paper we consider the recently-proposed, long wave-short wave YON (Yajima-
Oikawa-Newell) model (see [11]), an integrable model featuring two arbitrary
parameters, and unifying and generalising the Yajima-Oikawa model and the
Newell model.

We studied some relevant families of periodic and solitary wave solutions,
displaying the generation of very long waves. Among others, we also display
the expression of solutions that we term, with some abuse of language,
“rational”. Differently from the NLS equation, where the amplitude $|S|$ is
indeed rational (cf., the Peregrine soliton), in the present case it is rather
the function $|S|^{2}$ that comes to be rational. This is due to the quintic
nonlinearity appearing in the YON system for the short wave amplitude $|S|$,
rather then the usual cubic one as in the NLS equation. An analytical study of
the stability of the solutions presented in this paper is left to future
investigation. In this respect, we limit ourselves to report here that we
carried out a preliminary numerical study, solving the initial value problem
for initial conditions obtained by computing our solutions at $t=0$, using the
method of lines with pseudospectral, Fourier discretisation in space and an
adaptive Dormand-Prince embedded Runge-Kutta method for the time stepping: the
numerical results seem to suggest the existence of regions of stability and
regions where different forms of instability are observed, similarly to what
is predicted for plane wave solutions of the YON system [11, 29].

The families of explicit solutions presented in this paper have been obtained
by choosing a suitable Ansatz. A systematic derivation of soliton solutions of
bright, gray and dark type, as well as of breathers and rogue-waves,
exploiting the integrable character of the YON system, is currently in
progress.

In this paper we have also derived a few conservation laws, which are of
interest in view of numerical and analytical studies of this system. An
argument based on the effect of the surface tension on the dispersion relation
for short waves, allowing for short-long wave resonance, is presented to
justify the physical relevance of the YON model in a fluid-dynamical context
(and in particular for experimental set-ups of capillarity-gravity wave
propagation on the scales of the centimetres), where the value of the wave
number strongly depends on the water depth. In spite of the expected physical
relevance, a derivation, via multiscale techniques, of the full YON system –
similarly to the Newell model, which is contained within the YON system – as
an (integrable) reduction from a known physical PDE has not yet been achieved
and remains an intriguing open problem. It is worth observing that a subcase
of the YON model, namely the Yajima-Oikawa model, has been indeed derived via
the multiscale technique in more than one way [9], suggesting that this should
be possible also for the more general YON model.

Acknowledgments: SL and PLdS acknowledge support by the Royal Society; PLdS is
supported by the Royal Society under a Newton International Fellowship
(reference number 201625) hosted by SL. The work of MS has been carried out
under the auspices of the Italian GNFM (Gruppo Nazionale Fisica Matematica),
INdAM (Istituto Nazionale di Alta Matematica).

##  Appendix A Appendix: Conserved vectors depending on derivatives up to the
second order

In this Appendix we list the remaining densities and fluxes (in addition to
(41) and (42)) depending on derivatives up to the second order.These have been
computed using GeM [23, 24, 25, 26, 27], for $\alpha,\beta\in\mathbb{R}$.

$\displaystyle\rho_{3}=$ $\displaystyle
2\alpha^{2}tL|S|^{2}+\frac{\alpha^{2}}{2}xL^{2}-\frac{\beta}{2}xL-\beta
t|S|^{2}-\alpha x|S|^{2}+2\alpha t\,Im(S^{\ast}S_{x}),$ (43a) $\displaystyle
f_{3}=$ $\displaystyle-2\alpha^{2}xL|S|^{2}+\beta
x|S|^{2}+4\alpha^{2}tL\,Im(S^{\ast}S_{x})+2\alpha t|S_{x}|^{2}-2\alpha
x\,Im(S^{\ast}S_{x})$ (43b) $\displaystyle-2\beta t\,Im(S^{\ast}S_{x})-2\alpha
t\,Re(S^{\ast}S_{xx});$

$\displaystyle\rho_{4}=$
$\displaystyle\frac{1}{2}\alpha^{5}L^{4}-2\alpha^{4}L^{2}|S|^{2}-\alpha^{3}\beta
L^{3}+2\alpha^{3}|S|^{4}+4\alpha^{3}L\,Im(S^{\ast}S_{x})-\frac{1}{2}\alpha^{3}L_{x}^{2}+2\beta^{2}|S|^{2}-4\alpha\beta\,Im(S^{\ast}S_{x})+4\alpha^{2}|S_{x}|^{2},$
(44a) $\displaystyle f_{4}=$ $\displaystyle 6\alpha^{3}\beta
L^{2}|S|^{2}+8\alpha^{4}L|S|^{4}-4\alpha^{2}\beta|S|^{4}-4\alpha\beta|S_{x}|^{2}+4\alpha^{3}L|S_{x}|^{2}-4\alpha^{5}L^{3}|S|^{2}+8\alpha^{3}|S|^{2}\,Im(S^{\ast}S_{x})$
(44b)
$\displaystyle+4\beta^{2}\,Im(S^{\ast}S_{x})-4\alpha^{4}L^{2}\,Im(S^{\ast}S_{x})+8\alpha^{3}L_{x}\,Re(S^{\ast}S_{x})+4\alpha\beta\,Re(S^{\ast}S_{xx})-4\alpha^{3}L\,Re(S^{\ast}S_{xx})+8\alpha^{2}\,Im(S^{\ast}_{x}S_{xx}).$

## References

  * [1] Whitham, G. B., Nonlinear dispersion of water waves, J. Fluid Mech., 1966, 27, 399–412.
  * [2] Ablowitz, M. J., Kaup, D. J., Newell, A. C., Segur, H., The inverse scattering transform-Fourier analysis for nonlinear problems, Stud. Appl. Math., 1974, 53(4), 249–315.
  * [3] Calogero, F., Degasperis, A., Spectral Transform and Solitons. Amsterdam: North-Holland, 1982.
  * [4] Zakharov, V. E., Kuznetsov, E. A., Multi-scale expansions in the theory of systems integrable by the inverse scattering transform, Physica D, 1986, 18D, 455–463.
  * [5] Calogero, F., Why are certain nonlinear PDEs both widely applicable and integrable?, in What is integrability? (ed. V. E. Zakharov). Berlin: Springer-Verlag, 1991.
  * [6] Degasperis, A., Multiscale expansion and integrability of dispersive wave equations, in Integrability, vol. 767 (ed. A. V. Mikhailov). Lecture Notes in Physics. Berlin: Springer-Verlag, 2009, 215–244.
  * [7] Benney, D. J., A general theory for interactions between short and long waves, Stud. Appl. Math., 1976, 56(1), 81–94.
  * [8] Yajima, N., Oikawa, M., Formation and interaction of sonic-Langmuir solitons: Inverse scattering method, Prog. Theor. Phys., 1976, 56, 1719–1739.
  * [9] Calogero, F., Degasperis, A., Xiaoda, J., Nonlinear Schrödinger-type equations from multiscale reduction of PDEs. I. Systematic derivation, J. Math. Phys., 2000, 41, 6399–6443.
  * [10] Newell, A. C., Long waves-short waves: A solvable model, SIAM J. Appl. Math., 1978, 35, 650–664.
  * [11] Caso-Huerta, M., Degasperis, A., Lombardo, S.,Sommacal, M., A new integrable model of long wave-short wave interaction and linear stability spectra, Proc. R. Soc. A., 2021, 477, 20210408.
  * [12] Wright, O. C., Homoclinic connections of unstable plane waves of the long-wave–short-wave equations, Stud. Appl. Math., 2006, 117, 71–93 [in particular, see system (53)].
  * [13] Chen, J., Chen, Y., Feng, B-F., Maruno, K., Ohta, Y., General high-order rogue waves of the (1+ 1)-dimensional Yajima-Oikawa system, J. Phys. Soc. Jpn., 2018, 87(9), 094007.
  * [14] Li, R., Geng, X., Periodic-background solutions for the Yajima-Oikawa long-wave–short-wave equation, Nonlinear Dyn., 2022, 94.
  * [15] Chowdhury, A., Tataronis, J. A., Long-wave short-wave resonance in nonlinear negative refractive index media, Phys. Rev. Lett., 2008, 100(15), 153905.
  * [16] Djordjevic, V. D., Redekopp, L. G., On two-dimensional packets of capillary-gravity waves, J. Fluid Mech., 1977, 79(4), 703–714.
  * [17] Lannes, D., The Water Waves Problem: Mathematical Analysis and Asymptotics. Mathematical Surveys and Monographs, vol. 188. Providence: American Mathematical Society, 2013.
  * [18] Grimshaw, R. H. J., The modulation of an internal gravity-wave packet, and the resonance with the mean motion, Stud. Appl. Math., 1977, 56, 241–266.
  * [19] Koop, C. G. and Redekopp, L. G., The interaction of long and short internal gravity waves: Theory and experiment, J. Fluid. Mech. , 1981, 367–409.
  * [20] Anco, S. C., Bluman, G. W., Direct construction of conservation laws from field equations, Phys. Rev. Lett., 1997, 78, 2869–2873.
  * [21] Anco, S. C., Bluman, G. W., Direct construction method for conservation laws of partial differential equations. I: Examples of conservation law classifications, Eur. J. Appl. Math., 2002, 13, 545–566.
  * [22] Anco, S. C., Bluman, G. W., Direct construction method for conservation laws of partial differential equations. II: General treatment, Eur. J. Appl. Math., 2002, 13, 567–585.
  * [23] Cheviakov, A., GeM software package for computation of symmetries and conservation laws of differential equations, Comput. Phys. Commun., 2007, 176, 48–61.
  * [24] Cheviakov, A., Symbolic computation of local symmetries of nonlinear and linear partial and ordinary differential equations, Math. Comput. Sci., 2010, 4, 203–222.
  * [25] Cheviakov, A., Computation of fluxes of conservation laws, J. Eng. Math., 2010, 66, 153–173.
  * [26] Cheviakov, A., Symbolic computation of nonlocal symmetries and nonlocal conservation laws of partial differential equations using the GeM package for Maple, in Similarity and Symmetry Methods (eds J.-F. Ganghoffer, I. Mladenov). Lecture Notes in Applied and Computational Mechanics, vol. 73. Cham: Springer, 2014, 165–184.
  * [27] Cheviakov, A., Symbolic computation of equivalence transformations and parameter reduction for nonlinear physical models, Comput. Phys. Commun., 2017, 220, 56–73.
  * [28] Olver, P. J., Applications of Lie Groups to Differential Equations, 2nd edition. Graduate Texts in Mathematics, vol. 107. New York: Springer, 1993.
  * [29] Degasperis, A., Lombardo, S., Sommacal, M., Integrability and Linear Stability of Nonlinear Waves, J. Nonlinear Sci., 2018, 28(4), 1251–1291.


††thanks: corresponding author, co-first author††thanks: corresponding author,
co-first author

# Multi-angle Quantum Approximate Optimization Algorithm

Rebekah Herrman<EMAIL_ADDRESS>Department of Industrial and Systems
Engineering
University of Tennessee at Knoxville
Knoxville, TN 37996  Phillip C. Lotshaw<EMAIL_ADDRESS>Quantum Computing
Institute
Oak Ridge National Laboratory
Oak Ridge, TN 37830  James Ostrowski  Department of Industrial and Systems
Engineering
University of Tennessee at Knoxville
Knoxville, TN 37996  Travis S. Humble  Quantum Computing Institute
Oak Ridge National Laboratory
Oak Ridge, TN 37830  George Siopsis  Department of Physics and Astronomy
University of Tennessee at Knoxville
Knoxville, TN 37996

###### Abstract

The quantum approximate optimization algorithm (QAOA) generates an approximate
solution to combinatorial optimization problems using a variational ansatz
circuit defined by parameterized layers of quantum evolution. In theory, the
approximation improves with increasing ansatz depth but gate noise and circuit
complexity undermine performance in practice. Here, we introduce a multi-angle
ansatz for QAOA that reduces circuit depth and improves the approximation
ratio by increasing the number of classical parameters. Even though the number
of parameters increases, our results indicate that good parameters can be
found in polynomial time. This new ansatz gives a 33% increase in the
approximation ratio for an infinite family of MaxCut instances over QAOA. The
optimal performance is lower bounded by the conventional ansatz, and we
present empirical results for graphs on eight vertices that one layer of the
multi-angle anstaz is comparable to three layers of the traditional ansatz on
MaxCut problems. Similarly, multi-angle QAOA yields a higher approximation
ratio than QAOA at the same depth on a collection of MaxCut instances on fifty
and one-hundred vertex graphs. Many of the optimized parameters are found to
be zero, so their associated gates can be removed from the circuit, further
decreasing the circuit depth. These results indicate that multi-angle QAOA
requires shallower circuits to solve problems than QAOA, making it more viable
for near-term intermediate-scale quantum devices.

††This manuscript has been authored by UT-Battelle, LLC under Contract No. DE-
AC05-00OR22725 with the U.S. Department of Energy. The United States
Government retains and the publisher, by accepting the article for
publication, acknowledges that the United States Government retains a non-
exclusive, paid-up, irrevocable, world-wide license to publish or reproduce
the published form of this manuscript, or allow others to do so, for United
States Government purposes. The Department of Energy will provide public
access to these results of federally sponsored research in accordance with the
DOE Public Access Plan. (http://energy.gov/downloads/doe-public-access-plan).

##  I Introduction

Among several quantum algorithms implemented on noisy intermediate-scale
quantum (NISQ) devices cruz2019efficient ; zhang2020efficient ;
godfrin2017operating ; borle2020quantum ; karamlou2020analyzing ;
hempel2018quantum ; alderete2020quantum ; xue2009quantum ;
linke2017experimental ; pagano2020quantum ; bengtsson2019quantum ;
harrigan2021quantum , the quantum approximate optimization algorithm (QAOA)
offers an opportunity to approximately solve combinatorial optimization
problems such as MaxCut, Max Independent Set, and Max k-cover farhi2014quantum
; saleem2020 ; cook2020quantum ; crooks2018performance ; pichler2018quantum ;
farhi2020quantumwholegraph ; farhi2020quantum ; wurtz2020bounds ; ward2018qaoa
; shaydulin2019evaluating . QAOA tunes a set of classical parameters to
optimize the cost function expectation value for a quantum state prepared by
well-defined sequence of operators acting on a known initial state. Variations
to the original algorithm include alternative operators and initial states
bartschi2020grover ; hadfield2019quantum ; wurtz2021classically ; wang2019xy ;
egger2021warm ; zhu2021improving ; tate2020bridging ; rieffel2020xy ;
golden2021threshold  while purely classical aspects such as the parameter
optimization and problem structure have been tested as well herrman2021impact
; shaydulin2020classical ; alexeev2020reinforcement ; wauters2020reinforcement
; lotshaw2021bfgs . However, an outstanding concern is that practical
implementations of QAOA require large numbers of qubits and deep circuits
guerreschi2019qaoa ; herrman2021lower . Noise grows rapidly with circuit depth
and affects the fidelity of the prepared quantum state xue2021effects ;
wang2020noise ; murali2019noise ; sun2021mitigating ;
marshall2020characterizing ; alam2020design ; alam2019analysis ;
streif2021quantum ; maciejewski2021modeling .

One approach to reduce the circuit depth of QAOA is to increase the number of
classical parameters introduced in each layer, a variation that we term multi-
angle QAOA (ma-QAOA). Increasing the number of classical parameters allows for
finer-grain control over the optimization of the cost function and the
approximation ratio, which measures optimality relative to the known best
solution. While introducing more classical parameters can lead to a more
challenging optimization, a corresponding reduction in circuit depth preserves
the critical resource of the quantum state. In addition, finding the absolute
optimal angles is not necessary in order to see an improvement over QAOA.

Here, we quantify the advantages of using multiple parameters for each layer
of QAOA. First, we prove that the approximation ratio converges to one as the
number of iterations of ma-QAOA tends to infinity, a property that ensures the
optimal solution is the most likely. We next demonstrate that one iteration of
ma-QAOA gives an approximation ratio that is at least that of the
approximation ratio after one iteration of QAOA. This shows that ma-QAOA
performs at least as well as QAOA. We also show that ma-QAOA used to solve the
MaxCut problem on star graphs achieves an approximation ratio of one after one
iteration, while single-iteration QAOA tends to an approximation ratio of 0.75
as the number of vertices goes to infinity. This result gives a concrete
example where ma-QAOA gives a strictly larger approximation ratio than QAOA.
We simulate solving MaxCut using ma-QAOA and QAOA on all connected, non-
isomorphic eight vertex graphs and compare the performance of the two
ansatzes. In doing so, we find that the average approximation ratio for ma-
QAOA after one iteration is larger than the average approximation ratio of
QAOA after three iterations. In looking at larger, 50- and 100-vertex graphs,
we see that ma-QAOA retains its advantage over QAOA, giving approximation
ratios that are on average six percentage points higher after the first
iteration.

##  II Results

### Quantum approximate optimization algorithm

The quantum approximate optimization algorithm (QAOA) relies on a combination
of classical parameter optimization and applying cost and mixing operators to
a quantum state in order to approximately solve combinatorial optimization
(CO) problems farhi2014quantum . The algorithm is implemented by applying two
operators,

$U(C,\gamma)=e^{-iC\gamma}$

and

$U(B,\beta)=e^{-iB\beta}$

in succession on an initial state,

${\left|s\right\rangle}=\frac{1}{\sqrt{2^{n}}}\sum_{z}{\left|z\right\rangle}.$

Here the sum is over the computational basis ${\left|z\right\rangle}$. QAOA
applied $p$ times to ${\left|s\right\rangle}$ is denoted $p$-QAOA. The state
for $p$-QAOA is

${\left|\gamma,\beta\right\rangle}=U(B,\beta_{p})U(C,\gamma_{p})...U(B,\beta_{1})U(C,\gamma_{1}){\left|s\right\rangle}.$

$C$ encodes the problem to be solved and $B$ drives transitions between
computational basis states. Often, $C$ is the sum over a collection of
clauses,

$C=\sum_{a}C_{a},$

and $B$ is typically

$B=\sum_{v\in V(G)}B_{v}$

where $B_{v}=\sigma_{v}^{x}$ is the Pauli-x operator acting on the $v^{th}$
qubit. Thus, we may write

$U(C,\gamma)=e^{-i\gamma\sum_{a}C_{a}}$

and

$U(B,\beta)=e^{-i\beta\sum_{v\in V(G)}B_{v}}.$

Instead of focusing on minimizing the classical optimization efforts in QAOA,
we modify QAOA such that it requires more classical parameters. The new
classical parameters are introduced to QAOA by allowing each summand of the
cost and mixing operators to have its own angle instead of a single angle for
the cost operator and a second angle for the mixing operator. In this
modification,

$U(C,\vec{\gamma_{l}})=e^{-i\sum_{a}C_{a}\gamma_{l,a}}=\prod_{a}e^{-i\gamma_{l,a}C_{a}}$

and

$U(B,\vec{\beta_{l}})=e^{-i\sum_{v\in V(G)}B_{v}\beta_{l,v}}=\prod_{v\in
V(G)}e^{-i\beta_{l,v}B_{v}}$

where $\vec{\gamma_{l}}=(\gamma_{l,a_{1}},\gamma_{l,a_{2}},...)$ and
$\vec{\beta_{l}}=(\beta_{l,v_{1}},\beta_{l,v_{2}},...)$. Here, $l$ denotes the
layer, $a_{i}$ denotes an edges in the graph, and $v_{j}$ refers to a specific
vertex. We call this modification multi-angle QAOA and abbreviate it ma-QAOA.

### Convergence of ma-QAOA

For QAOA, the expected value of $C$ after $p$ iterations is $\langle
C\rangle_{p}={\left\langle\gamma,\beta\right|}C{\left|\gamma,\beta\right\rangle}$.
Let $M_{p}$ be the maximum of $\langle C\rangle_{p}$ over all angles. Then,
$M_{p}\geq M_{p-1}$. Farhi, Goldstone, and Gutmann showed that $M_{p}$ tends
to the maximum of the objective function, $C_{\mathrm{max}}$, for the CO
problem being solved as $p$ tends to infinity farhi2014quantum .

We similarly define the expected value of $C$ after $p$ iterations of ma-QAOA
as $\langle
C\rangle_{p}^{\mathrm{ma}}={\left\langle\vec{\gamma}_{\mathrm{ma}},\vec{\beta}_{\mathrm{ma}}\right|}C{\left|\vec{\gamma}_{\mathrm{ma}},\vec{\beta}_{\mathrm{ma}}\right\rangle}$
where
$\vec{\gamma}_{\mathrm{ma}}=(\vec{\gamma}_{1},\vec{\gamma}_{2},...\vec{\gamma}_{p})$
and
$\vec{\beta}_{\mathrm{ma}}=(\vec{\beta}_{1},\vec{\beta}_{2},...\vec{\beta}_{p})$.
We also define $M_{p}^{\mathrm{ma}}$ to be the maximum of $\langle
C\rangle_{p}^{\mathrm{ma}}$ over all angles. Clearly, $M_{p}^{\mathrm{ma}}\geq
M_{p}$ because QAOA is the special case of ma-QAOA where
$\beta_{p,u}=\beta_{p,v}$ for all $u\neq v$ and
$\gamma_{p,a_{i}}=\gamma_{p,a_{j}}$ for edges $a_{i}\neq a_{j}$.

In order to show ma-QAOA gives the optimal solution to a combinatorial
optimization problem, we must show $\langle C\rangle_{p}^{\mathrm{ma}}$
converges to $C_{\mathrm{max}}$ as $p$ tends to infinity. Convergence is the
first main result of this work.

######  Theorem II.1.

The multi-angle quantum approximate optimization algorithm converges to the
optimal solution of a combinatorial optimization problem as
$p\rightarrow\infty$.

The proof of convergence is given in Sec. IV.

### MaxCut on star graphs

Recall that the approximation ratio for $p$-QAOA is equal to $\langle
C\rangle_{p}/C_{\mathrm{max}}$, where $C_{\mathrm{max}}$ is the optimal
solution to the CO problem. Similarly, we define the approximation ratio for
$p$ iterations of ma-QAOA as $\langle
C\rangle_{p}^{\mathrm{ma}}/C_{\mathrm{max}}$. The approximation ratio is often
used to measure the success of QAOA.

In graph theory, a star graph on $n$ vertices is a graph that consists of one
vertex of degree $n-1$, called the center. All other vertices of the graph
have degree one, meaning each vertex is connected to the center and only the
center. An example can be seen in Fig. 1. All stars are trees, and are thus
bipartite, so the optimal MaxCut solution includes all edges of the graph. In
order to show ma-QAOA outperforms QAOA when solving MaxCut on star graphs, we
show that $\langle C\rangle_{1}^{\mathrm{ma}}=1$ and $\langle C\rangle_{1}$
tends to $0.75$ as $n$ tends to infinity. The proof is found in Sec. IV.

Figure 1: The star graph on five vertices.

### Computational results

In order to test how ma-QAOA performs, we simulated the algorithm on a
collection of one-hundred triangle-free 3-regular graphs with fifty vertices
and one-hundred triangle-free 3-regular graphs with 100 vertices and compared
the approximation ratios calculated with ma-QAOA to those of 1-QAOA. We also
performed the same calculations with fifty modified $G_{n,p}$ random graphs
with fifty and one-hundred vertices each. In the $G_{n,p}$ model, $n$ sets the
number of nodes, and $p$ is the probability that an edge exists. In
particular, we examined $G_{50,0.08}$ and $G_{100,0.035}$ in order to create
random graphs that have average degree approximately three. After randomly
generating the graphs, triangles were removed by randomly removing edges from
each triangle. For these sets of triangle-free graphs we can compute $\langle
C\rangle_{1}^{\mathrm{ma}}$ for large $n$ using the analytical result of Thm.
IV.1. Table 1 shows the average approximation ratios for each collection of
graphs with ma-QAOA and 1-QAOA, as well as the changes in the approximation
ratio and percent change in the approximation ratio gap. This approximation
ratio gap is the percent difference between one minus the approximation ratio
for 1-QAOA and one minus the approximation ratio for ma-QAOA. The ma-QAOA has
a higher average approximation ratio and gives a significant percent increase
in approximation ratio gap for each collection of graphs. These simulations
only compare 1-QAOA to ma-QAOA, however, the next set of computational results
compares ma-QAOA to $p$-QAOA for $p\leq 3$ on all connected, non-isomorphic
graphs.

In previous work, we determined $\langle C\rangle_{1}$, $\langle
C\rangle_{2}$, and $\langle C\rangle_{3}$ for all connected, non-isomorphic
eight vertex graphs and compiled them into an online data set lotshaw2021bfgs
; lotshawdataset . For this work, we calculated the angles that maximize
$\langle C\rangle_{1}^{\mathrm{ma}}$ and compared $\langle C\rangle_{p}$ to
$\langle C\rangle^{\mathrm{ma}}_{1}$. On average, the performance of ma-QAOA
is comparable to 3-QAOA on these graphs. Table 2 shows that ma-QAOA has a
higher average approximation ratio than 1-QAOA and 2-QAOA on all eight vertex
graphs. However, the average approximation ratio for one iteration of ma-QAOA
is larger than the average approximation ratio for 3-QAOA.

Figure 2 shows how the distribution of approximation ratios for ma-QAOA
compares to the approximation ratios for up to three iterations of QAOA for
all connected, non-isomorphic eight vertex graphs. The percentage of graphs
with approximation ratio at least 0.95 is significantly higher with ma-QAOA
than up to three levels of QAOA. The fraction of graphs with approximation
ratio at least 0.85 and 0.9 is higher for 3-QAOA than ma-QAOA, however
significantly more graphs have an approximation ratio of at least 0.95 with
ma-QAOA.

Figure 2:  Fractions of non-isomorphic, eight vertex graphs with approximation ratios (A.R.) at least $x$ for ma-QAOA and $p$-QAOA. The lines are included in order to outline the shape of each distribution. Graph type | Average A.R. for 1-QAOA | Average A.R. for ma-QAOA | Change in A.R. | Percent change in Gap (1-AR)
---|---|---|---|---
50 vertex 3-regular | 0.7617 | 0.8123 | 0.0506 | 21.26%
100 vertex 3-regular | 0.7562 | 0.8000 | 0.0438 | 17.98%
Modified $G_{50,.08}$ | 0.7554 | 0.8156 | 0.0602 | 24.65&
Modified $G_{100,.035}$ | 0.7497 | 0.8098 | 0.0602 | 24.04%
Table 1: The average approximation ratio (A.R.) for a collection of one-hundred 3-regular graphs with fifty vertices, one-hundred 3-regular graphs with 100 vertices, fifty modified $G_{50,0.08}$ random graphs, and fifty modified $G_{100,0.035}$ random graphs.  QAOA type | Average approximation ratio for all eight vertex graphs
---|---
ma-QAOA | .9257
1-QAOA | .8061
2-QAOA | .8767
3-QAOA | .9192
Table 2: The average approximation ratio for all connected, non-isomorphic
graphs eight vertices.

### Circuit Fidelity

On fully connected hardware, the number of one-qubit unitary operators and
two-qubit unitary operators per iteration of QAOA for MaxCut corresponds to
the number of vertices and edges in the graph respectively. In order to solve
MaxCut on all eight vertex graphs, eight one-qubit unitary operators and
between seven and twenty-eight two-qubit unitary operators are needed per
iteration of the algorithm, depending on the number of edges in the graph.
Assuming a simple noise model based on Kraus-operator error channels acting
after each unitary, the circuit produces a final state
$\rho=F\rho_{\mathrm{ideal}}+(1-F)\rho_{\mathrm{noise}}$, where $F$ is the
probability associated with the ideal noiseless evolution component
$\rho_{\mathrm{ideal}}$ lotshaw2021resource . The expected number of
measurements to sample from the ideal distribution is $1/F$ in the worst-case,
when $\mathrm{Tr}\rho_{\mathrm{ideal}}\rho_{\mathrm{noise}}=0$; more generally
the number of measurements can decrease depending on the specific state and
noise process, but to keep the discussion general we take the expected number
of measurements as $1/F$. With this model, for $p$-QAOA or $p$-ma-QAOA,
$F=(1-\epsilon_{n})^{np}(1-\epsilon_{m})^{mp}$, where $\epsilon_{n}$ is the
expected error per single-qubit unitary, $\epsilon_{m}$ is the expected error
per two-qubit unitary, $n$ is the number of qubits, and $m$ is the number of
edges in the MaxCut graph. Assuming $p=1$ and an error rate of 1% for each
unitary operator, the expected number of measurements to obtain a sample from
the noiseless distribution is 1.25 when $n=8$ and $m=14.4$. We choose $m=14.4$
since this is the average number of edges per graph in the collection of
connected, non-isomorphic eight vertex graphs, though note each specific graph
has an integer number of edges.

We find that parameter optimization with ma-QAOA yields angles of zero for a
subset of the edge and vertex unitary operators. Since
$\exp(-i\gamma_{p,a}C_{a})=\mathbb{I}=\exp(-i\beta_{p,v}B_{v})$ when
$\gamma_{p,a}=0$ and $\beta_{p,v}=0$, all unitary operators with an angle of
zero may be excluded from the optimized circuit. The presence of fewer
operators reduces the amount of noise in the circuit by decreasing the
exponent of the first and second terms in $F$ by the number of vertex and edge
operators that have zero angles, respectively. When maximizing the expected
value of the cost function using ma-QAOA, a sizeable percent of angles
associated with vertices and edges were zero. Table 3 gives the percent of
zero angles, rounded to three decimal places, for each collection of graphs
that were studied. These zero angles impact the fidelity of each solution, as
the operators associated with them need not be implemented on hardware. Table
4 shows the ratio of the expected number of measurements needed to sample from
the noiseless distribution for $p$-QAOA relative to ma-QAOA for each
collection of graphs with varying values of $\epsilon_{m}$, using the average
reduction in gates for ma-QAOA from Table 3. Note that if the
$\epsilon_{m}=0.05$, the number of samples increases rapidly with $p$. Since
one iteration of ma-QAOA is comparable to three iterations of QAOA on eight
vertex graphs, if the trend holds for larger graphs, ma-QAOA has the potential
to require significantly fewer samples than QAOA.

$n$ |  Percent of $v$ with $\beta_{v}=0$ |  Percent of $a$ with $\gamma_{a}=0$
---|---|---
8 | 15.030 | 25.449
50 (3-reg.) | 13.000 | 18.147
50 (E.R.) | 11.440 | 14.381
100 (3-reg.) | 14.690 | 19.973
100 (E.R.) | 12.900 | 16.541
Table 3: The percent of $\beta_{v}$ and $\gamma_{a}$, rounded to three decimal places, that are zero when optimizing ma-QAOA on the family of graphs found in the first column.  $n$ | $m$ | $\epsilon_{n}=\epsilon_{m}=0.01$ | $\epsilon_{n}=0.01,\epsilon_{m}=0.05$
---|---|---|---
|  | $p=1$ | $p=2$ | $p=3$ | $p=1$ | $p=2$ | $p=3$
8 | 14.4 | 1.05 | 1.32 | 1.65 | 1.22 | 2.77 | 6.28
50 (3-reg.) | 75 | 1.22 | 4.30 | 15.10 | 2.15 | 166.16 | $1\times 10^{4}$
50 (E.R.) | 87.2 | 1.20 | 4.77 | 18.94 | 2.02 | 291.78 | $4\times 10^{4}$
100 (3-reg.) | 150 | 1.57 | 19.32 | 238.39 | 5.39 | $3\times 10^{4}$ | $2\times 10^{8}$
100 (E.R.) | 167.34 | 1.50 | 22.08 | 324.26 | 4.71 | $7\times 10^{4}$ | $1\times 10^{9}$
Table 4: The ratio of the expected number of measurements to obtain a sample
from the noiseless distribution for $p$-QAOA relative to $1$-ma-QAOA on an $n$
vertex graph with $m$ edges.

### Computing Angles

With a larger number of variables to optimize, the ma-QAOA method requires
more classical effort to find angles that optimize the approximation ratio.
However, it is not necessary to identify exact optimal angles, only to find
angles that are better than QAOA angles.

We used the BFGS algorithm to compute angles for the 8-vertex graphs; details
can be found in Methods. Figure 3 shows how the approximation ratio improves
on average across all iterations of BFGS for each ansatz studied for a random
sample of eight vertex graphs. Note that after approximately ten iterations,
ma-QAOA tends to achieve a higher approximation ratio than any of the
$p$-QAOA. We do note that the time required to perform each iteration of BFGS
is slower for ma-QAOA, as the number of gradient components is linearly
dependent on the number of variables being optimized.

Figure 3: Typical behavior of the BFGS search algorithm for ma-QAOA at $p=1$
and regular QAOA at $p=1,2,3$ on a set of 100 random eight vertex graphs. Each
curve is computed as an average over the graphs and over 100 random angle
seeds for each graph. Each curve terminates at the average final iteration of
the BFGS algorithm for the dataset.

### Scaling

We assess the scalability of ma-QAOA using computed optimized $\langle
C\rangle$ for sets of triangle-free Erdős-Rényi and 3-regular graphs of sizes
$n=50$ and $n=100$. The computational details are given in the Methods. We
compare the run times for typical graph optimizations to assess how the ma-
QAOA parameter optimization time increases with graph size.

For the Erdős-Rényi graphs, the time for a single optimization for $n=50$ was
0.10 seconds, for $n=100$ it was 0.46 seconds. We attribute the difference
primarily to the scaling in the calculation of the gradient, which is the most
expensive calculation in the optimization. Our approach computes each of the
$n+m$ derivatives $\partial\langle
C_{p,uv}\rangle^{\mathrm{ma}}/\partial\beta_{p,w}$ and $\partial\langle
C_{p,uv}\rangle^{\mathrm{ma}}/\partial\gamma_{p,jk}$ for each of the $m$ terms
$\langle C_{p,uv}\rangle^{\mathrm{ma}}$ in the cost function, giving a total
number of terms $\sim(n+m)m$. The time to compute each term will vary with the
degree of the graph, as this determines the number of cosine terms in Thm.
IV.1; however, for our graphs the degree is approximately constant hence does
not affect the scaling. For our graphs $m\sim n$ on average, so the overall
scaling is $\sim n^{2}$, which is consistent with the $\approx 4\times$
increase in time when $n$ is doubled from $n=50$ to $n=100$. We attribute the
remainder of the time difference to variations in the number of iterations as
$n$ and $m$ increase.

##  III Discussion

We have shown that multi-angle QAOA converges to an optimal solution, and
furthermore that $\langle C\rangle_{1}^{\mathrm{ma}}\geq\langle C\rangle_{1}$,
as QAOA is a special case of ma-QAOA. We also show that the circuit depth for
one iteration of ma-QAOA can be less than the depth for one iteration of QAOA.
Additionally, the analysis of star graphs shows that there is a family of
graphs that gives larger $\langle C\rangle$ for MaxCut when solved with ma-
QAOA than when solved with QAOA. We find significant increases in the
approximation ratio in numerical optimizations for large triangle-free graphs
and over the set of all non-isomorphic graphs with eight vertices, hence fewer
layers are required to reach the same performance as QAOA. In the presence of
noise, the reduction in number of layers and in the number of unitary
operators per layer can significantly decrease the expected number of
measurements needed to sample from the noiseless distribution. This could be a
significant advantage for computations on noisy quantum hardware.

Interestingly, some graphs do not have a significantly higher $\langle
C\rangle$ when solving MaxCut with ma-QAOA versus QAOA. It would be useful to
characterize for which graphs the increase in $\langle C\rangle$ from QAOA to
ma-QAOA is insignificant. This would help determine the appropriate ansatz to
use in order to solve MaxCut on the graph.

One drawback to ma-QAOA is that the number of classically optimized parameters
is $n+m$ per layer, where $n$ is the number of vertices of $G$ and $m$ is the
number of edges. From a practical standpoint, one way to solve optimal ma-QAOA
angles would be to calculate $\beta$ and $\gamma$ that optimize QAOA. We can
use those angles as the initial point of a BFGS search for the optimal
$\beta_{p,v}$ and $\gamma_{p,a_{i}}$ for all vertices $v$ and edges $a_{i}$.
Overall, however, the results seem to indicate that good paramaters can be
found in polynomial time. As many combinatorial optimization problems, like
MaxCut, are NP-Hard, any polynomially-bounded effort that improves performance
is likely to improve performance at large scale.

##  IV Methods

### Proof of convergence

###### Proof.

Recall that QAOA converges to the optimal solution for a combinatorial
optimization problem, which is the maximum over the objective function
farhi2014quantum . Thus, in order to show convergence of ma-QAOA, we need only
bound ma-QAOA from below by the value of QAOA. However, it is clear that the
optimal expected value of the cost function for ma-QAOA can be no lower than
that of QAOA, since QAOA is a special case of ma-QAOA when all
$\gamma_{p,ij}=\gamma_{p,kl}$ and all $\beta_{p,a}=\beta_{p,b}$ for all edges
$ij,kl$ and all vertices $a,b$. ∎

###  $\langle C\rangle$ for solving MaxCut on triangle-free graphs using ma-
QAOA

In order to prove that $\langle C\rangle_{1}^{\mathrm{ma}}=1$ for MaxCut on
star graphs, we derive a formula that calculates $\langle
C\rangle_{1}^{\mathrm{ma}}$ for MaxCut on triangle-free graphs.

######  Theorem IV.1.

Let $\beta_{p,u}^{\prime}=2\beta_{p,u}$ and
$\beta_{p,v}^{\prime}=2\beta_{p,v}$ The expected value of $C$ after one
iteration of ma-QAOA applied to MaxCut for triangle free graphs $G$ is

$\displaystyle\begin{split}&{\left\langle\vec{\gamma_{1}}\vec{\beta_{1}}\right|}C_{uv}{\left|\vec{\gamma_{1}}\vec{\beta_{1}}\right\rangle}=\frac{1}{2}+\\\
&\frac{1}{2}\sin{\gamma_{1,uv}}(\cos{\beta_{1,v}^{\prime}}\sin{\beta_{1,u}^{\prime}}\prod_{w}\cos{\gamma_{1,uw}}+\\\
&\cos{\beta_{1,u}^{\prime}}\sin{\beta_{1,v}^{\prime}}\prod_{x}\cos{\gamma_{1,vx}})\end{split}$

where $w\in Nbhd(u)\setminus v$ and $x\in Nbhd(v)\setminus u$.

The neighborhood of a vertex $x$, denoted $Nbhd(x)$, is the set of vertices
$y$ such that $xy\in E(G)$.

###### Proof.

The proof of this result relies on the Pauli-solver algorithm, which is
explained in detail in hadfield2018quantum . The proof of the result is
virtually identical to that for QAOA on triangle-free graphs, but we include
the proof here for completeness.

Consider edge $uv$ and consider acting on
$C_{uv}=(1/2)(\mathbb{I}-Z_{u}Z_{v})$ by conjugation of the mixing operator,
$\prod_{i\in V}e^{-i\beta_{1,i}B_{i}},$ followed by conjugation of the phase
operator, $\prod_{uv\in E}e^{-i\gamma_{1,uv}C_{uv}}$. We have that

$\begin{split}&\prod_{i\in V}e^{i\beta_{1,i}B_{i}}Z_{u}Z_{v}\prod_{i\in
V}e^{-i\beta_{1,i}B_{i}}=e^{2i\beta_{1,u}X_{u}}e^{2i\beta_{1,v}X_{v}}Z_{u}Z_{v}=\\\
&\cos{2\beta_{1,u}}\cos{2\beta_{1,v}}Z_{u}Z_{v}+\cos{2\beta_{1,v}}\sin{2\beta_{1,u}}Y_{v}Z_{u}\\\
&+\cos{2\beta_{1,u}}\sin{2\beta_{1,v}}Z_{v}Y_{u}+\sin{2\beta_{1,u}}\sin{2\beta_{1,v}}Y_{u}Y_{v}.\end{split}$
(1)

Note that the first term commutes with $\prod_{uv\in
E}e^{-i\gamma_{1,uv}C_{uv}}$, so does not contribute to the expected value.
Let $V_{u}$ be the neighborhood of $u$ in $V(G)$. Conjugating the third term
of Eqn. (1) by $\prod_{uv\in E}e^{-i\gamma_{1,uv}C_{uv}}$, we get

$\displaystyle{\left\langle
s\right|}\Upsilon^{\dagger}Y_{u}Z_{v}\Upsilon{\left|s\right\rangle}=$
$\displaystyle{\left\langle s\right|}e^{2i\gamma_{1,uv}C_{uv}}e^{2i\sum_{a\in
V_{u}\setminus{v}}\gamma_{1,ua}C_{ua}}Y_{u}Z_{v}{\left|s\right\rangle}=$
$\displaystyle{\left\langle
s\right|}e^{-i\gamma_{1,uv}Z_{u}Z_{v}}e^{-i\sum_{a\in
V_{u}\setminus{v}}\gamma_{1,ua}Z_{u}Z_{a}}Y_{u}Z_{v}{\left|s\right\rangle}=$
$\displaystyle{\left\langle
s\right|}(\mathbb{I}\cos{\gamma_{1,uv}}-i\sin{\gamma_{1,uv}})Z_{u}Z_{v}$
$\displaystyle\prod_{a\in
V_{u}}(\mathbb{I}\cos{\gamma_{1,ua}}-i\sin{\gamma_{1,ua}}Z_{u}Z_{a})Y_{u}Z_{v}{\left|s\right\rangle}=$
$\displaystyle-\sin{\gamma_{1,uv}}\prod_{a\in V_{u}}\cos{\gamma_{1,ua}},$

where $\Upsilon=e^{-i\gamma_{1,uv}C_{uv}}e^{-i\sum_{a\in
V_{u}\setminus{v}}\gamma_{1,ua}C_{ua}}$, and $\Upsilon^{\dagger}$ is its
Hermitian conjugate. By symmetry, the term for $Z_{u}Y_{v}$ is
$-\sin{\gamma_{1,uv}}\prod_{b\in V_{v}\setminus{u}}\cos{\gamma_{1,vb}}$, where
$V_{v}$ is the neighborhood of $v$ in $V$. Factoring in the coefficient $-1/2$
of $Z_{u}Z_{v}$ in $C_{uv}$ gives the final two terms in the theorem.

Now, let us conjugate the last term of Eqn. (1). Doing so, we get

$\displaystyle{\left\langle s\right|}e^{i\sum_{gh\in
E}\gamma_{1,gh}C_{gh}}Y_{u}Y_{v}e^{-i\sum_{gh\in
E}\gamma_{1,gh}C_{gh}}{\left|s\right\rangle}=$ $\displaystyle{\left\langle
s\right|}\prod_{a\in V_{u}\setminus
v}(\mathbb{I}\cos{\gamma_{1,ua}}-i\sin{\gamma_{1,ua}}Z_{u}Z_{a})\times$
$\displaystyle\prod_{b\in V_{v}\setminus
u}(\mathbb{I}\cos{\gamma_{1,vb}}-i\sin{\gamma_{1,vb}}Z_{v}Z_{b})Y_{u}Y_{v}{\left|s\right\rangle}$

The simplest terms that contribute to the expected value are of the form

$\displaystyle\sin{\gamma_{1,uc}}\sin{\gamma_{1,vc}}\prod_{x\neq
y}\cos{\gamma_{1,ux}}\cos{\gamma_{1,vy}}$

and there are $f$ of these where $f$ is the number of triangles containing
$uv$. The higher order terms only contribute to the expected value if there
are triangles in the graph. Thus, the last term of Eqn. (1) contributes
nothing to the expected value of triangle-free graphs.

Combining these expressions gives the theorem. ∎

### Star graphs

First, we will show that $\langle C_{ij}\rangle$ approaches $0.75$ as $n$
tends to infinity for QAOA. Since there are $n-1$ edges in a star on $n$
vertices, this implies $\langle C\rangle$ tends to $0.75(n-1)$. Additionally,
$n-1$ is the size of the optimal MaxCut solution, so $\langle
C\rangle_{1}/C_{\mathrm{max}}=0.75$.

Wang, Hadfield, Jiang, and Rieffel showed that wang2018quantum

$\begin{split}&\langle
C_{ij}\rangle_{1}=\frac{1}{2}+\frac{1}{4}(\sin{4\beta}\sin{\gamma})(\cos^{d}{\gamma}+\cos^{e}{\gamma})\\\
&-\frac{1}{4}(\sin^{2}{2\beta}\cos^{d+e-f}{\gamma})(1-\cos^{f}{2\gamma})\end{split}$
(2)

where $d$ is the $\deg(i)-1$, $e$ is the $\deg(j)-1$ and $f$ is the number of
triangles containing edge $ij$ wang2018quantum ; hadfield2018quantum .

Let us consider the above formula applied to a star graph. Without loss of
generality, let $j$ be the center of each star. Then $d=0$, $e=n-2$, and
$f=0$, since star graphs are trees. For each edge of the star, Eq. (2) reduces
to

$\langle
C_{ij}\rangle_{1}=\frac{1}{2}+\frac{1}{4}(\sin{4\beta}\sin{\gamma})(1+\cos^{n-2}{\gamma}).$

We set $\beta=\pi/8$, which implies $\sin{4\beta}=1$, since only one
trigonometric function has $\beta$ as an argument. As $n$ tends to infinity,
note $\cos^{n-2}{\gamma}$ tends to zero unless $\gamma=k\pi$ for some
$k\in\mathbf{N}$. However, if $\gamma=k\pi$, $\sin{\gamma}=0$. Thus, this
quantity is maximized when $\gamma\neq k\pi$, which implies $\langle
C_{ij}\rangle_{1}$ tends to $0.75$ for star graphs.

In order to prove $\langle C\rangle^{\mathrm{ma}}=n-1$ for ma-QAOA on star
graphs, we examine Thm. IV.1. Without loss of generality, let $u$ be a leaf
vertex and $v$ be the center. Note that the first product is empty, since the
leaf vertices have no neighbors except the center. Thus, Thm. IV.1 reduces to

$\displaystyle{\left\langle\vec{\gamma_{1}}\vec{\beta_{1}}\right|}C_{uv}{\left|\vec{\gamma_{1}}\vec{\beta_{1}}\right\rangle}=\frac{1}{2}+$
$\displaystyle\frac{1}{2}\sin{\gamma_{1,uv}}(\cos{\beta_{1,v}^{\prime}}\sin{\beta_{1,u}^{\prime}}+\cos{\beta_{1,u}^{\prime}}\sin{\beta_{1,v}^{\prime}}\prod_{x}\cos{\gamma_{1,vx}})$

Now, recall
${\left\langle\vec{\gamma_{1}}\vec{\beta_{1}}\right|}C_{uv}{\left|\vec{\gamma_{1}}\vec{\beta_{1}}\right\rangle}\leq
1$, as two vertices that have an edge between them add one to the objective
function if they are in different sets. In order to obtain equality, we can
set $\gamma_{1,uv}=\pi/2$, as it is an argument for only a single sine term.
Next, note that if either term in the parenthesis is one, the other must be
zero. Also, setting one term equal to one allows gives an expected value of
one for the edge. Let $\beta_{1,u}^{\prime}=\pi/2$ and
$\beta_{1,v}^{\prime}=0$. Then
$\cos{\beta_{1,v}^{\prime}}=\sin{\beta_{1,u}^{\prime}}=1$ while
$\cos{\beta_{1,u}^{\prime}}=\sin{\beta_{1,v}^{\prime}}=0$. Thus, the first
term in the parenthesis is one and the second is zero. This allows us to set
$\gamma_{1,vx}=\pi/2$ for all $x\in Nbhd(v)$. Since each of the $n-1$ edges in
the star are described similarly, $\langle C\rangle_{1}^{\mathrm{ma}}=n-1$ for
all $n$. The size of the optimal cut on a star graph is $n-1$, so $\langle
C\rangle_{1}^{\mathrm{ma}}/C_{\mathrm{max}}=1$.

### Setup for Computational Results

In order to calculate the angles that maximize $\langle C\rangle_{p}$ and
$\langle C\rangle_{1}^{\mathrm{ma}}$ for the eight vertex graphs, we used the
Broyden-Fletcher-Goldfarb-Shanno (BFGS) algorithm NumericalRecipesBFGS . The
algorithm inputs an initial collection of angles and then uses a numerical
gradient and second order approximate Hessian to find angles that converge to
local maxima of $\langle C\rangle_{p}$ and $\langle
C\rangle_{1}^{\mathrm{ma}}$. For the eight vertex graphs, 100 random seeds
were used to optimize $\langle C\rangle^{\mathrm{ma}}_{1}$. The results for
the $\langle C\rangle_{p}$ were taken from the online dataset lotshawdataset
of Ref. lotshaw2021bfgs , where we performed an exhaustive analysis of QAOA
performance on small graphs. These used 50 seeds for $p=1$, 100 seeds for
$p=2$, and 1,000 seeds for $p=3$.

For the 50 and 100 node graphs, we used the method of moving asymptotes (MMA)
algorithm nlopt ; nloptmma , but note that calculations with BFGS gave similar
results. The $\langle C\rangle_{1}$ were computed using Eq. (2) and the
reported results were taken as the best from 1000 initial seeds in MMA
optimizations. The $\langle C\rangle_{1}^{\mathrm{ma}}$ were computed with
Theorem IV.1 and MMA optimization. We report results as the best optimized
values from 1,000 seeds at $n=50$ and from 100 seeds at $n=100$.

## Acknowledgements

This work was supported by DARPA ONISQ program under award W911NF-20-2-0051.
J. Ostrowski acknowledges the Air Force Office of Scientific Research award,
AF-FA9550-19-1-0147. G. Siopsis acknowledges the Army Research Office award
W911NF-19-1-0397. J. Ostrowski and G. Siopsis acknowledge the National Science
Foundation award OMA-1937008.

## References

  * [1] Diogo Cruz, Romain Fournier, Fabien Gremion, Alix Jeannerot, Kenichi Komagata, Tara Tosic, Jarla Thiesbrummel, Chun Lam Chan, Nicolas Macris, Marc-André Dupertuis, et al.  Efficient quantum algorithms for ghz and w states, and implementation on the ibm quantum computer.  Advanced Quantum Technologies, 2(5-6):1900015, 2019.
  * [2] Jingfu Zhang, Swathi S Hegde, and Dieter Suter.  Efficient implementation of a quantum algorithm in a single nitrogen-vacancy center of diamond.  Physical Review Letters, 125(3):030501, 2020.
  * [3] Clément Godfrin, Abdelkarim Ferhat, Rafik Ballou, Svetlana Klyatskaya, Mario Ruben, Wolfgang Wernsdorfer, and Franck Balestro.  Operating quantum states in single magnetic molecules: implementation of grover’s quantum algorithm.  Physical review letters, 119(18):187702, 2017.
  * [4] Ajinkya Borle, Vincent E Elfving, and Samuel J Lomonaco.  Quantum approximate optimization for hard problems in linear algebra.  arXiv preprint arXiv:2006.15438, 2020.
  * [5] Amir H Karamlou, William A Simon, Amara Katabarwa, Travis L Scholten, Borja Peropadre, and Yudong Cao.  Analyzing the performance of variational quantum factoring on a superconducting quantum processor.  arXiv preprint arXiv:2012.07825, 2020.
  * [6] Cornelius Hempel, Christine Maier, Jonathan Romero, Jarrod McClean, Thomas Monz, Heng Shen, Petar Jurcevic, Ben P Lanyon, Peter Love, Ryan Babbush, et al.  Quantum chemistry calculations on a trapped-ion quantum simulator.  Physical Review X, 8(3):031022, 2018.
  * [7] C Huerta Alderete, Shivani Singh, Nhung H Nguyen, Daiwei Zhu, Radhakrishnan Balu, Christopher Monroe, CM Chandrashekar, and Norbert M Linke.  Quantum walks and dirac cellular automata on a programmable trapped-ion quantum computer.  Nature communications, 11(1):1–7, 2020.
  * [8] Peng Xue, Barry C Sanders, and Dietrich Leibfried.  Quantum walk on a line for a trapped ion.  Physical review letters, 103(18):183602, 2009.
  * [9] Norbert M Linke, Dmitri Maslov, Martin Roetteler, Shantanu Debnath, Caroline Figgatt, Kevin A Landsman, Kenneth Wright, and Christopher Monroe.  Experimental comparison of two quantum computing architectures.  Proceedings of the National Academy of Sciences, 114(13):3305–3310, 2017.
  * [10] Guido Pagano, Aniruddha Bapat, Patrick Becker, Katherine S Collins, Arinjoy De, Paul W Hess, Harvey B Kaplan, Antonis Kyprianidis, Wen Lin Tan, Christopher Baldwin, et al.  Quantum approximate optimization of the long-range ising model with a trapped-ion quantum simulator.  Proceedings of the National Academy of Sciences, 117(41):25396–25401, 2020.
  * [11] Andreas Bengtsson, Pontus Vikstål, Christopher Warren, Marika Svensson, Xiu Gu, Anton Frisk Kockum, Philip Krantz, Christian Križan, Daryoush Shiri, I Svensson, et al.  Quantum approximate optimization of the exact-cover problem on a superconducting quantum processor.  arXiv preprint arXiv:1912.10495, 2019.
  * [12] Matthew P Harrigan, Kevin J Sung, Matthew Neeley, Kevin J Satzinger, Frank Arute, Kunal Arya, Juan Atalaya, Joseph C Bardin, Rami Barends, Sergio Boixo, et al.  Quantum approximate optimization of non-planar graph problems on a planar superconducting processor.  Nature Physics, 17(3):332–336, 2021.
  * [13] Edward Farhi, Jeffrey Goldstone, and Sam Gutmann.  A quantum approximate optimization algorithm.  arXiv preprint arXiv:1411.4028, 2014.
  * [14] Zain Hamid Saleem.  Max-independent set and the quantum alternating operator ansatz.  International Journal of Quantum Information, 18(04):2050011, 2020\.
  * [15] Jeremy Cook, Stephan Eidenbenz, and Andreas Bärtschi.  The quantum alternating operator ansatz on max-k vertex cover.  Bulletin of the American Physical Society, 65, 2020.
  * [16] Gavin E Crooks.  Performance of the quantum approximate optimization algorithm on the maximum cut problem.  arXiv preprint arXiv:1811.08419, 2018.
  * [17] Hannes Pichler, Sheng-Tao Wang, Leo Zhou, Soonwon Choi, and Mikhail D Lukin.  Quantum optimization for maximum independent set using rydberg atom arrays.  arXiv preprint arXiv:1808.10816, 2018.
  * [18] Edward Farhi, David Gamarnik, and Sam Gutmann.  The quantum approximate optimization algorithm needs to see the whole graph: Worst case examples.  arXiv preprint arXiv:2005.08747, 2020.
  * [19] Edward Farhi, David Gamarnik, and Sam Gutmann.  The quantum approximate optimization algorithm needs to see the whole graph: A typical case.  arXiv preprint arXiv:2004.09002, 2020.
  * [20] Jonathan Wurtz and Peter J Love.  Bounds on maxcut qaoa performance for p¿ 1.  arXiv e-prints, pages arXiv–2010, 2020.
  * [21] Jonathan Ward, Johannes Otterbach, Gavin Crooks, Nicholas Rubin, and Marcus da Silva.  Qaoa performance benchmarks using max-cut.  In APS March Meeting Abstracts, volume 2018, pages R15–007, 2018\.
  * [22] Ruslan Shaydulin and Yuri Alexeev.  Evaluating quantum approximate optimization algorithm: A case study.  In 2019 tenth international green and sustainable computing conference (IGSC), pages 1–6. IEEE, 2019.
  * [23] Andreas Bärtschi and Stephan Eidenbenz.  Grover mixers for qaoa: Shifting complexity from mixer design to state preparation.  In 2020 IEEE International Conference on Quantum Computing and Engineering (QCE), pages 72–82. IEEE, 2020.
  * [24] Stuart Hadfield, Zhihui Wang, Bryan O’Gorman, Eleanor G Rieffel, Davide Venturelli, and Rupak Biswas.  From the quantum approximate optimization algorithm to a quantum alternating operator ansatz.  Algorithms, 12(2):34, 2019.
  * [25] Jonathan Wurtz and Peter Love.  Classically optimal variational quantum algorithms.  arXiv preprint arXiv:2103.17065, 2021.
  * [26] Zhihui Wang, Nicholas C Rubin, Jason M Dominy, and Eleanor G Rieffel.  $xy$-mixers: analytical and numerical results for qaoa.  arXiv preprint arXiv:1904.09314, 2019.
  * [27] Daniel J Egger, Jakub Mareček, and Stefan Woerner.  Warm-starting quantum optimization.  Quantum, 5:479, 2021.
  * [28] Linghua Zhu, Ho Lun Tang, George Barron, Fernando Calderon-Vargas, Nicholas Mayhall, Edwin Barnes, and Sophia Economou.  Improving the performance of quantum approximate optimization algorithm through an adaptive, problem-tailored ansatz.  Bulletin of the American Physical Society, 2021.
  * [29] Reuben Tate, Majid Farhadi, Creston Herold, Greg Mohler, and Swati Gupta.  Bridging classical and quantum with sdp initialized warm-starts for qaoa.  arXiv preprint arXiv:2010.14021, 2020.
  * [30] Eleanor Rieffel, Jason M Dominy, Nicholas Rubin, and Zhihui Wang.  Xy-mixers: analytical and numerical results for qaoa.  2020\.
  * [31] John Golden, Andreas Bärtschi, Daniel O’Malley, and Stephan Eidenbenz.  Threshold-based quantum optimization.  arXiv preprint arXiv:2106.13860, 2021.
  * [32] Rebekah Herrman, Lorna Treffert, James Ostrowski, Phillip C. Lotshaw, Travis S. Humble, and George Siopsis.  Impact of graph structures for qaoa on maxcut.  arXiv preprint arXiv:2102.05997, 2021.
  * [33] Ruslan Shaydulin, Stuart Hadfield, Tad Hogg, and Ilya Safro.  Classical symmetries and qaoa.  arXiv preprint arXiv:2012.04713, 2020.
  * [34] Yuri Alexeev, Sami Khairy, Ruslan Shaydulin, Lukasz Cincio, and Prasanna Balaprakash.  Reinforcement learning for finding qaoa parameters.  Bulletin of the American Physical Society, 65, 2020.
  * [35] Matteo M Wauters, Emanuele Panizon, Glen B Mbeng, and Giuseppe E Santoro.  Reinforcement-learning-assisted quantum optimization.  Physical Review Research, 2(3):033446, 2020.
  * [36] Phillip C. Lotshaw, Travis S. Humble, Rebekah Herrman, James Ostrowski, and George Siopsis.  Empirical performance bounds for quantum approximate optimization.  arXiv preprint arXiv:2102.06813, 2021.
  * [37] Gian Giacomo Guerreschi and Anne Y Matsuura.  Qaoa for max-cut requires hundreds of qubits for quantum speed-up.  Scientific reports, 9, 2019.
  * [38] Rebekah Herrman, James Ostrowski, Travis S Humble, and George Siopsis.  Lower bounds on circuit depth of the quantum approximate optimization algorithm.  Quantum Information Processing, 20(2):1–17, 2021.
  * [39] Cheng Xue, Zhao-Yun Chen, Yu-Chun Wu, and Guo-Ping Guo.  Effects of quantum noise on quantum approximate optimization algorithm.  Chinese Physics Letters, 38(3):030302, 2021.
  * [40] Samson Wang, Enrico Fontana, M Cerezo, Kunal Sharma, Akira Sone, Lukasz Cincio, and Patrick J Coles.  Noise-induced barren plateaus in variational quantum algorithms.  arXiv preprint arXiv:2007.14384, 2020.
  * [41] Prakash Murali, Jonathan M Baker, Ali Javadi-Abhari, Frederic T Chong, and Margaret Martonosi.  Noise-adaptive compiler mappings for noisy intermediate-scale quantum computers.  In Proceedings of the Twenty-Fourth International Conference on Architectural Support for Programming Languages and Operating Systems, pages 1015–1029, 2019.
  * [42] Jinzhao Sun, Xiao Yuan, Takahiro Tsunoda, Vlatko Vedral, Simon C Benjamin, and Suguru Endo.  Mitigating realistic noise in practical noisy intermediate-scale quantum devices.  Physical Review Applied, 15(3):034026, 2021.
  * [43] Jeffrey Marshall, Filip Wudarski, Stuart Hadfield, and Tad Hogg.  Characterizing local noise in qaoa circuits.  IOP SciNotes, 1(2):025208, 2020.
  * [44] Mahabubul Alam, Abdullah Ash-Saki, and Swaroop Ghosh.  Design-space exploration of quantum approximate optimization algorithm under noise.  In 2020 IEEE Custom Integrated Circuits Conference (CICC), pages 1–4. IEEE, 2020.
  * [45] Mahabubul Alam, Abdullah Ash-Saki, and Swaroop Ghosh.  Analysis of quantum approximate optimization algorithm under realistic noise in superconducting qubits.  arXiv preprint arXiv:1907.09631, 2019.
  * [46] Michael Streif, Martin Leib, Filip Wudarski, Eleanor Rieffel, and Zhihui Wang.  Quantum algorithms with local particle-number conservation: Noise effects and error correction.  Physical Review A, 103(4):042412, 2021.
  * [47] Filip B Maciejewski, Flavio Baccari, Zoltán Zimborás, and Michał Oszmaniec.  Modeling and mitigation of realistic readout noise with applications to the quantum approximate optimization algorithm.  arXiv preprint arXiv:2101.02331, 2021.
  * [48] Phillip C. Lotshaw and Travis S. Humble.  QAOA dataset.  Found at https://code.ornl.gov/qci/qaoa-dataset-version1.
  * [49] Phillip C. Lotshaw, Thien Nguyen, Anthony Santana, Alexander McCaskey, Rebekah Herrman, James Ostrowski, George Siopsis, and Travis S. Humble.  Scaling quantum approximate optimization on near-term hardware.  forthcoming, 2021.
  * [50] Stuart Hadfield.  Quantum algorithms for scientific computing and approximate optimization.  arXiv preprint arXiv:1805.03265, 2018.
  * [51] Zhihui Wang, Stuart Hadfield, Zhang Jiang, and Eleanor G Rieffel.  Quantum approximate optimization algorithm for maxcut: A fermionic view.  Physical Review A, 97(2):022304, 2018.
  * [52] William H. Press, Brian P. Flannery, and Saul A. Teukolsky.  Numerical Recipes in Fortran 77: The Art of Scientific Computing.  Cambridge University Press, second edition, 1993.  https://people.sc.fsu.edu/$\sim$inavon/5420a/DFP.pdf.
  * [53] Steven G. Johnson.  The nlopt nonlinear-optimization package.  http://github.com/stevengj/nlopt.
  * [54] Krister Svanberg.  A class of globally convergent optimization methods based on conservative convex separable approximations.  SIAM J. Optim., 12:555–573, 2002.