Modulation instability in dispersive parity-broken systems

Sudheesh Srivastava1,2, Gustavo M. Monteiro 3, Sriram Ganeshan1,2, *

1 Department of Physics, City College, City University of New York, New York,
NY 10031, USA

2 CUNY Graduate Center, New York, NY 10031

3 Department of Physics and Astronomy, College of Staten Island, CUNY, Staten
Island, NY 10314, USA

<EMAIL_ADDRESS>
## Abstract

This work explores the interplay between dispersive parity breaking and non-
linearity in two contrasting continuous dynamical systems that exhibit
Modulation Instability (MI). We begin by examining deep water odd surface
gravity waves and derive the non-linear Schrödinger equation (NLSE) for the
modulated envelope dynamics using the method of multiple scales. The parity
breaking in the odd gravity waves results in distinct NLSEs for the right and
the left mover, leading to chirality-dependent stability properties for the
envelope dynamics. Moreover, the resonant interaction of gravity waves and odd
viscosity-induced capillary dynamics creates a window of wave numbers in one
of the chiral sectors where the envelope propagation remains stable. Following
the odd gravity results, we design a one-dimensional non-reciprocal PT-
symmetric dielectric model that exhibits parity-breaking effects analogous to
the odd viscosity term in 2D hydrodynamics. With cubic non-linearity in the
polarization dynamics, we derive the corresponding NLSE. Once again, parity
breaking stabilizes the modulated envelope dynamics in the lower polariton
bands. We then compare the similarities and differences between this case and
that of odd gravity waves.

###### Contents

  1. 1 Introduction
  2. 2 Review of a non-linear potential theory for odd gravity waves
    1. 2.1  Irrotational limit of the free surface dynamics and boundary layer approximation
  3. 3 Method of multiple scales for envelope dynamics and non-Linear Schrödinger equation
    1. 3.1 Linear order
    2. 3.2 Higher order non-linear terms and envelope dynamics
    3. 3.3 Stability of the NLSE and soliton formation
  4. 4 Nonreciprocal 1D dielectric
    1. 4.1 Dispersive parity breaking in a dielectric
    2. 4.2 Heuristics of the Chiral Term
    3. 4.3 Kerr Nonlinearity
  5. 5 Discussion
  6. 6 Acknowledgements
  7. 7 Appendix
    1. 7.1 Details for envelope solutions in parity-odd fluid
    2. 7.2 Envelope dynamics in the parity-odd dielectric medium

##  1 Introduction

Parity breaking is a ubiquitous phenomenon that occurs at a microscopic scale
in both synthetic and natural entities with handedness. How these microscale
chiral effects translate to the collective dynamics at the macro and mesoscale
systems has been an active area of research across several disciplines [1, 2,
3, 4].

One approach to studying parity-breaking effects at a macroscopic scale is to
write down classical field theories such as hydrodynamics with symmetry-
allowed terms in the constitutive equations encoding the chiral effects [5,
4]. The macroscopic approach based on symmetry has allowed the identification
of equivalent physical effects across different systems. For example, odd
viscosity [6, 7] is a unique term in two-dimensional fluids that preserves
isotropy while breaking parity without causing dissipation. In certain broken
parity two-dimensional quantum fluids, it has been investigated extensively as
Hall viscosity [8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23].
Subsequently, this concept has been extended to various fields such as
elasticity theory [24], acoustics [25], and active matter[26, 27, 3, 28, 29],
with each system exhibiting distinct physics related to this term. In
elasticity, for instance, energy must be locally supplied to achieve such
parity-odd responses. The diverse class of systems where odd viscosity
manifests have been extensively reviewed in Ref. [4] and the references
therein.

Despite these advances, the cross-fertilization of ideas has mostly been
confined to the linear-response regime, with some works investigating weakly
non-linear effects in incompressible fluid systems with odd viscosity [30, 31,
32, 33]. However, certain non-linear effects are also ubiquitous across
various platforms and disciplines. One such phenomenon is modulation
instability (MI), a fundamental concept in nonlinear wave theory, with
applications across many branches of physics and engineering. MI arises due to
the interplay between wave dispersion and nonlinear interactions leading to
soliton formation, sideband instabilities, and in some cases Fermi-Pasta-Ulam
recurrences [34].

Historically, modulation instability (MI) was discovered almost simultaneously
in surface gravity waves and optics, and it has since become a widely studied
nonlinear effect [34]. In fluid dynamics, it is known as the Benjamin-Feir
instability, named after its discoverers [35, 36]. Benjamin and Feir
demonstrated that gravity waves on deep water bodies are unstable to harmonic
sideband perturbations, leading to the breakdown of wave trains due to
nonlinear interactions leading to the formation of more stable solitons. In
optics, MI was experimentally observed as a spatial modulation instability of
high-power lasers in organic solvents by Pilipetskii and Rustamov [37], and
Bespalov and Talanov published the corresponding theory [38].

In this work, we consider the phenomenon of modulation instability (MI) in the
presence of a specific dispersive parity-breaking term for fluid dynamics and
electromagnetic waves propagating in a non-reciprocal medium. For fluid
dynamics, our starting point is a nonlinear irrotational odd surface wave
model in the presence of gravity (a.k.a odd gravity waves) developed by two of
us [30, 39]. We then systematically derive the Non-Linear Schrödinger Equation
(NLSE) using the method of multiple scales. The parity-breaking effects result
in completely different stability criteria for the left and right movers.

We also observe a resonant interaction between the gravity waves and odd
surface waves in one of the chiral sectors, leading to a narrow stability
window. Such resonant interaction-induced stability windows have been
understood for gravity-surface tension interactions, but they symmetrically
manifest in both sectors [40, 41]. This chirality-dependent stability is the
hallmark of parity breaking. The consequence of chirality-dependent
instability leads to the formation of envelope solitons that propagate in only
one direction starting from an arbitrary initial condition. A similar NLSE has
been derived starting from a three-dimensional fluid with vertical odd
viscosity [33]; however, this setup is investigating three-dimensional odd
viscosity and is different from our two-dimensional starting point.

Following our analysis of odd gravity waves, a natural question arises whether
light propagating in a nonreciprocal dielectric medium can exhibit an odd-
viscosity-like dispersive parity-breaking phenomenon. Specifically, we aim to
explore if mechanisms of parity-breaking and the resulting chiral dispersion
observed in fluid dynamics can be mirrored in the electromagnetic domain,
potentially revealing new insights and applications in optical materials and
waveguides. To this end, we first introduce a minimal model for breaking
parity symmetry in a 1D dielectric. The term that breaks parity couples the
velocity of one of the oscillators to the position of the next oscillator.
This term in the continuum limit contains a spatial and a time derivative
which is odd under parity (P) and time reversal (T) but invariant under the
combined PT action. The continuum equations of electromagnetic waves passing
through this medium contain a dispersive parity-breaking term similar to the
odd viscosity term arising in fluids.

We show that the dispersion of the propagating waves becomes chiral in the
presence of this term. Next, we investigate the effect of Kerr nonlinearity on
the polarization dynamics in the context of modulation instability (MI) in
such nonreciprocal systems. To analyze the MI phenomena, we derive a parity-
broken nonlinear Schrödinger equation (NLSE) using the method of multiple
scales.

There are key differences between the odd gravity wave and the non-reciprocal
dielectric system. Firstly, the former is a 2D system with a 1D free surface,
while the latter is a purely 1D system. In the odd gravity case, the two
scalar fields are the dynamical surface profile and the scalar potential,
involving only first-order temporal derivatives. In the optical case, the two
fields are polarization and electric field, with dynamics governed by the
second-order coupled matter Maxwell’s equations in the continuum limit. The
nonlinearity in this case arises from the Kerr-type cubic nonlinearity in the
polarization equations.

Despite the fundamentally different underlying dynamics compared to the odd
gravity fluid, the form of the NLSE is retained, and the parity-breaking term
leads to an asymmetric stability criterion for left-moving and right-moving
chiral sectors. Unlike the odd gravity case, the stability criterion here is
entirely determined by the form of the dispersion relation, rather than the
explicit form of the interaction coefficient.

This paper is organized as follows. We begin with the review of the potential
theory for the odd gravity waves in Sec. 2 and derive the odd viscosity-
modified free surface dynamics. We then introduce the method of multiple
scales applied to this system and derive the NLSE for the modulated envelope
dynamics in Sec. 3. In Sec. 3.3, we investigate the stability of the NLSE.
Following the insights from the odd gravity wave phenomena, in Sec. 4, we then
develop a non-reciprocal 1D dielectric model with Kerr non-linearity. We
derive the NLSE for the optical system and investigate how the interplay
between parity breaking and MI manifests in the optical system.

##  2 Review of a non-linear potential theory for odd gravity waves

We begin by reviewing our recent work on deriving non-linear potential theory
for the free surface dynamics of 2D incompressible fluids with odd viscosity
and gravity[30, 39]. The hydrodynamic equations for incompressible fluids with
odd viscosity consist of the Euler equation, together with the
incompressibility condition, that is,

$\displaystyle\partial_{i}v_{i}$ $\displaystyle=0\,.$ (1)
$\displaystyle\partial_{t}v_{i}+v_{j}\partial_{j}v_{i}$
$\displaystyle=\frac{1}{\rho}\partial_{j}T_{ij}-\partial_{i}(gy)\,,$ (2)

Here, $v_{i}$ are the components of the velocity field, $\rho$ is the constant
and uniform fluid density, and the summation over repeated indices is assumed
($i,j=1,2$). The term $gy$ is the external gravitational potential and
$T_{ij}$ is the the stress tensor. For a fluid with only odd viscosity, we
have

$\displaystyle
T_{ij}=-p\delta_{ij}+\nu_{o}\rho\,(\partial_{i}^{*}v_{j}+\partial_{i}v_{j}^{*})\,.$
(3)

The first term of the stress tensor (3) is the standard isotropic pressure
term. The second term in Eq. (3) is the odd viscosity term. The coefficient
$\nu_{o}$ is known as kinematic odd viscosity (or Hall viscosity). In writing
this term we introduced the notation $a_{i}^{*}\equiv\epsilon_{ij}a_{j}$ so
that the “starred” vector $\mathbf{a}^{*}$ is just a vector $\mathbf{a}$
rotated by 90 degrees clockwise.

Under the incompressibility condition (1), we have $\Delta
v_{i}^{*}=\partial_{i}\omega$ and using this identity Eq. (2) becomes

$\partial_{t}v_{i}+v_{j}\partial_{j}v_{i}=-\partial_{i}\left(\frac{\tilde{p}}{\rho}+gy\right)\,.$
(4)

This differs from the ordinary Euler equation in the definition of the
modified pressure [6, 7] $\tilde{p}=p-\nu_{o}\rho\,\omega\,$ where
$\omega=\partial_{i}v_{i}^{*}$ is the fluid vorticity.

###  2.1  Irrotational limit of the free surface dynamics and boundary layer
approximation

We now review how one can use boundary layer theory to write effective free
surface dynamics of the odd gravity waves only in terms of the scalar
potential $\theta$ and the free surface profile $\eta$. We consider the fluid
domain as shown in Fig. 1. The free surface is a dynamical interface
$y=\eta(t,x)$ between two fluids (water and air) where we impose one kinematic
and two dynamical boundary conditions. The kinematic boundary condition states
that the velocity of the fluid normal to the boundary is equal to the rate of
change of the boundary shape.

$\displaystyle\partial_{t}\eta=\left(v_{y}-v_{x}\partial_{x}\eta\right)\big{|}_{y=\eta(t,x)}\,.$
(5)

The pair of dynamical boundary conditions $n_{i}T_{ij}\big{|}_{y=\eta(t,x)}=0$
impose that there are no normal and tangent forces acting on an element of the
fluid surface, where, $n_{i}$ are the components of the normal vector to the
surface $y=\eta(t,x)$. In our previous works, we have shown that the presence
of two dynamical boundary conditions, together with incompressibility,
requires a singular boundary layer where the divergent vorticity is confined
[30, 42]. The role of this singular boundary layer is to ensure that there are
no tangential forces on this interface.

We have shown that regardless of the boundary layer mechanism — dissipative
[30] or compressible [42] — the normal component of the dynamical boundary
condition (DBC) is universal and geometric [43]. Assuming that the boundary
layer is stable and is confined to short length scales, the free surface
problem can be written as an effective description of the fluid with the
effects of the boundary layer encoded in the odd viscosity modified pressure
term,

$\displaystyle\tilde{p}\big{|}_{y=\eta(t,x)}=(p-\nu_{o}\rho\,\omega)\big{|}_{y=\eta(t,x)}=\frac{2\nu_{o}\rho}{\sqrt{1+(\partial_{x}\eta)^{2}}}\partial_{x}v_{n}\big{|}_{y=\eta(t,x)}\,.$
(6)

Here $v_{n}$ is the velocity component which is normal to the boundary, taken
at $y=\eta(t,x)$. Only after taking care of the tangent stress boundary
condition and taking the thin boundary layer limit, we can define velocity
field as a gradient of a scalar function $v_{i}=\partial_{i}\theta$. In the
following, we set the density $\rho=1$, for simplicity.

Figure 1: Schematics of the deep water system with free surface. Here, $a$
denotes the amplitude of the edge profile, $\delta$ is the boundary layer
thickness and $\lambda$ is the wavelength.

The bulk equation for the irrotational system is simply the Laplace equation
for the potential

$\displaystyle\Delta\theta=0,\quad-\infty<y<\eta(t,x)$ (7)

The boundary conditions at the top and bottom of the fluid domain can be
written as,

$\displaystyle\theta_{y}=0\,,\qquad\qquad\qquad\qquad y=-\infty,$ (8)
$\displaystyle\eta_{t}+\eta_{x}\theta_{x}=\theta_{y}\,,\qquad\qquad
y=\eta(t,x),$ (9)
$\displaystyle\theta_{t}+\frac{\theta_{x}^{2}+\theta_{y}^{2}}{2}+g\eta=\frac{2\nu_{o}}{\sqrt{1+\eta_{x}^{2}}}\,\partial_{x}\left[\frac{\eta_{t}}{\sqrt{1+\eta_{x}^{2}}}\right],\;\,y=\eta(t,x).$
(10)

Here, Eq. (9) is the kinematic boundary condition and Eq. (10) is the odd
viscosity modified dynamic boundary condition on the top surface, whereas Eq.
(8) simply states that $v_{y}=0$ in the bottom surface. Note that the boundary
condition at the bottom surface is a slip boundary condition. For the deep
water problem, a more convenient set of variables is given by
$(x,y,t)\rightarrow(x,\chi,t)$ where $\chi=y-\eta(x,t)$. The advantage of
these variables is that it removes the implicit dependence of the $\theta$ on
the free surface profile $\eta$ in the boundary conditions. On the other hand,
we pay the price by making the bulk Laplace equation non-linear. In terms of
these variables, the bulk equation is given by

$\displaystyle\theta_{xx}+(1+\eta_{x}^{2})\theta_{\chi\chi}-2\partial_{x}(\eta_{x}\theta_{\chi})+\eta_{xx}\theta_{\chi}=0,\qquad\qquad-\infty\leq\chi\leq
0$ (11)

Figure 2: Schematics of the deep water system with free surface in the the new
coordinate system, where the vertical coordinate $y$ is transformed to
$\chi=y-\eta(t,x)$, resulting in a flat top surface.

The schematic of the free surface system in the $(x,\chi,t)$ coordinates is
shown in Fig. 2. The equations at the top and bottom boundary can be written
as

$\displaystyle\theta_{\chi}=0\,,\qquad\qquad\qquad\qquad\qquad\,\,\,\,\chi=-\infty,$
(12)
$\displaystyle\eta_{t}+\eta_{x}\theta_{x}=(1+\eta_{x}^{2})\theta_{\chi}\,,\qquad\qquad\chi=0,$
(13)
$\displaystyle\theta_{t}-\eta_{t}\theta_{\chi}+\frac{\theta_{x}^{2}+\theta_{\chi}^{2}-2\eta_{x}\theta_{\chi}\theta_{x}+\eta_{x}^{2}\theta_{\chi}^{2}}{2}+g\eta=\frac{2\nu_{o}}{\sqrt{1+\eta_{x}^{2}}}\,\partial_{x}\left[\frac{\eta_{t}}{\sqrt{1+\eta_{x}^{2}}}\right],\;\,\chi=0.$
(14)

where the subscript notation corresponds to the action of partial derivative
(i.e. $\theta_{x}:=\partial_{x}\theta$). In the following, we will focus on
the weakly non-linear analysis using the method of multiple scales of the
above system of equations.

##  3 Method of multiple scales for envelope dynamics and non-Linear
Schrödinger equation

The effect of modulation instability results from the interaction between a
strong carrier harmonic wave and small sideband perturbations. This is a
particular case of four-wave interaction, where the weak modulation imposed on
a harmonic wave leads to the dynamics of the modulated envelope at much slower
time scales. The slow dynamics of the envelope can be accessed by the method
of multiple scales [44], and due to the non-linear interactions of the waves,
it can lead to the unstable growth of the sideband perturbations. These
instabilities can be a precursor to the formation of more stable entities,
such as solitons. To find envelope dynamics for the system, we introduce
multi-scale variables

$\displaystyle\sigma=kx-\omega
t\,,\quad\xi=\epsilon(x-c_{g}t)\,,\quad\tau=\epsilon^{2}t$ (15)

where $k$ is the arbitrary wavenumber, $\omega$ is the dispersion,
$c_{g}=d\omega/dk$ is the group velocity of the wave-packet and $\epsilon\ll
1$ is the wave steepness parameter. Furthermore, we seek asymptotic solutions
for the system of equations (11-14) that will lead to the envelope dynamics
for the slow degrees of freedom, that is,

$\displaystyle\theta(\sigma,\xi,\chi,\tau)$
$\displaystyle=\epsilon\,\theta_{1}(\sigma,\xi,\chi,\tau)+\epsilon^{2}\,\theta_{2}(\sigma,\xi,\chi,\tau)+\epsilon^{3}\,\theta_{3}(\sigma,\xi,\chi,\tau)+...,$
(16) $\displaystyle\eta(\sigma,\xi,\tau)$
$\displaystyle=\epsilon\,\eta_{1}(\sigma,\xi,\tau)+\epsilon^{2}\,\eta_{2}(\sigma,\xi,\tau)+\epsilon^{3}\,\eta_{3}(\sigma,\xi,\tau)+...$
(17)

At every order, the fields are further expanded to the relevant number of
harmonics that are expected to be generated due to the non-linear
interactions,

$\displaystyle\theta_{j}(\sigma,\xi,\chi,\tau)=\sum_{n=-j}^{j}\theta_{jn}(\xi,\chi,\tau)q^{n},\quad\eta_{j}(\sigma,\xi,\tau)=\sum_{n=-j}^{j}\eta_{jn}(\xi,\tau)q^{n},\quad
q=e^{i\sigma}.$ (18)

Variables $\theta_{j}$ and $\eta_{j}$ are the perturbative expansion of fields
which contributes to the solution at different orders in nonlinearity. They
are further expanded in ($\theta_{jn}$ and $\eta_{jn}$) corresponding to the
various possible harmonics generated at that order. This form is motivated by
the fact that higher harmonics become relevant to dynamics processes at higher
orders. A systematic calculation involves comparing terms at every order of
$\mathcal{O}(\epsilon^{n}q^{m})$ for both bulk and boundary equations. All the
higher harmonics are written in terms of the slow envelope amplitude at the
first order $\eta_{11}(\xi,\tau)$ and at the order
$\mathcal{O}(\epsilon^{3}q^{1})$, we obtain the NLSE as a solvability
condition on $\eta_{11}(\xi,\tau)$. Below we outline the order-by-order
results relegating the details of this calculation to the Appendix 7.1.

###  3.1 Linear order

Plugging the multiscale ansatz defined in Eqs. (17) and (18) into Eqs. (11-14)
and comparing terms at the order $\mathcal{O}(\epsilon\,q)$, we obtain,

$\displaystyle\partial_{\chi}^{2}\theta_{11}$
$\displaystyle=k^{2}\theta_{11},\quad-\infty\leq\chi\leq 0$ (19)
$\displaystyle-i\omega\eta_{11}$
$\displaystyle=\partial_{\chi}\theta_{11}\,\,\,\qquad\chi=0$ (20)
$\displaystyle-i\omega\theta_{11}+g\eta_{11}$ $\displaystyle=2\nu_{o}\omega
k\eta_{11},\quad\chi=0$ (21) $\displaystyle\partial_{\chi}\theta_{11}$
$\displaystyle=0,\qquad\quad\quad\chi=-\infty$ (22)

We consider $k>0$ from here on. For the $\mathcal{O}(\epsilon)$ order, we
additionally have $\theta_{10}=A_{0}(\xi,\tau)$ and $\eta_{10}=0$.
$A_{0}(\xi,\tau)$ is an arbitrary function that is fixed at higher order.
Solving these equations and boundary conditions, in the linear order we
obtain:

$\displaystyle\theta_{11}=A(\xi,\tau)e^{k\chi}\qquad\eta_{11}=i(k/\omega)A(\xi,\tau)$
(23)

where $A(\xi,\tau)$ is the slowly varying amplitude of the envelope whose
dynamics in these slow variables will be fixed to be NLSE at the
$\mathcal{O}(\epsilon^{3}q)$. Solving the equation above, we obtain the
dispersion relation for the odd gravity waves [30],

$\displaystyle\omega_{\pm}=-\nu_{o}k^{2}\pm\sqrt{gk+\nu_{o}^{2}k^{4}}\,.$ (24)

This dispersion has a different form for the two sectors indicating parity
breaking symmetry. In the zero gravity limit ($g\rightarrow 0$) the effects of
the parity breaking are maximal, where we observe that the right-movers tend
towards $\omega_{+}\rightarrow 0$, while left-movers have a parabolic
dispersion $\omega_{-}\rightarrow-2\nu_{o}k^{2}$ as shown in Fig. 3.

Figure 3: Dispersion for the gravity waves and odd-gravity waves
$\omega_{\pm}$. The dashed curve corresponds to the gravity waves
$(\nu_{o}=0)$. The characteristic scales $k_{c}$ and $\omega_{c}$ are defined
in terms of $\nu_{o}$ and $g$, such that,
$k_{c}=\left(\tfrac{g}{\nu_{o}^{2}}\right)^{1/3}$ and
$\omega_{c}=\nu_{o}k_{c}^{2}$.

###  3.2 Higher order non-linear terms and envelope dynamics

The higher-order calculations are a bit unwieldy so we only present the final
answers with the details relegated to the Appendix 7.1. At second order and
first harmonic $\mathcal{O}(\epsilon^{2}q)$ we get:

$\displaystyle\theta_{21}=A_{2}(\xi,\tau)e^{k\chi}-iA_{\xi}\chi
e^{k\chi},\qquad\eta_{21}=i(k/\omega)A_{2}+\left(\frac{\omega-
c_{g}k}{\omega^{2}}\right)A_{\xi}\,,$ (25)

where $A_{2}(\xi,\tau)$ is arbitrary and remains undetermined since this term
always appears in a particular combination with the $\eta_{21}$. This also
generates second harmonics due to non-linear interactions present at this
order, which are given by

$\displaystyle\theta_{22}=\left(\frac{-6i\nu_{o}k^{4}}{\omega(\omega+6\nu_{o}k^{2})}\right)A^{2}e^{2k\chi}+i\left(\frac{k^{2}}{\omega}\right)A^{2}e^{k\chi}\qquad\eta_{22}=\left(\frac{-k^{3}}{\omega(\omega+6\nu_{o}k^{2})}\right)A^{2}.$
(26)

Note that the second-harmonic solution has a resonance in one of the sectors
when the condition $\omega_{-}+6\nu_{o}k^{2}=0$ is satisfied. At this value of
$k$, our expansion breaks down and one needs a different counting scheme in
$\epsilon$ near this value of $k$, which is beyond the scope of this paper.
This resonance also manifests in the gravity-capillary waves for a particular
value of $k$ and has been addressed by McGoldrick [41, 45].

At the third order, that is, $\mathcal{O}(\epsilon^{3})$, the bulk equations
are solved for $\theta_{30}$ and $\theta_{31}$, which gives us

$\displaystyle\theta_{30}(\xi,\chi,\tau)=\left(\frac{-\partial_{\xi\xi}^{2}A_{0}}{2}\right)\chi^{2}+\left(\frac{k^{2}c_{g}}{\omega^{2}}+\frac{k^{2}\chi}{\omega}\right)e^{k\chi}\partial_{\xi\xi}^{2}|A_{1}|^{2}.$
(27)

The boundary condition imposes ($\partial_{\chi}\theta_{30}=0$ as
$\chi\rightarrow-\infty$) leading to the constraint $A_{0}=0$ 111Note that
this is not the case if we use the finite depth fluid and we need to express
$A_{0}$ in terms of $A$ and its derivatives in a consistent way. We further
need the bulk expression for the $\theta_{31}$ fields to substitute them into
the two surface equations at the order $\mathcal{O}(\epsilon^{3}q^{1})$. The
two surface equations in this order can be viewed as a system of two linear
equations for the fields $\theta_{31}(\xi,0,\tau)$ and $\eta_{31}(\xi,\tau)$.
Solving this system at the third order, we finally arrive at a Non-Linear
Schrödinger Equation (NLSE) for the envelope amplitude $A(\xi,\tau)$. This
NLSE turns out to be the consistency condition for the perturbative solution
to hold and is typical of how the method of multiple scales works. The
coefficients of NLSE now depend on parity-breaking parameter $\nu_{o}$,
suggesting different stability criteria of the envelope dynamics for the right
propagating NLSE and the left propagating NLSE. The general form of the NLSE
can be written as

$\displaystyle iA_{\tau}+LA_{\xi\xi}-MA|A|^{2}=0\,,$ (28)

where we have defined

$\displaystyle L=\frac{1}{2}\omega^{\prime\prime}(k),\qquad
M=\frac{2k^{4}}{\omega}\left(1-\frac{3\nu_{o}k^{2}\omega}{(\omega+6\nu_{o}k^{2})(\omega+\nu_{o}k^{2})}\right).$
(29)

The parity-breaking effect of odd viscosity enters in both the dispersion
relation $\omega(k)$ and in the interaction coefficient $M$. Before we discuss
the stability analysis of the NLSE for the general case, let us take the limit
when $\nu_{o}\rightarrow 0$ and compare it with the NLSE obtained for the
gravity waves. In this regime, the dispersion relation becomes
$\omega=\omega_{\text{gravity}}=\pm\sqrt{gk}$, the interaction term becomes
$M=2k^{2}/\omega_{\text{gravity}}$ and the NLSE takes the form

$\displaystyle
iA_{\tau}+\frac{1}{2}\omega^{\prime\prime}_{\text{gravity}}A_{\xi\xi}-\frac{2k^{2}}{\omega_{\text{gravity}}}A|A|^{2}=0,\qquad\nu_{o}\rightarrow
0.$ (30)

This equation is consistent with the NLSE for the deep water gravity waves
[44, 46].

In the opposite limit of $g\rightarrow 0$, we still have propagating odd
surface waves. They fully break parity symmetry, which manifests in completely
different envelope dynamics depending on whether the choice of sector, i.e.,
$\omega=\omega_{\pm}$. First, we will consider the $\omega_{-}$ case and take
the $g\rightarrow 0$ limit. In this case the dispersion relation becomes
$\omega_{-}=-2\nu_{o}k^{2}$ and consequently $L=-2\nu_{o}$ and
$M=k^{2}/(2\nu_{o})$ leading to the following NLSE,

$\displaystyle
iA_{\tau}-2\nu_{o}A_{\xi\xi}-\frac{k^{2}}{2\nu_{o}}A|A|^{2}=0\qquad
g\rightarrow 0,\,\,\omega=\omega_{-}.$ (31)

On the other hand, taking the $g\rightarrow 0$ limit is a bit pathological for
the $\omega_{+}$ sector since the dispersion $\omega_{+}=0$ becomes a zero
mode. This causes the interaction term $M$ to diverge and the dispersion term
$L=0$. We can consider the leading order behavior in the small $g$ limit as a
regulator and we obtain

$\displaystyle
iA_{\tau}+\frac{1}{2}\left(\frac{g}{\nu_{o}k^{3}}\right)A_{\xi\xi}-\frac{4\nu_{o}k^{5}}{g}A|A|^{2}=0\qquad
g\rightarrow 0,\,\,\omega=\omega_{+}\,.$ (32)

In fact, this forces the amplitude to vanish. However, scaling the amplitude
by $\sqrt{g}$, that is, $A=\sqrt{g}\tilde{A}$, and taking the $g\rightarrow
0$, we find that

$\displaystyle i\tilde{A}_{\tau}-4\nu_{o}k^{5}|\tilde{A}|^{2}\tilde{A}=0\,.$
(33)

This shows that the envelope amplitude becomes asymptotically proportional to
$\sqrt{g}$, with the proportionality constant $\tilde{A}$ satisfying Eq. (33).
The absence of the dispersion term keeps $\tilde{A}$ bounded, since Eq. (33)
can be solve analytically, giving us

$\displaystyle\tilde{A}=\tilde{A}_{0}\,e^{-i4\nu_{o}k^{5}|\tilde{A}_{0}|^{2}\tau}\,,$
(34)

for some constant $\tilde{A}_{0}$.

###  3.3 Stability of the NLSE and soliton formation

We now consider the linear stability of the envelope dynamics for the NLSE
derived in Eqs. (28) and (29). The stability of the envelope dynamics is
deduced by expanding the envelope around Stokes solution
$A=A_{0}e^{-iM|A_{0}|^{2}\tau}$. Let us consider infinitesimal perturbations
to the solution:

$A=A_{0}(1+a)e^{i(\mu-M|A_{0}|^{2}\tau)}$ (35)

Substituting the ansatz above in the NLSE, and keeping first-order terms in
$a$ and $\mu$, we get two equations, which correspond to the real and
imaginary parts

$\displaystyle\partial_{\tau}a+L\partial_{\xi}^{2}\mu=0,\qquad\partial_{\tau}\mu-L\partial_{\xi}^{2}a+2M|A_{0}|^{2}a=0.$
(36)

This linear system of equations can be expanded in modes
$(a,\mu)=(\kappa_{a},\kappa_{\mu})e^{i(\delta\xi-\lambda\tau)}$ and the
resulting system can be cast into a matrix form

$\begin{pmatrix}i\lambda&\delta^{2}L\\\
(2M|A_{0}|^{2}+\delta^{2}L)&-i\lambda\end{pmatrix}\begin{pmatrix}\kappa_{a}\\\
\kappa_{\mu}\end{pmatrix}=0\,.$ (37)

This expression admits real solutions for $\lambda$ when

$L(2M|A_{0}|^{2}+\delta^{2}L)\geq 0\,.$ (38)

The perturbation remains bounded when $\lambda$ is real and the Stokes’
solution remains stable. On the other hand, when the condition (38) is not
satisfied the perturbation grows in an unbounded fashion. In this case, the
harmonic perturbation on the mode grows with growth rate $\gamma$, which is
given by

$\displaystyle\gamma=|\delta|\sqrt{-2ML|A_{0}|^{2}-\delta^{2}L^{2}}.$ (39)

Figure 4: Time evolution of initial wave packet in two regimes of the NLSE.
The final state at $\tau=2$ for the two cases has been contrasted. (Left)
Initial state consisting of a primary harmonic with $\bar{k}=\pi/5$ and
sideband perturbations in bandwidth $\delta=0.001$,
$A(\xi,0)=\sin{\bar{k}\xi}+0.01\left(\sin{(\bar{k}+\delta)\xi}+\sin{(\bar{k}-\delta)\xi}\right)$.
(Right) Initial Gaussian wave packet, $A(\xi,0)=e^{-\xi^{2}}$.

In the zero bandwidth regime of the perturbation, that is, $\delta\rightarrow
0$ limit, the criterion for the envelope stability under modulations in
amplitude (or phase) simplifies to

$\displaystyle LM\geq 0.$ (40)

In Fig. 4, we investigate the initial value problem for the NLSE with $LM>0$
and $LM<0$, using a perturbed plane wave and a Gaussian packet as the initial
conditions. For the case where $LM>0$, we observe a “superlinear” spreading of
the initial profile, while for the unstable case where $LM<0$, we notice the
onset of soliton formation.

In the vanishing bandwidth limit, the stability of the NLSE only depends on
the relative signs of the $L$ and $M$ terms. We now analyze the relative signs
of these for the expressions given in Eq. (29). A compact way of studying the
stability criterion for different parameter regimes is in terms of
dimensionless variables $\bar{\omega}=\omega/\omega_{c}$, $\bar{k}=k/k_{c}$
and $\bar{A}=\nu_{o}A$, defined in terms of the characteristic scales
$k_{c}=\left(\frac{g}{\nu_{o}^{2}}\right)^{1/3}$,
$\omega_{c}=\nu_{o}k_{c}^{2}$. In terms of the dimensionless variables, we see
that the small wavenumber limit $(\bar{k}\ll 1)$ corresponds to the gravity-
dominated regime and the large wavenumber one $(\bar{k}\gg 1)$ is the odd
viscosity dominated regime, where parity breaking effects are most pronounced.
The interpolating regime $\bar{k}\sim 1$ is where the gravity effects are
comparable to that of the odd viscosity effects. There is also a resonance
present due to the interaction between the gravity waves and odd surface waves
in the intermediate regime $\bar{k}\sim 1$ albeit only in the $\omega_{-}$
sector.

It is straightforward to see that the gravity-dominated regime ($\bar{k}\ll
1$) with $LM=-\bar{k}^{2}/4<0$ is unstable for both $\omega_{\pm}$ sectors. On
the other hand, we observe a sector-dependent stability in the odd viscosity-
dominated regime of $\bar{k}\gg 1$. For the $\omega_{-}$ sector, we have
$LM=-\bar{k}^{2}<0$ which is unstable and leads to soliton formation. However,
the $\omega_{+}$ sector in this regime is fully stable with
$LM=2\bar{k}^{2}>0$ (see Figs. 5 and 6). The physical consequence of such a
strong parity breaking is that soliton formation post instability will only
manifest for the $\omega_{-}$ and the $\omega_{+}$ may become “superlinear”
spreading of an initially localized perturbation.

The intermediate regime of $\bar{k}\sim 1$ shows non-trivial interaction
between the gravity waves and odd surface waves and it is in this regime where
the stability criterion changes.

Figure 5: Left: Stability for the $\omega_{+}$ sector of the odd-gravity
waves, calculated in vanishing bandwidth $\delta\rightarrow 0$. Right:
Coefficients of NLSE. In this sector, the band curvature changes the sign for
$\bar{k}>0.965$ and solely determines the onset of stability since the
nonlinearity doesn’t change the sign.

Figure 6: Left: Stability for the $\omega_{-}$ sector of the odd-gravity
waves, calculated in vanishing bandwidth $\delta\rightarrow 0$. Right:
Coefficients of NLSE. In this sector, the band curvature changes sign for
$\bar{k}>0.218$. The nonlinearity changes sign for $\bar{k}>0.346$. The window
$0.218<\bar{k}<0.346$ is stable to MI.

For the $\omega_{+}$ sector, the stability arises because of the change in the
sign of the band curvature at $\bar{k}=0.965$. We calculate this by finding
the real root of the equation $\bar{\omega}^{{}^{\prime\prime}}(\bar{k})=0$
for $\omega_{+}$. The interaction term $M$ does not change its sign in this
sector as shown in Fig. 5.

The $\omega_{-}$ sector is interesting because there is a window of stability
that emerges in this sector due to the sign change in the band curvature at
$\bar{k}=0.218$ and also a resonance-induced sign change in the interaction
term $M$ at $\bar{k}=0.346$ (see Fig. 6). As discussed before, near this
resonance our original counting scheme breaks down and a different counting
scheme needs to be used. Similar resonance also manifests in the gravity
capillary interaction, however, it happens for both sectors. These two
consecutive sign changes in $L$ and $M$ for nearby values of $\bar{k}$ lead to
a narrow window of stability in the range ($0.218<\bar{k}<0.346$) as shown in
Fig. 6.

Moreover, in the unstable regions, the growth rate of MI also shows asymmetric
behavior amongst the two sectors, as shown in Fig. 7. We notice that for the
$\omega_{+}$ sector the growth rate is enhanced relative to the gravity case,
while for the $\omega_{-}$ sector the growth rate is suppressed relative to
the gravity case. Also, for $\omega_{-}$, at intermediate and large $\bar{k}$,
the growth rate is closer to the chiral limit.

Figure 7: (Left) Dimensionless growth rate of instability in sector
$\omega_{+}$: shows the enhanced growth rate of MI, relative to gravity waves.
(Right) The dimensionless growth rate of instability in $\omega_{-}$ sector.
In both cases, we set $\delta=0.0001\,k_{c}$ and $|A_{0}|=\nu_{o}$.

Building on our analysis of modulation instability (MI) in odd gravity waves,
a natural question arises: can light propagating in a nonreciprocal dielectric
medium exhibit a dispersive parity-breaking phenomenon similar to odd
viscosity? Furthermore, in the following sections, we seek to study how this
parity-breaking effect interacts with nonlinear MI in such optical systems.

##  4 Nonreciprocal 1D dielectric

Charges in a dielectric medium are bound in molecules and cannot move freely
throughout the whole sample. Hence, their response to an external electric
field causes a change in the charge distribution within each molecule, which
induces an electric dipole moment in it.

In particular, one-dimensional dielectric systems can be modeled by a group of
damped harmonic oscillators that imitate molecules in a medium. The molecules
have heavy nuclei that form a one-dimensional lattice along the $x$-axis, and
each molecule has a mobile degree of freedom that enables it to move
perpendicular to the lattice direction. When the spacing between adjacent
atoms becomes extremely small, the model can be approximated by a continuous
dielectric. In this case, the collection of molecular electric dipoles gives
rise to a polarization field, which modifies the propagation of the
electromagnetic field inside the medium.

When an electromagnetic wave propagates in this dielectric wire, the
electromagnetic pulse moves forward, it interacts with the medium by
initiating motion in the individual oscillators, which mimics the movement of
electrons bound to nuclei in the dielectric. This links the evolution of the
electromagnetic pulse in the medium to its effective polarization.

For the sake of simplicity, we will restrict the following analysis to
linearly polarized electromagnetic waves. In this case, electric and magnetic
fields are on the form $\vec{E}(\vec{r},t)=E(x,t)\hat{y}$ and
$\vec{B}(\vec{r},t)=B(x,t)\hat{z}$ such that, they are connected through
Faraday’s law:

$\displaystyle E_{x}+\frac{1}{c}B_{t}$ $\displaystyle=0\,.$ (41)

The electric field drives the oscillators such that their dynamics are
described by

$\displaystyle
m_{n}\left(\ddot{y}_{n}+\gamma_{n}\dot{y}_{n}+\omega_{n}^{2}y_{n}\right)=q_{n}E(x_{n},t)\,,$
(42)

where $m_{n}$ is the oscillator effective mass, $\omega_{n}$ is the oscillator
frequency, $q_{n}$ is the charge of the movable particle and $\gamma_{n}$ is
the oscillator damping coefficient. Assuming that the oscillator are identical
$\omega_{n}=\omega_{0}$ and separated by lattice-spacing $a$ along the
$x$-direction, it is possible to express Eq. (42) in terms of the molecular
electric dipole moment $p_{n}=qy_{n}$:

$\displaystyle\ddot{p}_{n}+\gamma\dot{p}_{n}+\omega_{0}^{2}p_{n}=\frac{q^{2}}{m}E(na,t)\,.$
(43)

The polarization field is nothing but the electric dipole density. Thus, it
can be defined by dividing the molecular electric dipole by the wire cross-
section area and by the spacing between molecules and taking the continuum
limit, that is,

$\displaystyle P(x,t)=\lim_{a\rightarrow 0}\frac{q}{Aa}y_{n}(t)\,.$ (44)

This limit is taken by setting $x=\lim_{a\rightarrow 0}na$ as the location of
the infinitesimal electric dipole.

Therefore, for a continuum one-dimensional dielectric medium, the governing
equations are given by

$\displaystyle P_{tt}+\gamma
P_{t}+\omega_{0}^{2}P=\frac{\omega_{P}^{2}}{4\pi}E\,,$ (45) $\displaystyle
B_{x}+\frac{1}{c}E_{t}=-\frac{4\pi}{c}P_{t}\,,$ (46)

together with Faraday’s law (41). The coefficient $\omega_{P}$ is often called
plasmon frequency. The magnetic field can be eliminated by combining Eqs. (41,
46) which leads to

$\displaystyle c^{2}E_{xx}-E_{tt}=4\pi P_{tt}\,.$ (47)

However, the presence of the damping coefficient $\gamma$ breaks the time-
reversal symmetry, which leads to the absorption of the electromagnetic wave
and a finite penetration depth in the wire. In this scenario, nonlinearities
become irrelevant and the linear theory completely describes the system.
Therefore, to study nonlinear media, one must ensure that $\gamma\rightarrow
0$, which preserves the time-reversal symmetry.

###  4.1 Dispersive parity breaking in a dielectric

Following the intuition drawn from the hydrodynamic case, it should be
possible to break the time-reversal symmetry keeping the system non-
dissipative; as long as the combined PT symmetry is preserved. Eqs. (45) and
(47) are symmetric under the parity transformation $x\rightarrow-x$.
Therefore, terms such as $P_{tx}$, $P_{t,xxx}$, etc., violate parity and time-
reversal symmetry, but preserve the combined PT symmetry 222In fact, any term
with an odd order of time derivatives and an odd order of spatial derivatives
will break T and P, but preserves PT.. Modifying the dynamical equation for
the polarization to accommodate for the simplest of these chiral terms gives
us the following set of equations

$\displaystyle
P_{tt}+\beta_{\text{odd}}P_{tx}+\omega_{0}^{2}P=\frac{\omega_{P}^{2}}{4\pi}E\,,$
(48) $\displaystyle c^{2}E_{xx}-E_{tt}=4\pi P_{tt}\,.$

Note that equation (47) since it is derived directly from Maxwell’s equations.

It is not hard to see that Eqs. (47-48) lead to wave propagation without any
absorption, which satisfies the following dispersion relation:

$(c^{2}k^{2}-\omega^{2})(\omega^{2}-\beta_{\text{odd}}\,\omega
k-\omega_{0}^{2})+\omega^{2}\omega_{P}^{2}=0\,.$ (49)

For $\beta_{\text{odd}}\neq 0$, the dispersion relation is not invariant under
$\omega\rightarrow-\omega$ (time-reversal symmetry) or $k\rightarrow-k$
(parity symmetry) alone, nevertheless it remains the same if we combine both.
Furthermore, note that parity and time-reversal become symmetries of the
system if $\beta_{\text{odd}}$ change sign under them, that is,
$\beta_{\text{odd}}\rightarrow-\beta_{\text{odd}}$ for both P and T
symmetries.

In this medium, the chiral nature is evident, since the presence of the
coefficient $\beta_{\text{odd}}$ skews the dispersion relation and creates an
asymmetry between left and right moving modes, as seen in Fig. 8. This effect
is more pronounced in the two middle bands at large $k$. In this limit,
polarization and electromagnetic waves decouple from each other and the middle
bands describe the polarization waves, which asymptotically satisfy the the
dispersion relation $\omega^{2}-\beta_{\text{odd}}\,\omega
k-\omega_{0}^{2}=0$, or

$\displaystyle\omega\sim\tfrac{1}{2}\left(\beta_{\text{odd}}k\pm\sqrt{\beta_{\text{odd}}^{2}k^{2}+4\omega_{c}^{2}}\right).$
(50)

Equivalently, the top and the bottom bands describe the uncoupled
electromagnetic waves at large $k$, which satisfy the dispersion relation
$\omega=\pm c|k|$. For these bands, the parity breaking shows up for small
$k$, when the coupling between polarization and electromagnetic waves is
substantial.

Figure 8: Dispersion relation for the parity broken dielectric media. The
dashed curve is for the parity-symmetric dielectric. Scaling is such that
$(\omega_{P}/\omega_{0})=0.5$.

###  4.2 Heuristics of the Chiral Term

The term $\beta_{\text{odd}}P_{tx}$ was phenomenologically introduced in the
setup, however, it can also be obtained from a microscopic model. In the
oscillator picture, we find that

$\displaystyle\beta_{\text{odd}}P_{tx}\approx\frac{\beta_{\text{odd}}\,q}{Aa}\left(\frac{\dot{y}_{n+1}-\dot{y}_{n-1}}{2a}\right).$
(51)

For $\beta q>0$, this corresponds to a force acting on the oscillator $n$,
which is along the direction of motion of its right neighbor and opposite to
the motion of its left neighbor. This automatically breaks parity, since each
molecule feels a drag force coming from the adjacent right neighbor and a push
to oppose the motion of the adjacent left neighbor.

Even though these forces could be derived directly from a Lagrangian model,
obtained from the term

$\displaystyle\mathcal{L}_{PT}=\frac{m\beta_{\text{odd}}}{2a}y_{n}\dot{y}_{n+1}\,,$
(52)

it is not clear to us what underlying mechanism allows for such a term.
Identifying this can lead to interesting non-reciprocal effects in 1D
polariton systems since this model can be viewed as a classical version of the
polariton systems.

###  4.3 Kerr Nonlinearity

Nonlinearities can be introduced into the model by adding anharmonic terms to
the restoring force. The simplest term that preserves the oscillator’s
symmetry is a cubic restoring force: $F_{Kerr}=-\kappa y_{n}^{3}$. This is
generally referred to as the Kerr medium and allows for the generation of
higher harmonics. The polarization dynamics now becomes

$P_{tt}+\beta_{\text{odd}}P_{tx}+\omega_{0}^{2}P+\alpha
P^{3}=\frac{\omega_{P}^{2}}{4\pi}E\,.$ (53)

In the optics case, the Maxwell equation (47) remains linear, while
nonlinearities emerge in the polarization dynamics. By applying the method of
multiple scales, as previously used in the hydrodynamic case, we can analyze
modulation instability (MI) in this classical polariton system.

After introducing the same variables defined in Eq. (15), and applying the
method of multiple scales to Eqs. (47) and (53), we obtain the NLSE
corresponding to the polariton wave envelope to be

$i\partial_{\tau}P_{11}+\tfrac{1}{2}\,\omega^{\prime\prime}(k)\partial_{\xi}^{2}P_{11}-M|P_{11}|^{2}P_{11}=0$
(54)

where

$\displaystyle
M=\frac{3\,\alpha\,\omega_{P}^{2}\,\omega^{3}}{\omega_{P}^{2}\,\omega^{2}(\omega^{2}+\omega_{0}^{2})+(\omega^{2}+c^{2}k^{2})(\omega_{0}^{2}+\beta_{\text{odd}}\,\omega
k-\omega^{2})^{2}}\,.$ (55)

For calculation details, we refer to Appendix 7.2. Since $\alpha>0$, the sign
of $M$ is determined by the sign of $\omega$. This means that the condition
$LM\geq 0$ can be recast as

$\omega\frac{d^{2}\omega}{dk^{2}}\geq 0.$ (56)

The stability criterion in this case is simpler than in the hydrodynamic case,
as it only depends on the nature of the dispersion relation. When
$\beta_{\text{odd}}=0$, the top and bottom bands are stable under modulation
instability (MI) as they satisfy Eq.( 56). The two central bands, on the
contrary, are unstable for all $k\neq 0$ modes in this parity-preserving case.
When $\beta_{\text{odd}}\neq 0$, while the top and bottom bands remain stable,
new stability pockets appear in the middle bands.

In particular, we focus on the lower polariton band (the positive-frequency
middle band) in the parity-broken case with $\beta_{\text{odd}}=-0.2c$ (see
Fig. 9 left panel). For this case, the dispersion relation satisfies
$\omega>0$ while $\omega^{\prime\prime}(k)$ is positive for
$k<-3.579\,\omega_{0}/c$ and $k>2.989\,\omega_{0}/c$, resulting in the
stabilization of modes in these regions. In the region where modes are
unstable, parity breaking modifies the growth rate of instability in an
asymmetric way (differently for $k>0$ and $k<0$) as is evident in Fig. 10.
Furthermore, in the ($\beta_{\text{odd}}$ vs $k$) plot, shown in Fig. 9 right
panel, the regions of stability appear only when $\beta_{\text{odd}}\neq 0$ at
certain critical values of $k$, which are determined by solving
$\omega^{\prime\prime}(k)=0$. Note that $k=0$ is trivially stable since the
nonlinear term vanishes for this case.

Figure 9: (Left) Lower Polariton band for parity broken case, with
$\beta_{\text{odd}}=-0.2c$. For the two sectors, the band curvature changes
sign in the regions $ck/\omega_{0}<-3.579$ and $ck/\omega_{0}>2.989$, which
leads to the stabilization of modes in the respective sectors. (Right)
Stability in the $\beta_{\text{odd}}$-k space of the lower polariton mode, in
$\delta\rightarrow 0$ limit. Dark regions are unstable and light regions are
stable. Stability is calculated in vanishing bandwidth limit ($\delta=0$). The
parameters of the model considered are the same as Fig. 8.  Figure 10:
Dimensionless growth rate for MI in parity broken dielectric, with
$\beta_{\text{odd}}=-0.2c$. Parameters considered are the same as Fig. 8, with
$\delta=0.0001\,\omega_{0}/c$ and $|A_{0}|=\omega_{0}/\sqrt{\alpha}\,$.

##  5 Discussion

In this work, we studied the interplay of parity breaking symmetry and the
modulation instability (MI) in two distinct physical systems, fluid dynamics
and optics, with characteristic wavelengths differing in principle by orders
of magnitude. For the hydrodynamics system, we considered the nonlinear odd
gravity waves: a 2D free surface problem with gravity $g$ and odd viscosity
$\nu_{o}$. Using the method of multiple scales, we derived the nonlinear
Schrödinger equation (NLSE) for the slowly evolving envelope dynamics. We
performed a linear stability analysis for this NLSE and showed that the
stability criteria for the left-moving and right-moving NLSE are very
different due to the parity breaking in the system. In the intermediate
regime, where gravity’s effects are comparable to parity breaking, gravity and
odd viscosity-induced capillary waves interact resonantly, resulting in a
narrow window of stability around this resonance. At resonance, our counting
scheme breaks down, and a different counting scheme is required. This
situation is reminiscent of the resonant interaction between gravity waves and
surface tension-induced capillary waves. However, in the case of gravity and
surface tension, both left and right propagating waves resonate, unlike the
present case where the resonance occurs in only one sector. One interesting
consequence of such chirality-dependent stability is that MI-triggered
solitons propagate in only one direction in these nonlinear systems.

Following the intuition drawn from the hydrodynamical case, we construct a
model for a non-reciprocal dielectric in 1D that realizes analogous parity-
broken effects. The term that breaks parity couples the velocity of one
oscillator to the position of the next oscillator. In the continuum limit,
this term includes a spatial and a time derivative, which is odd under parity
(P) and time reversal (T) but invariant under the combined PT action. The
continuum equations for electromagnetic waves passing through this medium
contain a dispersive parity-breaking term similar to the odd viscosity term
found in fluids.

We emphasize that, unlike usual non-reciprocal media, which are not
Hamiltonian, our model system is energy-conserving and Hamiltonian. The
presence of the second-order temporal derivatives of the coupled equations for
the polarization and electric field results in four branches for the
dispersion relation. With a cubic nonlinearity in the polarization equation
combined with the parity-breaking effects, we derive the nonlinear Schrödinger
equation (NLSE) using the method of multiple scales and analyze the stability
criterion as a function of the parameters and the wavenumber. In contrast to
the NLSE for odd gravity waves, the stability criterion in the optics case
depends only on the form of the dispersion relation. The non-reciprocity leads
to stability for large wavenumbers, where the effects of parity breaking are
most dominant (see Fig. 9).

In conclusion, we believe the phenomena discussed in the work are likely
observable in other fields of physics with the interplay of parity breaking
and MI being natural, such as astrophysics, plasma physics, non-reciprocal
metamaterials, and Bose-Einstein condensates of dipolar gases. As such, our
work could be the start of a direction in the investigation of chirality-
dependent instabilities and their consequences in a wide class of systems.

##  6 Acknowledgements

We thank Vinod Menon for motivating this work. This work is supported by NSF
CAREER Grant No. DMR-1944967 (SG and SS). SS and SG were also supported in
part by NSF OMA-1936351.

##  7 Appendix

###  7.1 Details for envelope solutions in parity-odd fluid

The deep-water surface system Eq.(11)-Eq.(14), can be expressed in terms of
envelope variables, defined in Eq.(15), that is,

$\sigma=kx-\omega
t\,,\qquad\xi=\epsilon(x-c_{g}t)\,,\qquad\tau=\epsilon^{2}t\,.$

In these variables, the bulk equation $(-\infty\leq\chi\leq 0)$ becomes

$\displaystyle
k^{2}\theta_{\sigma\sigma}+(1+k^{2}\eta_{\sigma}^{2})\theta_{\chi\chi}-2k^{2}\partial_{\sigma}(\eta_{\sigma}\theta_{\chi})+k^{2}\eta_{\sigma\sigma}\theta_{\chi}+2k\epsilon\big{[}\theta_{\sigma\xi}+\eta_{\sigma}\eta_{\xi}\theta_{\chi\chi}$
$\displaystyle-\partial_{\xi}(\eta_{\sigma}\theta_{\chi})-\eta_{\xi}\theta_{\sigma\chi}\big{]}+\epsilon^{2}\left[\theta_{\xi\xi}+\eta_{\xi}^{2}\theta_{\chi\chi}-2\partial_{\xi}(\eta_{\xi}\theta_{\chi})+\eta_{\xi\xi}\theta_{\chi}\right]=0,$
(57)

whereas the kinematic boundary condition can be written as

$\displaystyle
k^{2}\eta_{\sigma}\theta_{\sigma}-\omega\eta_{\sigma}-(1+k^{2}\eta_{\sigma}^{2})\theta_{\chi}+\epsilon\big{[}k(\eta_{\sigma}\theta_{\xi}+\eta_{\xi}\theta_{\sigma})$
$\displaystyle-
c_{g}\eta_{\xi}-2k\eta_{\xi}\eta_{\sigma}\theta_{\chi}\big{]}+\epsilon^{2}\left[\eta_{\tau}+\eta_{\xi}\theta_{\xi}-\eta_{\xi}^{2}\theta_{\chi}\right]=0\,,\qquad\qquad\chi=0.$
(58)

On the other hand, the no-penetration condition at the infinity does not get
modified, and it remains

$\displaystyle\theta_{\chi}=0\,,\qquad\qquad\qquad\qquad\chi=-\infty.$ (59)

In addition to these, we have the Bernoulli equation taken at the surface,
where the effects of odd viscosity are encoded as the modified dynamical
pressure term. This is in contrast to the ideal fluids case where the pressure
is simply zero. In terms of the envelope variables, this becomes

$\displaystyle-\omega(\theta_{\sigma}-\eta_{\sigma}\theta_{\chi})+\tfrac{1}{2}\big{[}k^{2}\theta_{\sigma}^{2}+(1+k^{2}\eta_{\sigma}^{2})\theta_{\chi}^{2}-2k^{2}\eta_{\sigma}\theta_{\sigma}\theta_{\chi}\big{]}+g\eta+\epsilon\big{[}c_{g}(\eta_{\xi}\theta_{\chi}\\-\theta_{\xi})+k\theta_{\sigma}\theta_{\xi}$
$\displaystyle-k\theta_{\chi}\left(\eta_{\sigma}\theta_{\xi}+\eta_{\xi}\theta_{\sigma}-\eta_{\sigma}\eta_{\xi}\theta_{\chi}\right)\big{]}+\epsilon^{2}\left[\theta_{\tau}-\eta_{\tau}\theta_{\chi}+\tfrac{1}{2}\left(\theta^{2}_{\xi}+\eta-\xi^{2}\theta_{\chi}^{2}-2\eta_{\xi}\theta_{\xi}\theta_{\chi}\right)\right]=$
$\displaystyle-\frac{2\nu_{o}}{\sqrt{1+(k\eta_{\sigma}+\epsilon\eta_{\xi})^{2}}}(k\partial_{\sigma}+\epsilon\partial_{\xi})\left[\frac{\omega\eta_{\sigma}+c_{g}\eta_{\xi}-\epsilon^{2}\eta_{\tau}}{\sqrt{1+(k\eta_{\sigma}+\epsilon\eta_{\xi})^{2}}}\right],\qquad\qquad\quad\chi=0\,.$
(60)

Even though the odd viscosity term in the last line looks horrendous, it
simplifies quite substantially when we consider that
$\eta,\theta\sim\mathcal{O}(\epsilon)$. We seek a perturbative expansion for
this system of equations of the type:

$\eta=\sum_{j=1}^{\infty}\,\epsilon^{j}\sum_{n=-j}^{j}\eta_{jn}q^{n}\,,\qquad\qquad\theta=\sum_{j=1}^{\infty}\,\epsilon^{j}\sum_{n=-j}^{j}\theta_{jn}q^{n}\,,$
(61)

which will lead to the envelope dynamics for the slow degrees of freedom, as
prescribed by Eq.(18). The equations at order $\mathcal{O}(\epsilon q)$ can be
written as

$\displaystyle\partial_{\chi}^{2}\theta_{11}$
$\displaystyle=k^{2}\theta_{11},\quad-\infty\leq\chi\leq 0$ (62)
$\displaystyle\partial_{\chi}\theta_{11}$
$\displaystyle=0,\qquad\quad\quad\chi=-\infty$ (63)
$\displaystyle-i\omega\eta_{11}$
$\displaystyle=\partial_{\chi}\theta_{11}\,\,\,\qquad\chi=0$ (64)
$\displaystyle-i\omega\theta_{11}+g\eta_{11}$ $\displaystyle=2\nu_{o}\omega
k\eta_{11},\quad\chi=0$ (65) $\displaystyle\theta_{10}$
$\displaystyle=A_{0}(\xi,\tau)\quad-\infty\leq\chi\leq 0$ (66)
$\displaystyle\eta_{10}$ $\displaystyle=0$ (67)

From here on, we will consider $k>0$. In fact, the analysis for $k<0$, can be
recovered under the replacement $K\rightarrow-k$ and
$\omega\rightarrow-\omega$. The first two equations can be satisfied by the
solution

$\displaystyle\theta_{11}=A(\xi,\tau)e^{k\chi},\qquad\eta_{11}=\frac{i\omega}{g-2\nu_{o}\omega
k}A(\xi,\tau)=i(k/\omega)A(\xi,\tau)\,$ (68)

where $A$ is the slowly-varying unknown carrier wave-train amplitude. The
dispersion of the general system is also obtained in this order and gives us

$\displaystyle\omega^{2}+2\nu_{o}k^{2}\omega-gk=0\,.$ (69)

Since we have fixed $k>0$, this leads to the two sectors $\omega_{\pm}$ (right
and left movers), be defined by the two branches of the above dispersion.

$\displaystyle\omega_{\pm}=-\nu_{o}k^{2}\pm\sqrt{gk+\nu_{o}^{2}k^{4}}$ (70)

In the following, we use the general $\omega$ which has two possible solutions
$\omega_{\pm}$ for $k>0$. In the final NLSE, we can analyze the dynamics for
both solutions. We solve for $\mathcal{O}(\epsilon^{2})$ equations at the free
surface. The higher order bulk involves more unknown functions. Comparing the
harmonics, we obtain the following equations for the bulk solution:

$\displaystyle\mathcal{O}(\epsilon^{2}q^{0}):\quad\partial_{\chi}^{2}\theta_{20}=0\,,$
(71)
$\displaystyle\mathcal{O}(\epsilon^{2}q^{1}):\quad\partial_{\chi}^{2}\theta_{21}-k^{2}\theta_{21}+2ik\partial_{\xi}Ae^{k\chi}=0\,,$
(72)
$\displaystyle\mathcal{O}(\epsilon^{2}q^{2}):\quad\partial_{\chi}^{2}\theta_{22}-4k^{2}\theta_{22}+\frac{3ik^{4}A^{2}e^{k\chi}}{\omega}=0\,.$
(73)

The associated solutions for the bulk subject to the boundary condition
$\partial_{\chi}\theta_{jn}\Big{|}_{\chi=-\infty}=0$ at this order is given by

$\displaystyle\theta_{21}$
$\displaystyle=\left[A_{2}(\xi,\tau)e^{k\chi}-i(\partial_{\xi}A)\chi
e^{k\chi}\right],$ (74) $\displaystyle\theta_{22}$
$\displaystyle=A_{3}(\xi,\tau)e^{2k\chi}+i\left(\frac{k^{2}A^{2}}{\omega}\right)e^{k\chi}\,,$
(75) $\displaystyle\theta_{20}$ $\displaystyle=A_{4}(\xi,\tau)\,,$ (76)

where $A_{2},A_{3},A_{4}$ are unknown functions of $\xi$ and $\tau$. Using the
boundary conditions we can express them in terms of the linear order envelope
$A(\xi,\tau)$.

The zeroth harmonic trivially satisfies the kinematic boundary at this order.
Comparing the first harmonic and second harmonic at this order, we obtain

$\displaystyle\mathcal{O}(\epsilon^{2}q^{1}):-i\omega\eta_{21}=i\left(\frac{c_{g}k-\omega}{\omega}\right)\partial_{\xi}A(\xi,\tau)+kA_{2}\,,$
(77)
$\displaystyle\mathcal{O}(\epsilon^{2}q^{2}):-2i\omega\eta_{22}=\left(\frac{2ik^{3}}{\omega}\right)A^{2}+2kA_{3}(\xi,\tau)\,.$
(78)

From the dynamical boundary condition, we obtain

$\displaystyle\mathcal{O}(\epsilon^{2}q^{0}):\quad
c_{g}\partial_{\xi}A_{0}=g\,\eta_{20}\,,$ (79)
$\displaystyle\mathcal{O}(\epsilon^{2}q^{1}):\quad(g-2\nu_{o}\omega
k)\eta_{21}=i\omega
A_{2}+\left(\frac{2c_{g}\nu_{o}k^{2}}{\omega}+2\nu_{o}k+c_{g}\right)\partial_{\xi}A\,,$
(80) $\displaystyle\mathcal{O}(\epsilon^{2}q^{2}):\quad(g-8\nu_{o}\omega
k)\eta_{22}=2i\omega A_{3}-k^{2}A^{2}\,.$ (81)

The equation involving $\eta_{21}$ and $A_{2}$ are not independent and
consequently, we only obtain the following combination as a condition

$\displaystyle\omega\,\eta_{21}-ikA_{2}=\left(1-\frac{c_{g}k}{\omega}\right)\partial_{\xi}A\,.$
(82)

The higher-order equations depend only on the aforementioned combination or
its derivatives. Solving the equations appearing in the second harmonic, we
obtain,

$\displaystyle\eta_{22}=\frac{-k^{3}}{\omega(\omega+6\nu_{o}k^{2})}A^{2},\qquad
A_{3}=\frac{-6i\nu_{o}k^{4}}{\omega(\omega+6\nu_{o}k^{2})}A^{2}.$ (83)

For one of the branches, we observe the second-harmonic resonance in the
solution $\omega_{-}(k)+6\nu_{o}k^{2}=0$ due to interactions between the
gravity waves and odd surface waves. At the resonance value of $k$, the
present counting scheme breaks down and one needs to choose a different
counting scheme. The remainder of this analysis assumes that we are away from
this resonance point. At $\mathcal{O}(\epsilon^{3}q^{0})$ we get for the bulk

$\displaystyle\mathcal{O}(\epsilon^{3}q^{0}):\partial_{\chi\chi}^{2}\theta_{30}+\partial_{\xi\xi}^{2}\theta_{10}+\left(2ik\eta_{1-1}\partial_{\chi\xi}^{2}\theta_{11}+c.c.\right)-k^{2}\left(\eta_{2-1}\partial_{\chi}\theta_{11}+\eta_{11}\partial_{\chi}\theta_{2-1}+c.c.\right)=0,$

where $c.c.$ refers to complex conjugation. Solving the inhomogeneous bulk
equation, we get

$\displaystyle\theta_{30}=\left(\frac{-\partial_{\xi\xi}^{2}A_{0}}{2}\right)\chi^{2}+\left(\frac{k^{2}c_{g}}{\omega^{2}}+\frac{k^{2}\chi}{\omega}\right)e^{k\chi}\partial_{\xi\xi}^{2}|A|^{2}\,.$
(84)

However, the boundary condition imposes ($\partial_{\chi}\theta_{30}=0$ as
$\chi=-\infty$) leading to the constraint $A_{0}=0$. We further need the bulk
equation at $\mathcal{O}(\epsilon^{3}q^{1})$, which will be used to eliminate
$\theta_{31}$ fields in the surface equations

$\displaystyle\mathcal{O}(\epsilon^{3}q^{1}):\partial_{\chi\chi}^{2}\theta_{31}-k^{2}\theta_{31}=\frac{-6ik^{4}}{\omega}A^{*}A_{3}e^{2k\chi}-2ik\partial_{\xi}A_{2}e^{k\chi}-\partial_{\xi\xi}^{2}A(2k\chi+1)e^{k\chi}\,.$
(85)

Solving this bulk equation gives us another unknown function $A_{5}$ and the
full solution is given by

$\displaystyle\theta_{31}=A_{5}e^{k\chi}-\frac{12\nu_{o}k^{6}}{\omega^{2}(\omega+6\nu_{o}k^{2})}A|A|^{2}e^{2k\chi}-i\partial_{\xi}A_{2}\chi
e^{k\chi}-\chi^{2}e^{k\chi}\partial_{\xi\xi}^{2}A/2.$ (86)

The last set of equations we need to consider is the two boundary conditions
at $\mathcal{O}(\epsilon^{3}q^{1})$. After simplifying these equations using
the lower order solutions, we obtain

$\displaystyle-\frac{\omega^{2}}{k}\eta_{31}+i\omega
A_{5}+\partial_{\tau}A+i\frac{c_{g}\omega}{k}\partial_{\xi}\eta_{21}+\frac{\omega}{k}\partial_{\xi}A_{2}-\left(\frac{ik^{4}}{\omega}-\frac{4ik^{4}}{6k^{2}\nu_{o}+\omega}\right)A|A|^{2}=0\,,$
(87) $\displaystyle\frac{\omega^{2}}{k}\eta_{31}-i\omega
A_{5}+\left(1+\frac{2\nu_{o}k^{2}}{\omega}\right)\partial_{\tau}A+2i\nu_{o}(\omega+kc_{g})\partial_{\xi}\eta_{21}-c_{g}\partial_{\xi}A_{2}+\frac{2i\nu_{o}kc_{g}}{\omega}\partial_{\xi\xi}^{2}A$
$\displaystyle+ik^{4}\left(\frac{4\nu_{o}k^{2}}{\omega^{2}}+\frac{12\nu_{o}k^{2}}{\omega(\omega+6k^{2}\nu_{o})}+\frac{1}{\omega}\right)A|A|^{2}=0\,.$
(88)

We can add the equations above to simultaneously eliminate $\eta_{3,1}$ and
$A_{5}$. Further using the linear combination Eq.(82), and the dispersion, we
find a consistency condition in the form of the non-linear Schrödinger
equation for the envelope as discussed in the main text,

$\displaystyle i\partial_{\tau}A+L\partial_{\xi\xi}^{2}A-MA|A|^{2}=0.$ (89)

where

$\displaystyle L=\frac{1}{2}\frac{d^{2}\omega}{dk^{2}},\qquad
M=\frac{2k^{4}}{\omega}\left(1-\frac{3\nu_{o}k^{2}\omega}{(\omega+\nu_{o}k^{2})(\omega+6\nu_{o}k^{2})}\right).$
(90)

This is the final result that is presented in the main text.

###  7.2 Envelope dynamics in the parity-odd dielectric medium

In this section, we present the calculation details for the envelope dynamics
in the 1D light-matter dielectric system with broken parity. The starting
point is the set of following equations of motion:

$\displaystyle P_{tt}+\beta_{\text{odd}}P_{tx}+\omega_{0}^{2}P+\alpha
P^{3}=\frac{\omega_{P}^{2}}{4\pi}E\,$ $\displaystyle c^{2}E_{xx}-E_{tt}=4\pi
P_{tt}\,.$

In terms of the multiple scales variables introduce in Eq. (15), that is,

$\sigma=kx-\omega
t\,,\quad\xi=\epsilon(x-c_{g}t)\,,\quad\tau=\epsilon^{2}t\,,$

these dynamical equations become

$\displaystyle(\omega^{2}-c^{2}k^{2})E_{\sigma\sigma}+4\pi\omega^{2}P_{\sigma\sigma}+\epsilon\left[2(c_{g}\,\omega-c^{2}k)E_{\sigma\xi}+8\pi
c_{g}\omega P_{\sigma\xi}\right]+\epsilon^{2}[-2\omega E_{\sigma\tau}$
$\displaystyle+(c_{g}^{2}-c^{2})E_{\xi\xi}-8\pi\omega P_{\sigma\tau}+4\pi
c_{g}^{2}P_{\xi\xi}]+O(\epsilon^{3})=0\,,$ (91)
$\displaystyle(\omega^{2}-\beta_{\text{odd}}\,\omega
k)P_{\sigma\sigma}+\omega_{0}^{2}P+\alpha
P^{3}-\frac{\omega_{P}^{2}}{4\pi}E+\epsilon\,[2\omega
c_{g}-\beta_{\text{odd}}(\omega+kc_{g})]P_{\sigma\xi}$
$\displaystyle+\epsilon^{2}[(\beta_{\text{odd}}k-2\omega)P_{\sigma\tau}+c_{g}(c_{g}-\beta_{\text{odd}})P_{\xi\xi}]+O(\epsilon^{3})=0\,.$
(92)

Using the intuition from the hydrodynamic case, we only need terms up to
$\mathcal{O}(\epsilon^{2})$ in this expansion, since this will produce all the
terms at $\mathcal{O}(\epsilon^{3})$ when we impose that
$P,E\sim\mathcal{O}(\epsilon)$. Expanding $E$ and $P$ perturbatively as

$\displaystyle E=\sum_{j=1}^{\infty}\epsilon^{j}E_{j}\,,\qquad
P=\sum_{j=1}^{\infty}\epsilon^{j}P_{j},$

and arranging the system order by order in $\epsilon$ we get

$\displaystyle\mathcal{O}(\epsilon):\qquad$
$\displaystyle(\omega^{2}-c^{2}k^{2})E_{1\sigma\sigma}+4\pi\omega^{2}P_{1\sigma\sigma}=0,$
(93) $\displaystyle(\omega^{2}-\beta_{\text{odd}}\omega
k)P_{1\sigma\sigma}+\omega_{0}^{2}P_{1}-\frac{\omega_{p}^{2}}{4\pi}E_{1}=0\,,$
(94) $\displaystyle\mathcal{O}(\epsilon^{2}):\qquad$
$\displaystyle(\omega^{2}-c^{2}k^{2})E_{2\sigma\sigma}+4\pi\omega^{2}P_{2\sigma\sigma}+2(\omega
c_{g}-c^{2}k)E_{1\sigma\xi}+8\pi\omega c_{g}P_{1\sigma\xi}=0,$ (95)
$\displaystyle(\omega^{2}-\beta_{\text{odd}}\,\omega
k)P_{2\sigma\sigma}+\omega_{0}^{2}P_{2}-\frac{\omega_{p}^{2}}{4\pi}E_{2}+(2\omega
c_{g}-\beta_{\text{odd}}(kc_{g}+\omega))P_{1\sigma\xi}=0\,,$ (96)
$\displaystyle\mathcal{O}(\epsilon^{3}):\qquad$
$\displaystyle(\omega^{2}-c^{2}k^{2})E_{3\sigma\sigma}+4\pi\omega^{2}P_{3\sigma\sigma}+2(\omega
c_{g}-c^{2}k)E_{2\sigma\xi}+8\pi\omega c_{g}P_{2\sigma\xi}$
$\displaystyle-2\omega E_{1\sigma\tau}+(c_{g}^{2}-c^{2})E_{1\xi\xi}-8\pi\omega
P_{1\sigma\tau}+4\pi c_{g}^{2}P_{1\xi\xi}=0,$ (97)
$\displaystyle(\omega^{2}-\beta_{\text{odd}}\,\omega
k)P_{3\sigma\sigma}-\frac{\omega_{p}^{2}}{4\pi}E_{3}+(2\omega
c_{g}-\beta_{\text{odd}}\,(kc_{g}+\omega))P_{2\sigma\xi}$
$\displaystyle+(\beta_{\text{odd}}\,k-2\omega)P_{1\sigma\tau}+c_{g}(c_{g}-\beta_{\text{odd}})P_{1\xi\xi}+\omega_{0}^{2}P_{3}+\alpha
P_{1}^{3}=0.$ (98)

Let us we further expanding each $\epsilon$-order into harmonic modes as

$\displaystyle E_{j}=\sum_{n=-j}^{j}E_{jn}(\xi,\tau)e^{in\sigma}\,,\qquad
P_{j}=\sum_{n=-j}^{j}P_{jn}(\xi,\tau)e^{in\sigma}.$

Due to reality condition of $E$ and $P$, we have $E_{j,-n}=E_{j,n}^{*}$ and
$P_{j,-n}=P_{j,n}^{*}$. Note that under this scheme, higher harmonics can only
be generated perturbatively. Therefore, the leading order equations, that is,
Eqs. (93-94) can be written as

$\displaystyle\sum_{n=-1}^{1}\left[n^{2}(k^{2}c^{2}-\omega^{2})E_{1n}-4\pi
n^{2}\omega^{2}P_{1n}\right]e^{in\sigma}=0\,,$ (99)
$\displaystyle\sum_{n=-1}^{1}\left[(n^{2}\omega(\beta_{\text{odd}}k-\omega)+\omega_{0}^{2})P_{1n}-\frac{\omega_{p}^{2}}{4\pi}E_{1n}\right]e^{in\sigma}=0\,.$
(100)

The zeroth harmonic mode gives us
$\omega_{0}^{2}P_{10}=\frac{\omega_{p}^{2}}{4\pi}E_{10}$. However, this
corresponds to a uniform and constant electric (polarization) field in the
medium. Since we are interested in the field propagation, we can set
$E_{10}=P_{10}=0$.

The rest of the leading order equations give us the dispersion relation inside
the medium and how the electric field relates to the polarization field:

$\displaystyle(\omega^{2}-c^{2}k^{2})(\omega_{0}^{2}+\beta_{\text{odd}}\,\omega
k-\omega^{2})+\omega^{2}\omega_{p}^{2}=0\,,$ $\displaystyle
E_{11}=\frac{4\pi}{\omega_{p}^{2}}(\omega_{0}^{2}+\beta_{\text{odd}}\,\omega
k-\omega^{2})P_{11}=\frac{4\pi\omega^{2}}{c^{2}k^{2}-\omega^{2}}P_{11}\,.$
(101)

Plugging the solutions of $\mathcal{O}(\epsilon)$ into the equation for order
$\mathcal{O}(\epsilon^{2})$, we find that the equations do not depend on
$E_{11}$, $P_{11}$. Since nonlinear terms only appear in the third-order
equations, the terms $P_{2}$ and $E_{2}$ do not contribute for the NLSE and
can be set to zero without any loss of generality. From the same argument, we
can also set $E_{32}=E_{30}=P_{32}=P_{30}=0$, leaving us with the following
equations to consider

$\displaystyle(c^{2}k^{2}-\omega^{2})E_{33}-4\pi\omega^{2}P_{33}=0\,,$ (102)
$\displaystyle
9(\omega\beta_{\text{odd}}\,k-\omega^{2})P_{33}+\omega_{0}^{2}P_{33}-\frac{\omega_{p}^{2}}{4\pi}E_{33}+\alpha
P_{11}^{3}=0\,,$ (103)

which determines the generation of the third harmonics in the medium in terms
of the incident wave envelope, and

$\displaystyle(c^{2}k^{2}-\omega^{2})E_{31}-4\pi\omega^{2}P_{31}-2i\omega
E_{11\tau}+(c_{g}^{2}-c^{2})E_{11\xi\xi}-i8\pi\omega P_{11\tau}+4\pi
c_{g}^{2}P_{11\xi\xi}=0\,,$ (104)
$\displaystyle\frac{\omega^{2}\omega_{P}^{2}}{c^{2}k^{2}-\omega^{2}}P_{31}-\frac{\omega_{p}^{2}}{4\pi}E_{31}+3\alpha
P_{11}^{2}P_{1,-1}+i(\beta_{\text{odd}}k-2\omega)P_{11\tau}+c_{g}(c_{g}-\beta_{\text{odd}})P_{11\xi\xi}=0\,.$
(105)

After some algebra, the set of Eqs. (104-105) can be expressed in form of the
NLSE:

$\displaystyle\left(i\partial_{\tau}+\frac{1}{2}\frac{d^{2}\omega}{dk^{2}}\partial_{\xi}^{2}\right)P_{11}-\frac{3\alpha\omega^{3}\omega_{p}^{2}}{(\omega^{2}\omega_{p}^{2})(\omega^{2}+\omega_{0}^{2})+(\omega_{0}^{2}+\beta_{\text{odd}}\omega
k-\omega^{2})^{2}(\omega^{2}+c^{2}k^{2})}|P_{11}|^{2}P_{11}=0\,.$ (106)

Since $\alpha>0$, the envelope solution is stable to MI, in zero bandwidth
limit, when:

$\displaystyle\omega\frac{d^{2}\omega}{dk^{2}}\geq 0\,.$

## References

  * [1] PC Martin, O Parodi, and Peter S Pershan.  Unified hydrodynamic theory for crystals, liquid crystals, and normal fluids.  Physical Review A, 6(6):2401, 1972.
  * [2] Michel Fruchart, Ming Han, Colin Scheibner, and Vincenzo Vitelli.  The odd ideal gas: Hall viscosity and thermal conductivity from non-hermitian kinetic theory.  2022\.
  * [3] Tomer Markovich and Tom C. Lubensky.  Odd viscosity in active matter: microscopic origin and 3d effects, 2020\.
  * [4] Michel Fruchart, Colin Scheibner, and Vincenzo Vitelli.  Odd viscosity and odd elasticity.  Annual Review of Condensed Matter Physics, 14:471–510, 2023.
  * [5] LD Landau.  Em lifshitz fluid dynamics, 1987.
  * [6] JE Avron.  Odd viscosity.  Journal of statistical physics, 92(3-4):543–557, 1998.
  * [7] Sriram Ganeshan and Alexander G Abanov.  Odd viscosity in two-dimensional incompressible fluids.  Physical Review Fluids, 2(9):094101, 2017.
  * [8] JE Avron, Ruedi Seiler, and Petr G Zograf.  Viscosity of quantum hall fluids.  Physical review letters, 75(4):697, 1995.
  * [9] IV Tokatly.  Magnetoelasticity theory of incompressible quantum hall liquids.  Physical Review B, 73(20):205340, 2006.
  * [10] IV Tokatly and Giovanni Vignale.  New collective mode in the fractional quantum hall liquid.  Physical review letters, 98(2):026805, 2007.
  * [11] FDM Haldane.  Geometrical description of the fractional quantum hall effect.  Physical review letters, 107(11):116801, 2011.
  * [12] Carlos Hoyos and Dam Thanh Son.  Hall viscosity and electromagnetic response.  Physical review letters, 108(6):066805, 2012.
  * [13] Barry Bradlyn, Moshe Goldstein, and N Read.  Kubo formulas for viscosity: Hall viscosity, ward identities, and the relation with conductivity.  Physical Review B, 86(24):245309, 2012.
  * [14] Alexander G Abanov.  On the effective hydrodynamics of the fractional quantum hall effect.  Journal of Physics A: Mathematical and Theoretical, 46(29):292001, 2013.
  * [15] Carlos Hoyos.  Hall viscosity, topological states and effective theories.  International Journal of Modern Physics B, 28(15):1430007, 2014\.
  * [16] M Laskin, T Can, and P Wiegmann.  Collective field theory for quantum hall states.  Physical Review B, 92(23):235141, 2015.
  * [17] Tankut Can, Michael Laskin, and Paul B Wiegmann.  Geometry of quantum hall states: Gravitational anomaly and transport coefficients.  Annals of Physics, 362:752–794, 2015.
  * [18] Semyon Klevtsov and Paul Wiegmann.  Geometric adiabatic transport in quantum hall states.  Physical review letters, 115(8):086801, 2015.
  * [19] Thomas Scaffidi, Nabhanila Nandi, Burkhard Schmidt, Andrew P Mackenzie, and Joel E Moore.  Hydrodynamic electron flow and hall viscosity.  Physical review letters, 118(22):226601, 2017.
  * [20] Francesco MD Pellegrino, Iacopo Torre, and Marco Polini.  Nonlocal transport and the hall viscosity of two-dimensional hydrodynamic electron liquids.  Physical Review B, 96(19):195401, 2017.
  * [21] AI Berdyugin, SG Xu, FMD Pellegrino, R Krishna Kumar, Alessandro Principi, Iacopo Torre, M Ben Shalom, Takashi Taniguchi, Kenji Watanabe, IV Grigorieva, et al.  Measuring hall viscosity of graphene’s electron fluid.  Science, page eaau0685, 2019.
  * [22] Gabriel S. Denicol, Xu-Guang Huang, Etele Molnár, Gustavo M. Monteiro, Harri Niemi, Jorge Noronha, Dirk H. Rischke, and Qun Wang.  Nonresistive dissipative magnetohydrodynamics from the boltzmann equation in the 14-moment approximation.  Phys. Rev. D, 98:076009, Oct 2018.
  * [23] J Korving, H Hulsman, HFP Knaap, and JJM Beenakker.  Transverse momentum transport in viscous flow of diatomic gases in a magnetic field.  Physics Letters, 21(1):5–7, 1966.
  * [24] Colin Scheibner, Anton Souslov, Debarghya Banerjee, Piotr Surówka, William TM Irvine, and Vincenzo Vitelli.  Odd elasticity.  Nature Physics, 16(4):475–480, 2020.
  * [25] Hussein Nassar, Behrooz Yousefzadeh, Romain Fleury, Massimo Ruzzene, Andrea Alù, Chiara Daraio, Andrew N Norris, Guoliang Huang, and Michael R Haberman.  Nonreciprocity in acoustic and elastic materials.  Nature Reviews Materials, 5(9):667–685, 2020.
  * [26] Debarghya Banerjee, Anton Souslov, Alexander G Abanov, and Vincenzo Vitelli.  Odd viscosity in chiral active fluids.  Nature Communications, 8(1):1573, 2017.
  * [27] Vishal Soni, Ephraim S Bililign, Sofia Magkiriadou, Stefano Sacanna, Denis Bartolo, Michael J Shelley, and William TM Irvine.  The odd free surface flows of a colloidal chiral fluid.  Nature Physics, 15(11):1188–1194, 2019.
  * [28] Gustavo M. Monteiro, Alexander G. Abanov, and Sriram Ganeshan.  Hamiltonian structure of 2d fluid dynamics with broken parity, 2021.
  * [29] Dylan Reynolds, Gustavo M. Monteiro, and Sriram Ganeshan.  Hele-shaw flow for parity odd three-dimensional fluids.  Phys. Rev. Fluids, 7:114201, Nov 2022.
  * [30] Alexander Abanov, Tankut Can, and Sriram Ganeshan.  Odd surface waves in two-dimensional incompressible fluids.  SciPost Physics, 5(1):010, 2018.
  * [31] E Kirkinis and AV Andreev.  Odd-viscosity-induced stabilization of viscous thin liquid films.  Journal of Fluid Mechanics, 878:169–189, 2019.
  * [32] Gustavo M. Monteiro and Sriram Ganeshan.  Non-linear shallow water dynamics with odd viscosity.  arXiv:2009.00118, 2020.
  * [33] Alex Doak, Guido Baardink, Paul A Milewski, and Anton Souslov.  Nonlinear shallow-water waves with vertical odd viscosity.  SIAM Journal on Applied Mathematics, 83(3):938–965, 2023.
  * [34] Vladimir E Zakharov and Lev A Ostrovsky.  Modulation instability: The beginning.  Physica D: Nonlinear Phenomena, 238(5):540–548, 2009.
  * [35] Thomas Brooke Benjamin.  Instability of periodic wavetrains in nonlinear dispersive systems.  Proceedings of the Royal Society of London. Series A. Mathematical and Physical Sciences, 299(1456):59–76, 1967.
  * [36] T Brooke Benjamin and James E Feir.  The disintegration of wave trains on deep water part 1. theory.  Journal of Fluid Mechanics, 27(3):417–430, 1967.
  * [37] NF Pilipetskii and AR Rustamov.  Observation of self-focusing of light in liquids.  JETP LETTERS-USSR, 2(2):55–+, 1965.
  * [38] VI Bespalov and VI Talanov.  Filamentary structure of light beams in nonlinear liquids.  Soviet Journal of Experimental and Theoretical Physics Letters, 3:307, 1966.
  * [39] Alexander G. Abanov, Tankut Can, Sriram Ganeshan, and Gustavo M. Monteiro.  Hydrodynamics of two-dimensional compressible fluid with broken parity: Variational principle and free surface dynamics in the absence of dissipation.  Physical Review Fluids, 5(10).
  * [40] V Dj Djordjevic and Larry G Redekopp.  On two-dimensional packets of capillary-gravity waves.  Journal of Fluid Mechanics, 79(4):703–714, 1977.
  * [41] Lawrence F McGoldrick.  Resonant interactions among capillary-gravity waves.  Journal of Fluid Mechanics, 21(2):305–331, 1965.
  * [42] Alexander G. Abanov, Tankut Can, Sriram Ganeshan, and Gustavo M. Monteiro.  Hydrodynamics of two-dimensional compressible fluid with broken parity: variational principle and free surface dynamics in the absence of dissipation.  arXiv preprint arXiv:1907.11196, 2019.
  * [43] Alexander G Abanov and Gustavo M Monteiro.  Free surface variational principle for an incompressible fluid with odd viscosity.  Phys Rev Lett, 122(15):154501, Apr 2019.
  * [44] Hidenori Hasimoto and Hiroaki Ono.  Nonlinear modulation of gravity waves.  Journal of the Physical Society of Japan, 33(3):805–811, 1972.
  * [45] Lawrence F McGoldrick.  An experiment on second-order capillary gravity resonant wave interactions.  Journal of Fluid Mechanics, 40(2):251–271, 1970.
  * [46] A Davey and Keith Stewartson.  On three-dimensional packets of surface waves.  Proceedings of the Royal Society of London. A. Mathematical and Physical Sciences, 338(1613):101–110, 1974.