#  Square Pattern Formation as Stable Fixed Point in Driven Two-Dimensional
Bose-Einstein Condensates

Keisuke Fujii<EMAIL_ADDRESS>Institut für Theoretische
Physik, Universität Heidelberg, Philosophenweg 19, 69120 Heidelberg, Germany
Sarah L. Görlitz  Institut für Theoretische Physik, Universität Heidelberg,
Philosophenweg 19, 69120 Heidelberg, Germany  Nikolas Liebster  pattern-
<EMAIL_ADDRESS>Kirchhoff-Institute für Physik, Universität
Heidelberg, Im Neuenheimer Feld 227, 69120 Heidelberg, Germany  Marius Sparn
Kirchhoff-Institute für Physik, Universität Heidelberg, Im Neuenheimer Feld
227, 69120 Heidelberg, Germany  Elinor Kath  Kirchhoff-Institute für Physik,
Universität Heidelberg, Im Neuenheimer Feld 227, 69120 Heidelberg, Germany
Helmut Strobel  Kirchhoff-Institute für Physik, Universität Heidelberg, Im
Neuenheimer Feld 227, 69120 Heidelberg, Germany  Markus K. Oberthaler
Kirchhoff-Institute für Physik, Universität Heidelberg, Im Neuenheimer Feld
227, 69120 Heidelberg, Germany  Tilman Enss  Institut für Theoretische Physik,
Universität Heidelberg, Philosophenweg 19, 69120 Heidelberg, Germany

###### Abstract

We investigate pattern formation in two-dimensional Bose-Einstein condensates
(BECs) caused by temporal periodic modulation of the interatomic interaction.
Temporal modulation of the interaction causes the so-called Faraday
instability in the condensate, which we show generically leads to a stable
square grid density pattern. We take the amplitudes in each of the two
directions spanning the two-dimensional density pattern as order parameters in
pattern formation and derive a set of simultaneous time evolution equations
for those order parameters from the Gross–Pitaevskii (GP) equation with a
time-periodic interaction. We identify the fixed points of the time evolution
and show by stability analysis that the inhomogeneous density exhibits a
square grid pattern as a stable fixed point.

Introduction.—Spontaneous pattern formation is a phenomenon in which a uniform
state loses its stability and becomes inhomogeneous as external parameters are
varied. As a spontaneous translational symmetry breaking in nonequilibrium
dynamics, pattern formation appears in nature at diverse scales, not only in
physics [1, 2] but also in chemical reactions [3] and biology [4].
Understanding what kinds of patterns are formed and how they are formed is a
fundamental interdisciplinary question common to these studies.

In BECs, temporal modulations of system parameters, such as the magnitudes of
trapping potentials and interactions, cause parametric instability. This
instability, called the Faraday instability as an analogy to a similar
phenomenon in classical fluids [5, 6], results in spontaneous pattern
formation [7]. Since ultracold atomic systems provide BECs with experimentally
controllable parameters, one-dimensional patterns were observed in one-
dimensional BECs [8, 9] and as surface waves in two-dimensional BECs [10].
Theoretically, parametric instabilities appearing in driven quantum gases,
including simple BECs, have been intensively studied [11, 12, 13, 14, 15, 16,
17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35,
36, 37, 38, 39, 40, 41, 42, 43, 44, 45].

For two-dimensional Faraday patterns, it was observed that the symmetry of
selected patterns can be engineered by combining multiple temporal modulations
of the interactions [46]. Furthermore, in our companion work [47], we have
observed two-dimensional Faraday patterns, clearly exhibiting a square grid,
with a single-frequency temporal modulation of the interaction. The
realization of two-dimensional Faraday patterns opens the stage for
investigating nonlinear effects, such as the correlations that emerge between
instability-induced excitations in different directions in the plane.

In this Letter, we derive time-evolution equations for the amplitudes of
density waves in two directions spanning two-dimensional density patterns in
two-dimensional BECs due to the Faraday instability. We show that
inhomogeneous condensates exhibit a square grid as a stable pattern. In our
model, we consider density modulations in two directions $\bm{k}$ and $\bm{p}$
in the plane with amplitudes $R_{k}$ and $R_{p}$, respectively. When the angle
between the two directions is near $\pi/2$, both amplitudes grow to finite
values and we find that the system exhibits a stable grid pattern. Conversely,
for small angles, only one of the amplitudes grows to a finite value while the
other is suppressed, and the BEC exhibits a stripe pattern. This result is
clearly seen in Fig. 1, where the global stability of the patterns differs
significantly depending on the angle between the two directions. In the
experiment, many modes at different angles will initially grow due to the
instability, and our analysis shows that two of these modes with an angle
close to $\pi/2$ between them will reinforce each other and grow into a grid
pattern, while other modes at small angles are suppressed.

Figure 1: Global stability of the patterns formed by two standing waves in a
planar BEC in directions $\bm{k}$ and $\bm{p}$. Schematic figures at the fixed
points represent corresponding stationary solutions, i.e., the grid-pattern,
stripe-pattern, and uniform solutions. The angle $\theta$ between the two
excited modes is $\pi/6$ (left) and $\pi/2$ (right). In the latter case, the
square grid pattern emerges as a stable fixed point. Parameters are
$\omega/\mu=2$ and $A=0.6$ with small dissipation $\Gamma=0.1\alpha$, where
$\omega,\mu$ and $\alpha$ are the driving frequency, the chemical potential,
and the drive amplitude for the Bogoliubov modes, respectively.

BEC with a time-periodically modulated interaction.—We consider a two-
dimensional BEC with an interaction strength $g(t)=\bar{g}[1-A\sin(\omega
t)]$, which is modulated periodically in time with drive amplitude $|A|<1$
around its mean value $\bar{g}$. The dynamics of a BEC with wave function
$\Psi(t,\bm{x})$ is described by the Gross–Pitaevskii (GP) equation [48]:

$i\frac{\partial}{\partial
t}\Psi(t,\bm{x})=\biggl{[}-\frac{\nabla^{2}}{2m}+g(t)|\Psi(t,\bm{x})|^{2}\biggr{]}\Psi(t,\bm{x}).$
(1)

We assume that the BEC is confined to a sufficiently large flat potential and
dropped the trapping potential term; in experiment [47] the potential has
absorptive boundaries, which mimics an infinitely extended system.

Equation (1) has a uniform solution,
$\Psi_{\textrm{uni}}(t)=\Psi_{0}\exp[-i\mu t-i(\mu/\omega)A\cos(\omega t)]$
with the chemical potential $\mu=\bar{g}|\Psi_{0}|^{2}$, but this solution
becomes unstable due to the Faraday instability induced by the oscillating
interaction [7]. This instability can be understood to be caused by an
anomalous amplification of excited modes with wavevector $\bm{k}$ satisfying
the resonance condition $n\omega/2=E_{\bm{k}}$ for $n\in\mathbb{N}$. Here,
$E_{\bm{k}}=\sqrt{\epsilon_{\bm{k}}(\epsilon_{\bm{k}}+2\mu)}$ and
$\epsilon_{\bm{k}}=\bm{k}^{2}/(2m)$ represent the Bogoliubov quasiparticle and
single-particle dispersions, respectively. The resonance condition
$n\omega/2=E_{\bm{k}}$ comes from the fact that the energy quantum $n\omega$,
injected into the system by the oscillation, is split between two
quasiparticle excitations with wavevectors $\pm\bm{k}$. Within a linear
analysis, one can indeed derive Mathieu’s differential equation from Eq. (1),
which shows the amplification of modes with wavenumbers around the resonance
condition [49].

The amplitude equation.—In two-dimensional systems, the density distribution
resulting from the instability exhibits a grid pattern spanned by two non-
parallel wavevectors. The realized grid pattern is determined by the
competition between the linear instability due to the parametric resonance and
the nonlinear suppression from the interaction between the excited modes. We
determine the magnitude of the amplitude and the angle of the realized grid,
assuming that the drive amplitude $A$ is small. The small amplitude $A$
weakens the instability and slows the amplitude growth of the density pattern.
As a result, the time evolution of the pattern amplitude is systematically
obtained as slow-timescale dynamics. Using the multiple-scale method [1, 2],
we derive the time-evolution equation for the pattern amplitude from Eq. (1).

To investigate two-dimensional patterns spanned by two wavevectors, we expand
the wave function as 111 As discussed in the last paragraph of the discussion
section, there is no relevant three-mode scattering in our system due to the
constraints imposed by the rotating-wave approximation, and the ansatz with
two excited modes, as in Eq. (2), is sufficient.

$\Psi(t,\bm{x})=\Psi_{\textrm{uni}}(t)\biggl{[}1+\phi_{k}(t)\cos(\bm{k}\cdot\bm{x})+\phi_{p}(t)\cos(\bm{p}\cdot\bm{x})\biggr{]}.$
(2)

The small drive amplitude $A$ maintains the excitation $\phi_{k/p}$ in the
form of the Bogoliubov basis:

$\displaystyle\phi_{k/p}(t)$
$\displaystyle=\biggl{(}1-\frac{\epsilon_{\bm{k}/\bm{p}}+2\mu}{E_{\bm{k}/\bm{p}}}\biggr{)}R_{k/p}(t)e^{i\omega
t/2}$
$\displaystyle\quad+\biggl{(}1+\frac{\epsilon_{\bm{k}/\bm{p}}+2\mu}{E_{\bm{k}/\bm{p}}}\biggr{)}R^{\ast}_{k/p}(t)e^{-i\omega
t/2},$ (3)

where the complex amplitudes $R_{k/p}(t)$ obey the complex Ginzburg-Landau
equation (for details, see the supplement [51]):

$\displaystyle i\frac{\mathrm{d}}{\mathrm{d}t}R_{k}(t)=\Delta R_{k}(t)-i\alpha
R_{k}^{\ast}(t)+\lambda\Bigl{(}|R_{k}(t)|^{2}R_{k}(t)$
$\displaystyle\quad+c_{1}(\theta)|R_{p}(t)|^{2}R_{k}(t)+c_{2}(\theta)R_{p}(t)^{2}R^{\ast}_{k}(t)\Bigr{)},$
(4)

with detuning $\Delta=\omega/2-E$, drive amplitude for the Bogoliubov mode
$\alpha=\mu A\epsilon/(2E)$, and nonlinearity $\lambda=\mu(5\epsilon+3\mu)/E$.
The same equation holds for $R_{p}(t)$ after exchange of the $\bm{k}$ and
$\bm{p}$ labels. The $c_{2}(\theta)$ term would violate momentum conservation
for traveling wave patterns in the Ginzburg-Landau equation [2], but in our
case of standing waves (2), which are agnostic to the sign of $\bm{p}$,
$c_{2}(\theta)$ arises as a new coupling. In order to focus on the angle of
the realized pattern, we assume that the absolute values of the wavevectors
$\bm{k}$ and $\bm{p}$ are equal, as determined by the resonance condition
$\Delta=0$ for $n=1$, and set $\epsilon_{\bm{k}}=\epsilon_{\bm{p}}=\epsilon$
and $E_{\bm{k}}=E_{\bm{p}}=E$. The coupling coefficients $c_{1}(\theta)$ and
$c_{2}(\theta)$ between modes in different directions are then given as
functions of the angle $\theta\in[0,\pi/2]$ between $\bm{k}$ and $\bm{p}$,

$\displaystyle
c_{1}(\theta)=\frac{\mu}{5\epsilon+3\mu}\biggl{[}4\frac{\epsilon^{2}-\mu^{2}}{\mu\epsilon}+\biggl{(}\frac{2\epsilon+\mu}{\epsilon}\frac{2\epsilon+\mu}{\epsilon_{+}/2+\mu}$
$\displaystyle-\frac{(2\epsilon-\mu)(\epsilon+2\mu)+(2\epsilon^{2}+\mu^{2})\epsilon_{+}/(2\epsilon)}{E^{2}-E^{2}_{+}/4}+(\epsilon_{+}\rightarrow\epsilon_{-})\biggr{)}\biggr{]},$
(5a) $\displaystyle
c_{2}(\theta)=\frac{\mu}{5\epsilon+3\mu}\biggl{[}-2\frac{\epsilon^{2}+3\mu\epsilon+\mu^{2}}{\mu\epsilon}$
$\displaystyle+\frac{2\epsilon+\mu}{\epsilon}\biggl{(}\frac{2\epsilon+\mu}{\epsilon_{+}/2+\mu}+(\epsilon_{+}\to\epsilon_{-})\biggr{)}\biggr{]},$
(5b)

where we introduced $E_{\pm}=\sqrt{\epsilon_{\pm}(\epsilon_{\pm}+2\mu)}$ with
$\epsilon_{+}=\epsilon_{\bm{k}+\bm{p}}=4\epsilon\cos^{2}\frac{\theta}{2}$ and
$\epsilon_{-}=\epsilon_{\bm{k}-\bm{p}}=4\epsilon\sin^{2}\frac{\theta}{2}$
222Since we have introduced the wavevectors $\bm{k}$ and $\bm{p}$ via the
cosines in Eq. (2), the signs of them are irrelevant. Thus, we can take the
angle $\theta$ being in $[0,\pi/2]$, and the case with $\theta\in[\pi/2,\pi]$
is given by replacing $\theta\in[0,\pi/2]$ with $\pi-\theta$..

The coefficient $c_{1}(\theta)$ diverges at the singular angle satisfying
$2E=E_{+}$, but it is an artifact of the ansatz (2) considering only two
modes. We present the complete theory without divergence later and find that
no other (e.g., triangular) patterns appear at the singular angle. But first,
we focus on angles away from the singular angle where this model provides
analytical and quantitatively reliable results.

The fixed points and their stability analysis.—The time-dependent solutions of
the evolution (amplitude) equation trace out trajectories in the four-
dimensional space of the two complex amplitudes. In the following, we analyze
their fixed points and stability. We focus only on the excited modes
satisfying the resonance condition at zero detuning $\Delta=0$, and introduce
dissipation $\Gamma>0$ to capture the suppression from interaction with other
modes besides the $\bm{k}$\- and $\bm{p}$-modes:

$\displaystyle i\frac{\mathrm{d}}{\mathrm{d}t}R_{k}(t)=-i\Gamma
R_{k}(t)-i\alpha R_{k}^{\ast}(t)+\lambda\Bigl{(}|R_{k}(t)|^{2}R_{k}(t)$
$\displaystyle\quad+c_{1}(\theta)|R_{p}(t)|^{2}R_{k}(t)+c_{2}(\theta)R_{p}(t)^{2}R^{\ast}_{k}(t)\Bigr{)}.$
(6)

Setting $\mathrm{d}R_{k}(t)/\mathrm{d}t=0$ in Eq. (6) and similarly for
$R_{p}(t)$, we find four possible fixed-point values of $R_{k}$ and $R_{p}$:

$\displaystyle(\textrm{U})\quad$ $\displaystyle R_{k}=R_{p}=0,$ (7a)
$\displaystyle(\textrm{S}_{k})\quad$ $\displaystyle
R_{k}=\bar{R}e^{i\bar{\eta}},\quad R_{p}=0,$ (7b)
$\displaystyle(\textrm{S}_{p})\quad$ $\displaystyle R_{k}=0,\quad
R_{p}=\bar{R}e^{i\bar{\eta}},$ (7c) $\displaystyle(\textrm{G})\quad$
$\displaystyle
R_{k}=R_{p}=\frac{\bar{R}e^{i\bar{\eta}}}{\sqrt{1+c_{1}(\theta)+c_{2}(\theta)}},$
(7d)

with $\bar{R}^{2}=\sqrt{\alpha^{2}-\Gamma^{2}}/\lambda$ and
$\exp(i\bar{\eta})=(\sqrt{\alpha-\Gamma}+i\sqrt{\alpha+\Gamma})/\sqrt{2\alpha}$.
The fixed points correspond to the following density patterns: (U) a uniform
pattern, (Sk) and (Sp) stripe patterns for each direction, and (G) a grid
pattern (see Fig. 1).

We first investigate the stability of the uniform fixed point (U). For small
$R_{k}(t)$ and $R_{p}(t)$ around (U), only the linear terms of Eq. (6) remain,
and the $\bm{k}$\- and $\bm{p}$-directions become independent. By dividing the
linearized equation into its real and imaginary parts, we directly find the
eigenvalues of the Jacobian (scaling dimensions) $-\Gamma-\alpha$ and
$-\Gamma+\alpha$. Since both $\alpha$ and $\Gamma$ are positive, the fixed
point (U) is unstable for $\alpha>\Gamma$. This corresponds to the Faraday
instability, in which the uniform solution becomes unstable when the drive is
stronger than the dissipation.

We next study the stability of the grid fixed point (G). The four eigenvalues
of the amplitude equation linearized around (G) are found to be

$\displaystyle\Lambda_{1}^{\pm}=-\Gamma\pm i\sqrt{4\alpha^{2}-5\Gamma^{2}},$
(8a)
$\displaystyle\Lambda_{2}^{\pm}=-\Gamma\pm\frac{\sqrt{4(\alpha^{2}-\Gamma^{2})D(\theta)+\Gamma^{2}(1+c_{1}+c_{2})^{2}}}{1+c_{1}+c_{2}}$
(8b)

with

$D(\theta)\equiv-1+c_{1}(\theta)^{2}+2c_{2}(\theta)-c_{2}(\theta)^{2}.$ (9)

While the real part of the first eigenvalue $\Lambda_{1}$ is negative for
$\alpha>\Gamma$, the second eigenvalue $\Lambda_{2}$ has a negative real part
only when

$D(\theta)<0.$ (10)

The inequality (10) thus provides the condition for the grid pattern to be
stable, regardless of the magnitude of the dissipation $\Gamma$. Note that the
stability analysis around the fixed points (Sk) and (Sp) leads to the same
inequality (10) as the condition for the stripe patterns to be _unstable_. As
seen in Fig. 2, the grid pattern is stable around an angle of $\theta=\pi/2$,
which is consistent with the experimental results [47].

Figure 2:  Stability criterion for the grid pattern: $D(\theta)$ in Eq. (9) as
a function of $\theta=\angle(\bm{k},\bm{p})\in[0,\pi/2]$ with fixed $A=0.6$
for different values of $\omega/\mu$. The grid pattern for angles
$\theta\approx\pi/2$ is stable where $D(\theta)<0$.

In the absence of dissipation, the eigenvalues always appear in positive and
negative pairs, such as Eq. (8) with $\Gamma=0$, because the amplitude
equation without dissipation enjoys the time-reversal symmetry inherited from
the GP equation. The behavior of the solution in the four-dimensional space
near the fixed points can be understood separately for each two-dimensional
subspace corresponding to each pair of positive and negative eigenvalues.
Without dissipation, the pair of eigenvalues with a real part makes the fixed
point in the corresponding two-dimensional subspace a saddle point, while the
pure imaginary pair makes the fixed point a center. A small dissipation
$\Gamma<\alpha$ keeps the saddle fixed point as a saddle while turning the
center into a stable focus (in-spiral).

Let us investigate the global behavior of the solutions of the amplitude
equation beyond the local behavior around the fixed points. When we introduce
the real and imaginary parts of the phase-rotated amplitudes as
$\rho_{k/p}=\mathop{\mathrm{Re}}[R_{k/p}e^{-i\bar{\eta}}]$ and
$\nu_{k/p}=\mathop{\mathrm{Im}}[R_{k/p}e^{-i\bar{\eta}}]$, all four fixed
points lie in the two-dimensional subspace spanned by $\rho_{k}$ and
$\rho_{p}$ with $\nu_{k}=\nu_{p}=0$. The four-dimensional flow trajectories
still depart from an unstable fixed point (saddle) and approach an attractive
one (in-spiral). This global behavior can be visualized by utilizing the
square eigenvalue $\Lambda^{2}$, which is positive (repulsive) for real
$\Lambda$ at the saddle, while it is negative (attractive) for imaginary
$\Lambda$ at the in-spiral. We can efficiently obtain the square eigenvalues
via the second-order differential equation derived from the amplitude equation
(6), which is given by

$\displaystyle\frac{\mathrm{d}^{2}}{\mathrm{d}t^{2}}\rho_{k}(t)\biggl{.}\biggr{|}_{\nu_{k}=\nu_{p}=0}=\lambda^{2}\biggl{[}\bar{R}^{2}+\rho_{k}(t)^{2}+(c_{1}-c_{2})\rho_{p}(t)^{2}\biggr{]}$
$\displaystyle\times\biggl{[}\bar{R}^{2}-\rho_{k}(t)^{2}-(c_{1}+c_{2})\rho_{p}(t)^{2}\biggr{]}\rho_{k}(t)$
(11)
$\displaystyle+2\lambda^{2}c_{2}\biggl{[}\bar{R}^{2}-\rho_{p}(t)^{2}-(c_{1}+c_{2})\rho_{k}(t)^{2}\biggr{]}\rho_{p}(t)^{2}\rho_{k}(t)$

and likewise for $\rho_{p}(t)$ after exchanging the $\rho_{k}(t)$ and
$\rho_{p}(t)$ variables. The force field described by the right-hand side of
Eq. (11) captures the global behavior of the solution of the original
amplitude equation, although it does not correspond to the trajectories of the
solutions themselves. This global behavior is shown in Fig. 1, and it changes
drastically depending on the angle between $\bm{k}$ and $\bm{p}$ 333Figure 1
shows the direction of the second derivative given by Eq. (11), which
represents the force. Also, since $\rho_{k/p}$ is equal to the absolute value
of the amplitude $R_{k/p}$ on the plane with fixed $\nu_{k}=\nu_{p}=0$, we
simply denote $\rho_{k/p}$ as $R_{k/p}$..

The amplitude equation with no divergence.—Finally, we discuss the divergence
of the coefficient $c_{1}(\theta)$ at the angles satisfying $2E=E_{+}$. At
this angle, two Bogoliubov modes with wavevectors $\bm{k}$ and $\bm{p}$ can
resonantly scatter into a Bogoliubov mode with wavevector $\bm{k}+\bm{p}$
without violating energy conservation, which enhances the contribution of this
collision process. Therefore, the amplitude of the wavevector $\bm{k}+\bm{p}$,
which grows proportionally to both the amplitudes $R_{k}$ and $R_{p}$, cannot
be neglected, and its omission in the two-mode ansatz (3) causes the
divergence of the coefficient $c_{1}(\theta)$. By additionally including the
$\bm{k}+\bm{p}$ mode described by the complex amplitude $R_{+}(t)$ in the
ansatz (3), we can derive the coupled amplitude equations for the three modes
$R_{k/p}(t)$ and $R_{+}(t)$ as (for details, see [51])

$\displaystyle i\frac{\mathrm{d}}{\mathrm{d}t}R_{k}(t)=-i\Gamma
R_{k}(t)-i\alpha R_{k}^{\ast}(t)+\beta(\theta)R_{+}(t)R_{p}^{\ast}(t)$
$\displaystyle+\lambda\Bigl{(}|R_{k}(t)|^{2}R_{k}(t)+\tilde{c}_{1}(\theta)|R_{p}(t)|^{2}R_{k}(t)$
$\displaystyle+c_{2}(\theta)R_{p}(t)^{2}R^{\ast}_{k}(t)\Bigr{)},$ (12a)
$\displaystyle
i\frac{\mathrm{d}}{\mathrm{d}t}R_{+}(t)=-i\Gamma_{+}R_{+}(t)+\Delta_{+}(\theta)R_{+}(t)$
$\displaystyle+\beta_{+}(\theta)R_{k}(t)R_{p}(t)+\lambda_{+}(\theta)|R_{+}(t)|^{2}R_{+}(t),$
(12b)

where $\Gamma$ and $\Gamma_{+}$ are dissipation coefficients. The parameters
are given by the detuning $\Delta_{+}(\theta)=2E-E_{+}$ and the nonlinearities
$\lambda_{+}(\theta)=\mu(5\epsilon_{+}+3\mu)/E_{+}$,
$\beta(\theta)=\mu[(\epsilon-\mu)/E+(\epsilon_{+}+2\mu)/E_{+}]$,
$\beta_{+}(\theta)=\mu[(\epsilon+2\mu)/E+\epsilon_{+}(\epsilon-\mu)/(\epsilon
E_{+})]$, and

$\displaystyle\tilde{c}_{1}(\theta)=c_{1}(\theta)$
$\displaystyle+\frac{\mu}{5\epsilon+3\mu}\frac{(2\epsilon-\mu)(\epsilon+2\mu)+(2\epsilon^{2}+\mu^{2})\epsilon_{+}/(2\epsilon)}{E^{2}-E^{2}_{+}/4}.$
(13)

We find that the coefficient $\tilde{c}_{1}(\theta)$ is now regular and the
divergence of the coefficient $c_{1}(\theta)$ at the singular angle satisfying
$2E-E_{+}$ is removed. We note that the $\bm{k}$ and $\bm{p}$ modes now
interact with the $\bm{k}+\bm{p}$ mode via quadratic terms in the amplitude
equation (three-mode scattering).

We first analyze the stability of the stripe-pattern fixed point, which can be
found analytically, and later confirm this result by numerical analysis of the
grid-pattern fixed point. Linearizing the amplitude equation (12) around the
stripe-pattern fixed point given by $R_{+}=0$ in addition to (Sk) in Eq. (7b),
we obtain six eigenvalues; two of them are given by $-\Gamma\pm
i\sqrt{4\alpha^{2}-5\Gamma^{2}}$, which must have negative real parts for
$\alpha>\Gamma$, and the real parts of the other four are plotted in Fig. 3.
In the case of Fig. 3, i.e., when the dissipation coefficients are
$\Gamma=0.5\alpha$ and $\Gamma_{+}=\alpha$, the stability (instability) of the
stripe pattern corresponds to the instability (stability) of the grid pattern,
as confirmed numerically 444The stability of the stripe pattern does not fully
correspond to the instability of the grid pattern unless the two dissipation
coefficients are large.. Thus, we are able to identify the angular region
where the largest eigenvalue in Fig. 3 has a positive real part (unstable
stripes) as the region where the grid pattern is stable. Also in this model,
the stable region in Fig. 3 is insensitive to changes in the magnitude of the
dissipation, and this insensitivity is consistent with the dissipation-
independent inequality (10) in the previous two-mode model. Figure 3 agrees
closely with Fig. 2 in that the grid pattern is stable around $\theta=\pi/2$;
moreover, it describes the stability for all angles reliably, including the
angle satisfying $2E=E_{+}$.

Figure 3:  Real parts of the four eigenvalues of the amplitude equation (12)
linearized around the stripe-pattern fixed point for $\omega/\mu=2$ and
$A=0.6$ (Each line is doubly degenerate for small $\theta$). The two
$\theta$-independent eigenvalues are excluded. The dissipation coefficients
are set as $\Gamma=0.5\alpha$ and $\Gamma_{+}=\alpha$, and the dashed and
dotted lines represent $-\Gamma$ and $-\Gamma_{+}$, respectively. The stripe
pattern at angle $\theta$ is unstable whenever $\mathop{\mathrm{Re}}\Lambda$
is positive. Notably, there is no singularity when $2E=E_{+}$, which occurs at
$\theta\approx 0.34\pi$ for our parameters in this figure.

Discussion and outlook.—In this Letter, we have derived the amplitude equation
(4) (and (12)) for pattern formation in two-dimensional BECs caused by the
Faraday instability. The amplitude equation can be considered as a complex
Ginzburg-Landau equation for pattern formation with the amplitudes in the two
directions as order parameters, so that it provides a simple description of
the system dynamics. Our method to derive the amplitude equation is equivalent
to the renormalization group theory for asymptotic analysis [55, 56, 57].
Accordingly, the amplitude equation describes the order parameter dynamics as
an effective model for the only two relevant modes that remain at long times,
while incorporating in a renormalization group sense the multitude of
irrelevant modes of the full GP solution. The dissipation that we introduced
in Eq. (6) appears effectively as a result of the renormalization onto the
model for only the two relevant modes. While the microscopic derivation of
dissipation from interactions and thermal fluctuations is still an open
question, one may also attempt to extract its value from experiments.

Using the obtained amplitude equation, we have analyzed the stability between
the uniform, stripe-pattern, and grid-pattern solutions. For $\alpha>\Gamma$,
where the drive amplitude is stronger than the dissipation, the uniform
solution becomes unstable, resulting in an inhomogeneous density pattern.
Figures 2 and 3 show that the grid pattern becomes stable around the angle
$\pi/2$ between the two excitation directions. The global stability of the
amplitude is shown in Fig. 1. Our results explain the experimental data
presented in our companion paper [47]. Furthermore, the amplitude equation has
been experimentally validated under various initial conditions and has been
confirmed to give a good description.

Finally, we briefly discuss the scattering between the three modes satisfying
$\bm{k}_{1}\pm\bm{k}_{2}\pm\bm{k}_{3}=\bm{0}$, which usually leads to a
triangular pattern [2]. In the Faraday pattern formation, the leading three-
mode scattering occurs at the frequency $\omega/2\pm\omega/2\pm\omega/2$
because the excited modes have the energy $\omega/2$ from the first resonance
condition $n=1$. Because of $\omega/2\pm\omega/2\pm\omega/2\neq 0$, the three-
mode scattering is a fast-rotating contribution in the rotation wave basis and
can be neglected. If the third mode is from the second resonance condition
$n=2$ with frequency $\omega$ instead of $\omega/2$, the three-mode scattering
becomes relevant because it has a slow-rotating contribution with
$\omega/2+\omega/2-\omega=0$. In fact, the amplitude equation (12) includes
three-mode scattering. However, we found that the stability of the patterns
does not change significantly between Eq. (4) and Eq. (12). Indeed, the second
resonant mode is less amplified by the Faraday instability than the first one
in the perturbative regime with a small drive amplitude $|A|<1$.

###### Acknowledgements.

The authors thank F. Ziebert for useful discussions. This work is supported by
the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under
project-ID 273811115 (SFB1225 ISOQUANT) and under Germany’s Excellence
Strategy EXC2181/1-390900948 (the Heidelberg STRUCTURES Excellence Cluster),
and QUANTERA DYNAMITE PCI2022-132919. This project was funded within the
QuantERA II Programme that has received funding from the European Union’s
Horizon 2020 research and innovation programme under Grant Agreement No
101017733. N.L. acknowledges support by the Studienstiftung des Deutschen
Volkes.

## References

  * Cross and Hohenberg [1993] M. C. Cross and P. C. Hohenberg, Pattern formation outside of equilibrium, Rev. Mod. Phys. 65, 851 (1993).
  * Cross and Greenside [2009] M. C. Cross and H. Greenside,  _Pattern Formation and Dynamics in Nonequilibrium Systems_ (Cambridge University Press, 2009).
  * Turing [1952] A. M. Turing, The chemical basis of morphogenesis, Phil. Trans. R. Soc. Lond. B 237, 37 (1952).
  * Koch and Meinhardt [1994] A. J. Koch and H. Meinhardt, Biological pattern formation: from basic mechanisms to complex structures, Rev. Mod. Phys. 66, 1481 (1994).
  * Faraday [1837] M. Faraday, On a peculiar class of acoustical figures; and on certain forms assumed by groups of particles upon vibrating elastic surfaces, Phil. Trans. R. Soc. Lond. 3, 49 (1837).
  * Wagner  _et al._ [2000] C. Wagner, H. W. Müller, and K. Knorr, Crossover from a square to a hexagonal pattern in Faraday surface waves, Phys. Rev. E 62, R33(R) (2000).
  * Staliunas  _et al._ [2002] K. Staliunas, S. Longhi, and G. J. de Valcárcel, Faraday Patterns in Bose-Einstein Condensates, Phys. Rev. Lett. 89, 210406 (2002).
  * Engels  _et al._ [2007] P. Engels, C. Atherton, and M. A. Hoefer, Observation of Faraday Waves in a Bose-Einstein Condensate, Phys. Rev. Lett. 98, 095301 (2007).
  * Nguyen  _et al._ [2019] J. H. V. Nguyen, M. C. Tsatsos, D. Luo, A. U. J. Lode, G. D. Telles, V. S. Bagnato, and R. G. Hulet, Parametric Excitation of a Bose-Einstein Condensate: From Faraday Waves to Granulation, Phys. Rev. X 9, 011052 (2019).
  * Kwon  _et al._ [2021] K. Kwon, K. Mukherjee, S. J. Huh, K. Kim, S. I. Mistakidis, D. K. Maity, P. G. Kevrekidis, S. Majumder, P. Schmelcher, and J.-y. Choi, Spontaneous Formation of Star-Shaped Surface Patterns in a Driven Bose-Einstein Condensate, Phys. Rev. Lett. 127, 113001 (2021).
  * Staliunas  _et al._ [2004] K. Staliunas, S. Longhi, and G. J. de Valcárcel, Faraday patterns in low-dimensional Bose-Einstein condensates, Phys. Rev. A 70, 011601(R) (2004).
  * Kevrekidis and Frantzeskakis [2004] P. G. Kevrekidis and D. J. Frantzeskakis, Pattern Forming Dynamical Instabilities of Bose-Einstein Condensates, Mod. Phys. Lett. B 18, 173 (2004).
  * Nicolin  _et al._ [2007] A. I. Nicolin, R. Carretero-González, and P. G. Kevrekidis, Faraday waves in Bose-Einstein condensates, Phys. Rev. A 76, 063609 (2007).
  * Bhattacherjee [2008] A. B. Bhattacherjee, Faraday instability in a two-component BoseâEinstein condensate, Phys. Scripta 78, 045009 (2008).
  * Capuzzi and Vignolo [2008] P. Capuzzi and P. Vignolo, Faraday waves in elongated superfluid fermionic clouds, Phys. Rev. A 78, 043613 (2008).
  * Nath and Santos [2010] R. Nath and L. Santos, Faraday patterns in two-dimensional dipolar Bose-Einstein condensates, Phys. Rev. A 81, 033626 (2010).
  * Nicolin and Raportaru [2010] A. I. Nicolin and M. C. Raportaru, Faraday waves in high-density cigar-shaped BoseâEinstein condensates, Physica A 389, 4663 (2010).
  * Capuzzi  _et al._ [2011] P. Capuzzi, M. Gattobigio, and P. Vignolo, Suppression of Faraday waves in a Bose-Einstein condensate in the presence of an optical lattice, Phys. Rev. A 83, 013603 (2011).
  * Tang  _et al._ [2011] R.-A. Tang, H.-C. Li, and J.-K. Xue, Faraday instability and Faraday patterns in a superfluid Fermi gas, J. Phys. B: At. Mol. Opt. Phys. 44, 115303 (2011).
  * Nicolin [2011] A. I. Nicolin, Resonant wave formation in Bose-Einstein condensates, Phys. Rev. E 84, 056202 (2011).
  * Sabbatini  _et al._ [2011] J. Sabbatini, W. H. Zurek, and M. J. Davis, Phase Separation and Pattern Formation in a Binary Bose-Einstein Condensate, Phys. Rev. Lett. 107, 230402 (2011).
  * Balaž and Nicolin [2012] A. Balaž and A. I. Nicolin, Faraday waves in binary nonmiscible Bose-Einstein condensates, Phys. Rev. A 85, 023613 (2012).
  * Nicolin [2012] A. I. Nicolin, Variational treatment of Faraday waves in inhomogeneous BoseâEinstein condensates, Physica A 391, 1062 (2012).
  * Łakomy  _et al._ [2012] K. Łakomy, R. Nath, and L. Santos, Faraday patterns in coupled one-dimensional dipolar condensates, Phys. Rev. A 86, 023620 (2012).
  * Abdullaev  _et al._ [2013] F. K. Abdullaev, M. Ögren, and M. P. Sørensen, Faraday waves in quasi-one-dimensional superfluid Fermi-Bose mixtures, Phys. Rev. A 87, 023616 (2013).
  * Balaž  _et al._ [2014] A. Balaž, R. Paun, A. I. Nicolin, S. Balasubramanian, and R. Ramaswamy, Faraday waves in collisionally inhomogeneous Bose-Einstein condensates, Phys. Rev. A 89, 023609 (2014).
  * Cairncross and Pelster [2014] W. Cairncross and A. Pelster, Parametric resonance in Bose-Einstein condensates with periodic modulation of attractive interaction, Eur. Phys. J. D 68, 106 (2014).
  * Abdullaev  _et al._ [2015] F. K. Abdullaev, A. Gammal, and L. Tomio, Faraday waves in BoseâEinstein condensates with engineering three-body interactions, J. Phys. B: At. Mol. Opt. Phys. 49, 025302 (2015).
  * Conforti  _et al._ [2016] M. Conforti, A. Mussot, A. Kudlinski, S. Rota Nodari, G. Dujardin, S. De Biévre, A. Armaroli, and S. Trillo, Heteroclinic Structure of Parametric Resonance in the Nonlinear Schrödinger Equation, Phys. Rev. Lett. 117, 013901 (2016).
  * Sudharsan  _et al._ [2016] J. B. Sudharsan, R. Radha, M. C. Raportaru, A. I. Nicolin, and A. BalaÅ¾, Faraday and resonant waves in binary collisionally-inhomogeneous BoseâEinstein condensates, J. Phys. B: At. Mol. Opt. Phys. 49, 165303 (2016).
  * Tomio  _et al._ [2017] L. Tomio, A. Gammal, and F. K. Abdullaev, Faraday Waves in Cold-Atom Systems with Two- and Three-Body Interactions, Few-Body Syst. 58 (2017).
  * Zhu  _et al._ [2019] C.-X. Zhu, W. Yi, G.-C. Guo, and Z.-W. Zhou, Parametric resonance of a Bose-Einstein condensate in a ring trap with periodically driven interactions, Phys. Rev. A 99, 023619 (2019).
  * VudragoviÄ and BalaÅ¾ [2019] D. VudragoviÄ and A. BalaÅ¾, Faraday and Resonant Waves in Dipolar Cigar-Shaped Bose-Einstein Condensates, Symmetry 11, 1090 (2019).
  * Abdullaev  _et al._ [2019] F. K. Abdullaev, A. Gammal, R. K. Kumar, and L. Tomio, Faraday waves and droplets in quasi-one-dimensional Bose gas mixtures, J. Phys. B: At. Mol. Opt. Phys. 52, 195301 (2019).
  * Chen  _et al._ [2019] T. Chen, K. Shibata, Y. Eto, T. Hirano, and H. Saito, Faraday patterns generated by Rabi oscillation in a binary Bose-Einstein condensate, Phys. Rev. A 100, 063610 (2019).
  * Turmanov  _et al._ [2020] B. K. Turmanov, B. B. Baizakov, and F. K. Abdullaev, Generation of density waves in dipolar quantum gases by time-periodic modulation of atomic interactions, Phys. Rev. A 101, 053616 (2020).
  * Okazaki  _et al._ [2021] K. Okazaki, J. Han, and M. Tsubota, Faraday waves in Bose–Einstein condensate: From instability to nonlinear dynamics (2021), arXiv:2012.02391 [cond-mat.quant-gas] .
  * Bougas  _et al._ [2021] G. Bougas, S. I. Mistakidis, and P. Schmelcher, Pattern formation of correlated impurities subjected to an impurity-medium interaction pulse, Phys. Rev. A 103, 023313 (2021).
  * Otlaadisa  _et al._ [2021] P. Otlaadisa, C. B. Tabi, and T. C. Kofané, Modulation instability in helicoidal spin-orbit coupled open Bose-Einstein condensates, Phys. Rev. E 103, 052206 (2021).
  * Cheng and Shi [2021] Y. Cheng and Z.-Y. Shi, Many-body dynamics with time-dependent interaction, Phys. Rev. A 104, 023307 (2021).
  * HernÃ¡ndez-Rajkov  _et al._ [2021] D. HernÃ¡ndez-Rajkov, J. E. Padilla-Castillo, A. del RÃo-Lima, A. GutiÃ©rrez-ValdÃ©s, F. J. Poveda-Cuevas, and J. A. Seman, Faraday waves in strongly interacting superfluids, New J. Phys. 23, 103038 (2021).
  * DÃaz  _et al._ [2021] P. DÃaz, L. PÃ©rez, L. Reyes, D. Laroze, and J. Bragard, Taming Faraday waves in binary fermionic clouds: The effect of Zeeman interaction, Chaos, Soliton. Fract. 153, 111416 (2021).
  * Zhang  _et al._ [2022] H. Zhang, S. Liu, and Y.-S. Zhang, Faraday patterns in spin-orbit-coupled Bose-Einstein condensates, Phys. Rev. A 105, 063319 (2022).
  * Jose  _et al._ [2023] S. M. Jose, K. Sah, and R. Nath, Patterns, spin-spin correlations and competing instabilities in driven quasi-two-dimensional spin-1 Bose-Einstein condensates, Phys. Rev. A 108, 023308 (2023).
  * Shukuno  _et al._ [2023] N. Shukuno, Y. Sano, and M. Tsubota, Faraday Waves in BoseâEinstein Condensates—The Excitation by the Modulation of the Interaction and the Potential, J. Phys. Soc. Jpn. 92, 064602 (2023).
  * Zhang  _et al._ [2020] Z. Zhang, K.-X. Yao, L. Feng, J. Hu, and C. Chin, Pattern formation in a driven Bose–Einstein condensate, Nat. Phys. 16, 652 (2020).
  * Liebster  _et al._ [2023] N. Liebster, M. Sparn, E. Kath, K. Fujii, S. L. Görlitz, T. Enss, H. Strobel, and M. K. Oberthaler, Spontaneous formation of persistent square pattern in a driven superfluid (2023), arXiv:2309.03792 [cond-mat.quant-gas] .
  * Pethick and Smith [2008] C. J. Pethick and H. Smith,  _BoseâEinstein Condensation in Dilute Gases_, 2nd ed. (Cambridge University Press, Cambridge, 2008).
  * Landau and Lifshitz [1976] L. D. Landau and E. M. Lifshitz,  _Mechanics_ , Vol. 1 (Butterworth-Heinemann, Oxford, 1976).
  * Note [1] As discussed in the last paragraph of the discussion section, there is no relevant three-mode scattering in our system due to the constraints imposed by the rotating-wave approximation, and the ansatz with two excited modes, as in Eq. (2), is sufficient.
  * [51] See Supplemental Material for the derivation of the amplitude equation from the GP equation .
  * Note [2] Since we have introduced the wavevectors $\bm{k}$ and $\bm{p}$ via the cosines in Eq. (2), the signs of them are irrelevant. Thus, we can take the angle $\theta$ being in $[0,\pi/2]$, and the case with $\theta\in[\pi/2,\pi]$ is given by replacing $\theta\in[0,\pi/2]$ with $\pi-\theta$.
  * Note [3] Figure 1 shows the direction of the second derivative given by Eq. (11), which represents the force. Also, since $\rho_{k/p}$ is equal to the absolute value of the amplitude $R_{k/p}$ on the plane with fixed $\nu_{k}=\nu_{p}=0$, we simply denote $\rho_{k/p}$ as $R_{k/p}$.
  * Note [4] The stability of the stripe pattern does not fully correspond to the instability of the grid pattern unless the two dissipation coefficients are large.
  * Chen  _et al._ [1994] L.-Y. Chen, N. Goldenfeld, and Y. Oono, Renormalization Group Theory for Global Asymptotic Analysis, Phys. Rev. Lett. 73, 1311 (1994).
  * Kunihiro [1995] T. Kunihiro, A Geometrical Formulation of the Renormalization Group Method for Global Analysis, Prog. Theo. Phys. 94, 503 (1995).
  * Chen  _et al._ [1996] L.-Y. Chen, N. Goldenfeld, and Y. Oono, Renormalization group and singular perturbations: Multiple scales, boundary layers, and reductive perturbation theory, Phys. Rev. E 54, 376 (1996).
  * de Valcárcel [2002] G. J. de Valcárcel, Faraday patterns in Bose-Einstein condensates. Amplitude equation for rolls in the parametrically driven, damped Gross-Pitaevskii equation (2002), arXiv:cond-mat/0204406 [cond-mat.soft] .

Supplemental Materials:
Square Pattern Formation as Stable Fixed Point in Driven Two-Dimensional Bose-
Einstein Condensates

##  Appendix A Derivation of the amplitude equation (4)

We perform a perturbative expansion of the wave function $\Psi(t,\bm{x})$
obeying the GP equation,

$\displaystyle i\frac{\partial}{\partial
t}\Psi(t,\bm{x})=\biggl{[}-\frac{\nabla^{2}}{2m}+g(t)|\Psi(t,\bm{x})|^{2}\biggr{]}\Psi(t,\bm{x}),\quad\textrm{with}\quad
g(t)=\bar{g}[1-A\sin(\omega t)],$ (S1)

using the multiple-scale method. This section is partly based on Ref. [58].

##  A.0.1 Order counting and perturbative equations

Introducing a small bookkeeping parameter $\varepsilon$, we denote quantities
as

$\displaystyle
A=\varepsilon^{2}\tilde{A},\qquad\frac{\omega}{2}-E_{\bm{k}}=\varepsilon^{2}\tilde{\omega}_{k},\qquad\frac{\omega}{2}-E_{\bm{p}}=\varepsilon^{2}\tilde{\omega}_{p}.$
(S2)

According to the multiple-scale method, we introduce a slow time
$\tau=\varepsilon^{2}t$ and expand the wave function as

$\displaystyle\Psi(t,\bm{x})=\Psi_{\textrm{uni}}(t)\biggl{[}1+\sum_{n=1}^{\infty}\varepsilon^{n}\phi_{n}(t,\tau,\bm{x})\biggr{]},$
(S3)

where we regarded $\phi_{n}$ as a function not only of $t$ and $\bm{x}$ but
also of $\tau$. Here, $\Psi_{\textrm{uni}}(t)=\Psi_{0}\exp[-i\mu
t-i(\mu/\omega)A\cos(\omega t)]$ is a uniform solution of Eq. (S1) with the
chemical potential $\mu=\bar{g}|\Psi_{0}|^{2}$. For later convenience, we
expand the wave function around $\Psi_{\textrm{uni}}(t)$ instead of the
zeroth-order solution, $\Psi_{0}e^{-i\mu t}$, of Eq. (S1). To focus on the two
excited modes spanning the two-dimensional pattern, we suppose
$\phi_{1}(t,\tau,\bm{x})$ to have the form of

$\displaystyle\phi_{1}(t,\tau,\bm{x})=\phi^{k}_{1}(t,\tau)\cos(\bm{k}\cdot\bm{x})+\phi^{p}_{1}(t,\tau)\cos(\bm{p}\cdot\bm{x}).$
(S4)

Substituting the above reparameterizations and the expansion, together with
$\partial_{t}\to\partial_{t}+\varepsilon^{2}\partial_{\tau}$, into Eq. (S1),
we find

$\displaystyle i\frac{\partial}{\partial
t}\phi_{1}(t,\tau,\bm{x})+\frac{\nabla^{2}}{2m}\phi_{1}(t,\tau,\bm{x})-2\mu\mathop{\mathrm{Re}}[\phi_{1}(t,\tau,\bm{x})]=0,$
(S5a) $\displaystyle i\frac{\partial}{\partial
t}\phi_{2}(t,\tau,\bm{x})+\frac{\nabla^{2}}{2m}\phi_{2}(t,\tau,\bm{x})-2\mu\mathop{\mathrm{Re}}[\phi_{2}(t,\tau,\bm{x})]=\mu\Bigl{(}|\phi_{1}(t,\tau,\bm{x})|^{2}+2\mathop{\mathrm{Re}}[\phi_{1}(t,\tau,\bm{x})]\phi_{1}(t,\tau,\bm{x})\Bigr{)},$
(S5b) $\displaystyle i\frac{\partial}{\partial
t}\phi_{3}(t,\tau,\bm{x})+\frac{\nabla^{2}}{2m}\phi_{3}(t,\tau,\bm{x})-2\mu\mathop{\mathrm{Re}}[\phi_{3}(t,\tau,\bm{x})]$
$\displaystyle=-i\frac{\partial}{\partial\tau}\phi_{1}(t,\tau,\bm{x})+\mu\Bigl{(}-2\tilde{A}\sin(\omega
t)\mathop{\mathrm{Re}}[\phi_{1}(t,\tau,\bm{x})]+|\phi_{1}(t,\tau,\bm{x})|^{2}\phi_{1}(t,\tau,\bm{x})+2\mathcal{K}(\phi_{1}(t,\tau,\bm{x}),\phi_{2}(t,\tau,\bm{x}))\Bigr{)},$
(S5c)

for increasing orders of $\varepsilon$, where we introduced

$\displaystyle\mathcal{K}(\phi_{1},\phi_{2})\equiv\mathop{\mathrm{Re}}[\phi_{1}]\phi_{2}+\mathop{\mathrm{Re}}[\phi_{2}]\phi_{1}+\mathop{\mathrm{Re}}[\phi_{1}\phi^{\ast}_{2}].$
(S6)

The solution of Eq. (S5a) describes the usual Bogoliubov quasiparticles and,
with the use of Eq. (S4), is found as

$\displaystyle\phi_{1}^{k}(t,\tau)$
$\displaystyle=\biggl{(}1-\frac{\epsilon_{\bm{k}}+2\mu}{E_{\bm{k}}}\biggr{)}r_{k}(\tau)e^{iE_{\bm{k}}t}+\biggl{(}1+\frac{\epsilon_{\bm{k}}+2\mu}{E_{\bm{k}}}\biggr{)}r^{\ast}_{k}(\tau)e^{-iE_{\bm{k}}t},$
(S7)

and the same for $\phi^{p}_{1}(t,\tau)$ with $\bm{k}$ replaced by $\bm{p}$.
Our remaining task in the multiple-scale method is to derive the equation
determining the complex amplitude $r_{k/p}(\tau)$ from the solvability
condition that Eqs. (S5b) and (S5c) have no secular terms.

##  A.0.2 Solution of the second-order equation

The general solution of Eq. (S5b) is irrelevant for our purposes; it is
sufficient to consider only modes resulting from scattering of the two modes
with wavevectors $\bm{k}$ and $\bm{p}$ in $\phi_{1}(t,\tau)$. This allow us to
take $\phi_{2}(t,\tau,\bm{x})$ in the form of

$\displaystyle\phi_{2}(t,\tau,\bm{x})=\phi_{20}(t,\tau)+\phi_{22}^{k}(t,\tau)\cos(2\bm{k}\cdot\bm{x})+\phi_{22}^{p}(t,\tau)\cos(2\bm{p}\cdot\bm{x})+\phi_{2}^{+}(t,\tau)\cos((\bm{k}+\bm{p})\cdot\bm{x})+\phi_{2}^{-}(t,\tau)\cos((\bm{k}-\bm{p})\cdot\bm{x}).$
(S8)

Plugging this form into Eq. (S5b) leads to

$\displaystyle i\frac{\partial}{\partial
t}\phi_{20}(t,\tau)-2\mu\mathop{\mathrm{Re}}[\phi_{20}(t,\tau)]=\frac{\mu}{2}\Bigl{(}|\phi^{k}_{1}(t,\tau)|^{2}+|\phi^{p}_{1}(t,\tau)|^{2}+2\mathop{\mathrm{Re}}[\phi^{k}_{1}(t,\tau)]\phi^{k}_{1}(t,\tau)+2\mathop{\mathrm{Re}}[\phi^{p}_{1}(t,\tau)]\phi^{p}_{1}(t,\tau)\Bigr{)}$
(S9)

for $\phi_{20}(t,\tau)$,

$\displaystyle i\frac{\partial}{\partial
t}\phi^{k}_{22}(t,\tau)-4\epsilon_{\bm{k}}\phi^{k}_{22}(t,\tau)-2\mu\mathop{\mathrm{Re}}[\phi^{k}_{22}(t,\tau)]=\frac{\mu}{2}\Bigl{(}|\phi^{k}_{1}(t,\tau)|^{2}+2\mathop{\mathrm{Re}}[\phi^{k}_{1}(t,\tau)]\phi^{k}_{1}(t,\tau)\Bigr{)}$
(S10)

for $\phi^{k}_{22}(t,\tau)$ and the same with $\bm{k}$ replaced by $\bm{p}$
for $\phi_{22}^{p}(t,\tau)$, and

$\displaystyle i\frac{\partial}{\partial
t}\phi^{+}_{2}(t,\tau)-\epsilon_{+}\phi^{+}_{2}(t,\tau)-2\mu\mathop{\mathrm{Re}}[\phi^{+}_{2}(t,\tau)]=\mu\mathcal{K}(\phi^{k}_{1}(t,\tau),\phi^{p}_{1}(t,\tau))$
(S11)

for $\phi^{+}_{2}(t,\tau)$ and the same with
$\epsilon_{+}=\epsilon_{\bm{k}+\bm{p}}$ replaced by
$\epsilon_{-}=\epsilon_{\bm{k}-\bm{p}}$ for $\phi^{-}_{2}(t,\tau)$.

By dividing Eqs. (S9), (S10), and (S11) into their real and imaginary parts,
and by using Eq. (S7), we can find the solutions to these equations. First,
the solution of Eq. (S9) is found as

$\displaystyle\phi_{20}(t,\tau)=\phi_{20}^{k}(t,\tau)+\phi^{p}_{20}(t,\tau)$
(S12)

with

$\displaystyle\phi_{20}^{k}(t,\tau)=-\frac{2\epsilon_{\bm{k}}+\mu}{\epsilon_{\bm{k}}}|r_{k}(\tau)|^{2}+\frac{\mu}{2\epsilon_{\bm{k}}}\biggl{(}1-\frac{\epsilon_{\bm{k}}}{E_{\bm{k}}}\biggr{)}r_{k}(\tau)^{2}e^{2iE_{\bm{k}}t}+\frac{\mu}{2\epsilon_{\bm{k}}}\biggl{(}1+\frac{\epsilon_{\bm{k}}}{E_{\bm{k}}}\biggr{)}r^{\ast}_{k}(\tau)^{2}e^{-2iE_{\bm{k}}t}+i\tilde{v}_{20}^{k}(\tau)$
(S13)

and the same with $\bm{k}$ replaced by $\bm{p}$ for $\phi^{p}_{20}(t,\tau)$.
Here, $\tilde{v}_{20}^{k}(\tau)$ is an arbitrary real function of $\tau$
independent of $t$, which does not contribute to the results.

Second, the solution of Eq. (S10) is found as

$\displaystyle\phi_{22}^{k}(t,\tau)=-\frac{\mu}{\epsilon_{\bm{k}}}|r_{k}(\tau)|^{2}-\frac{\mu}{2\epsilon_{\bm{k}}}\biggl{(}1-\frac{E_{\bm{k}}}{\epsilon_{\bm{k}}}\biggr{)}r_{k}(\tau)^{2}e^{2iE_{\bm{k}}t}-\frac{\mu}{2\epsilon_{\bm{k}}}\biggl{(}1+\frac{E_{\bm{k}}}{\epsilon_{\bm{k}}}\biggr{)}r^{\ast}_{k}(\tau)^{2}e^{-2iE_{\bm{k}}t},$
(S14)

and the same with $\bm{k}$ replaced by $\bm{p}$ for $\phi_{22}^{p}(t,\tau)$.
Although $\phi_{22}^{k}(t,\tau)$ has additional terms proportional to $e^{\pm
iE_{2\bm{k}}t}$ in addition to the particular solution (S14), we dropped them
because those additional terms do not contribute to the results.

Finally, the solution of Eq. (S11) is found in the form of

$\begin{split}\mathop{\mathrm{Re}}[\phi_{2}^{+}(t,\tau)]&=\varphi_{1}^{+}r_{k}(\tau)r_{p}(\tau)e^{i(E_{\bm{k}}+E_{\bm{p}})t}+\varphi_{2}^{+}r_{k}(\tau)r^{\ast}_{p}(\tau)e^{i(E_{\bm{k}}-E_{\bm{p}})t}+(\textrm{complex
conjugate}),\\\
\mathop{\mathrm{Im}}[\phi_{2}^{+}(t,\tau)]&=i\varphi_{3}^{+}r_{k}(\tau)r_{p}(\tau)e^{i(E_{\bm{k}}+E_{\bm{p}})t}+i\varphi_{4}^{+}r_{k}(\tau)r^{\ast}_{p}(\tau)e^{i(E_{\bm{k}}-E_{\bm{p}})t}+(\textrm{complex
conjugate}),\end{split}$ (S15)

where the coefficients $\varphi_{n}^{+}$ for $n=1,2,3,4$ satisfy the following
equations and are consequently real:

$\displaystyle\left\\{\begin{array}[]{l}\displaystyle(E_{\bm{k}}+E_{\bm{p}})\varphi^{+}_{1}-\epsilon_{+}\varphi^{+}_{3}=\mu\biggl{(}\frac{E_{\bm{k}}}{\epsilon_{\bm{k}}}+\frac{E_{\bm{p}}}{\epsilon_{\bm{p}}}\biggr{)},\\\\[8.61108pt]
\displaystyle(E_{\bm{k}}+E_{\bm{p}})\varphi^{+}_{3}-(\epsilon_{+}+2\mu)\varphi^{+}_{1}=\mu\biggl{(}3-\frac{E_{\bm{k}}}{\epsilon_{\bm{k}}}\frac{E_{\bm{p}}}{\epsilon_{\bm{p}}}\biggr{)},\end{array}\right.\quad\left\\{\begin{array}[]{l}\displaystyle(E_{\bm{k}}-E_{\bm{p}})\varphi^{+}_{2}-\epsilon_{+}\varphi^{+}_{4}=\mu\biggl{(}\frac{E_{\bm{k}}}{\epsilon_{\bm{k}}}-\frac{E_{\bm{p}}}{\epsilon_{\bm{p}}}\biggr{)},\\\\[8.61108pt]
\displaystyle(E_{\bm{k}}-E_{\bm{p}})\varphi^{+}_{4}-(\epsilon_{+}+2\mu)\varphi^{+}_{2}=\mu\biggl{(}3+\frac{E_{\bm{k}}}{\epsilon_{\bm{k}}}\frac{E_{\bm{p}}}{\epsilon_{\bm{p}}}\biggr{)}.\end{array}\right.$
(S20)

We again dropped the additional terms in the solution for
$\phi^{+}_{2}(t,\tau)$ because they do not contribute to the results. The
solution of $\phi^{-}_{2}(t,\tau)$ is given in the same form with
$\varphi^{+}_{n}$ replaced by $\varphi^{-}_{n}$ in Eq. (S15), and its
coefficients $\varphi^{-}_{n}$ satisfy the same equations (S20) with
$\varphi^{+}_{n}$ and $\epsilon_{+}$ replaced by $\varphi^{-}_{n}$ and
$\epsilon_{-}$, respectively.

##  A.0.3 The solvability condition

Similarly to $\phi_{2}(t,\tau,\bm{x})$, we expand $\phi_{3}(t,\tau,\bm{x})$ in
the Fourier cosine basis as

$\displaystyle\phi_{3}(t,\tau,\bm{x})=\phi_{31}^{k}(t,\tau)\cos(\bm{k}\cdot\bm{x})+\phi_{31}^{p}(t,\tau)\cos(\bm{p}\cdot\bm{x})+\cdots,$
(S21)

where the ellipsis represents terms with other bases. We need only
$\phi_{31}^{k/p}(t,\tau)$ for our purpose. From Eq. (S5c), the differential
equation for $\phi^{k}_{31}(t,\tau)$ is obtained as

$\displaystyle i\frac{\partial}{\partial
t}\phi_{31}^{k}(t,\tau)-\epsilon_{\bm{k}}\phi^{k}_{31}(t,\tau)-2\mu\mathop{\mathrm{Re}}[\phi_{31}^{k}(t,\tau)]=f(t,\tau)$
(S22)

with

$\displaystyle f(t,\tau)$
$\displaystyle=-i\frac{\partial}{\partial\tau}\phi_{1}^{k}(t,\tau)+\mu\biggl{\\{}-2\tilde{A}\sin(\omega
t)\mathop{\mathrm{Re}}[\phi^{k}_{1}(t,\tau)]+\frac{3}{4}|\phi^{k}_{1}(t,\tau)|^{2}\phi^{k}_{1}(t,\tau)+|\phi^{p}_{1}(t,\tau)|^{2}\phi^{k}_{1}(t,\tau)+\frac{1}{2}\phi^{k\ast}_{1}(t,\tau)\phi^{p}_{1}(t,\tau)^{2}$
$\displaystyle\quad+\mathcal{K}(\phi^{k}_{1}(t,\tau),2\phi_{20}(t,\tau)+\phi^{k}_{22}(t,\tau))+\mathcal{K}(\phi^{p}_{1}(t,\tau),\phi_{2}^{+}(t,\tau)+\phi^{-}_{2}(t,\tau))\biggr{\\}}.$
(S23)

Here, we write Eq. (S22) in the matrix form as

$\displaystyle\frac{\partial}{\partial
t}\left(\begin{array}[]{c}\mathop{\mathrm{Re}}[\phi_{31}^{k}(t,\tau)]\\\\[5.16663pt]
\mathop{\mathrm{Im}}[\phi_{31}^{k}(t,\tau)]\end{array}\right)+\left(\begin{array}[]{cc}0&-\epsilon_{\bm{k}}\\\\[5.16663pt]
\epsilon_{\bm{k}}+2\mu&0\end{array}\right)\left(\begin{array}[]{c}\mathop{\mathrm{Re}}[\phi_{31}^{k}(t,\tau)]\\\\[5.16663pt]
\mathop{\mathrm{Im}}[\phi_{31}^{k}(t,\tau)]\end{array}\right)=\left(\begin{array}[]{c}\mathop{\mathrm{Im}}[f(t,\tau)]\\\\[5.16663pt]
-\mathop{\mathrm{Re}}[f(t,\tau)]\end{array}\right).$ (S32)

Diagonalizing this equation, we obtain

$\displaystyle\frac{\partial}{\partial
t}\chi(t,\tau)-iE_{\bm{k}}\chi(t,\tau)=g(t,\tau)$ (S33)

and its complex conjugate with

$\displaystyle\chi(t,\tau)$
$\displaystyle=\frac{1}{\sqrt{2}}\biggl{(}\frac{\epsilon_{\bm{k}}+2\mu}{E_{\bm{k}}}\mathop{\mathrm{Re}}[\phi^{k}_{31}(t,\tau)]-i\mathop{\mathrm{Im}}[\phi^{k}_{31}(t,\tau)]\biggr{)},$
(S34) $\displaystyle g(t,\tau)$
$\displaystyle=\frac{1}{\sqrt{2}}\biggl{(}\frac{\epsilon_{\bm{k}}+2\mu}{E_{\bm{k}}}\mathop{\mathrm{Im}}[f(t,\tau)]+i\mathop{\mathrm{Re}}[f(t,\tau)]\biggr{)}.$
(S35)

Therefore, since the homogeneous version of Eq. (S33), i.e., Eq. (S33) with
$g(t,\tau)=0$, has an oscillating solution $e^{iE_{\bm{k}}t}$, the condition
for Eq. (S33) to have no secular term is equivalent to that $g(t,\tau)$ has no
term proportional to $e^{iE_{\bm{k}}t}$.

Substituting the above results into $g(t,\tau)$ straightforwardly yields

$\displaystyle\frac{i}{\sqrt{2}}\frac{E_{\bm{k}}}{\epsilon_{\bm{k}}+2\mu}g(t,\tau)=-i\frac{\mathrm{d}r_{k}(\tau)}{\mathrm{d}\tau}e^{iE_{\bm{k}}t}-i\mu\frac{\tilde{A}}{2}\frac{\epsilon_{\bm{k}}}{E_{\bm{k}}}r^{\ast}_{k}(\tau)e^{i\tilde{\omega}_{k}\tau}e^{iE_{\bm{k}}t}$
$\displaystyle\quad+\mu\frac{5\epsilon_{\bm{k}}+3\mu}{E_{\bm{k}}}\Bigl{(}|r_{k}(\tau)|^{2}r_{k}(\tau)+c_{1}(\bm{k},\bm{p})|r_{p}(\tau)|^{2}r_{k}(\tau)+c_{2}(\bm{k},\bm{p})r_{p}(\tau)^{2}r_{k}^{\ast}(\tau)e^{2i(\tilde{\omega}_{k}-\tilde{\omega}_{p})\tau}\Bigr{)}e^{iE_{\bm{k}}t}+\cdots$
(S36)

with

$\displaystyle c_{1}(\bm{k},\bm{p})$
$\displaystyle=\frac{E_{\bm{k}}}{5\epsilon_{\bm{k}}+3\mu}\biggl{[}4\frac{\epsilon_{\bm{k}}\epsilon_{\bm{p}}-\mu^{2}}{\epsilon_{\bm{p}}E_{\bm{k}}}-\frac{1}{2}(\varphi^{+}_{1}+\varphi^{-}_{1})\biggl{(}-\frac{\epsilon_{\bm{p}}+2\mu}{E_{\bm{p}}}+3\frac{\epsilon_{\bm{k}}}{E_{\bm{k}}}\biggr{)}-\frac{1}{2}(\varphi^{+}_{2}+\varphi^{-}_{2})\biggl{(}\frac{\epsilon_{\bm{p}}+2\mu}{E_{\bm{p}}}+3\frac{\epsilon_{\bm{k}}}{E_{\bm{k}}}\biggr{)}$
$\displaystyle\quad-\frac{1}{2}(\varphi^{+}_{3}+\varphi^{-}_{3})\biggl{(}1+\frac{\epsilon_{\bm{p}}+2\mu}{E_{\bm{p}}}\frac{\epsilon_{\bm{k}}}{E_{\bm{k}}}\biggr{)}-\frac{1}{2}(\varphi^{+}_{4}+\varphi^{-}_{4})\biggl{(}1-\frac{\epsilon_{\bm{p}}+2\mu}{E_{\bm{p}}}\frac{\epsilon_{\bm{k}}}{E_{\bm{k}}}\biggr{)}\biggr{]},$
(S37a) $\displaystyle c_{2}(\bm{k},\bm{p})$
$\displaystyle=-\frac{E_{\bm{k}}}{5\epsilon_{\bm{k}}+3\mu}\biggl{[}\frac{\epsilon_{\bm{k}}+\mu}{E_{\bm{k}}}+\frac{\mu}{\epsilon_{\bm{p}}}\frac{2\epsilon_{\bm{k}}+2\mu}{E_{\bm{k}}}+\frac{\epsilon_{\bm{p}}+3\mu}{E_{\bm{p}}}$
$\displaystyle\quad+\frac{1}{2}(\varphi^{+}_{2}+\varphi^{-}_{2})\biggl{(}\frac{\epsilon_{\bm{p}}+2\mu}{E_{\bm{p}}}+3\frac{\epsilon_{\bm{k}}}{E_{\bm{k}}}\biggr{)}-\frac{1}{2}(\varphi^{+}_{4}+\varphi^{-}_{4})\biggl{(}1-\frac{\epsilon_{\bm{p}}+2\mu}{E_{\bm{p}}}\frac{\epsilon_{\bm{k}}}{E_{\bm{k}}}\biggr{)}\biggr{]},$
(S37b)

where the ellipsis represents terms not proportional to $e^{iE_{\bm{k}}t}$.
Note that
$e^{iE_{\bm{p}}t}=e^{iE_{\bm{k}}t}e^{i(\tilde{\omega}_{k}-\tilde{\omega}_{p})\tau}$
is considered proportional to $e^{iE_{\bm{k}}t}$ as a function of $t$ because
of
$E_{\bm{p}}=\frac{\omega}{2}-\varepsilon^{2}\tilde{\omega}_{p}=E_{\bm{k}}+\varepsilon^{2}(\tilde{\omega}_{k}-\tilde{\omega}_{p})$.
Thus, the solvability condition that the coefficient of $e^{iE_{\bm{k}}t}$ in
$g(t,\tau)$ is equal to zero is found to be

$\displaystyle i\frac{\mathrm{d}}{\mathrm{d}\tau}r_{k}(\tau)$
$\displaystyle=-i\mu\frac{\tilde{A}}{2}\frac{\epsilon_{\bm{k}}}{E_{\bm{k}}}r^{\ast}_{k}(\tau)e^{i\tilde{\omega}_{k}\tau}+\mu\frac{5\epsilon_{\bm{k}}+3\mu}{E_{\bm{k}}}\Bigl{(}|r_{k}(\tau)|^{2}r_{k}(\tau)+c_{1}(\bm{k},\bm{p})|r_{p}(\tau)|^{2}r_{k}(\tau)+c_{2}(\bm{k},\bm{p})r_{p}(\tau)^{2}r_{k}^{\ast}(\tau)e^{2i(\tilde{\omega}_{k}-\tilde{\omega}_{p})\tau}\Bigr{)}.$
(S38)

##  A.0.4 Resulting amplitude equation

Assuming the absolute values of $\bm{k}$ and $\bm{p}$ to be equal, we
introduce

$\displaystyle\epsilon=\epsilon_{\bm{k}}=\epsilon_{\bm{p}},\qquad
E=E_{\bm{k}}=E_{\bm{p}}=\sqrt{\epsilon(\epsilon+2\mu)}.$ (S39)

Under this assumption, we find the coefficients $\varphi^{+}_{n}$ from Eq.
(S20) as

$\displaystyle\varphi^{+}_{1}=\frac{\mu\left[\epsilon+2\mu+(\epsilon-\mu)\epsilon_{+}/(2\epsilon)\right]}{E^{2}-E^{2}_{+}/4},\quad\varphi^{+}_{2}=-\frac{\mu}{\epsilon}\frac{2\epsilon+\mu}{\epsilon_{+}/2+\mu},\quad\varphi^{+}_{3}=\frac{\mu
E\left[1+\epsilon_{+}/(2\epsilon)\right]}{E^{2}-E^{2}_{+}/4},\quad\varphi^{+}_{4}=0,$
(S40)

where we introduced $E_{+}=\sqrt{\epsilon_{+}(\epsilon_{+}+2\mu)}$. The
coefficients $\varphi_{n}^{-}$ are obtained by replacing $\epsilon_{+}$ with
$\epsilon_{-}$ in Eq. (S40). Under the assumption (S39),
$\epsilon_{+}=4\epsilon\cos^{2}(\theta/2)$ and
$\epsilon_{-}=4\epsilon\sin^{2}(\theta/2)$ turn into functions of the angle
between $\bm{k}$ and $\bm{p}$, and thus the coefficients
$c_{1}(\bm{k},\bm{p})$ and $c_{2}(\bm{k},\bm{p})$ are also simplified as
functions of the angle,

$\displaystyle c_{1}(\theta)$
$\displaystyle=\frac{\mu}{5\epsilon+3\mu}\biggl{[}4\frac{\epsilon^{2}-\mu^{2}}{\mu\epsilon}+\biggl{(}\frac{2\epsilon+\mu}{\epsilon}\frac{2\epsilon+\mu}{\epsilon_{+}/2+\mu}-\frac{(2\epsilon-\mu)(\epsilon+2\mu)+(2\epsilon^{2}+\mu^{2})\epsilon_{+}/(2\epsilon)}{E^{2}-E^{2}_{+}/4}+(\epsilon_{+}\rightarrow\epsilon_{-})\biggr{)}\biggr{]},$
(S41a) $\displaystyle c_{2}(\theta)$
$\displaystyle=\frac{\mu}{5\epsilon+3\mu}\biggl{[}-2\frac{\epsilon^{2}+3\mu\epsilon+\mu^{2}}{\mu\epsilon}+\frac{2\epsilon+\mu}{\epsilon}\biggl{(}\frac{2\epsilon+\mu}{\epsilon_{+}/2+\mu}+(\epsilon_{+}\to\epsilon_{-})\biggr{)}\biggr{]},$
(S41b)

which are Eqs. (5a) and (5b) in the main text.

To get back the original parameters, let us introduce

$\displaystyle R_{k}(t)=\varepsilon
r_{k}(\varepsilon^{2}t)e^{-i\varepsilon^{2}\tilde{\omega}_{k}t}.$ (S42)

Then, the solvability condition (S38) turns into

$\displaystyle i\frac{\mathrm{d}}{\mathrm{d}t}R_{k}(t)=\Delta R_{k}(t)-i\alpha
R^{\ast}_{k}(t)+\lambda\Bigl{(}|R_{k}(t)|^{2}R_{k}(t)+c_{1}(\theta)|R_{p}(t)|^{2}R_{k}(t)+c_{2}(\theta)R_{p}(t)^{2}R^{\ast}_{k}(t)\Bigr{)},$
(S43)

with

$\displaystyle\Delta=\biggl{(}\frac{\omega}{2}-E\biggr{)},\qquad\alpha=\mu\frac{A}{2}\frac{\epsilon}{E},\qquad\lambda=\mu\frac{5\epsilon+3\mu}{E},$
(S44)

which is the amplitude equation (4) in the main text. Also,
$\varepsilon\phi^{k/p}_{1}(t,\tau)$ with Eq. (S7) expressed in terms of
$R_{k}(t)$ using Eq. (S42) corresponds to Eq. (3) in the main text.

##  Appendix B Derivation of the amplitude equation (12)

We again start with the GP equation (S1). Considering the growth of the mode
with the wavevector $\bm{k}+\bm{p}$ associated with the growth of the
$\bm{k}$\- and $\bm{p}$-modes, we take the following order counting with
respect to the small bookkeeping parameter $\varepsilon$:

$\displaystyle
A=\varepsilon^{2}\tilde{A},\qquad\frac{\omega}{2}-E_{\bm{k}}=\varepsilon\tilde{\omega}_{k},\qquad\frac{\omega}{2}-E_{\bm{p}}=\varepsilon\tilde{\omega}_{p},\qquad\omega-
E_{+}=\varepsilon\tilde{\omega}_{+}.$ (S45)

Introducing two slow timescales $\tau_{1}=\varepsilon t$ and
$\tau_{2}=\varepsilon^{2}t$, we expand the wave function as

$\displaystyle\Psi(t,\bm{x})=\Psi_{\textrm{uni}}(t)\Biggl{[}1+\sum_{n=1}^{\infty}\varepsilon^{n}\phi_{n}(t,\tau,\bm{x})\Biggr{]},$
(S46)

where we regard $\phi_{n}$ as a function not only of $t$ and $\bm{x}$ but also
of $\tau_{1}$ and $\tau_{2}$, and labeled the two slow times collectively as
$\tau=(\tau_{1},\tau_{2})$ for brevity. In this section, we suppose
$\phi_{1}(t,\tau,\bm{x})$ to have the form of

$\displaystyle\phi_{1}(t,\tau,\bm{x})=\phi^{k}_{1}(t,\tau)\cos(\bm{k}\cdot\bm{x})+\phi^{p}_{1}(t,\tau)\cos(\bm{p}\cdot\bm{x})+\phi^{+}_{1}(t,\tau)\cos((\bm{k}+\bm{p})\cdot\bm{x}).$
(S47)

Substituting the above reparameterizations and the expansion, together with
$\partial_{t}\to\partial_{t}+\varepsilon\partial_{\tau_{1}}+\varepsilon^{2}\partial_{\tau_{2}}$,
into Eq. (S1), we find

$\displaystyle i\frac{\partial}{\partial
t}\phi_{1}(t,\tau,\bm{x})+\frac{\nabla^{2}}{2m}\phi_{1}(t,\tau,\bm{x})-2\mu\mathop{\mathrm{Re}}[\phi_{1}(t,\tau,\bm{x})]=0,$
(S48a) $\displaystyle i\frac{\partial}{\partial
t}\phi_{2}(t,\tau,\bm{x})+\frac{\nabla^{2}}{2m}\phi_{2}(t,\tau,\bm{x})-2\mu\mathop{\mathrm{Re}}[\phi_{2}(t,\tau,\bm{x})]=-i\frac{\partial}{\partial\tau_{1}}\phi_{1}(t,\tau,\bm{x})+\mu\Bigl{(}|\phi_{1}(t,\tau,\bm{x})|^{2}+2\mathop{\mathrm{Re}}[\phi_{1}(t,\tau,\bm{x})]\phi_{1}(t,\tau,\bm{x})\Bigr{)},$
(S48b) $\displaystyle i\frac{\partial}{\partial
t}\phi_{3}(t,\tau,\bm{x})+\frac{\nabla^{2}}{2m}\phi_{3}(t,\tau,\bm{x})-2\mu\mathop{\mathrm{Re}}[\phi_{3}(t,\tau,\bm{x})]=-i\frac{\partial}{\partial\tau_{1}}\phi_{2}(t,\tau,\bm{x})-i\frac{\partial}{\partial\tau_{2}}\phi_{1}(t,\tau,\bm{x})$
$\displaystyle+\mu\Bigl{(}-2\tilde{A}\sin(\omega
t)\mathop{\mathrm{Re}}[\phi_{1}(t,\tau,\bm{x})]+|\phi_{1}(t,\tau,\bm{x})|^{2}\phi_{1}(t,\tau,\bm{x})+2\mathcal{K}(\phi_{1}(t,\tau,\bm{x}),\phi_{2}(t,\tau,\bm{x}))\Bigr{)},$
(S48c)

for increasing orders of $\varepsilon$. The solution of Eq. (S48a) again
describes the usual Bogoliubov quasiparticles and, with the use of Eq. (S47),
is found as

$\displaystyle\phi_{1}^{+}(t,\tau)$
$\displaystyle=\biggl{(}1-\frac{\epsilon_{+}+2\mu}{E_{+}}\biggr{)}r_{+}(\tau)e^{iE_{+}t}+\biggl{(}1+\frac{\epsilon_{+}+2\mu}{E_{+}}\biggr{)}r^{\ast}_{+}(\tau)e^{-iE_{+}t},$
(S49)

and the same as before for $\phi_{1}^{k/p}(t,\tau)$, except that $\tau$
represents the two variables $\tau_{1}$ and $\tau_{2}$.

Our task, as in the previous section, is to find the solvability conditions
that Eqs. (S48b) and (S48c) have no secular terms, although it is complicated
by the presence of the new timescale $\tau_{1}$ and the additional mode with
the wavevector $\bm{k}+\bm{p}$. In the following, we consider only the minimal
contributions necessary to remove the divergence from the previous amplitude
equation as interactions between the $\bm{k}+\bm{p}$ mode and the $\bm{k}$\-
or $\bm{p}$-mode.

##  B.0.1 Solvability conditions from the second-order equation

For $\phi_{2}(t,\tau,\bm{x})$, it is sufficient to consider only modes
resulting from scattering of the two modes in $\phi_{1}(t,\tau)$, as in the
previous model. This allows us to take $\phi_{2}(t,\tau,\bm{x})$ in the form
of

$\displaystyle\phi_{2}(t,\tau,\bm{x})$
$\displaystyle=\phi_{20}(t,\tau)+\phi^{k}_{22}(t,\tau)\cos(2\bm{k}\cdot\bm{x})+\phi^{p}_{22}(t,\tau)\cos(2\bm{p}\cdot\bm{x})+\phi^{+}_{22}(t,\tau)\cos(2(\bm{k}+\bm{p})\cdot\bm{x})$
$\displaystyle\quad+\phi^{+}_{2}(t,\tau)\cos((\bm{k}+\bm{p})\cdot\bm{x})+\phi^{-}_{2}(t,\tau)\cos((\bm{k}-\bm{p})\cdot\bm{x})+\phi^{k}_{21}(t,\tau)\cos(\bm{k}\cdot\bm{x})+\phi^{p}_{21}(t,\tau)\cos(\bm{p}\cdot\bm{x}).$
(S50)

Here, we dropped the contributions of the modes with the wavevectors
$2\bm{k}+\bm{p}$ and $\bm{k}+2\bm{p}$ because they are irrelevant to the
divergence in the previous model.

Plugging the expansion (S50) into Eq. (S48b), one obtains the differential
equation for each mode. We first consider the differential equation for
$\phi^{k}_{21}(t,\tau)$:

$\displaystyle i\frac{\partial}{\partial
t}\phi^{k}_{21}(t,\tau)-\epsilon_{\bm{k}}\phi^{k}_{21}(t,\tau)-2\mu\mathop{\mathrm{Re}}[\phi^{k}_{21}(t,\tau)]=f_{21}^{k}(t,\tau)$
(S51)

with

$\displaystyle f^{k}_{21}(t,\tau)$
$\displaystyle=-i\frac{\partial}{\partial\tau_{1}}\phi_{1}^{k}(t,\tau)+\mu\mathcal{K}(\phi^{p}_{1}(t,\tau),\phi^{+}_{1}(t,\tau)).$
(S52)

Taking the real and imaginary parts of the differential equation and
diagonalizing them, we get

$\displaystyle\frac{\partial}{\partial
t}\chi^{k}_{21}(t,\tau)-iE_{\bm{k}}\chi^{k}_{21}(t,\tau)=g^{k}_{21}(t,\tau)$
(S53)

and its complex conjugate with

$\displaystyle\chi^{k}_{21}(t,\tau)=\frac{1}{\sqrt{2}}\biggl{(}\frac{\epsilon_{\bm{k}}+2\mu}{E_{\bm{k}}}\mathop{\mathrm{Re}}[\phi^{k}_{21}(t,\tau)]-i\mathop{\mathrm{Im}}[\phi^{k}_{21}(t,\tau)]\biggr{)},$
(S54) $\displaystyle
g^{k}_{21}(t,\tau)=\frac{1}{\sqrt{2}}\biggl{(}\frac{\epsilon_{\bm{k}}+2\mu}{E_{\bm{k}}}\mathop{\mathrm{Im}}[f^{k}_{21}(t,\tau)]+i\mathop{\mathrm{Re}}[f^{k}_{21}(t,\tau)]\biggr{)}.$
(S55)

Since Eq. (S53) with $g^{k}_{21}(t,\tau)=0$ has an oscillating solution
$e^{iE_{\bm{k}}t}$, the condition for Eq. (S53) to have no secular term is
equivalent to that $g^{k}_{21}(t,\tau)$ has no term proportional to
$e^{iE_{\bm{k}}t}$. The solvability condition is obtained as

$\displaystyle
i\frac{\partial}{\partial\tau_{1}}r_{k}(\tau)=-\beta(\bm{k},\bm{p})r_{+}(\tau)r^{\ast}_{p}(\tau)e^{-i(\tilde{\omega}_{+}-\tilde{\omega}_{k}-\tilde{\omega}_{p})\tau_{1}},$
(S56)

with

$\displaystyle\beta(\bm{k},\bm{p})=\frac{\mu}{2}\biggl{[}\frac{\epsilon_{+}+2\mu}{E_{+}}-\frac{\epsilon_{\bm{p}}+2\mu}{E_{\bm{p}}}+\frac{\epsilon_{\bm{k}}}{E_{\bm{k}}}\biggl{(}3+\frac{\epsilon_{\bm{p}}+2\mu}{E_{\bm{p}}}\frac{\epsilon_{+}+2\mu}{E_{+}}\biggr{)}\biggr{]}$
(S57)

Under this solvability condition, $\phi^{k}_{21}(t,\tau)$ is found as

$\displaystyle\phi^{k}_{21}(t,\tau)$
$\displaystyle=\biggl{(}1-\frac{\epsilon_{\bm{k}}+2\mu}{E_{\bm{k}}}\biggr{)}\tilde{r}_{k}(\tau)e^{iE_{\bm{k}}t}+\biggl{(}1+\frac{\epsilon_{\bm{k}}+2\mu}{E_{\bm{k}}}\biggr{)}\tilde{r}^{\ast}_{k}(\tau)e^{-iE_{\bm{k}}t}$
$\displaystyle\quad+(\varphi^{k}_{1}-\varphi^{k}_{2})r_{p}(\tau)r_{+}(\tau)e^{i(E_{\bm{p}}+E_{+})t}+(\varphi^{k}_{1}+\varphi^{k}_{2})r^{\ast}_{p}(\tau)r^{\ast}_{+}(\tau)e^{-i(E_{\bm{p}}+E_{+})t},$
(S58)

where $\tilde{r}_{k}(\tau)$ is an arbitrary complex amplitude. Although the
coefficient $\varphi^{k}_{1/2}$ can be found by substituting this solution
into the differential equation, it is irrelevant to our purpose. For
$\phi^{p}_{21}(t,\tau)$, we can find the same solvability condition (S56) and
the solution (S58) with $\bm{k}$ and $\bm{p}$ replaced with each other.

We next consider the differential equation for $\phi^{+}_{2}(t,\tau)$:

$\displaystyle i\frac{\partial}{\partial
t}\phi^{+}_{2}(t,\tau)-\epsilon_{+}\phi^{+}_{2}(t,\tau)-2\mu\mathop{\mathrm{Re}}[\phi^{+}_{2}(t,\tau)]=f_{2}^{+}(t,\tau)$
(S59)

with

$\displaystyle f^{+}_{2}(t,\tau)$
$\displaystyle=-i\frac{\partial}{\partial\tau_{1}}\phi_{1}^{+}(t,\tau)+\mu\mathcal{K}(\phi^{k}_{1}(t,\tau),\phi^{p}_{1}(t,\tau)).$
(S60)

Taking the real and imaginary parts of the differential equation and
diagonalizing them, we obtain

$\displaystyle\frac{\partial}{\partial
t}\chi^{+}_{2}(t,\tau)-iE_{+}\chi^{+}_{2}(t,\tau)=g^{+}_{2}(t,\tau)$ (S61)

and its complex conjugate with

$\displaystyle\chi^{+}_{2}(t,\tau)=\frac{1}{\sqrt{2}}\biggl{(}\frac{\epsilon_{+}+2\mu}{E_{+}}\mathop{\mathrm{Re}}[\phi^{+}_{2}(t,\tau)]-i\mathop{\mathrm{Im}}[\phi^{+}_{2}(t,\tau)]\biggr{)},$
(S62) $\displaystyle
g^{+}_{2}(t,\tau)=\frac{1}{\sqrt{2}}\biggl{(}\frac{\epsilon_{+}+2\mu}{E_{+}}\mathop{\mathrm{Im}}[f^{+}_{2}(t,\tau)]+i\mathop{\mathrm{Re}}[f^{+}_{2}(t,\tau)]\biggr{)}.$
(S63)

Since Eq. (S61) with $g^{+}_{2}(t,\tau)=0$ has an oscillating solution
$e^{iE_{+}t}$, the condition for Eq. (S61) to have no secular term is
equivalent to that $g^{+}_{2}(t,\tau)$ has no term proportional to
$e^{iE_{+}t}$. Using
$e^{i(E_{\bm{k}}+E_{\bm{p}})t}=e^{iE_{+}t}e^{i(\tilde{\omega}_{+}-\tilde{\omega}_{k}-\tilde{\omega}_{p})\tau_{1}}$,
we get the solvability condition as

$\displaystyle
i\frac{\partial}{\partial\tau_{1}}r_{+}(\tau)=-\beta_{+}(\bm{k},\bm{p})r_{k}(\tau)r_{p}(\tau)e^{i(\tilde{\omega}_{+}-\tilde{\omega}_{k}-\tilde{\omega}_{p})\tau_{1}},$
(S64)

with

$\displaystyle\beta_{+}(\bm{k},\bm{p})=\frac{\mu}{2}\biggl{[}\frac{\epsilon_{\bm{k}}+2\mu}{E_{\bm{k}}}+\frac{\epsilon_{\bm{p}}+2\mu}{E_{\bm{p}}}+\frac{\epsilon_{+}}{E_{+}}\biggl{(}3-\frac{\epsilon_{\bm{k}}+2\mu}{E_{\bm{k}}}\frac{\epsilon_{\bm{p}}+2\mu}{E_{\bm{p}}}\biggr{)}\biggr{]}.$
(S65)

Under this solvability condition, $\phi^{+}_{2}(t,\tau)$ is found as

$\displaystyle\phi^{+}_{2}(t,\tau)$
$\displaystyle=\biggl{(}1-\frac{\epsilon_{+}+2\mu}{E_{+}}\biggr{)}\tilde{r}_{+}(\tau)e^{iE_{+}t}+\biggl{(}1+\frac{\epsilon_{+}+2\mu}{E_{+}}\biggr{)}\tilde{r}^{\ast}_{+}(\tau)e^{-iE_{+}t}$
$\displaystyle\quad+(\varphi^{+}_{2}-\varphi^{+}_{4})r_{k}(\tau)r^{\ast}_{p}(\tau)e^{i(E_{\bm{k}}-E_{\bm{p}})t}+(\varphi^{+}_{2}+\varphi^{+}_{4})r^{\ast}_{k}(\tau)r_{p}(\tau)e^{-i(E_{\bm{k}}-E_{\bm{p}})t},$
(S66)

where $\tilde{r}_{+}(\tau)$ is an arbitrary complex amplitude, and
$\varphi^{+}_{2}$ and $\varphi^{+}_{4}$ are given as solutions of Eq. (S20).
Comparing Eq. (S66) with the previous solution (S15), one can see that the
divergent coefficients $\varphi^{+}_{1}$ and $\varphi^{+}_{3}$ are removed
once the solvability condition (S64) is imposed. In other words, the
solvability condition (S64) provides the interaction between the
$(\bm{k}+\bm{p})$-mode and the $\bm{k}$\- and $\bm{p}$-modes to remove the
divergence in the previous model.

For the coefficient functions of the other modes in the expansion (S50), only
the particular solution is sufficient for our purposes and can be obtained as
before.

##  B.0.2 Solvability conditions from the third-order equation

We expand $\phi_{3}(t,\tau,\bm{x})$ as

$\displaystyle\phi_{3}(t,\tau,\bm{x})=\phi^{k}_{31}(t,\tau)\cos(\bm{k}\cdot\bm{x})+\phi^{p}_{31}(t,\tau)\cos(\bm{p}\cdot\bm{x})+\phi^{+}_{31}(t,\tau)\cos((\bm{k}+\bm{p})\cdot\bm{x})+\cdots,$
(S67)

where the ellipsis represents terms with other irrelevant bases. From the
third-order differential equation (S48c), the differential equations for
$\phi^{k}_{31}(t,\tau)$ and $\phi^{+}_{31}(t,\tau)$ are found in the form of

$\displaystyle i\frac{\partial}{\partial
t}\phi^{k}_{31}(t,\tau)-\epsilon_{\bm{k}}\phi^{k}_{31}(t,\tau)-2\mu\mathop{\mathrm{Re}}[\phi^{k}_{31}(t,\tau)]=f^{k}_{31}(t,\tau),$
(S68a) $\displaystyle i\frac{\partial}{\partial
t}\phi^{+}_{31}(t,\tau)-\epsilon_{+}\phi^{+}_{31}(t,\tau)-2\mu\mathop{\mathrm{Re}}[\phi^{+}_{31}(t,\tau)]=f^{+}_{31}(t,\tau),$
(S68b)

with

$\displaystyle
f^{k}_{31}(t,\tau)=-i\frac{\partial}{\partial\tau_{1}}\phi^{k}_{21}(t,\tau)-i\frac{\partial}{\partial\tau_{2}}\phi^{k}_{1}(t,\tau)+\mu\biggl{\\{}-2\tilde{A}\sin(\omega
t)\mathop{\mathrm{Re}}[\phi^{k}_{1}(t,\tau)]+\frac{3}{4}|\phi^{k}_{1}(t,\tau)|^{2}\phi^{k}_{1}(t,\tau)$
$\displaystyle+|\phi^{p}_{1}(t,\tau)|^{2}\phi^{k}_{1}(t,\tau)+\frac{1}{2}\phi^{k\ast}_{1}(t,\tau)\phi^{p}(t,\tau)^{2}+|\phi^{+}_{1}(t,\tau)|^{2}\phi^{k}_{1}(t,\tau)+\frac{1}{2}\phi^{k\ast}_{1}(t,\tau)\phi^{+}(t,\tau)^{2}$
$\displaystyle+\mathcal{K}(\phi^{k}_{1}(t,\tau),2\phi_{20}(t,\tau)+\phi^{k}_{22}(t,\tau))+\mathcal{K}(\phi^{p}_{1}(t,\tau),\phi_{2}^{+}(t,\tau)+\phi^{-}_{2}(t,\tau))+\mathcal{K}(\phi^{+}_{1}(t,\tau),\phi^{p}_{21}(t,\tau))\biggr{\\}},$
(S69a) $\displaystyle
f^{+}_{31}(t,\tau)=-i\frac{\partial}{\partial\tau_{1}}\phi^{+}_{2}(t,\tau)-i\frac{\partial}{\partial\tau_{2}}\phi^{+}_{1}(t,\tau)+\mu\biggl{\\{}-2\tilde{A}\sin(\omega
t)\mathop{\mathrm{Re}}[\phi^{+}_{1}(t,\tau)]+\frac{3}{4}|\phi^{+}_{1}(t,\tau)|^{2}\phi^{+}_{1}(t,\tau)$
$\displaystyle+|\phi^{k}_{1}(t,\tau)|^{2}\phi^{+}_{1}(t,\tau)+\frac{1}{2}\phi^{+\ast}_{1}(t,\tau)\phi^{k}(t,\tau)^{2}+|\phi^{p}_{1}(t,\tau)|^{2}\phi^{+}_{1}(t,\tau)+\frac{1}{2}\phi^{+\ast}_{1}(t,\tau)\phi^{p}(t,\tau)^{2}$
$\displaystyle+\mathcal{K}(\phi^{+}_{1}(t,\tau),2\phi_{20}(t,\tau)+\phi^{+}_{22}(t,\tau))+\mathcal{K}(\phi^{p}_{1}(t,\tau),\phi_{21}^{k}(t,\tau))+\mathcal{K}(\phi^{k}_{1}(t,\tau),\phi_{21}^{p}(t,\tau))\biggr{\\}}.$
(S69b)

To find the solvability conditions from these differential equations, we take
their real and imaginary parts and diagonalize them as

$\displaystyle\frac{\partial}{\partial
t}\chi^{k}_{31}(t,\tau)-iE_{\bm{k}}\chi^{k}_{31}(t,\tau)$
$\displaystyle=g^{k}_{31}(t,\tau),$ (S70a)
$\displaystyle\frac{\partial}{\partial
t}\chi^{+}_{31}(t,\tau)-iE_{+}\chi^{+}_{31}(t,\tau)$
$\displaystyle=g^{+}_{31}(t,\tau),$ (S70b)

and their complex conjugates. Here, $\chi^{k/+}_{31}(t,\tau)$ and
$g^{k/+}_{31}(t,\tau)$ are introduced as

$\displaystyle\chi^{k}_{31}(t,\tau)=\frac{1}{\sqrt{2}}\biggl{(}\frac{\epsilon_{\bm{k}}+2\mu}{E_{\bm{k}}}\mathop{\mathrm{Re}}[\phi^{k}_{31}(t,\tau)]-i\mathop{\mathrm{Im}}[\phi^{k}_{31}(t,\tau)]\biggr{)},$
(S71a)
$\displaystyle\chi^{+}_{31}(t,\tau)=\frac{1}{\sqrt{2}}\biggl{(}\frac{\epsilon_{+}+2\mu}{E_{+}}\mathop{\mathrm{Re}}[\phi^{+}_{31}(t,\tau)]-i\mathop{\mathrm{Im}}[\phi^{+}_{31}(t,\tau)]\biggr{)},$
(S71b)

and

$\displaystyle g^{k}_{31}(t,\tau)$
$\displaystyle=\frac{1}{\sqrt{2}}\biggl{(}\frac{\epsilon_{\bm{k}}+2\mu}{E_{\bm{k}}}\mathop{\mathrm{Im}}[f^{k}_{31}(t,\tau)]+i\mathop{\mathrm{Re}}[f^{k}_{31}(t,\tau)]\biggr{)},$
(S72a) $\displaystyle g^{+}_{31}(t,\tau)$
$\displaystyle=\frac{1}{\sqrt{2}}\biggl{(}\frac{\epsilon_{+}+2\mu}{E_{+}}\mathop{\mathrm{Im}}[f^{+}_{31}(t,\tau)]+i\mathop{\mathrm{Re}}[f^{+}_{31}(t,\tau)]\biggr{)}.$
(S72b)

Eqs. (S70a) and (S70b) with $g^{k}_{31}(t,\tau)=0$ and $g^{+}_{31}(t,\tau)=0$
have oscillating solutions $e^{iE_{\bm{k}}t}$ and $e^{iE_{+}t}$, respectively.
Accordingly, the conditions for Eqs. (S70a) and (S70b) to have no secular
terms is equivalent to that $g^{k}_{31}(t,\tau)$ and $g^{+}_{31}$ have no
terms proportional to $e^{iE_{\bm{k}}t}$ and $e^{iE_{+}t}$, respectively. The
solvability conditions are found as

$\displaystyle
i\frac{\partial}{\partial\tau_{1}}\tilde{r}_{k}(\tau)+i\frac{\partial}{\partial\tau_{2}}r_{k}(\tau)=-i\mu\frac{\tilde{A}}{2}\frac{\epsilon_{\bm{k}}}{E_{\bm{k}}}r^{\ast}_{k}(\tau)e^{i\tilde{\omega}_{k}\tau_{1}}-\beta(\bm{k},\bm{p})\Bigl{(}\tilde{r}_{+}(\tau)r_{p}^{\ast}(\tau)+r_{+}(\tau)\tilde{r}_{p}^{\ast}(\tau)\Bigr{)}e^{i(\tilde{\omega}_{k}+\tilde{\omega}_{p}-\tilde{\omega}_{+})\tau_{1}}$
$\displaystyle+\mu\frac{5\epsilon_{\bm{k}}+3\mu}{E_{\bm{k}}}\Bigl{(}|r_{k}(\tau)|^{2}r_{k}(\tau)+\tilde{c}_{1}(\bm{k},\bm{p})|r_{p}(\tau)|^{2}r_{k}(\tau)+c_{2}(\bm{k},\bm{p})r_{p}(\tau)^{2}r^{\ast}_{k}(\tau)e^{2i(\tilde{\omega}_{k}-\tilde{\omega}_{p})\tau_{1}}\Bigr{)},$
(S73a) $\displaystyle
i\frac{\partial}{\partial\tau_{1}}\tilde{r}_{+}(\tau)+i\frac{\partial}{\partial\tau_{2}}r_{+}(\tau)=-\beta_{+}(\bm{k},\bm{p})\Bigl{(}\tilde{r}_{k}(\tau)r_{p}(\tau)+r_{k}(\tau)\tilde{r}_{p}(\tau)\Bigr{)}e^{i(\tilde{\omega}_{+}-\tilde{\omega}_{k}-\tilde{\omega}_{p})\tau_{1}}+\mu\frac{5\epsilon_{+}+3\mu}{E_{+}}|r_{+}(\tau)|^{2}r_{+}(\tau).$
(S73b)

Here, we dropped the cubic interaction terms between the
$(\bm{k}+\bm{p})$-mode and the $\bm{k}$\- and $\bm{p}$-modes, i.e.,
$|r_{+}(\tau)|^{2}r_{k/p}(\tau)$ and $|r_{k/p}(\tau)|^{2}r_{+}(\tau)$, as they
are not essential to remove the divergence in the previous model. The
coefficients $\beta(\bm{k},\bm{p})$, $\beta_{+}(\bm{k},\bm{p})$ and
$c_{2}(\bm{k},\bm{p})$ are given by Eqs. (S57), (S65) and (S37b),
respectively, and the coefficient $\tilde{c}_{1}(\bm{k},\bm{p})$ is given by

$\displaystyle\tilde{c}_{1}(\bm{k},\bm{p})$
$\displaystyle=\frac{E_{\bm{k}}}{5\epsilon_{\bm{k}}+3\mu}\biggl{[}4\frac{\epsilon_{\bm{k}}\epsilon_{\bm{p}}-\mu^{2}}{\epsilon_{\bm{p}}E_{\bm{k}}}-\frac{1}{2}\varphi^{-}_{1}\biggl{(}-\frac{\epsilon_{\bm{p}}+2\mu}{E_{\bm{p}}}+3\frac{\epsilon_{\bm{k}}}{E_{\bm{k}}}\biggr{)}-\frac{1}{2}(\varphi^{+}_{2}+\varphi^{-}_{2})\biggl{(}\frac{\epsilon_{\bm{p}}+2\mu}{E_{\bm{p}}}+3\frac{\epsilon_{\bm{k}}}{E_{\bm{k}}}\biggr{)}$
$\displaystyle\quad-\frac{1}{2}\varphi^{-}_{3}\biggl{(}1+\frac{\epsilon_{\bm{p}}+2\mu}{E_{\bm{p}}}\frac{\epsilon_{\bm{k}}}{E_{\bm{k}}}\biggr{)}-\frac{1}{2}(\varphi^{+}_{4}+\varphi^{-}_{4})\biggl{(}1-\frac{\epsilon_{\bm{p}}+2\mu}{E_{\bm{p}}}\frac{\epsilon_{\bm{k}}}{E_{\bm{k}}}\biggr{)}\biggr{]}.$
(S74)

##  B.0.3 The resulting amplitude equation

To combine the obtained solvability conditions and go back to the original
parameters, we introduce

$\displaystyle R_{k}(t)$ $\displaystyle=\Bigl{(}\varepsilon r_{k}(\varepsilon
t,\varepsilon^{2}t)+\varepsilon^{2}\tilde{r}_{k}(\varepsilon
t,\varepsilon^{2}t)\Bigr{)}e^{-i\varepsilon\tilde{\omega}_{k}t},$ (S75a)
$\displaystyle R_{p}(t)$ $\displaystyle=\Bigl{(}\varepsilon r_{p}(\varepsilon
t,\varepsilon^{2}t)+\varepsilon^{2}\tilde{r}_{p}(\varepsilon
t,\varepsilon^{2}t)\Bigr{)}e^{-i\varepsilon\tilde{\omega}_{p}t},$ (S75b)
$\displaystyle R_{+}(t)$ $\displaystyle=\Bigl{(}\varepsilon r_{+}(\varepsilon
t,\varepsilon^{2}t)+\varepsilon^{2}\tilde{r}_{+}(\varepsilon
t,\varepsilon^{2}t)\Bigr{)}e^{-i\varepsilon\tilde{\omega}_{+}t}.$ (S75c)

Then, with the use of the solvability conditions (S56), (S64), and (S73), the
time derivatives of these amplitudes are computed as

$\displaystyle i\frac{\mathrm{d}}{\mathrm{d}t}R_{k}(t)$ $\displaystyle=\Delta
R_{k}(t)-i\alpha
R_{k}^{\ast}(t)-\beta(\theta)R_{+}(t)R^{\ast}_{p}(t)+\lambda\Bigl{[}|R_{k}(t)|^{2}R_{k}(t)+\tilde{c}_{1}(\theta)|R_{p}(t)|^{2}R_{k}(t)+c_{2}(\theta)R_{p}(t)^{2}R^{\ast}_{k}(t)\Bigr{]},$
(S76a) $\displaystyle i\frac{\mathrm{d}}{\mathrm{d}t}R_{+}(t)$
$\displaystyle=\Delta_{+}(\theta)R_{k}(t)-\beta_{+}(\theta)R_{k}(t)R_{p}(t)+\lambda_{+}(\theta)|R_{+}(t)|^{2}R_{+}(t),$
(S76b)

where we assumed the absolute values of $\bm{k}$ and $\bm{p}$ to be equal to
focus on the angle between them. Then, with the use of
$\epsilon=\epsilon_{\bm{k}}=\epsilon_{\bm{p}}$ and $E=E_{\bm{k}}=E_{\bm{p}}$,
the coefficients become functions of the angle
$\theta=\angle(\bm{k},\bm{p})\in[0,\pi/2]$ given by

$\displaystyle\Delta_{+}(\theta)=\omega-
E_{+},\quad\beta(\theta)=\mu\biggl{(}\frac{\epsilon-\mu}{E}+\frac{\epsilon_{+}+2\mu}{E_{+}}\biggr{)},\quad\beta_{+}(\theta)=\mu\biggl{(}\frac{\epsilon+2\mu}{E}+\frac{\epsilon_{+}}{E_{+}}\frac{\epsilon-\mu}{\epsilon}\biggr{)},\quad\lambda_{+}(\theta)=\mu\frac{5\epsilon_{+}+3\mu}{E_{+}},$
(S77)

and

$\displaystyle\tilde{c}_{1}(\theta)$
$\displaystyle=\frac{\mu}{5\epsilon+3\mu}\biggl{[}4\frac{\epsilon^{2}-\mu^{2}}{\mu\epsilon}-\frac{(2\epsilon-\mu)(\epsilon+2\mu)+(2\epsilon^{2}+\mu^{2})\epsilon_{-}/(2\epsilon)}{E^{2}-E^{2}_{-}/4}+\biggl{(}\frac{2\epsilon+\mu}{\epsilon}\frac{2\epsilon+\mu}{\epsilon_{+}/2+\mu}+(\epsilon_{+}\rightarrow\epsilon_{-})\biggr{)}\biggr{]}.$
(S78)

The other coefficients $\Delta,\alpha,\lambda$, and $c_{2}(\theta)$ are the
same as before and given by Eqs. (S44) and (S41b). The coefficient
$\tilde{c}_{1}(\theta)$ is the coefficient $c_{1}(\theta)$ with the divergent
term at the angle satisfying $2E=E_{+}$ removed, and thus all coefficients
have no divergence in $\theta\in[0,\pi/2]$. Here, $2E=E_{-}$ cannot be
satisfied for $\theta\in[0,\pi/2]$ and therefore is not resonant; note that it
is sufficient to consider only the case with $\theta\in[0,\pi/2]$ because the
case with $\theta\in[\pi/2,\pi]$ can be taken into account with the
replacement $\theta\in[0,\pi/2]$ by $\pi-\theta$.

Finally, assuming that $\bm{k}$ and $\bm{p}$ satisfy the resonance condition
$\Delta=0$, and including phenomenological dissipation coefficients $\Gamma$
and $\Gamma_{+}$, we arrive at

$\displaystyle i\frac{\mathrm{d}}{\mathrm{d}t}R_{k}(t)$
$\displaystyle=-i\Gamma R_{k}(t)-i\alpha
R_{k}^{\ast}(t)-\beta(\theta)R_{+}(t)R^{\ast}_{p}(t)+\lambda\Bigl{[}|R_{k}(t)|^{2}R_{k}(t)+\tilde{c}_{1}(\theta)|R_{p}(t)|^{2}R_{k}(t)+c_{2}(\theta)R_{p}(t)^{2}R^{\ast}_{k}(t)\Bigr{]},$
(S79a) $\displaystyle i\frac{\mathrm{d}}{\mathrm{d}t}R_{+}(t)$
$\displaystyle=-i\Gamma_{+}R_{k}(t)+\Delta_{+}(\theta)R_{k}(t)-\beta_{+}(\theta)R_{k}(t)R_{p}(t)+\lambda_{+}(\theta)|R_{+}(t)|^{2}R_{+}(t),$
(S79b)

with $\Delta_{+}(\theta)=2E-E_{+}$, which are the amplitude equations of Eq.
(12) in the main text.