# Bright and dark breathers of the Benjamin–Ono equation on the traveling
periodic background

Jinbing Chen  School of Mathematics, Southeast University, Nanjing, Jiangsu
210096, PR China<EMAIL_ADDRESS>and  Dmitry E. Pelinovsky  Department of
Mathematics and Statistics, McMaster University, Hamilton, Ontario, Canada,
L8S 4K1<EMAIL_ADDRESS>
###### Abstract.

The Benjamin–Ono (BO) equation describes long internal waves of small
amplitude in deep fluids. Compared to its counterpart for shallow fluids, the
Korteweg-de Vries (KdV) equation, the BO equation admits exact solutions for
the traveling periodic and solitary waves as well as their interactions
expressed in elementary (trigonometric and polynomial) functions. Motivated by
a recent progress for the KdV equation, we discover here two scenarios of the
soliton-periodic wave interactions which result in the propagation of either
elevation (bright) or depression (dark) breathers (periodic in time coherent
structures). The existence of two different breathers is related to the band-
gap spectrum of the Lax operator associated with the traveling periodic wave.
Given a simple structure of the exact solutions in the BO equation, we obtain
a closed-form expression for multi-solitons interacting with the traveling
periodic wave.

This article is written in memory of N. F. Smyth for his many contributions to
studies of the modulation theory and, in particular, periodic waves of the
Benjamin–Ono equation.

##  1\. Introduction

Two limiting cases of long internal waves of small amplitude can be described
within the integrable evolution equations: the Korteweg–de Vries (KdV)
equation for shallow fluids and the Benjamin–Ono (BO) equation for deep fluids
[16]. These integrable equations have a large class of exact solutions
describing spatially periodic and solitary traveling waves as well as their
interactions [26].

Exact solutions for the time-periodic interactions of a single solitary wave
with the spatially periodic traveling wave were constructed recently for the
KdV equation in [5, 18] after much earlier work [15, 21]. Such periodic in
time coherent structures were referred to as breathers. Two families of
breathers were constructed as elevation (bright) breathers and depression
(dark) breathers. The two families correspond to two scenarios of a more
complicated interaction of a solitary wave with dispersive shock waves (DSWs),
see a recent review in [4]. Bright breathers describe a transmission of a
solitary wave overtaking DSWs and dark breathers describe a trapping of a
solitary wave in DSWs. The exact bright and dark breather solutions of the KdV
equation expressed in terms of the Jacobi theta functions were used for
comparison with numerical and laboratory experiments of soliton-DSW
interactions in shallow fluids [24, 31, 43].

Similar breathers of the defocusing nonlinear Schrödinger equation were
obtained in different solution forms in [22, 23, 41, 44]. Recently, breathers
of the defocusing modified KdV equation were elaborated in [33] by analyzing
properties of the Jacobi theta functions.

The purpose of this work is to obtain breather solutions for the BO equation
which covers applications of soliton-DSW interactions in deep fluids. See [11,
39] for two particular examples where the DSWs of the BO equation arise in
various physical applications.

We take the BO equation in the standard normalized form:

$u_{t}+2uu_{x}+H(u_{xx})=0,\quad H(f):=\frac{1}{\pi}{\rm
p.v.}\int_{-\infty}^{\infty}\frac{f(y)dy}{y-x}.$ (1)

The periodic traveling wave is expressed in terms of elementary functions:

$u(x,t)=\frac{k\sinh\phi}{\cos(k\xi)+\cosh\phi},$ (2)

where $\xi=x-c_{0}t-\xi_{0}$ is the traveling wave coordinate, $k,\phi>0$ are
arbitrary parameters, $\xi_{0}\in\mathbb{R}$ is an arbitrary phase, and
$c_{0}=k\coth\phi$. The solitary wave corresponds to the limit $\phi\to 0$
with $k=c_{0}\tanh\phi\to 0$ uniquely determined by fixed $c_{0}>0$. After
phase transformation $\xi\to\xi+\pi k^{-1}$, the periodic wave solution (2)
reduces to the solitary wave solution written in the form:

$u(x,t)=\frac{2c_{0}}{1+c_{0}^{2}\xi^{2}}.$ (3)

The main result of this work is that the BO equation admits two families of
breather solutions similar to the KdV equation. However, the breather
solutions of the BO equation are expressed in terms of elementary
(trigonometric and power) functions instead of the Jacobi theta functions. The
new exact solutions are given by

$\hat{u}(x,t)=\frac{2(c+k\beta)\cosh\phi+[k(1+\beta^{2}+c^{2}\eta^{2})+2\beta
c]\sinh\phi+2c\cos(k\xi)}{(1+\beta^{2}+c^{2}\eta^{2})\cosh\phi+2\beta\sinh\phi+(1-\beta^{2}+c^{2}\eta^{2})\cos(k\xi)+2\beta
c\eta\sin(k\xi)},$ (4)

where $\eta=x-ct-\eta_{0}$ is the traveling wave coordinate,
$\eta_{0}\in\mathbb{R}$ is an arbitrary phase, $\beta$ is uniquely defined by

$\beta=\frac{2ck}{(c-c_{0})^{2}-k^{2}},$ (5)

and $c>0$ is an arbitrary speed of the solitary wave defined in two disjoint
regions:

$(\mbox{\rm BB})\;\;c-c_{0}>k\quad\mbox{\rm and}\quad(\mbox{\rm
DB})\;\;c-c_{0}<-k.$ (6)

The bright breathers (BB) correspond to $c>c_{0}+k$ and the dark breathers
(DB) correspond to $0<c<c_{0}-k$, hence the former overtakes the periodic wave
and the latter is trapped in the periodic wave, similarly to the KdV equation
[18]. As $\eta\to\pm\infty$, $\hat{u}(x,t)\to u(x,t)$ given by (2) so that no
phase shift arises in the interactions between the solitary wave and the
periodic wave, contrary to the KdV equation [18]. The localization of the
solitary wave is defined by its speed $c$: the larger $c$ is, the smaller is
the width of the solitary wave.

Figure 1 show the solution surfaces computed from the exact solution (4) with
two choices of $c$ in (6). The left panel shows propagation of a bright
(elevation) breather with the speed exceeding the speed of the background
periodic wave. The right panel shows propagation of a dark (depression)
breather with the speed being smaller that the speed of the background
periodic wave.

Figure 1. Left: bright breather with $k=1$, $\phi=3$, and $c=c_{0}+1.5k$.
Right: dark breather with $k=1$, $\phi=0.5$, and $c=c_{0}-1.5k$.

By generalizing the breather solution (4) in the determinant form [29], we
present the general multi-breather solution on the traveling periodic wave (2)
as the following theorem.

######  Theorem 1.

The BO equation (1) admits the following general multi-breather solution on
the traveling periodic wave (2) in the form

$u(x,t)=i\partial_{x}\log\frac{\det(F^{\prime})}{\det(F)},$ (7)

with

$F=\left[\begin{matrix}\displaystyle
1+e^{ik\xi+\phi}&\displaystyle\frac{2k}{k+c_{0}-c_{1}}&\dots&\displaystyle\frac{2k}{k+c_{0}-c_{N}}\\\
\displaystyle\frac{2c_{1}}{k+c_{1}-c_{0}}&\displaystyle-
ic_{1}\eta_{1}-1-\frac{2kc_{1}}{(c_{0}-c_{1})^{2}-k^{2}}&\dots&\displaystyle\frac{2c_{1}}{c_{1}-c_{N}}\\\
\vdots&\vdots&\ddots&\vdots\\\
\displaystyle\frac{2c_{N}}{k+c_{N}-c_{0}}&\displaystyle\frac{2c_{N}}{c_{N}-c_{1}}&\dots&\displaystyle-
ic_{N}\eta_{N}-1-\frac{2kc_{N}}{(c_{0}-c_{N})^{2}-k^{2}}\end{matrix}\right]$

and

$F^{\prime}=\left[\begin{matrix}\displaystyle
1+e^{ik\xi-\phi}&\displaystyle\frac{2k}{k+c_{0}-c_{1}}&\dots&\displaystyle\frac{2k}{k+c_{0}-c_{N}}\\\
\displaystyle\frac{2c_{1}}{k+c_{1}-c_{0}}&\displaystyle-
ic_{1}\eta_{1}+1-\frac{2kc_{1}}{(c_{0}-c_{1})^{2}-k^{2}}&\dots&\displaystyle\frac{2c_{1}}{c_{1}-c_{N}}\\\
\vdots&\vdots&\ddots&\vdots\\\
\displaystyle\frac{2c_{N}}{k+c_{N}-c_{0}}&\displaystyle\frac{2c_{N}}{c_{N}-c_{1}}&\dots&\displaystyle-
ic_{N}\eta_{N}+1-\frac{2kc_{N}}{(c_{0}-c_{N})^{2}-k^{2}}\end{matrix}\right]$

where for $1\leq j\leq N$, we have defined $\eta_{j}=x-c_{j}t-x_{j}$ with
arbitrary $x_{j}\in\mathbb{R}$ and arbitrary distinct $c_{j}>0$ satisfing
$|c_{j}-c_{0}|>k$. The solution (7) is real-valued and bounded for every
$(x,t)\in\mathbb{R}\times\mathbb{R}$.

The solution in Theorem 1 corresponds to $N$ solitary waves propagating on the
traveling periodic wave (2) with the speeds $c_{1},\dots,c_{N}$. For $N=2$,
Figure 2 shows the top views of three solution surfaces computed from the
exact solution (7) with $N=2$ for three choices of $c_{1}$ and $c_{2}$. Phase
parameters $\xi_{0}$, $x_{1}$, and $x_{2}$ are all set to zero. The left panel
shows propagation of two bright breathers for $c_{1},c_{2}>c_{0}+k$. The
middle panel shows propagation of one bright and one dark breather for
$c_{1}<c_{0}-k$ and $c_{2}>c_{0}+k$. The right panel shows propagation of two
dark breathers for $c_{1},c_{2}<c_{0}-k$.

Figure 2. Left: bright–bright breather with $k=1$, $\phi=3$,
$c_{1}=c_{0}+1.5k$, and $c_{2}=c_{0}+3k$. Middle: bright–dark breather with
$k=1$, $\phi=0.5$, $c_{1}=c_{0}-1.5k$, and $c_{2}=c_{0}+3k$. Right: dark–dark
breather with $k=1$, $\phi=0.5$, $c_{1}=c_{0}-1.75k$, and $c_{2}=c_{0}-1.25k$.

Our work relies on several recent advances in the literature devoted to the
periodic and solitary wave solutions of the BO equation (1).

  1. (1)

Dispersive shock waves (DSWs) in the BO equation were studied with the Whitham
modulation theory in [8, 19, 27, 28]. The Whitham equations for the BO
equation reduce to the Hopf (inviscid Burgers) equations. The breather
solutions equivalent to the solution form (4) were obtained in Appendix B of
[19] but it was missed that two families of bright and dark breathers coexist
on the background of the traveling periodic wave. These earlier results have
motivated many recent works [9, 17, 30, 37], where the DSWs have been studied
in the nonlocal generalizations of the BO equation.

  2. (2)

Lax spectrum associated with traveling periodic waves was rigorously justified
only recently in [12] and was generalized to other periodic potentials [13,
14]. In our notations, the Lax spectrum of the traveling periodic wave (2) is
located at $[\lambda_{0},\mu_{0}]\cup[0,\infty)$, where

$\lambda_{0}=-\frac{1}{2}(c_{0}+k),\qquad\mu_{0}=-\frac{1}{2}(c_{0}-k).$ (8)

The two families of the bright (dark) breathers (6) are associated with the
semi-infinite (finite) spectral gap in the Lax spectrum:

$(\mbox{\rm BB})\;\;-\frac{1}{2}c\in(-\infty,\lambda_{0})\quad\mbox{\rm
and}\quad(\mbox{\rm DB})\;\;-\frac{1}{2}c\in(\mu_{0},0),$ (9)

where $\lambda=-\frac{1}{2}c$ is the eigenvalue of the Lax spectrum associated
with the single solitary wave [10, 20]. The division of breathers into two
families related to the choice of eigenvalue $\lambda$ in the semi-infinite
versus finite spectral gaps of the Lax spectrum is very similar to the KdV
equation [18].

  3. (3)

Modulational stability of traveling periodic waves of the BO equation was
studied with the explicit computations of eigenfunctions of the linearized BO
equation [42] and the asymptotic reduction of the BO equation to the nonlocal
defocusing NLS equation for long perturbations of the small-amplitude waves
[38]. Modulational stability was later confirmed for the BO equation with a
rigorous analysis of the spectral stability problem [6]. Breather solutions
like the one given by (4) may only propagate on the modulationally stable
traveling periodic wave background in agreement with these stability results.

  4. (4)

Particle dynamics of the Calogero–Moser system has been studied in relation to
poles of the meromorphic solutions of the BO equation [7]. In this context and
in the presence of a harmonic potential, an interaction of a single solitary
wave with $N$ solitary waves for large $N$ was studied in [2] (see their
Figure 1). Interactions of solitary waves via the particle dynamics were also
considered for the BO equation and its modifications in [1, 3]. Poles of the
exact solution (4) describe the exact trajectories of particles of the
Calogero–Moser system for such interactions in the BO equation (1).

We note that the KdV equation and the BO equation do not admit conventional
breather solutions which are localized in space and periodic in time. The
proof of non-existence of breathers in these equations was recently given in
[34, 35]. The breather wave packets derived here and in [5, 18] are obtained
as a result of periodic interaction of solitary waves and the traveling
periodic wave. Both the KdV and BO equations are the limiting cases of the
integrable intermediate long-wave (ILW) equation [32]. Similar breathers are
expected to exist in the ILW equation but expressed in the Jacobi theta
functions. Only the BO equation describes such complicated wave interactions
in elementary functions.

The rest of this paper is organized as follows. The derivation of the single-
breather solutions (4) is given in Section 2. The relation between the
existence region (6) and the Lax spectrum of the traveling periodic wave is
explained in Section 3. The exact solution (7) for multi-breather solutions is
derived in Section 4.

##  2\. Derivation of the single-breather solutions

In order to obtain the breather solutions (4), we follow the multi-periodic
solutions obtained in the bilinear method [40]. Equivalently the same
solutions can be obtained by the Bäcklund–Darboux transformation of the BO
equation [36] with an eigenvalue $\lambda=-\frac{1}{2}c$ in either of the two
intervals (9).

The class of multi-periodic and multi-soliton solutions of the BO equation can
be obtained in the form:

$u(x,t)=i\partial_{x}\log\frac{f^{\prime}(x,t)}{f(x,t)},$ (10)

where $f$ ($f^{\prime}$) has zeros only in the upper (lower) half-plane of
$z:=x+iy$ and satisfies the bilinear equation

$i(f^{\prime}_{t}f-f^{\prime}f_{t})=f^{\prime}_{xx}f-2f^{\prime}_{x}f_{x}+f^{\prime}f_{xx}.$
(11)

Let us now consider the simplest solutions of the bilinear equation (11).

1-periodic solutions. The traveling periodic wave (2) is obtained from (10) by
using the following solution of the bilinear equation (11):

$\left\\{\begin{array}[]{l}f(x,t)=1+e^{ik\xi+\phi},\\\
f^{\prime}(x,t)=1+e^{ik\xi-\phi},\end{array}\right.$ (12)

where $\xi=x-c_{0}t-\xi_{0}$, $k,\phi>0$, $\xi_{0}\in\mathbb{R}$, and
$c_{0}=k\coth\phi>0$.

2-periodic solutions. The $2$-periodic solutions are generated from the
following solution of the bilinear equation (11):

$\left\\{\begin{array}[]{l}f(x,t)=1+\alpha e^{ik\xi+\phi}+\alpha
e^{i\kappa\eta+\psi}+e^{ik\xi+i\kappa\eta+\phi+\psi},\\\
f^{\prime}(x,t)=1+\alpha e^{ik\xi-\phi}+\alpha
e^{i\kappa\eta-\psi}+e^{ik\xi+i\kappa\eta-\phi-\psi},\end{array}\right.$ (13)

where $\xi=x-c_{0}t-\xi_{0}$, $\eta=x-ct-\eta_{0}$, $k\phi>0$, $\kappa\psi>0$,
$\xi_{0},\eta_{0}\in\mathbb{R}$, $c_{0}=k\coth\phi>0$, $c=\kappa\coth\psi>0$,
and

$\alpha=\frac{\sqrt{(c-c_{0})^{2}-(\kappa+k)^{2}}}{\sqrt{(c-c_{0})^{2}-(\kappa-k)^{2}}}.$

Since

$f^{\prime}(x,t)=e^{ik\xi+i\kappa\eta-\phi-\psi}\bar{f}(x,t),$

with bar being the complex conjugation, the transformation (10) with $f$,
$f^{\prime}$ in (13) generates a real solution $u(x,t)$ of the BO equation
(1). Moreover, it was shown in [40] that the solution $u(x,t)$ is bounded for
every $(x,t)\in\mathbb{R}\times\mathbb{R}$ and the assumption on zeros of $f$
($f^{\prime}$) in the upper (lower) half-plane is satisfied if parameters $c$,
$c_{0}$, $\kappa$, and $k$ satisfy the constraint

$(c-c_{0})^{2}>(|k|+|\kappa|)^{2}.$ (14)

Breather solutions. In the limit $\kappa\to 0$, there are two ways to satisfy
the constraint (14): $c-c_{0}>k$ or $c-c_{0}<-k$ with $k>0$ assumed here. This
choice leads to the two disjoint solution families (6). The breather solution
(4) arises in the limit $\psi\to 0$ with $\kappa=c\tanh\psi\to 0$ uniquely
defined for fixed $c>0$. Expanding $\alpha$ yields

$\displaystyle\alpha$
$\displaystyle=1-\frac{2k\kappa}{(c-c_{0})^{2}-k^{2}}+\frac{2k^{2}\kappa^{2}}{[(c-c_{0})^{2}-k^{2}]^{2}}+\mathcal{O}(\kappa^{3})$
$\displaystyle=1-\beta\psi+\frac{1}{2}\beta^{2}\psi^{2}+\mathcal{O}(\psi^{3}),$

where $\beta>0$ is given by (5). Using the transformation
$\eta\to\eta+\pi\kappa^{-1}$ and expanding as $\psi\to 0$ yield at the leading
order:

$\left\\{\begin{array}[]{l}\displaystyle
f(x,t)=\beta\psi(1-e^{ik\xi+\phi})-\psi(1+ic\eta)(1+e^{ik\xi+\phi})+\mathcal{O}(\psi^{2}),\\\
\displaystyle
f^{\prime}(x,t)=\beta\psi(1-e^{ik\xi-\phi})+\psi(1-ic\eta)(1+e^{ik\xi-\phi})+\mathcal{O}(\psi^{2}).\end{array}\right.$
(15)

Using (10) with the leading $\mathcal{O}(\psi)$ order of $f$ and $f^{\prime}$
yields the breather solutions (4) which are obviously real. To show that the
breather solutions (4) are bounded for every
$(x,t)\in\mathbb{R}\times\mathbb{R}$ and the assumption on zeros of $f$
($f^{\prime}$) in the upper (lower) half-plane is satisfied, we consider the
zeros of $f(x,t)$ at the leading $\mathcal{O}(\psi)$ order in (15), that is,

$\beta(1-e^{ik\xi+\phi})-(1+ic\eta)(1+e^{ik\xi+\phi})=0.$

Zeros are defined from

$e^{ik\xi+\phi}=\frac{\beta-1-ic\eta}{\beta+1+ic\eta}$

for which

${\rm Im}(\xi)=k^{-1}(\phi+q),\quad
q:=\frac{1}{2}\log\frac{(\beta+1)^{2}+c^{2}\eta^{2}}{(\beta-1)^{2}+c^{2}\eta^{2}}.$

Since $k,\beta,\phi,q>0$, zeros of $f$ are located in the upper half-plane
away from the real axis.

$2$-soliton solutions. The $2$-soliton solutions arise from the breather
solutions in the limit $\phi\to 0$ with $k=c_{0}\tanh\phi\to 0$ uniquely
defined for fixed $c_{0}>0$. After the transformation $\xi\to\xi+\pi k^{-1}$,
the expansion (15) is further reduced to the form:

$\left\\{\begin{array}[]{l}\displaystyle
f(x,t)=\left[(1+ic_{0}\xi)(1+ic\eta)+\frac{4cc_{0}}{(c-c_{0})^{2}}\right]\phi\psi+\mathcal{O}(\phi^{2}\psi^{2}),\\\
\displaystyle
f^{\prime}(x,t)=\left[(1-ic_{0}\xi)(1-ic\eta)+\frac{4cc_{0}}{(c-c_{0})^{2}}\right]\phi\psi+\mathcal{O}(\phi^{2}\psi^{2}).\end{array}\right.$

The $2$-soliton solutions are well-known from [25, 36, 40]. After the limiting
transition $k\to 0$ in (6), the bright breathers correspond to $c>c_{0}$ and
the dark breathers correspond to $0<c<c_{0}$, where $c_{0}$ and $c$ are wave
speeds of the two distinct solitary waves. This agrees with the property that
the wave amplitudes are proportional to the wave speeds. The new soliton’s
speed $c$ exceeds the background soliton’s speed $c_{0}$ for a bright breather
which degenerates into a larger soliton overtaking the background soliton. The
new soliton’s speed $c$ is smaller than the background soliton’s speed $c_{0}$
for a dark breather which degenerates into a smaller soliton staying behind
the background soliton.

##  3\. Existence of breathers and the Lax spectrum of the periodic wave

The BO equation (1) arises as a compatibility condition of the following
system of Lax equations [10, 20]:

$\left\\{\begin{array}[]{l}i\varphi^{+}_{x}+\lambda(\varphi^{+}-\varphi^{-})+u\varphi^{+}=0,\\\
i\varphi^{\pm}_{t}-2i\lambda\varphi^{\pm}_{x}+\varphi^{\pm}_{xx}-[\pm
iu_{x}+H(u_{x})]\varphi^{\pm}=0,\end{array}\right.$ (16)

where $\varphi^{\pm}$ are analytic functions in $\mathbb{C}^{\pm}$. If
$u(x+L)=u(x)$ is periodic in $x$ with period $L=\frac{2\pi}{k}$ as in (2),
then Lax spectrum can be defined from admissible values of the linear operator
acting on $\varphi^{+}$:

$L_{u}:=-i\partial_{x}-P^{+}(u\;\cdot),$ (17)

where $P^{+}=\frac{1}{2}(1-iH)$ is the projection operator acting on functions
defined on the real line and returning functions which can be analytically
continued in $\mathbb{C}^{+}$.

As is shown in Proposition 2.2 of [12], the Lax spectrum of $L_{u}$ on the
periodic domain consists of simple eigenvalues enumerated as
$\\{\lambda_{n}\\}_{n=0}^{\infty}$ such that $\lambda_{n+1}\geq\lambda_{n}+k$
for $n\geq 0$, where $k=\frac{2\pi}{L}$ ($k=1$ was used in [12] with
$L=2\pi$). On the other hand, as is shown in Proposition C.2 of [12], the Lax
spectrum of $L_{u}$ on the real line is absolutely continuous and consists of
a union of bands $\cup_{n=0}^{\infty}[\lambda_{n},\lambda_{n}+k]$.

The traveling periodic wave (2) corresponds to a particular case of the one-
gap potential for which $\lambda_{n+1}=\lambda_{n}+k$ for every
$n\in\mathbb{N}$. As follows from the computations in Appendix B of [12]
(generalized in our notations), that $\lambda_{0}$ is given by (8) and
$\lambda_{n}=n-1$ for $n\in\mathbb{N}$. Then, $\mu_{0}=\lambda_{0}+k$ is
defined in (8) in agreement with formula (1.19) in [8], as well as Proposition
C.2 in [12].

To summarize, the Lax spectrum for the traveling periodic wave (2) is located
in $[\lambda_{0},\mu_{0}]\cup[0,\infty)$. The existence region (6) of the
breather solutions (4) corresponds to the two choices in (9) for
$\lambda=-\frac{1}{2}c$, the eigenvalue of the point spectrum associated to
the breather.

The rest of this section contains computations of the explicit expressions for
eigenfunctions of the Lax system (16), which visualize the general results
above. As an interesting new feature, we discover the following properties:

  1. (P1)

Eigenvalues $\\{\lambda_{n}\\}_{n=0}^{\infty}$ are defined from the zero-mean
condition on $\varphi^{-}$ among the $L$-periodic and
$\mathbb{C}^{\pm}$-analytic eigenfunctions $\varphi^{\pm}$.

  2. (P2)

Eigenvalues $\\{\mu_{n}\\}_{n=0}^{\infty}$ with $\mu_{n}=\lambda_{n}+k$ is
defined from the zero-mean condition on $\varphi^{+}$ among the $L$-periodic
and $\mathbb{C}^{\pm}$-analytic eigenfunctions $\varphi^{\pm}$.

  3. (P3)

The $L$-periodic and $\mathbb{C}^{\pm}$-analytic eigenfunctions
$\varphi^{\pm}$ without the zero-mean conditions exist for every
$\lambda\in\mathbb{R}$.

There exist two general solutions of the linear system (16). The first
solution has $\varphi^{-}\equiv 0$ and the second solution has both components
$\varphi^{+}$ and $\varphi^{-}$ nonzero.

If $\varphi^{-}\equiv 0$, the exact solution of the linear system (16) exists
in the form:

$\varphi^{+}(x,t)=e^{i\lambda
x+i\lambda^{2}t}\frac{1+e^{ik\xi+\phi}}{1+e^{ik\xi-\phi}}=:e^{i\lambda
x+i\lambda^{2}t}p^{+}(\xi,t),$ (18)

where $p^{+}(\cdot,t)$ is an $L$-periodic and $\mathbb{C}^{+}$-analytic
function with period $L=2\pi/k$ for every $t\in\mathbb{R}$. If we require that
$\varphi^{+}(\cdot,t)$ be bounded as ${\rm Im}(x)\to+\infty$ for every
$t\in\mathbb{R}$, then we must take $\lambda\geq 0$. The function
$\varphi^{+}(\cdot,t)$ is $L$-periodic if and only if
$\lambda\in\\{\lambda_{n}\\}_{n=1}^{\infty}$ with $\lambda_{n}:=n$ for
$n\in\mathbb{N}$. It has zero mean if and only if
$\lambda\in\\{\mu_{n}\\}_{n=1}^{\infty}$ with $\mu_{n}:=n+1$ for
$n\in\mathbb{N}$. This proves (P1) and (P2) for
$\\{\lambda_{n}\\}_{n=1}^{\infty}$ and $\\{\mu_{n}\\}_{n=1}^{\infty}$.

If both components $\varphi^{+}$ and $\varphi^{-}$ are nonzero, we can
separate the variables as

$\varphi^{\pm}(x,t)=m^{\pm}(\xi)e^{-i\Omega t}.$

The system (16) is rewritten in ordinary derivatives:

$\left\\{\begin{array}[]{l}i(m^{+})^{\prime}+\lambda(m^{+}-m^{-})+um^{+}=0,\\\
-(m^{\pm})^{\prime\prime}+i(c+2\lambda)(m^{\pm})^{\prime}+[\pm
iu^{\prime}+H(u^{\prime})]m^{\pm}=\Omega m^{\pm},\end{array}\right.$ (19)

where $u=u(\xi)$ is given by (2). Since the second equation of the system (19)
is the second-order ODE, a general solution can be obtained as a superposition
of two linearly independent solutions:

$m^{+}(\xi)=e^{i(\theta+\lambda+\frac{1}{2}c)\xi}\frac{(2\theta+k)+(2\theta-k)e^{ik\xi-\phi}}{1+e^{ik\xi-\phi}}A_{1}+e^{i(-\theta+\lambda+\frac{1}{2}c)\xi}\frac{(2\theta-k)+(2\theta+k)e^{ik\xi-\phi}}{1+e^{ik\xi-\phi}}A_{2}$

and

$m^{-}(\xi)=e^{i(\theta+\lambda+\frac{1}{2}c)\xi}\frac{(2\theta+k)+(2\theta-k)e^{ik\xi+\phi}}{1+e^{ik\xi+\phi}}B_{1}+e^{i(-\theta+\lambda+\frac{1}{2}c)\xi}\frac{(2\theta-k)+(2\theta+k)e^{ik\xi+\phi}}{1+e^{ik\xi+\phi}}B_{2},$

where $\theta>0$ is found from $\theta^{2}=\Omega+(\lambda+\frac{1}{2}c)^{2}$
and the constants $A_{1,2}$, $B_{1,2}$ are arbitrary. Substituting the
expressions for $m^{\pm}$ into the first equation of the system (19) only
defines equation between the constants:

$2\lambda B_{1}=-(2\theta+c)A_{1},\quad 2\lambda B_{2}=(2\theta-c)A_{2},$

whereas $\lambda$ and $\theta$ are left arbitrary.

If we require $\varphi^{\pm}(\cdot,t)$ to be bounded as ${\rm
Im}(x)\to\pm\infty$, then we must set uniquely $\theta=|\lambda+\frac{1}{2}c|$
which yields $\Omega=0$ and $\varphi^{\pm}(x,t)=m^{\pm}(\xi)$. For this choice
of $\theta$, the linear superposition of two solutions yields only one
solution in the form (after taking $A_{1,2}=-1$):

$m^{+}(\xi)=\frac{(2\lambda+c-k)+(2\lambda+c+k)e^{ik\xi-\phi}}{1+e^{ik\xi-\phi}}$
(20)

and

$m^{-}(\xi)=\frac{(2\lambda+c-k)+(2\lambda+c+k)e^{ik\xi+\phi}}{1+e^{ik\xi+\phi}}.$
(21)

Parameter $\lambda$ is left arbitrary, whereas we can see that
$\varphi^{\pm}(\cdot,t)=m^{\pm}(\cdot-c_{0}t)$ are $L$-periodic and
$\mathbb{C}^{\pm}$-analytic for every $t\in\mathbb{R}$ and every
$\lambda\in\mathbb{R}$. This proves (P3).

It remains to confirm (P1) and (P2) for $\lambda_{0}$ and $\mu_{0}$. If
$\lambda=\lambda_{0}=-(c+k)/2$, then the eigenfunction (up to a norming
factor) is found from (20) and (21) as

$m^{+}(\xi)=\frac{1}{1+e^{ik\xi-\phi}},\qquad
m^{-}(\xi)=\frac{1}{1+e^{ik\xi+\phi}}=\frac{e^{-ik\xi-\phi}}{1+e^{-ik\xi-\phi}}.$
(22)

It follows from (22) that the mean-value of $m^{-}$ is zero. If
$\lambda=\mu_{0}=-(c-k)/2$, then the eigenfunction (up to a norming factor) is
found from (20) and (21) as

$m^{+}(\xi)=\frac{e^{ik\xi-\phi}}{1+e^{ik\xi-\phi}},\qquad
m^{-}(\xi)=\frac{e^{ik\xi+\phi}}{1+e^{ik\xi+\phi}}=\frac{1}{1+e^{-ik\xi-\phi}}.$
(23)

It follows from (23) that the mean-value of $m^{+}$ is zero. For other values
of $\lambda\in\mathbb{R}$, the mean values of $m^{+}$ and $m^{-}$ defined by
(20) and (21) are nonzero. This proves (P1) and (P2) for $\lambda_{0}$ and
$\mu_{0}$. Note that (20) and (21) agree with (1.11) and (1.22) in [8].

##  4\. Derivation of the multi-breather solutions

To obtain the multi-breather solutions of the BO equation (1) on the
background of the traveling wave (2), we define solutions in the form (10)
with $f=\det(F)$ and $f^{\prime}=\det(F^{\prime})$ being determinants of the
square matrices $F$ and $F^{\prime}$ with the following elements:

$\displaystyle F_{jm}$
$\displaystyle=\displaystyle\frac{1}{k_{j}}\exp\left(ik_{j}\xi_{j}+\phi_{j}+\frac{1}{2}\sum\limits_{m\neq
j}A_{jm}\right)\delta_{jm}+\frac{2}{k_{j}+k_{m}+c_{j}-c_{m}},$ (24)
$\displaystyle F_{jm}^{\prime}$
$\displaystyle=\displaystyle\frac{1}{k_{j}}\exp\left(ik_{j}\xi_{j}-\phi_{j}+\frac{1}{2}\sum\limits_{m\neq
j}A_{jm}\right)\delta_{jm}+\frac{2}{k_{j}+k_{m}+c_{j}-c_{m}},$ (25)

where $\xi_{j}=x-c_{j}t-x_{j}$, $k_{j}\phi_{j}>0$, $x_{j}\in\mathbb{R}$,
$c_{j}=k_{j}\coth\phi_{j}>0$, and

$\displaystyle
e^{A_{jm}}=\frac{(c_{j}-c_{m})^{2}-(k_{j}-k_{m})^{2}}{(c_{j}-c_{m})^{2}-(k_{j}+k_{m})^{2}}.$
(26)

These determinant representations were derived in [29, Eqs. (2.16)] based on
the determinant formulas used in [8]. They are particularly useful in the
derivation of multi-soliton solutions by degeneration of $k_{j}\to 0$ with
$c_{j}$ being fixed, see [29, Eqs. (2.48)]. It also follows from the
representation, see [29, Eq.(2.29), Appendix A],

$f^{\prime}=\det(F^{\prime})=\exp\left(\sum_{j}ik_{j}\xi_{j}-\phi_{j}\right)\det(\bar{F})=\exp\left(\sum_{j}ik_{j}\xi_{j}-\phi_{j}\right)\bar{f}$
(27)

that the solution $u(x,t)$ in the form (10) is a real quantity given by

$u=-\sum_{j}k_{j}+i\partial_{x}\log\frac{\det(\bar{F})}{\det(F)},$

where the bar denotes the complex conjugation. Moreover, it was shown in Lemma
1.1 of [8] that the solution $u(x,t)$ is bounded for every
$(x,t)\in\mathbb{R}\times\mathbb{R}$ and the assumption on zeros of $f$
($f^{\prime}$) in the upper (lower) half-plane is satisfied if parameters of
the solution satisfy the constraints

$\displaystyle(c_{j}-c_{m})^{2}>(|k_{j}|+|k_{m}|)^{2},\quad j\neq m.$ (28)

Constraints (28) ensure that the Lax spectrum of the multiperiodic solutions
is located on

$\cup_{j}[\lambda_{j},\mu_{j}]\cup[0,\infty),\quad\lambda_{j}=-\frac{c_{j}+|k_{j}|}{2},\quad\mu_{j}=-\frac{c_{j}-|k_{j}|}{2}$

with the spectral band $[\lambda_{j},\mu_{j}]$ being disjoint from the other
spectral bands [8, 12]. See the upper panel of Figure 3 for illustration.

Figure 3. The degeneration procedure used to obtain multi-breather solutions.
Top: three finite bands in the Lax spectrum associated with the three-periodic
solutions. Bottom: two eigenvalues and one finite band in the Lax spectrum
associated with the two-breather solution on the periodic wave.

Multi-breather solutions. For multi-breather solutions, we consider the square
matrices $F$ and $F^{\prime}$ of size $(N+1)\times(N+1)$ and reserve the first
element for the traveling periodic wave (2) with $\xi=x-c_{0}t-\xi_{0}$,
$k_{0}=k>0$, and $\phi_{0}=\phi>0$. For the remaining $N$ elements, we use
enumeration $1\leq j\leq N$, add the phase transformations $x_{j}\to x_{j}+\pi
k_{j}^{-1}$, denote $\eta_{j}=x-c_{j}t-x_{j}$, and take the limit $k_{j}\to
0$, $\phi_{j}\to 0$ with $c_{j}=k_{j}\coth\phi_{j}>0$ being fixed. This
degeneration procedure transforms elements of $F$ and $F^{\prime}$ to the
form:

$\displaystyle F_{00}$
$\displaystyle=\frac{1}{k}\left[1+e^{ik\xi+\phi}\right],\quad
F_{0m}=\frac{2}{k+c_{0}-c_{m}},\quad F_{m0}=\frac{2}{k+c_{m}-c_{0}},$
$\displaystyle F_{00}^{\prime}$
$\displaystyle=\frac{1}{k}\left[1+e^{ik\xi-\phi}\right],\quad
F_{0m}^{\prime}=\frac{2}{k+c_{0}-c_{m}},\quad
F_{m0}^{\prime}=\frac{2}{k+c_{m}-c_{0}},$

and for $1\leq j\neq m\leq N$:

$\displaystyle F_{jj}$
$\displaystyle=-i\eta_{j}-\frac{1}{c_{j}}-\frac{2k}{(c_{0}-c_{j})^{2}-k^{2}},\quad
F_{jm}=\frac{2}{c_{j}-c_{m}},$ $\displaystyle F_{jj}^{\prime}$
$\displaystyle=-i\eta_{j}+\frac{1}{c_{j}}-\frac{2k}{(c_{0}-c_{j})^{2}-k^{2}},\quad
F_{jm}^{\prime}=\frac{2}{c_{j}-c_{m}}.$

The terms with $e^{A_{jm}}$ with $1\leq j,m\leq N$ do not contribute to the
limit because of the power expansion:

$e^{A_{jm}}=1+\frac{4k_{j}k_{m}}{(c_{j}-c_{m})^{2}}+\mathcal{O}(k_{j}^{2}k_{m}^{2}).$

The terms with $e^{A_{0j}}$ with $1\leq j\leq N$ do not contribute to the
limit for $F_{00}$ but contribute to the limit for $F_{jj}$ because of the
power expansion:

$e^{A_{0j}}=1+\frac{4kk_{j}}{(c_{0}-c_{j})^{2}-k^{2}}+\mathcal{O}(k_{j}^{2}).$

The choice of $c_{j}$ for $1\leq j\leq N$ is arbitrary as long as
$|c_{j}-c_{0}|>k$ and all $c_{j}$ are distinct. These constraints follow from
(28). After the degeneration, the Lax spectrum of the breather solutions
consist on

$\cup_{j=1}^{N}\lambda_{j}\cup[\lambda_{0},\mu_{0}]\cup[0,\infty),\quad\lambda_{0}=-\frac{c_{0}+k}{2},\quad\mu_{0}=-\frac{c_{0}-k}{2},\quad\lambda_{j}=-\frac{c_{j}}{2},$

where simple eigenvalues $\\{\lambda_{j}\\}_{j=1}^{N}$ are located either in
$(-\infty,\lambda_{0})$ or $(\mu_{0},0)$. See the lower panel of Figure 3 for
illustration. The degeneration procedure has been used for the KdV equation in
[5]. It has been widely used to obtain multi-soliton solutions from the multi-
periodic solutions of the BO equation [29, 40].

Multiplying the rows of $F$ and $F^{\prime}$ by $k,c_{1},\dots,c_{N}$
respectively yields the solution formula (7) in Theorem 1. The solution
$u(x,t)$ is real which follows from (27) in the same degeneration limit as

$f^{\prime}=\det(F^{\prime})=(-1)^{N}e^{ik\xi-\phi}\det(\bar{F})=(-1)^{N}e^{ik\xi-\phi}\bar{f}.$
(29)

Since (24), (25), (26), and (28) yield roots of $f=\det(F)$
($f^{\prime}=\det(F^{\prime})$) in the upper (lower) half-plane for every
$k_{j}>0$ in Lemma 1.1 of [8], the roots can only approach the real line but
cannot cross the real line in the limit $k_{j}\to 0$. Since $f^{\prime}$ is a
related to the complex conjugation of $f$ by (29), the roots of $f$ and
$f^{\prime}$ on the real line coincide up to their multiplicities and the
quotient in the solution form (10) can be used to redefine $f$ and
$f^{\prime}$ free of roots on the real line. This proves that the solution
$u(x,t)$ arising in the limit $k_{j}\to 0$ is bounded for every
$(x,t)\in\mathbb{R}\times\mathbb{R}$.

Similar to Lemma 1.1 in [8], it should be possible to represent the solution
of the Lax equations (16) for the multi-breather solutions of Theorem 1 and
prove rigorously that the individual determinants of the matrices $F$ and
$F^{\prime}$ given below (7) do not vanish for every
$(x,t)\in\mathbb{R}\times\mathbb{R}$. This would rule out the possibility of
coalescence of roots of $f$ and $f^{\prime}$ at the real axis in the limit
$k_{j}\to 0$. The proof of the above fact is left out for further work.

Example with $N=1$. Expanding the determinants in Theorem 1 for $N=1$, we
obtain

$\displaystyle\det(F)$
$\displaystyle=-(1+e^{ik\xi+\phi})(ic\eta+1+\beta)+2\beta,$
$\displaystyle\det(F^{\prime})$
$\displaystyle=-(1+e^{ik\xi-\phi})(ic\eta-1+\beta)+2\beta,$

where $\beta$ is given by (5) and $\eta=x-ct-\eta_{0}$ is defined with
$c_{1}=c$. These expressions for $f=\det(F)$ and $f^{\prime}=\det(F^{\prime})$
are equivalent to those in (15) and lead to the exact solution (4) shown on
Figure 1 for either $c>c_{0}+k$ or $c<c_{0}-k$. The explicit expressions for
$\det(F)$ and $\det(F^{\prime})$ agree with (29) with $N=1$.

Example with $N=2$. Expanding the determinants in Theorem 1 for $N=2$ and
simplifying the answers, we obtain

$\displaystyle\det(F)$
$\displaystyle=(1+e^{ik\xi+\phi})\left[(1+ic_{1}\eta_{1})(1+ic_{2}\eta_{2})+\frac{4c_{1}c_{2}}{(c_{1}-c_{2})^{2}}+\beta_{1}\beta_{2}\right]$
$\displaystyle\quad+(e^{ik\xi+\phi}-1)\left[\beta_{1}(1+ic_{2}\eta_{2})+\beta_{2}(1+ic_{1}\eta_{1})\right],$
$\displaystyle\det(F^{\prime})$
$\displaystyle=(1+e^{ik\xi-\phi})\left[(1-ic_{1}\eta_{1})(1-ic_{2}\eta_{2})+\frac{4c_{1}c_{2}}{(c_{1}-c_{2})^{2}}+\beta_{1}\beta_{2}\right]$
$\displaystyle\quad+(1-e^{ik\xi-\phi})\left[\beta_{1}(1-ic_{2}\eta_{2})+\beta_{2}(1-ic_{1}\eta_{1})\right],$

where

$\beta_{j}:=\frac{2kc_{j}}{(c_{j}-c_{0})^{2}-k^{2}},\quad 1\leq j\leq 2.$

The explicit expressions for $\det(F)$ and $\det(F^{\prime})$ agree with (29)
with $N=2$. By using these formulas in (10), we obtain three different
solutions for three choices of $c_{1}$ and $c_{2}$: a bright–bright breather
if $c_{1},c_{2}>c_{0}+k$, a bright-dark breather if $c_{1}<c_{0}-k$ and
$c_{2}>c_{0}+k$, and a dark-dark breather if $c_{1},c_{2}<c_{0}-k$. The
solution surfaces for these breather solutions are shown on Figure 2.

Acknowledgement. This work was supported in part by the National Natural
Science Foundation of China (No. 12371248) and the Project “333” of Jiangsu
Province. The authors thank A. Abanov for asking questions if the breathers
similar to those in the KdV equation also exist in the BO equation. The
authors also thank P. Gérard and P. Miller for discussions of the Lax spectrum
related to the periodic wave solutions and Y. Matsuno for discussions of the
zeros of $f$ and $f^{\prime}$.

## References

  * [1] A. G. Abanov, E. Betterlheim, and P. B. Wiegmann, “Integrable hydrodynamics of Calogero–Sutherland model: bidirectional Benjamin–Ono equation”, J. Phys. A: Math. Theor. 42 (2009) 135201
  * [2] A. G. Abanov, A. Gromov, and M, Kulkarni, “Soliton solutions of a Calogero model in a harmonic potential”, J. Phys. A: Math. Theor. 44 (2011) 295203 (21 pages)
  * [3] A. G. Abanov and P. B. Wiegmann, “Quantum hydrodynamics, the quantum Benjamin–Ono equation, and the Calogero model”, Phys. Rev. Lett. 95 (2005) 076402
  * [4] M. J. Ablowitz, J. T. Cole, G. A. El, M. A. Hoefer, and X. Luo, “Soliton-mean field interaction in Korteweg-de Vries dispersive hydrodynamics,” Stud. Appl. Math. 151 (2023) 795–856
  * [5] M. Bertola, R. Jenkins, and A. Tovbis, “Partial degeneration of finite gap solutions to the Korteweg–de Vries equation: soliton gas and scattering on elliptic background”, Nonlinearity 36 (2023) 3622–3660
  * [6] J. C. Bronski and V. M. Hur, “Modulational instability and variational structure”, Stud. Appl. Math. 132 (2014) 285–331
  * [7] K. M. Case, “Meromorphic solutions of the Benjamin–Ono equation”, Phys. A 96 (1979) 173–182
  * [8] S. Yu. Dobrokhotov and I. M. Krichever, “Multiphase solutions of the Benjamin-Ono equation and their averaging”, Math. Notes 49 (1991) 583–594
  * [9] G. A. El, L. T. K. Ngueyn, and N. F. Smyth, “Dispersive shock waves in systems with nonlocal dispersion of Benjamin–Ono type”, Nonlinearity 31 (2018) 1392–1416
  * [10] A. S. Fokas and M. J. Ablowitz, “The inverse scattering transform for the Benjamin–Ono equation: a pivot to multidimensional problems”, Stud. Appl. Math. 68 (1983) 1–10
  * [11] J. Garnier, G. Xu, S. Trillo, and A. Picozzi, “Incoherent dispersive shocks in the spectral evolution of random waves”, Phys. Rev. Lett. 111 (2013) 113902 (5 pages)
  * [12] P. Gérard and T. Kappeler, “On the integrability of the Benjamin–Ono equation on the torus”, Comm. Pure Appl. Math. 74 (2021) 1685–1747
  * [13] P. Gérard, T. Kappeler, and P. Topalov, “On the spectrum of the Lax operator of the Benjamin–Ono equation on the torus”, J. Funct. Anal. 279 (2020) 108762 (75 pages)
  * [14] P. Gérard, T. Kappeler, and P. Topalov, “On the Benjamin–Ono equation on $\mathbb{T}$ and its periodic and quasiperiodic solutions”, J. Spectr. Theory 12 (2022) 169–193
  * [15] F. Gesztesy and R. Svirsky, “(m)KdV solitons on the background of quasi-periodic finite-gap solutions”, Mem. Amer. Math. Soc. 118 (1995) 563 (88 pages)
  * [16] R. Grimshaw, “Evolution equations for long nonlinear internal waves in stratified shear flows”, Stud. Appl. Math. 65 (1981) 159–188
  * [17] R. H. J. Grimshaw, N. F. Smyth, and Y. A. Stepanyants, “Interaction of internal solitary waves with long periodic waves within the rotation modified Benjamin–Ono equation”, Physica D 419 (2021) 132867 (10 pages)
  * [18] M. Hoefer, A. Mucalica, and D. E. Pelinovsky, “KdV breathers on cnoidal wave background”, J. Phys. A: Mathem. Theor. 56 (2023) 185701 (25 pages)
  * [19] M. C. Jorge, A. A. Minzoni, and N. F. Smyth, “Modulation solutions for the Benjamin–Ono equation”, Physica D 132 (1999) 1–18
  * [20] D. J. Kaup and Y. Matsuno, “‘The inverse scattering transform for the Benjamin–Ono equation”, Stud. Appl. Math. 101 (1998) 73–98
  * [21] E. A. Kuznetsov and A. V. Mikhailov, “Stability of stationary waves in nonlinear weakly dispersive media”, Sov. Phys. JETP 40 (1974) 855–859
  * [22] L. Ling and X. Sun, “Multi-elliptic-dark soliton solutions in the defocusing nonlinear Schrödinger equation”, Appl. Math. Lett. 148 (2023) 108866 (9 pages)
  * [23] S. Lou, X. Cheng, and X. Tang, “Dressed dark solitons of the defocusing nonlinear Schrödinger equation”, Chinese Phys. Lett. 31 (2014) 070201 (4 pages)
  * [24] M. D. Maiden, D. V. Anderson, A. A. Franco, G. A. El, and M. A. Hoefer, “Solitonic dispersive hydrodynamics: Theory and observation,” Phys. Rev. Lett. 120 (2018) 144101
  * [25] Y. Matsuno, “Exact multi-soliton solution of the Benjamin–Ono equation”, J. Phys. A : Math. Gen. 12 (1979) 619–621
  * [26] Y. Matsuno, Bilinear Transformation Method (Academic Press, New York, 1984)
  * [27] Y. Matsuno, “Nonlinear modulation of periodic waves in the small dispersion limit of the Benjamin–Ono equation”, Phys. Rev. E 58 (1998) 7934–7940
  * [28] Y. Matsuno, “The small dispersion limit of the Benjamin–Ono equation and the evolution of a step initial condition”, J. Phys. Soc. Japan 67 (1998) 1814–1817
  * [29] Y. Matsuno, “New representations of multiperiodic and multisoliton solutions for a class of nonlocal soliton equations”, J. Phys. Soc. Japan 73 (2004) 3285–3293
  * [30] Y. Matsuno, V. S. Shchenovich, A. M. Kamchatnov and R.A. Kraenkel, “Whitham method for the Benjamin–Ono–Burgers equation and dispersive shocks”, Phys. Rev. E 75 (2007) 016307 (5 pages)
  * [31] Y. Mao, S. Chandramouli, W. Xu, and M. Hoefer, “Observation of traveling breathers and their scattering in a two-fluid system”, Phys. Rev. Lett. 131 (2023) 147201 (7 pages)
  * [32] T. Miloh, “On periodic and solitary wavelike solutions of the intermediate long-wave equation”, J. Fluid Mech. 211 (1990) 617–627
  * [33] A. Mucalica and D. E. Pelinovsky, “Dark breathers on the snoidal wave background in the dofucusing mKdV equation”, arXiv: 2312.08969 (2023)
  * [34] C. Munoz and G. Ponce, “Breathers and the dynamics of solutions in KdV type equations”, Comm. Math. Phys. 367 (2019) 581–598
  * [35] C. Munoz and G. Ponce, “On the asymptotic behavior of solutions to the Benjamin–Ono equation”, Proceedings of the AMS 147 (2019) 5303–5312
  * [36] A. Nakamura, “Bäcklund transform and conservation laws of the Benjamin–Ono equation”, J. Phys. Soc. Japan 47 (1979) 1335–1340
  * [37] L. T. K. Ngueyn and N. F. Smyth, “Dispersive shock waves for the Boussinesq Benjamin–Ono equation”, Stud. Appl. Math. 147 (2021) 32–59
  * [38] D. Pelinovsky, “Intermediate nonlinear Schroödinger equation for internal waves in a fluid of finite depth”, Phys. Lett. A 197 (1995) 401–406
  * [39] A. Porter and N. F. Smyth, “Modelling the morning glory of the Gulf of Carpentaria”, J. Fluid Mech. 454 (2002) 1–20
  * [40] J. Satsuma and Y. Ishimori, “Periodic wave and rational soliton solutions of the Benjamin–Ono equation”, J. Phys. Soc. Japan 46 (1979) 681–687
  * [41] H. Shin, “The dark soliton on a cnoidal wave background”, J. Phys. A: Math. Gen. 38 (2005) 3307–3315
  * [42] M. D. Spector and T. Miloh, “Stability of nonlinear periodic internal waves in a deep stratified fluid”, SIAM J. Appl. Math. 54 (1994) 688–707
  * [43] P. Sprenger, M. A. Hoefer, and G. A. El, “Hydrodynamic optical soliton tunneling”, Phys. Rev. E 97 (2018) 032218
  * [44] D. A. Takahashi, “Integrable model for density-modulated quantum condensates: Solitons passing through a soliton lattice“, Phys. Rev E. 93 (2016) 062224 (20 pages)


Note that to derive the classical limit of the quantum “soft factor” from (<ref>), first we need to incorporate the sum over particles indexed by $i$. Subsequently, in their respective expressions, when the orbital momentum operator operates on $\mathcal{A}^{(N)}$, it must be substituted with the classical orbital momenta of the scattered objects. Furthermore, based on explicit classical computations of gravitational waveforms presented in <cit.>, it has been conjectured that the validity of the relation (<ref>) holds true solely when, in the classical limit, the soft graviton factor is calculated using the retarded propagator for the graviton field instead of Feynman propagator. In practice this prescription suggests that in the integral representations of (<ref>) and (<ref>), the term $(\ell^2-i\epsilon)^{-1}$ should be substituted with $-\left((\ell^0+i\epsilon)^2-\vec{\ell}^2\right)^{-1}$ in order to extract $\mathbb{S}^{gr}_{\text{classical}}$. With this substitution, $\mathbb{S}^{gr}_{\text{classical}}$ only receives contributions from the part of the loop integrals in which the virtual hard particle propagator goes on-shell (referred to as the potential region) and the graviton propagator with principal value. However, it does not receive contributions from the part of the loop integrals in which the virtual graviton goes on-shell (known as the radiation-reaction region). Unfortunately, we lack a fundamental understanding of why the classical limit of the quantum soft factor does not include contributions from radiation reactions.

For a $2\rightarrow 2$ scattering process with a large impact parameter or low momentum transfer, the contribution of radiation reaction to the quantum soft factor turns out to be suppressed compared to the contribution from the potential region (classical), as discussed in <cit.>. This has also been established in <cit.> using the KMOC formalism <cit.> after substituting the final momenta of the scattered particles in terms of initial momenta and the perturbatively computed momentum impulse in terms of initial scattering data and specified interaction. However, for hard scattering (small impact parameter scattering), the reason why the radiation reaction contribution to the quantum soft factor does not affect the classical waveform in the classical limit has not yet been resolved. We believe that a generalization of recent investigations into deriving classical gravitational waveforms from Eikonal exponentiation in <cit.> could potentially resolve this puzzle even for hard scattering. It would be interesting to explore whether there are any observable consequences resulting from the contribution of radiation reaction in the quantum soft factor.

There have been many applications of soft theorems both in the context of scattering amplitudes and in relation to gravitational memory. For instance, the universal characteristic of Weinberg's soft theorem imposes an infinite hierarchy of constraints on the linear momentum impulse within the KMOC formalism, as derived in <cit.>. It is also anticipated that the universality of the $\ln\omega$ soft theorem should impose an infinite hierarchy of constraints on the angular momentum impulse within the KMOC formalism.[BS acknowledges Alok Laddha for private communication.] As discussed in Section <ref>, while the soft theorems alone cannot impose non-trivial constraints on the quantum theory in the UV, combining the results of the soft theorems with certain physical assumptions about scattering amplitudes, such as analyticity, unitarity and crossing grants them the ability to constrain the UV quantum theory. As an illustration of this concept, the article <cit.> derived non-perturbative bounds on the $a$-anomaly coefficient of the UV conformal field theory (CFT), the deformation of which leads to a massive QFT along the renormalization group (RG) flow. These bounds were established by incorporating the constraints from the double soft dilaton theorem within the framework of non-perturbative S-matrix bootstrap. For numerous astrophysical scattering events, the order of magnitude of gravitational tail memory follows from the $\ln\omega$ soft theorem has been estimated in <cit.>. In certain classical scattering scenarios, the gravitational waveform resulting from the $\ln\omega$ soft theorem has also been derived in <cit.>, carrying observable consequences in the present era of gravitational wave physics.

We would like to thank Waël Aoun, P.V. Athira, Hofie Hannesdottir, Alok Laddha, Joao Penedones and Ashoke Sen for useful discussion.
The work of HK is supported by NSF grant PHY 2210533. The work of BS is supported by the Simons Foundation grant 488649 (Simons Collaboration on the Nonperturbative Bootstrap) and by the Swiss National Science Foundation through the National Centre of Competence in Research SwissMAP.

§ INTERMEDIATE STEPS IN DERIVING SOFT PHOTON THEOREM

Let us evaluate the following expression appears in (<ref>) after substituting the vertex (<ref>) and using the identities (<ref>), (<ref>) and (<ref>),
ζ_1μ   ≡    ϵ_i^T(-p_i)ε^ρ(k)Γ^(3)_ρ(p_i,-p_i-k,k)Ξ_i(-p_i-k)Γ^(3)_μ(p_i+k,-p_i-k+ℓ,-ℓ)Ξ_i(-p_i-k+ℓ)

=   e_i^2ϵ_i^T[ε^ρ𝒦_i(-p_i)p_i^ρΞ_i(-p_i)+ε^ρk^σ𝒦_i(-p_i)p_i^ρΞ_i(-p_i)p_i^σ+12ε^ρk^σ^2 𝒦_i(-p_i)p_i^ρp_i^σΞ_i(-p_i)

+2ie_iε^ρk^σℬ^i_ρσ(-p_i)Ξ_i(-p_i)+𝒪(kk)]×[ -𝒦_i(-p_i)p_i^μΞ_i(-p_i)-(k-ℓ)^κ𝒦_i(-p_i)p_i^μΞ_i(-p_i)p_i^κ

-12(2k-ℓ)^κ^2 𝒦_i(-p_i) p_i^μp_i^κΞ_i(-p_i)+2ie_iℓ^κℬ^i_μκ(-p_i)Ξ_i(-p_i)+𝒪(kk,kℓ,ℓℓ)] .
Now using the identities (<ref>),(<ref>) and (<ref>) we simplify the above expression and get
ζ_1μ   =   e_i^2 ϵ_i^T[2iε·p_i +12(ε^ρk^σ-ε^σk^ρ)𝒦_i(-p_i)p_i^ρΞ_i(-p_i)p_i^σ+ie_i(ε^ρk^σ-ε^σk^ρ) ℬ^i_ρσ(-p_i)Ξ_i(-p_i)

+ 𝒪(kk)]×[ -i(2p_i+2k-ℓ)_μ+𝒦_i(-p_i)Ξ_i(-p_i)p_i^μ+12(2k-ℓ)^κ𝒦_i(-p_i)^2 Ξ_i(-p_i)p_i^μp_i^κ

=   e_i^2 ϵ_i^T[ 2ε·p_i(2p_i+2k-ℓ)_μ+iε·p_i ℓ^ρ{𝒦_i(-p_i)p_i^μΞ_i(-p_i)p_i^ρ- 𝒦_i(-p_i)p_i^ρΞ_i(-p_i)p_i^μ

+4ie_iℬ^i_μρ(-p_i)Ξ_i(-p_i)}- i2p_iμ(ε^ρk^σ-ε^σk^ρ){𝒦_i(-p_i)p_i^ρΞ_i(-p_i)p_i^σ-𝒦_i(-p_i)p_i^σΞ_i(-p_i)p_i^ρ

+4ie_iℬ^i_ρσ(-p_i)Ξ_i(-p_i)}+2ip_i.k ε^ρ𝒦_i(-p_i)p_i^ρΞ_i(-p_i)p_i^μ +𝒪(kk,kℓ,ℓℓ)] .
Let us define a specific tensor structure which will appear together in all the computations in section-<ref> , which is also the non-universal contribution to the tree-level subleading soft photon theorem as derived in <cit.>
With the above definition, the expression of $\zeta_{1\mu}$ can be written in the following compact form
ζ_1μ   =    e_i^2 ϵ_i^T[ 2ε·p_i(2p_i+2k-ℓ)_μ-8ε·p_i ℓ^ρ𝒩^i_μρ(-p_i) +4 p_iμ(ε^ρk^σ-ε^σk^ρ) 𝒩^i_ρσ(-p_i)

+2i p_i.k ε^ρ𝒦_i(-p_i)p_i^ρΞ_i(-p_i)p_i^μ +𝒪(kk,kℓ,ℓℓ)] .
Let us evaluate the following expression appears in (<ref>) after substituting the vertex (<ref>) and using the identities (<ref>), (<ref>) and (<ref>),
ζ_2μ   ≡    ϵ_i^T(-p_i)ε^ρ(k)Γ^(4)_ρμ (p_i,-p_i-k+ℓ,k,-ℓ)Ξ_i(-p_i-k+ℓ)

=    ie_i^2 ϵ_i^T ε^ρ^2 𝒦_i(-p_i)p_i^μp_i^ρΞ_i(-p_i) +𝒪(ℓ, k) .
Using the identity in (<ref>) and on-shell condition (<ref>) the above expression reduces to
ζ_2μ   =    -e_i^2ϵ_i^T[2ε_μ+iε^ρ{𝒦_i(-p_i)p_i^μΞ_i(-p_i)p_i^ρ+𝒦_i(-p_i)p_i^ρΞ_i(-p_i)p_i^μ}+𝒪(ℓ, k)].
The expression below appears in (<ref>) and can be evaluated analogous to the evaluation of $\zeta_{1\mu}$. The final result reads
ζ_3μ   ≡   ϵ_i^T(-p_i)ε^ρ(k)Γ^(3)_μ(p_i,-p_i+ℓ,-ℓ)Ξ_i(-p_i+ℓ)Γ^(3)_ρ(p_i-ℓ,-p_i-k+ℓ,k)Ξ_i(-p_i-k+ℓ)

=    e_i^2 ϵ_i^T[ 2ε·p_i(2p_i-ℓ)_μ-4ε.ℓp_iμ-8ε·p_i ℓ^ρ𝒩^i_μρ(-p_i) +4 p_iμ(ε^ρk^σ-ε^σk^ρ) 𝒩^i_ρσ(-p_i)

-2i p_i.ℓ ε^ρ𝒦_i(-p_i)p_i^μΞ_i(-p_i)p_i^ρ +𝒪(kk,kℓ,ℓℓ)] .
Let us evaluate the following expression appears in (<ref>), (<ref>),(<ref>) and (<ref>). After substituting the vertex from (<ref>) and using the identities (<ref>), (<ref>) and (<ref>) we get
ζ_4ν   ≡    ϵ_j^T(-p_j)Γ^(3)_ν(p_j,-p_j-ℓ,ℓ)Ξ_j(-p_j-ℓ)

=    -ie_j ϵ_j^T[ 𝒦_j(-p_j)p_j^νΞ_j(-p_j)+ℓ^σ𝒦_j(-p_j)p_j^νΞ_j(-p_j)p_j^σ+12ℓ^σ^2 𝒦_j(-p_j)p_j^νp_j^σΞ_j(-p_j)

+2ie_jℓ^σℬ^j_νσ(-p_j)Ξ_j(-p_j)+𝒪(ℓℓ)] .
Using the identities (<ref>),(<ref>) and the on-shell condition (<ref>), the above expression reduces to
ζ_4ν   =   e_j ϵ_j^T[ (2p_j+ℓ)_ν+4ℓ^σ𝒩^j_νσ(-p_j)+𝒪(ℓℓ)] .

§ AMPUTATED GREEN'S FUNCTION INVOLVING SINGLE GRAVITON

Following the covariantization prescription described in section-<ref>, we compute the amputated Green's function involving $N$ number of massive spinning particles and one off-shell graviton, where the graviton is not attached to any external spinning particle leg. The resulting expression is given by
=   κ (2π)^4δ^(4)(p_1+⋯+p_N+ℓ)∑_i=1^N[δ_β_i^α_i h(ℓ)-δ_β_i^α_i  h^μν(ℓ)p_iμp_i^ν

+h^μν(ℓ)(Σ_iνb)_β_i^ α_i ℓ^b  p_i^μ -12δ_β_i^α_ih^μν(ℓ){ℓ^ρp_iμ^2p_i^ρp_i^ν+ℓ^ρp_iν^2p_i^ρp_i^μ

-p_i.ℓ^2p_i^μp_i^ν}+𝒪(ℓ^2) ]Γ^(N)α_1⋯α_i-1β_iα_i+1⋯α_N .
Above $\widehat{\Gamma}^{(N)}$ is defined after stripping out the momentum conserving delta function $(2\pi)^4\delta^{(4)}(p_1+\cdots+p_N)$ from the expression of $\Gamma^{(N)}$, i.e. $\Gamma^{(N)}\equiv (2\pi)^4\delta^{(4)}(p_1+\cdots+p_N)\widehat{\Gamma}^{(N)}$. Now starting from the above covariantized expression, the goal is to express $\widetilde{\Gamma}^{(N+1)}_{\mu\nu}(\ell)$ in terms of some operator operating on $\Gamma^{(N)}$ up to linear order in $\ell$. To do that we mostly follow the analysis of <cit.> with the only difference being $h^{\mu\nu}(\ell)$ is an off-shell graviton so we can not impose traceless or transverse condition. Instead in de Donder gauge, we use $\ell_\mu h^{\mu\nu}(\ell)=\f{1}{2}\ell^\nu h(\ell)$ in the intermediate stages of calculation.

Let us first analyze the first two terms within the square bracket in (<ref>) which is defined as
J_1    ≡    κ (2π)^4δ^(4)(p_1+⋯+p_N+ℓ)[ h(ℓ)- h^μν(ℓ)∑_i=1^N p_iμp_i^ν] Γ^(N)α_1⋯α_N

=   κ (2π)^4δ^(4)(p_1+⋯+p_N)[ h(ℓ)- h^μν(ℓ)∑_i=1^N p_iμp_i^ν] Γ^(N)α_1⋯α_N

+κ (2π)^4 ℓ^ρ{ℓ^ρδ^(4)(p_1+⋯+p_N+ℓ)}_ℓ=0 [ h(ℓ)- h^μν(ℓ)∑_i=1^N p_iμp_i^ν] Γ^(N)α_1⋯α_N

+𝒪(ℓ^2) .
Above we have just Taylor expanded the delta function and kept terms up to linear order in $\ell$. Now to evaluate the second line above we use the property involving momentum conserving delta function
∑_i=1^Np_iμ p_i^νδ^(4)(p_1+⋯+p_N) =∑_i=1^Np_iμ p_1^νδ^(4)(p_1+⋯+p_N)

=   p_1^ν∑_i=1^Np_iμ δ^(4)(p_1+⋯+p_N)-η_μν δ^(4)(p_1+⋯+p_N)

=    -η_μν δ^(4)(p_1+⋯+p_N) .
So using the above property and commuting the delta function through the momentum derivative the expression of $J_1$ becomes,
J_1    =    -κ (2π)^4 h^μν(ℓ)∑_i=1^N p_iμp_i^ν{δ^(4)(p_1+⋯+p_N) Γ^(N)α_1⋯α_N}

+κ (2π)^4 ℓ^ρ{ℓ^ρδ^(4)(p_1+⋯+p_N+ℓ)}_ℓ=0 [ h(ℓ)- h^μν(ℓ)∑_i=1^N p_iμp_i^ν] Γ^(N)α_1⋯α_N

+𝒪(ℓ^2) .
Last two terms within the square bracket in (<ref>) can be evaluated in the following way
J_2   ≡    κ (2π)^4δ^(4)(p_1+⋯+p_N+ℓ)h^μν(ℓ)∑_i=1^N[(Σ_iνb)_β_i^ α_i ℓ^b  p_i^μ-12δ_β_i^α_i{ℓ^ρp_iμ^2p_i^ρp_i^ν

=   κ (2π)^4 h^μν(ℓ)∑_i=1^N[(Σ_iνb)_β_i^ α_i ℓ^b  p_i^μ-12δ_β_i^α_i{ℓ^ρp_iμ^2p_i^ρp_i^ν+ℓ^ρp_iν^2p_i^ρp_i^μ-p_i.ℓ^2p_i^μp_i^ν}]

{δ^(4)(p_1+⋯+p_N+ℓ) Γ^(N)α_1⋯α_i-1β_iα_i+1⋯α_N} + ℒ ,
ℒ   =    -κ (2π)^4 Γ^(N)α_1⋯α_i-1β_iα_i+1⋯α_N h^μν(ℓ)∑_i=1^N[(Σ_iνb)_β_i^ α_i ℓ^b  p_i^μδ^(4)(p_1+⋯+p_N+ℓ)

+κ2(2π)^4 h^μν(ℓ)∑_i=1^N{p_i^ρδ^(4)(p_1+⋯+p_N+ℓ)}Γ^α_1⋯α_Np_i^σ{ℓ^ρp_iμδ_ν^σ+ℓ^σp_iμδ_ν^ρ

+ℓ^ρp_iνδ_μ^σ+ℓ^σp_iνδ_μ^ρ-p_i.ℓδ^ρ_μδ^σ_ν-p_i.ℓδ^σ_μδ^ρ_ν} .
To evaluate $\mathcal{L}$ we use the same trick as described earlier i.e. derivative w.r.t. $p_i$ on the delta function is same as derivative w.r.t. $\ell$ and then use the momentum conservation relation enforced by the delta function. After all these steps we get
ℒ   =   -κ (2π)^4 Γ^(N)α_1⋯α_i-1β_iα_i+1⋯α_Nh^μν(ℓ)[∑_i=1^N(Σ_iνb)_β_i^ α_i ℓ^b  ℓ^μδ^(4)(p_1+⋯+p_N+ℓ)

-p_i.ℓh^ρσ(ℓ)} .
To evaluate the first term within the square bracket above we use the conservation of total angular momenta follows from the Lorentz covariance of $\widehat{\Gamma}$
∑_i=1^N[(Σ_iνb)_β_i^ α_iΓ^(N)α_1⋯α_i-1β_iα_i+1⋯α_N-{p_iνp_i^b-p_i bp_i^ν}Γ^(N)α_1⋯α_N] =0 .
Using the above relation in the first line of the expression of $\mathcal{L}$ in (<ref>), operating the derivatives w.r.t. $\ell$ in the second term within the square bracket in (<ref>) and using de Donder gauge condition, we find
ℒ   =   -κ2  (2π)^4  Γ^α_1⋯α_N{h(ℓ)ℓ^νℓ^ρ^2ℓ^ρℓ^ν-ℓ^2 h^μν(ℓ)^2ℓ^μℓ^ν+2h(ℓ)ℓ^ρℓ^ρ}δ^(4)(p_1+⋯+p_N+ℓ)

+κ(2π)^4∑_i=1^N{ℓ^ρδ^(4)(p_1+⋯+p_N+ℓ)}Γ^α_1⋯α_Np_i^σ  p_iμℓ^ρh^μσ(ℓ) .
Now after substituting the expression of $\mathcal{L}$ from (<ref>) in (<ref>), then Taylor expanding the delta function and keeping terms up to linear in $\ell$ we get
J_2   =    κ (2π)^4 h^μν(ℓ)∑_i=1^N[(Σ_iνb)_β_i^ α_i ℓ^b  p_i^μ-12δ_β_i^α_i{ℓ^ρp_iμ^2p_i^ρp_i^ν+ℓ^ρp_iν^2p_i^ρp_i^μ-p_i.ℓ^2p_i^μp_i^ν}]

{δ^(4)(p_1+⋯+p_N) Γ^(N)α_1⋯α_i-1β_iα_i+1⋯α_N}

-κ (2π)^4  Γ^α_1⋯α_Nh(ℓ)ℓ^ρ{ℓ^ρ δ^(4)(p_1+⋯+p_N+ℓ)}_ℓ=0

+κ(2π)^4∑_i=1^N{ℓ^ρδ^(4)(p_1+⋯+p_N+ℓ)}_ℓ=0 Γ^α_1⋯α_Np_i^σ  p_iμℓ^ρh^μσ(ℓ) +𝒪(ℓ^2) .
Finally after adding the expressions in (<ref>) and (<ref>) we get

=    -κ (2π)^4 h^μν(ℓ)∑_i=1^N p_iμp_i^ν{δ^(4)(p_1+⋯+p_N) Γ^(N)α_1⋯α_N}

+κ (2π)^4 h^μν(ℓ)∑_i=1^N[(Σ_iνb)_β_i^ α_i ℓ^b  p_i^μ-12δ_β_i^α_i{ℓ^ρp_iμ^2p_i^ρp_i^ν+ℓ^ρp_iν^2p_i^ρp_i^μ-p_i.ℓ^2p_i^μp_i^ν}]

{δ^(4)(p_1+⋯+p_N) Γ^(N)α_1⋯α_i-1β_iα_i+1⋯α_N} +𝒪(ℓ^2) .

§ INTERMEDIATE STEPS IN DERIVING SOFT GRAVITON THEOREM

We want to compute the following expression involving the vertex in (<ref>) which takes the following form once we expanded in power of small $\ell$ and keep terms up to quadratic order

=    iκ[ η_μν𝒦(-q)Ξ(-q)+η_μν𝒦(-q)ℓ^ρΞ(-q)q^ρ+12η_μν𝒦(-q)ℓ^ρℓ^σ^2Ξ(-q)q^ρq^σ

+12η_μνℓ^ρ𝒦(-q)q^ρΞ(-q)+12η_μνℓ^ρℓ^σ𝒦(-q)q^ρΞ(-q)q^σ+14η_μνℓ^ρℓ^σ^2 𝒦(-q)q^ρq^σΞ(-q)

- q_(μ 𝒦(-q)q^ν)Ξ(-q)- q_(μℓ^ρ𝒦(-q)q^ν)Ξ(-q)q^ρ- 12q_(μℓ^ρℓ^σ𝒦(-q)q^ν)^2Ξ(-q)q^ρq^σ

-12ℓ_(μ 𝒦(-q)q^ν)Ξ(-q)-12ℓ_(μℓ^ρ𝒦(-q)q^ν)Ξ(-q)q^ρ-12q_(μℓ^ρ^2 𝒦(-q)q^ν) q^ρΞ(-q)

-12q_(μℓ^ρℓ^σ^2 𝒦(-q)q^ν) q^ρΞ(-q)q^σ -14q_(μℓ^ρℓ^σ^3 𝒦(-q)q^ν) q^ρq^σΞ(-q)

-12ℓ_(μ ℓ^ρ^2 𝒦(-q)q^ν) q^ρΞ(-q)-12ℓ^b𝒦(-q)q^(μΣ_ν)bΞ(-q)-12ℓ^bℓ^ρ𝒦(-q)q^(μΣ_ν)bΞ(-q)q^ρ

+12ℓ^b  Σ^T_(νb 𝒦(-q)q^μ)Ξ(-q)+12ℓ^b ℓ^ρ Σ^T_(νb 𝒦(-q)q^μ)Ξ(-q)q^ρ -14ℓ^b ℓ^ρ ^2 𝒦(-q)q^(μ q^ρΣ_ν) bΞ(-q)

+14ℓ^b ℓ^ρΣ^T_(νb^2 𝒦(-q)q^μ) q^ρΞ(-q)+14{ℓ_μ ℓ^ρ^2 𝒦(-q)q^ρq^ν+ℓ_ν ℓ^ρ^2 𝒦(-q)q^ρq^μ-ℓ^2   ^2 𝒦(-q)q^μq^ν}Ξ(-q)

+112ℓ^σ{q_μ ℓ^ρ^3 𝒦(-q)q^ρq^ν q^σ+q_ν ℓ^ρ^3 𝒦(-q)q^ρq^μ q^σ-q.ℓ^3 𝒦(-q)q^μq^νq^σ}Ξ(-q)

-ℓ^ρℓ^σ{𝒢_(μρσν)(-q)+𝒢_σ(νμ)ρ(-q)-𝒢_(μρν)σ(-q) -𝒢_σ(νρμ)(-q)}Ξ(-q)+𝒪(ℓ^3)].
Now to simplify the above expression we need to move the momenta derivatives from $\mathcal{K}$ to $\Xi$ as much possible and also move the spin operator to the extreme right using the identities (<ref>)-(<ref>). Following these steps the result up to order $\mathcal{O}(\ell^2)$ reads

=    iκ[ iη_μν(q^2+m^2+q.ℓ)-2iq_μq_ν-2iq_(μℓ_ν) +2iℓ^b q_(μΣ^T_ν)b+iℓ^b ℓ_(μΣ^T_ν)b

+ℓ_(μ 𝒦(-q)Ξ(-q)q^ν)-12ℓ^b 𝒦(-q)Ξ(-q)q^(μΣ^T_ν)b -14ℓ^ρℓ^b 𝒦(-q)^2 Ξ(-q)q^(μq^ρΣ^T_ν)b

-12ℓ^b   𝒦(-q)Ξ(-q)q^(μΣ^T_ν) b  +14ℓ^b ℓ^ρ𝒦(-q)Σ_(νb^2 Ξ(-q)q^μ) q^ρ+14ℓ^2 𝒦(-q)^2 Ξ(-q)q^μq^ν

+112ℓ^σ{q_(μ ℓ^ρ𝒦(-q) ^3 Ξ(-q)q^ρq^ν) q^σ+q.ℓ𝒦(-q)  ^3 Ξ(-q)q^μq^νq^σ}

+13q.ℓℓ^ρ𝒦(-q)q^ρ ^2 Ξ(-q)q^μq^ν +13ℓ^b ℓ^ρq_(μ𝒦(-q)q^ν)^2Ξ(-q)q^ρq^b-13q.ℓℓ^ρ𝒦(-q)q^(μ^2Ξ(-q)q^ν)q^ρ

-13ℓ^b ℓ^ρq_(μ𝒦(-q)q^ρ^2 Ξ(-q)q^ν)q^b   +16q.ℓℓ^ρ^2 𝒦(-q)q^(μq^ρΞ(-q)q^ν)  -16ℓ^bℓ^ρq_(μ^2 𝒦(-q)q^ρq^bΞ(-q)q^ν)

+16q_(μℓ^ρℓ^σ^2𝒦(-q)q^ρq^ν)Ξ(-q)q^σ-16q.ℓℓ^ρ^2 𝒦(-q)q^μq^νΞ(-q)q^ρ

-12ℓ^b ℓ^ρ  𝒦(-q)q^ρΞ(-q)q^(μΣ^T_ν) b +12ℓ^bℓ^ρ𝒦(-q)q^(μΞ(-q)q^ρΣ^T_ν)b

-ℓ^ρℓ^σ{𝒢_(μρσν)(-q)+𝒢_σ(νμ)ρ(-q)-𝒢_(μρν)σ(-q) -𝒢_σ(νρμ)(-q)}Ξ(-q)+𝒪(ℓ^3)].
The symmetrization in the exchange between $\mu$ and $\nu$ can be omitted since any contraction involving the above expression in any Feynman diagram calculation will always exhibit symmetry under the exchange of $\mu$ and $\nu$. Also the above expression can be written in a compact way by introducing the following tensor structure
Δ_μρνσ(-q)   ≡    13q_μ𝒦(-q)q^ν^2Ξ(-q)q^ρq^σ-16q_ρ^2𝒦(-q)q^μq^νΞ(-q)q^σ+14𝒦(-q)q^μΞ(-q)q^ρΣ_νσ^T

+ 𝒢_μρνσ(-q)Ξ(-q) .
Hence removing the $\mu\leftrightarrow\nu$ symmetrization, the expression in (<ref>) can be compactly written as

=    iκ[ iη_μν(q^2+m^2+q.ℓ)-2iq_μq_ν-2iq_μℓ_ν +2iℓ^b q_μΣ^T_νb+iℓ^b ℓ_μΣ^T_νb

+q_μ𝒦(-q)Ξ(-q)q^ν+14η_μν𝒦(-q)ℓ^ρℓ^σ^2Ξ(-q)q^ρq^σ +12q.ℓ𝒦(-q)^2Ξ(-q)q^μq^ν

+ℓ_μ 𝒦(-q)Ξ(-q)q^ν-ℓ^b 𝒦(-q)Ξ(-q)q^μΣ^T_νb -14ℓ^ρℓ^b 𝒦(-q)^2 Ξ(-q)q^μq^ρΣ^T_νb

+14ℓ^b ℓ^ρ𝒦(-q)Σ_νb^2 Ξ(-q)q^μ q^ρ+14ℓ^2 𝒦(-q)^2 Ξ(-q)q^μq^ν

+112ℓ^σ{q_μ ℓ^ρ𝒦(-q) ^3 Ξ(-q)q^ρq^ν q^σ+q.ℓ𝒦(-q)  ^3 Ξ(-q)q^μq^νq^σ}

Now we want to compute the following expression involving two vertices of kind (<ref>) which takes the following form once we expanded in power of small $\ell_1,\ell_2$ and keep terms up to quadratic order

=    -κ^2 [ iη_μν(q^2+m^2+q.ℓ_1)-2iq_μq_ν-2iq_μℓ_1ν +2iℓ_1^b q_μΣ^T_νb+iℓ_1^b ℓ_1μΣ^T_νb

+q_μ𝒦(-q)Ξ(-q)q^ν+14η_μν𝒦(-q)ℓ_1^λℓ_1^b ^2Ξ(-q)q^λq^b +12q.ℓ_1 𝒦(-q)^2Ξ(-q)q^μq^ν

+ℓ_1μ 𝒦(-q)Ξ(-q)q^ν-ℓ_1^b 𝒦(-q)Ξ(-q)q^μΣ^T_νb -14ℓ_1^λℓ_1^b 𝒦(-q)^2 Ξ(-q)q^μq^λΣ^T_νb

+14ℓ_1^b ℓ_1^λ𝒦(-q)Σ_νb^2 Ξ(-q)q^μ q^λ+14ℓ_1^2 𝒦(-q)^2 Ξ(-q)q^μq^ν

+112ℓ_1^b{q_μ ℓ_1^λ𝒦(-q) ^3 Ξ(-q)q^λq^ν q^b+q.ℓ_1 𝒦(-q)  ^3 Ξ(-q)q^μq^νq^b}

+ℓ_1^λℓ_1^b {Δ_μλνb(-q)+Δ_λμbν(-q)-Δ_λμνb(-q)-Δ_μλbν(-q)}]

×[ iη_ρσ(q^2+m^2+2q.ℓ_1+q.ℓ_2+ℓ_1^2+ℓ_1.ℓ_2)-2i(q+ℓ_1)_ρ(q+ℓ_1)_σ-2i(q+ℓ_1)_ρℓ_2σ

+2iℓ_2^a (q+ℓ_1)_ρΣ^T_σa+iℓ_2^a ℓ_2ρΣ^T_σa

+(q+ℓ_1)_ρ𝒦(-q)Ξ(-q)q^σ+(q+ℓ_1)_ρℓ_1^τ𝒦(-q)q^τΞ(-q)q^σ+(q+ℓ_1)_ρℓ_1^τ𝒦(-q) ^2Ξ(-q)q^σq^τ

+12q_ρℓ_1^τℓ_1^κ{^2 𝒦(-q)q^τq^κΞ(-q)q^σ+2𝒦(-q)q^τ^2Ξ(-q)q^σq^κ+ 𝒦(-q)^3Ξ(-q)q^σq^κq^τ}

+14η_ρσ𝒦(-q)ℓ_2^κℓ_2^τ^2Ξ(-q)q^κq^τ +12(q+ℓ_1).ℓ_2 𝒦(-q)^2Ξ(-q)q^ρq^σ

+12 q.ℓ_2 ℓ_1^τ𝒦(-q)q^τ^2Ξ(-q)q^ρq^σ+12q.ℓ_2 ℓ_1^τ𝒦(-q)^3Ξ(-q)q^ρq^σq^τ

+ℓ_2ρ 𝒦(-q)Ξ(-q)q^σ+ℓ_2ρ ℓ_1^τ𝒦(-q)q^τΞ(-q)q^σ+ℓ_2ρℓ_1^τ𝒦(-q)^2 Ξ(-q)q^σ q^τ

-ℓ_2^a 𝒦(-q)Ξ(-q)q^ρΣ^T_σa-ℓ_2^a ℓ_1^τ𝒦(-q)q^τΞ(-q)q^ρΣ^T_σa-ℓ_2^a ℓ_1^τ𝒦(-q)^2Ξ(-q)q^ρq^τΣ^T_σa

-14ℓ_2^τℓ_2^a 𝒦(-q)^2 Ξ(-q)q^ρq^τΣ^T_σa  +14ℓ_2^a ℓ_2^τ𝒦(-q)Σ_σa^2 Ξ(-q)q^ρ q^τ+14ℓ_2^2 𝒦(-q)^2 Ξ(-q)q^ρq^σ

+112ℓ_2^τ{q_ρ ℓ_2^κ𝒦(-q) ^3 Ξ(-q)q^κq^σ q^τ+q.ℓ_2 𝒦(-q)  ^3 Ξ(-q)q^ρq^σq^τ}

+ℓ_2^κℓ_2^τ{Δ_ρκστ(-q)+Δ_κρτσ(-q)-Δ_κρστ(-q)-Δ_ρκτσ(-q)}] .
In principle the above expression can be evaluated using the identities (<ref>)-(<ref>), but it is tedious and we don't need the full contribution. Instead we only evaluate the above expression by contracting $\epsilon(-q)^{T}$ from the left considering the particle with momentum $q$ being on-shell i.e. $q^2+m^2=0$ and $\epsilon(-q)^{T}\mathcal{K}(-q)=0$. Implementing these conditions and using the identities (<ref>)-(<ref>), the above expression simplifies to

≡    ϵ(-q)^TΓ^(3)_μν(q,-q-ℓ_1,ℓ_1)Ξ(-q-ℓ_1)Γ^(3)_ρσ(q+ℓ_1,-q-ℓ_1-ℓ_2,ℓ_2)Ξ(-q-ℓ_1-ℓ_2)

=     -κ^2 ϵ^T[ -η_μνη_ρσq.ℓ_1(2q.ℓ_1+q.ℓ_2)+2q_μq_νη_ρσ(2q.ℓ_1+q.ℓ_2+ℓ_1^2+ℓ_1.ℓ_2)

+2q_μℓ_1νη_ρσ(2q.ℓ_1+q.ℓ_2)-2η_ρσq_μℓ_1^bΣ_νb^T(2q.ℓ_1+q.ℓ_2)+2η_μνq.ℓ_1 (q_ρq_σ+q_ρℓ_1σ+q_σℓ_1ρ)

-4q_μq_ν(q+ℓ_1)_ρ(q+ℓ_1)_σ-4q_μℓ_1ν(q_ρq_σ+q_ρℓ_1σ+q_σℓ_1ρ)+4q_μℓ_1^b Σ_νb^T(q_ρq_σ+q_ρℓ_1σ+q_σℓ_1ρ)

+2ℓ_1^b ℓ_1μΣ_νb^T q_ρq_σ+2η_μνq.ℓ_1 q_ρℓ_2σ-4q_μq_ν(q+ℓ_1)_ρℓ_2σ-4q_μℓ_1νq_ρℓ_2σ+4q_μℓ_1^b Σ_νb^T q_ρℓ_2σ

-2η_μνq.ℓ_1 ℓ_2^a q_ρΣ_σa^T +4q_μq_νℓ_2^a (q+ℓ_1)_ρΣ_σa^T +4q_μℓ_1νℓ_2^a q_ρΣ_σa^T -4q_μℓ_1^b Σ_νb^T ℓ_2^a q_ρΣ_σa^T

+2q_μq_νℓ_2^a ℓ_2ρΣ_σa^T -2iq.ℓ_1 q_μq_ρ𝒦(-q)q^νΞ(-q)q^σ

-2iq_ρq_σℓ_1^λℓ_1^b {Δ_μλνb(-q)+Δ_λμbν(-q)-Δ_λμνb(-q)-Δ_μλbν(-q)}

-2i q_μq_νℓ_2^κℓ_2^τ{Δ_ρκστ(-q)+Δ_κρτσ(-q)-Δ_κρστ(-q)-Δ_ρκτσ(-q)}

+ℓ_1^λℓ_1^b {-16q_μq_ρ𝒦(-q)q^νΞ(-q)q^λ𝒦(-q)q^bΞ(-q)q^σ-16q_μq_ρ𝒦(-q)q^λΞ(-q)q^ν𝒦(-q)q^bΞ(-q)q^σ

-iℓ_1νℓ_1^b q_μq_ρ𝒦(-q)q^bΞ(-q)q^σ+2i η_νσq_ρq.ℓ_1 ℓ_1^b 𝒦(-q)q^μΞ(-q)q^b-2iℓ_1σq_ρq.ℓ_1 𝒦(-q)q^μΞ(-q)q^ν

-i(q.ℓ_1)^2 q_ρ^2 𝒦(-q)q^μq^νΞ(-q)q^σ +iq.ℓ_1 ℓ_1^b q_νq_ρ^2 𝒦(-q)q^μq^bΞ(-q)q^σ+2i q_ρq.ℓ_1 ℓ_1^b 𝒦(-q)q^μΞ(-q)q^σΣ_νb^T

+2i q_ρq.ℓ_1 ℓ_1^b q_ν𝒦(-q)q^μ^2 Ξ(-q)q^b q^σ-2i q_ρ(q.ℓ_1)^2  𝒦(-q)q^μ^2 Ξ(-q)q^νq^σ

-2iq.ℓ_1 q_μℓ_1ρ  𝒦(-q)q^νΞ(-q)q^σ+iℓ_1^b ℓ_1μ q_ρq_ν𝒦(-q)q^bΞ(-q)q^σ-2iq.ℓ_1 ℓ_1μ q_ρ𝒦(-q)q^νΞ(-q)q^σ

+2iη_μνq.ℓ_1 q_ρℓ_1^τ𝒦(-q)q^τΞ(-q)q^σ-iq.ℓ_1 q_μq_ρℓ_1^τ^2 𝒦(-q)q^νq^τΞ(-q)q^σ -iℓ_1^2 q_μq_ρ𝒦(-q)q^νΞ(-q)q^σ

-2iq.ℓ_1 q_μq_ρℓ_1^τ𝒦(-q)q^ν^2Ξ(-q)q^σq^τ+iℓ_1^b q_μ q.ℓ_2 q_ν𝒦(-q)q^b^2Ξ(-q)q^ρq^σ-iq.ℓ_1 q_μ q.ℓ_2  𝒦(-q)q^ν^2Ξ(-q)q^ρq^σ

-iq_μq_νq.ℓ_2 ℓ_1^τ𝒦(-q)q^τ^2Ξ(-q)q^ρq^σ+2iℓ_1^b q_μℓ_2ρ q_ν𝒦(-q)q^bΞ(-q)q^σ-2iq.ℓ_1 q_μℓ_2ρ  𝒦(-q)q^νΞ(-q)q^σ

-2iq_μq_νℓ_2ρ ℓ_1^τ𝒦(-q)q^τΞ(-q)q^σ-2iℓ_1^b q_μℓ_2^a q_ν𝒦(-q)q^bΞ(-q)q^ρΣ^T_σa+2iq.ℓ_1 q_μℓ_2^a  𝒦(-q)q^νΞ(-q)q^ρΣ^T_σa

+2iq_μq_νℓ_2^a ℓ_1^τ𝒦(-q)q^τΞ(-q)q^ρΣ^T_σa+𝒪(ℓ_1,2^3)] .
To write the above expression in a compact form we used our compact structure $\Delta$ defined in (<ref>).

We also need to compute the following expression involving the vertex in (<ref>) which takes the following form once we expanded in power of small $\ell_1,\ell_2$ and keep terms up to linear order

=   iκ^2 [ (η_μνη_ρσ-2η_μρη_νσ){𝒦(-q)Ξ(-q)+(ℓ_1+ℓ_2)^κ𝒦(-q)Ξ(-q)q^κ+12(ℓ_1+ℓ_2)^κ𝒦(-q)q^κΞ(-q)}

+12q_ρ(ℓ_1+ℓ_2)^κ^2 𝒦(-q)q^σq^κΞ(-q)+12ℓ_2^b 𝒦(-q)q^ρΣ_σbΞ(-q) -12ℓ_2^b Σ^T_σb𝒦(-q)q^ρΞ(-q)}

+12q_μ(ℓ_1+ℓ_2)^κ^2 𝒦(-q)q^νq^κΞ(-q)+12ℓ_1^b 𝒦(-q)q^μΣ_νbΞ(-q) -12ℓ_1^b Σ^T_νb𝒦(-q)q^μΞ(-q)}

+12q_σ(ℓ_1+ℓ_2)^κ^2 𝒦(-q)q^κq^νΞ(-q)+q_ν𝒦(-q)q^σΞ(-q)+q_ν(ℓ_1+ℓ_2)^κ𝒦(-q)q^σΞ(-q)q^κ

+12(ℓ_1+ℓ_2)_ν𝒦(-q)q^σΞ(-q)+12q_ν(ℓ_1+ℓ_2)^κ^2 𝒦(-q)q^κq^σΞ(-q)}+q_μq_ρ^2 𝒦(-q)q^νq^σΞ(-q)

+q_μq_ρ(ℓ_1+ℓ_2)^κ^2 𝒦(-q)q^νq^σΞ(-q)q^κ+12{q_μ(ℓ_1 +ℓ_2)_ρ+q_ρ(ℓ_1+ℓ_2)_μ}^2 𝒦(-q)q^νq^σΞ(-q)

+12q_μq_ρ(ℓ_1 +ℓ_2)^κ^3 𝒦(-q)q^νq^σq^κΞ(-q) -14η_νσ(ℓ_2-ℓ_1)^κ{𝒦(-q)q^κΣ_ρμ-Σ_ρμ^T𝒦(-q)q^κ}Ξ(-q)

+12η_μρ{ℓ_2^b 𝒦(-q)q^νΣ_σb-ℓ_2^b Σ^T_σb𝒦(-q)q^ν +ℓ_1^b 𝒦(-q)q^σΣ_νb -ℓ_1^b Σ^T_νb𝒦(-q)q^σ}Ξ(-q)

+14q_μℓ_2^b {^2 𝒦(-q)q^νq^ρΣ_σb -Σ_σb^T ^2 𝒦(-q)q^νq^ρ}Ξ(-q) +14q_ρℓ_1^b {^2 𝒦(-q)q^σq^μΣ_νb-Σ_νb^T ^2 𝒦(-q)q^σq^μ}Ξ(-q)

-12η_νσℓ_2^a {𝒦(-q)q^ρΣ_aμ -Σ_aμ^T 𝒦(-q)q^ρ}Ξ(-q) -12η_νσℓ_1^a {𝒦(-q)q^μΣ_aρ -Σ_aρ^T 𝒦(-q)q^μ}Ξ(-q)

-12ℓ_2ν{𝒦(-q)q^σΣ_μρ-Σ_μρ^T 𝒦(-q)q^σ}Ξ(-q) -12ℓ_1σ{𝒦(-q)q^νΣ_ρμ-Σ_ρμ^T 𝒦(-q)q^ν}Ξ(-q)

+𝒪(ℓ_1^2, ℓ_1ℓ_2,ℓ_2^2)] .
Now we contract the above expression with $\epsilon(-q)^{T}$ from the left considering the particle with momentum $q$ being on-shell i.e. $q^2+m^2=0$ and $\epsilon(-q)^{T}\mathcal{K}(-q)=0$. Then using the identities (<ref>)-(<ref>), the above expression up to order $\mathcal{O}(\ell_1,\ell_2)$ simplifies to

≡    ϵ(-q)^T Γ^(4)_μν,ρσ(q,-q-ℓ_1-ℓ_2,ℓ_1,ℓ_2)Ξ(-q-ℓ_1-ℓ_2)

=   iκ^2 ϵ^T[ (η_μνη_ρσ-2η_μρη_νσ){i(q^2+m^2)+iq.(ℓ_1 +ℓ_2)}-η_μν{2iq_ρq_σ+i(ℓ_1+ℓ_2)_ρq_σ

+i q_ρ(ℓ_1+ℓ_2)_σ}-η_ρσ {2iq_μq_ν+i(ℓ_1+ℓ_2)_μq_ν+i q_μ(ℓ_1+ℓ_2)_ν}+3η_μρ{2iq_νq_σ

+iq_ν(ℓ_1+ℓ_2)_σ+iq_σ(ℓ_1+ℓ_2)_ν}+η_νσ{2iq_μq_ν+i q_μ(ℓ_1 +ℓ_2)_ρ+iq_ρ(ℓ_1+ℓ_2)_μ}

-iη_νρq_μℓ_2^b Σ_σb^T-iη_σμq_ρℓ_1^b Σ_νb^T +2iη_μν ℓ_2^b q_ρΣ^T_σb+2iη_ρσℓ_1^b q_μΣ^T_νb

+2iη_νσq_ρℓ_2^a Σ_aμ^T +2iη_νσq_μℓ_1^a Σ_aρ^T +2iℓ_2νq_σΣ_μρ^T +2iℓ_1σq_νΣ_ρμ^T

-2iη_μρq_νℓ_2^b Σ_σb^T -2iη_μρq_σℓ_1^b Σ_νb^T +iη_νσ q.(ℓ_2-ℓ_1)Σ^T_ρμ

-η_μν{12q_ρℓ_1^κ𝒦(-q)q^σΞ(-q)q^κ-12q_ρℓ_1^κ𝒦(-q)q^κΞ(-q)q^σ }

-η_ρσ{12q_μℓ_2^κ𝒦(-q)q^νΞ(-q)q^κ-12q_μℓ_2^κ𝒦(-q)q^κΞ(-q)q^ν }

+32η_μρ{-12q_σ(ℓ_1+ℓ_2)^κ𝒦(-q)q^κΞ(-q)q^ν  +12q_σ(ℓ_1+ℓ_2)^κ𝒦(-q)q^νΞ(-q)q^κ

-12q_ν(ℓ_1+ℓ_2)^κ𝒦(-q)q^κΞ(-q)q^σ  +12q_ν(ℓ_1+ℓ_2)^κ𝒦(-q)q^σΞ(-q)q^κ}

-q_μq_ρ𝒦(-q)q^νΞ(-q)q^σ-q_μq_ρ𝒦(-q)q^σΞ(-q)q^ν  +q_μq_ρ(ℓ_1+ℓ_2)^κ^2 𝒦(-q)q^νq^σΞ(-q)q^κ

-12{q_μ(ℓ_1 +ℓ_2)_ρ+q_ρ(ℓ_1+ℓ_2)_μ}{𝒦(-q)q^νΞ(-q)q^σ+ 𝒦(-q)q^σΞ(-q)q^ν}

-12q_μq_ρ(ℓ_1 +ℓ_2)^κ{^2 𝒦(-q)q^νq^σΞ(-q)q^κ+^2 𝒦(-q)q^νq^κΞ(-q)q^σ+^2 𝒦(-q)q^κq^σΞ(-q)q^ν

+ 𝒦(-q)q^ν^2 Ξ(-q)q^κq^σ+𝒦(-q) q^κ^2 Ξ(-q)q^σq^ν+𝒦(-q) q^σ^2 Ξ(-q)q^νq^κ}

+12η_μρℓ_2^b{-q_σ𝒦(-q)q^νΞ(-q)q^b+q_σ𝒦(-q)q^bΞ(-q)q^ν+q_b 𝒦(-q)q^νΞ(-q)q^σ-q_b 𝒦(-q)q^σΞ(-q)q^ν}

+12η_μρℓ_1^b {-q_ν𝒦(-q)q^σΞ(-q)q^b+q_ν𝒦(-q)q^bΞ(-q)q^σ+q_b 𝒦(-q)q^σΞ(-q)q^ν-q_b 𝒦(-q)q^νΞ(-q)q^σ }

+14q_μℓ_2^b {2𝒦(-q)q^νΞ(-q)q^ρΣ_σb^T+2𝒦(-q)q^ρΞ(-q)q^νΣ_σb^T -q_σ^2 𝒦(-q)q^νq^ρΞ(-q)q^b

+q_σ^2 𝒦(-q)q^νq^bΞ(-q)q^ρ+q_σ^2 𝒦(-q)q^b q^ρΞ(-q)q^ν-q_b ^2 𝒦(-q)q^σq^ρΞ(-q)q^ν+η_σν𝒦(-q)q^bΞ(-q)q^ρ

-η_bν𝒦(-q)q^σΞ(-q)q^ρ +η_σρ𝒦(-q)q^bΞ(-q)q^ν-η_bρ𝒦(-q)q^σΞ(-q)q^ν+η_σρ𝒦(-q)q^νΞ(-q)q^b

-η_bρ𝒦(-q)q^νΞ(-q)q^σ+η_σν𝒦(-q)q^ρΞ(-q)q^b -η_bν𝒦(-q)q^ρΞ(-q)q^σ+q_σ𝒦(-q)q^b^2Ξ(-q)q^νq^ρ

-q_b 𝒦(-q)q^σ^2Ξ(-q)q^νq^ρ  +q_σ𝒦(-q)q^ν^2Ξ(-q)q^b q^ρ -q_b 𝒦(-q)q^ν^2Ξ(-q)q^σq^ρ +q_σ𝒦(-q)q^ρ^2Ξ(-q)q^νq^b

-q_b 𝒦(-q)q^ρ^2Ξ(-q)q^νq^σ}

+14q_ρℓ_1^b {2𝒦(-q)q^σΞ(-q)q^μΣ_νb^T+2𝒦(-q)q^μΞ(-q)q^σΣ_νb^T -q_ν^2 𝒦(-q)q^σq^μΞ(-q)q^b

+q_ν^2 𝒦(-q)q^σq^bΞ(-q)q^μ+q_ν^2 𝒦(-q)q^b q^μΞ(-q)q^σ-q_b ^2 𝒦(-q)q^νq^μΞ(-q)q^σ+η_σν𝒦(-q)q^bΞ(-q)q^μ

-η_bσ𝒦(-q)q^νΞ(-q)q^μ +η_νμ𝒦(-q)q^bΞ(-q)q^σ-η_bμ𝒦(-q)q^νΞ(-q)q^σ   +η_νμ𝒦(-q)q^σΞ(-q)q^b

-η_bμ𝒦(-q)q^σΞ(-q)q^ν+η_σν𝒦(-q)q^μΞ(-q)q^b  -η_bσ𝒦(-q)q^μΞ(-q)q^ν+q_ν𝒦(-q)q^b^2Ξ(-q)q^σq^μ

-q_b 𝒦(-q)q^ν^2Ξ(-q)q^σq^μ  +q_ν𝒦(-q)q^σ^2Ξ(-q)q^b q^μ -q_b 𝒦(-q)q^σ^2Ξ(-q)q^νq^μ +q_ν𝒦(-q)q^μ^2Ξ(-q)q^σq^b

-q_b 𝒦(-q)q^μ^2Ξ(-q)q^σq^ν}

-14η_νσ(ℓ_2+ℓ_1)^κ{-q_ρ𝒦(-q)q^κΞ(-q)q^μ+q_ρ𝒦(-q)q^μΞ(-q)q^κ -q_μ𝒦(-q)q^κΞ(-q)q^ρ

+q_μ𝒦(-q)q^ρΞ(-q)q^κ}-12η_νσq.(ℓ_1-ℓ_2) {- 𝒦(-q)q^μΞ(-q)q^ρ+ 𝒦(-q)q^ρΞ(-q)q^μ }

-12ℓ_2ν{q_ρ𝒦(-q)q^σΞ(-q)q^μ-q_ρ𝒦(-q)q^μΞ(-q)q^σ }-12ℓ_1σ{q_μ𝒦(-q)q^νΞ(-q)q^ρ

-q_μ𝒦(-q)q^ρΞ(-q)q^ν }+𝒪(ℓ_1^2, ℓ_1ℓ_2,ℓ_2^2)] .

Above we removed the terms which are anti-symmetric under $\mu\leftrightarrow\nu$ and/or $\rho\leftrightarrow\sigma$ exchanges.