# Resonant _Kushi_ -comb-like multi-frequency radiation of oscillating two-
color soliton molecules

O Melchert1,2, S Willms1,2, I Oreshnikov3, A Yulin4, U Morgner1,2, I
Babushkin1,2 and A Demircan1,2 1 Leibniz Universität Hannover, Institute of
Quantum Optics, Welfengarten 1, 30167 Hannover, Germany  2 Leibniz Universität
Hannover, Cluster of Excellence PhoenixD, Welfengarten 1A, 30167 Hannover,
Germany  3 Max Planck Institute for Intelligent Systems, Max-Planck-Ring 4,
72076 Tübingen, Germany  4 Department of Nanophotonics and Metamaterials, ITMO
University, Birzhevaya line 16, 199034 St. Petersburg, Russia

###### Abstract

Nonlinear waveguides with two distinct domains of anomalous dispersion can
support the formation of molecule-like two-color pulse compounds. They consist
of two tightly bound subpulses with frequency loci separated by a vast
frequency gap. Perturbing such a two-color pulse compound triggers periodic
amplitude and width variations, reminiscent of molecular vibrations. With
increasing strength of perturbation, the dynamics of the pulse compound
changes from harmonic to nonlinear oscillations. The periodic amplitude
variations enable coupling of the pulse compound to dispersive waves,
resulting in the resonant emission of multi-frequency radiation. We
demonstrate that the location of the resonances can be precisely predicted by
phase-matching conditions. If the pulse compound consists of a pair of
identical subpulses, inherent symmetries lead to degeneracies in the resonance
spectrum. Weak perturbations lift existing degeneracies and cause a splitting
of the resonance lines into multiple lines. Strong perturbations result in
more complex emission spectra, characterized by well separated spectral bands
caused by resonant Cherenkov radiation and additional four-wave mixing
processes.

##  1 Introduction

Emission of radiation by physical objects is omnipresent in physical systems.
For instance, in classical electrodynamics, oscillating and linear
accelerating point charges generate electric and magnetic fields. Charged
particles confined to circular orbits emit synchrotron radiation, decelerating
charged particles produce Bremsstrahlung, and oscillating electric dipoles
cause dipole-radiation characteristic of radiative transitions in atoms in
molecules [1].

In the field of nonlinear optics many processes can be described by wave
equations of the type of the nonlinear Schrödinger equation (NSE) [2, 3],
which, in the integrable case, supports localized field pulses given by
solitary waves [4]. Such localized field pulses resemble particles, and can,
under the right circumstances, sustain periodic variations of their amplitude
and width. The emission of resonance radiation caused by such oscillations has
been studied in different nonlinear-optics settings. For example, considering
a NSE in which a varying dispersion coefficient perturbs transmitted pulses
periodically, dispersion-managed solitons where found to radiate through a
resonance mechanism. The underlying parametric resonance and the radiative
decay of dispersion-managed solitons were studied for rapidly varying
dispersion maps [5], and for weakly varying dispersion maps [6]. In dispersion
oscillating fibers, such parametric resonances where investigated numerically
and also experimentally for solitons (in case of anomalous dispersion) and
shock fronts (in case of normal dispersion) [7], demonstrating the generation
of resonant radiation at multiple resonant frequencies. In a standard NSE with
added third-order disperion, higher-order solitons where shown to support a
quite similar radiation mechanism [8]. Therein, the periodic variation of the
solitons peak intensity enabled the coupling of the localized state to
continuum modes on the opposing side of the zero-dispersion point. As a
result, resonant radiation in form of multi-frequency comb-like spectral bands
is emitted. The generation of such resonant radiation was shown to also
proceed in presence of the Raman-effect and pulse self-steepening [8]. In the
case of the Lugiato-Lefever equation with third-order dispersion, i.e. a
variant of the NSE with added damping an driving, a Hopf-bifurcation can
result in oscillating dissipative solitons that shed resonant radiation in a
series of well resolved spectral peaks [9].

In analogy to quantum mechanics, where molecules are formed as bound-states of
multiple atoms, nonlinear optics supports various types of molecule-like bound
states of multiple field pulses. For instance, usual soliton molecules which
consist of two identical subpulses separated by a fixed time-delay, can be
realized in a dispersion managed NSE [10, 11]. Mutually bound solitons also
arise in models of coupled NSEs [12, 13, 14, 15, 16, 17, 18, 19].
Specifically, for twin-core fibers with higher order dispersion, bound
solitons near the zero-dispersion point where studied and shown to sustain
periodic center-of mass and amplitude variations with different oscillation
periods, giving rise to resonant radiation with different frequencies [20].
Further, dissipative optical soliton molecules can be generated in mode-locked
fiber lasers [21, 22]. Quite recently, a different kind of molecule-like pulse
compound has been demonstrated theoretically [23, 24]. It consists of two
subpulses at distinctly different center frequencies, which superimpose to
form a single complex in the time-domain. An attractive, cross-phase
modulation (XPM) induced potential holds these two-color pulses compounds
together [23]. The binding mechanism that enables their stable propagation is
different from that of usual soliton molecules. It requires a propagation
constant with (at least) two separate domains of anomalous dispersion within
which group-velocity matched co-propagation of pulses is possible. A strong
incoherent attractive interaction between the subpulses of such a two-color
soliton, mediated accross (at least) two zero-dispersion points, is ensured if
the group-velocity mismatch between the pulses is negligible [25]. The
emphasis on group-velocity matching is similar as for the repulsion between a
soliton and a dispersive wave mediated across a single zero-dispersion point
[26, 27, 28, 29, 30], where quasi group-velocity matching is vital to ensure a
strong and efficient interaction mechanism. In terms of a NSE with added
negative fourth-order dispersion, two-color soliton molecules were identified
as members of a large family of generalized dispersion Kerr solitons [24].
Such two-color solitons were recently verified experimentally in mode-locked
laser cavities [31, 32]. Two-color soliton microcomb states of similar
structure also exist in the framework of the Lugiato-Lefever equation [33,
34]. The strong mutual interaction of the subpulses of such molecule-like
bound states [35] renders them ideal for studying the effects of oscillations
exhibited by coupled localized states in nonlinear optics systems.

Here, we discuss the rich dynamical behavior of such two-color soliton
molecules, consisting of group-velocity matched pulses located in distinct
domains of anomalous dispersion, separated by a vast frequency gap. We show
that nonstationary dynamics of the subpulses results in emission of resonant
radiation. The location of the newly generated multi-peaked spectral bands can
be understood by extending existing approaches [36, 37, 7, 20, 9] to two-
frequency pulse compounds [38]. Due to the manifold of two-frequency pulse
compounds with differing substructure, their emission spectra manifest in
various complex forms. Viewing a two-color soliton molecule as a superposition
of two distinct subpulses immediately suggests the existence of two different
types of dynamics: first, oscillations of the subpulses with respect to a
common center in time, and second, pure amplitude oscillations of the
subpulses. In previous studies of two-color soliton molecules [23, 25, 35], a
combination of both types of oscillations has generally been found. Therein,
the first type of oscillation was observed to be heavily damped [25, 35, 38].
In the presented work we focus on pure amplitude oscillations, as they
represent the persistent oscillation mode of two-frequency soliton molecules.

The article is organized as follows. In Sec. 2 we detail the propagation model
used for our theoretical investigations and we specify initial conditions that
directly generate special soliton-like two-color soliton pairs [39], which we
below also refer to as two-color soliton _molecules_ (or soliton molecules).
This approach surpasses previous “seeding” procedures, which result in a
largely uncontrolled generation of localized two-color pulse compounds. In
Sec. 3 we will perturb such a two-color soliton molecule by systematically
increasing its initial amplitude. In analogy to usual single-pulse nonlinear
Schrödinger solitons, where increasing the amplitude results in an oscillating
pulse, this triggers amplitude oscillations of the entire pulse compound,
yielding a vibrating soliton molecule. We numerically demonstrate that the
mutual driving of the subpulses results in the generation of resonant multi-
frequency radiation and we compare the observed location of the resonances to
analytic predictions obtained from a phase-matching analysis. We demonstrate
that symmetries, inherent in the propagation setting, can lead to degenerate
resonance conditions, and we show how existing degeneracies can be lifted by a
perturbation of the propagation constant. Section 4 concludes with summary.

##  2 Methods

#### Propagation model.

For studying the propagation dynamics of oscillating two-color pulse
compounds, we consider the modified nonlinear Schrödinger equation (NSE)

$\displaystyle
i\partial_{z}A=\left(\frac{\beta_{2}}{2}\partial_{t}^{2}-\frac{\beta_{4}}{24}\partial_{t}^{4}\right)A-\gamma|A|^{2}A,$
(1)

for a complex-valued field $A\equiv A(z,t)$, on a periodic time-domain of
extent $T$ with boundary conditions $A(z,-T/2)=A(z,T/2)$. In Eq. (1),
$\beta_{2}>0$ (in units of $\mathrm{fs^{2}/\mu m}$) is a positive-valued
group-velocity dispersion coefficient, $\beta_{4}<0$ ($\mathrm{fs^{4}/\mu m}$)
is a negative-valued fourth-order dispersion coefficient, and $\gamma$
($\mathrm{W^{-1}/\mu m}$) is the nonlinear coefficient. For the discrete set
of angular frequency detunings $\Omega\in\frac{2\pi}{T}\mathbb{Z}$, the
expressions

$\displaystyle
A_{\Omega}(z)=\mathsf{F}[A(z,t)]\equiv\frac{1}{T}\int_{-T/2}^{T/2}A(z,t)\,e^{i\Omega
t}~{}{\rm{d}}t,$ (2) $\displaystyle
A(z,t)=\mathsf{F}^{-1}[A_{\Omega}(z)]\equiv\sum_{\Omega}A_{\Omega}(z)\,e^{-i\Omega
t},$ (3)

specify forward [Eq. (2)], and inverse [Eq. (3)] Fourier transforms that
relate the field envelope $A(z,t)$ to the spectral envelope $A_{\Omega}(z)$.
Using Parseval’s identity for these transforms, the total energy in both
domains is given by

$\displaystyle
E(z)=\int_{-T/2}^{T/2}|A(z,t)|^{2}~{}{\rm{d}}t=T\sum_{\Omega}|A_{\Omega}(z)|^{2},$
(4)

wherein $|A(z,t)|^{2}$ ($\mathrm{W=J/s}$) specifies the instantaneous power,
and $|A_{\Omega}(z)|^{2}$ ($\mathrm{W}$) the power spectrum.

#### Propagation constant.

Using the identity $\partial_{t}^{n}\,e^{-i\Omega
t}=(-i\Omega)^{n}\,e^{-i\Omega t}$, the frequency-domain representation of the
propagation constant is given by the polynomial expression
$\beta(\Omega)=\frac{\beta_{2}}{2}\Omega^{2}+\frac{\beta_{4}}{24}\Omega^{4}$.
The inverse group-velocity of a mode at detuning $\Omega$ is given by
$1/v_{g}(\Omega)=\beta_{1}(\Omega)\equiv\partial_{\Omega}\beta(\Omega)$, and
the group-velocity dispersion is
$\beta_{2}(\Omega)\equiv\partial_{\Omega}^{2}\beta(\Omega)$. Assuming
$\beta_{2}>0$, and $\beta_{4}<0$, the group-velocity dispersion exhibits a
downward parabolic shape with zero-dispersion points defined by
$\beta_{2}(\Omega)=0$ at
$\Omega_{\rm{Z1},\rm{Z2}}=\mp\sqrt{2\beta_{2}/|\beta_{4}|}$. Subsequently we
set $\beta_{2}=1~{}\mathrm{fs^{2}/\mu m}$, and
$\beta_{4}=-1~{}\mathrm{fs^{4}/\mu m}$, resulting in the dispersion
characteristics shown in Figs. 1(a-c). We further set
$\gamma=1~{}\mathrm{W^{-1}/\mu m}$.

#### Propagation algorithm.

The numerical simulations reported in Sec. 3 are performed using an algorithm
with adaptive stepsize control, implementing the “Conservation quantity error”
method (CQE) [40, 41], which employs a conservation law of the underlying
propagation equation to guide selection of the propagation stepsize $h$.
Specifically, we use the relative error

$\displaystyle\delta_{E}(z)=\frac{|E(z+h)-E(z)|}{E(z)},$ (5)

where $E$ is the energy (4), conserved by the NSE (1). The CQE method keeps
$\delta_{E}$ within the goal error range
$[0.1\,\delta_{\rm{G}},\delta_{\rm{G}}]$, specified by the local goal error
$\delta_{\rm{G}}$, by decreasing $h$ when necessary while increasing $h$ when
possible. For our numerical experiments we use $\delta_{\rm{G}}=10^{-10}$. To
advance the field from position $z$ to $z+h$, we use the fourth-order “Runge-
Kutta in the interaction picture” method (RK4IP) [42].

#### Spectrograms.

To assess the time-frequency interrelations within the field $A$ at a selected
propagation distance $z$, we use the spectrogram [43, 44]

$P_{S}(t,\Omega)=\frac{1}{2\pi}\left|\int_{-T/2}^{T/2}A(z,t^{\prime})h(t^{\prime}-t)e^{-i\Omega
t^{\prime}}~{}{\rm d}t^{\prime}\right|^{2}.$ (6)

Therein, $h(x)=\exp(-x^{2}/2\sigma^{2})$ specifies a Gaussian window function
with root-mean-square width $\sigma$, used to localize the field in time.

#### Initial conditions.

Given the modified NSE (1), two-color pulse compounds can be “seeded” by using
initial conditions consisting of two superimposed fundamental solitons in the
form

$\displaystyle
A_{0}(t)=a_{1}\,{\rm{sech}}(t/t_{1})\,e^{-i\Omega_{1}t}+a_{2}\,{\rm{sech}}(t/t_{2})\,e^{-i\Omega_{2}t},$
(7)

where $a_{n}^{2}=|\beta_{2}(\Omega_{n})|/(\gamma\,t_{n}^{2})$, with
$n\in(1,2)$, are the peak intensities of the two subpulses which are subject
to the group-velocity matching condition
$\beta_{1}(\Omega_{1})=\beta_{1}(\Omega_{2})$. Upon propagation, such an
initial condition will evolve into a narrow localized pulse which can exhibit
internal dynamics, reminiscent of molecular vibrations, and is accompanied by
free radiation [23]. Recently we demonstrated that if the spectral extent of
each subpulse is small compared to the frequency gap separating them, the
mutual interaction of both pulses is effectively incoherent, allowing to
construct special solutions in the form of two-color soliton pairs [39].
Considering the propagation constant used in Eq. (1), an entirely symmetric
two-color soliton pair can be constructed by choosing both subpulse loci at
the two group-velocity matched frequencies
$\Omega_{1,2}=\mp\sqrt{6\beta_{2}/|\beta_{4}|}$, highlighted in Figs. 1(a-c).
They exhibit $\beta_{1}(\Omega_{1})=\beta_{1}(\Omega_{2})=\beta_{1}(0)=0$
[Fig. 1(b)], and $\beta_{2}(\Omega_{1})=\beta_{2}(\Omega_{2})=-2\beta_{2}<0$
[Fig. 1(c)]. Assuming two identical subpulse envelopes in Eq. (7), realized by
$t_{1}=t_{2}\equiv t_{0}$ and $a_{1}^{2}=a_{2}^{2}=2\beta_{2}/(3\gamma
t_{0}^{2})$ [39], then yields

$A_{0}(t)=F(0,t)\,\cos(\sqrt{2\beta_{2}/|\beta_{4}|}t),$ (8)

wherein

$F(z,t)=\sqrt{P_{0}}\,{\rm{sech}}(t/t_{0})\,e^{i\kappa z}$ (9)

specifies a pulse-enclosing envelope with peak power
$P_{0}=8\beta_{2}/(3\gamma t_{0}^{2})$ and nonlinear wavenumber
$\kappa=\beta_{2}/t_{0}^{2}$ [39]. These two-color soliton pairs correspond to
the fundamental metasoliton describing generalized dispersion Kerr solitons,
recently obtained in terms of a multi-scales analysis [24].

#### Fundamental two-color soliton molecule.

In Figs. 1(d-g), we demonstrate such a two-color soliton pair, having
$\Omega_{1,2}=\mp\sqrt{6}\,\mathrm{rad/fs}\approx\mp 2.449~{}\mathrm{rad/fs}$,
and $t_{0}=6~{}\mathrm{fs}$. Specifically, Fig. 1(f) shows the meta-envelope,
enclosing the initial field. The vast frequency gap between the subpulse loci
is evident in Fig. 1(g). The above condition for incoherent coupling,
requiring the frequency separation of the pulses, given by
$\Omega_{\rm{sep}}\equiv|\Omega_{2}-\Omega_{1}|=\sqrt{24\beta_{2}/|\beta_{4}|}$,
to be large compared to the spectral extent of the individual subpulses, i.e.
$\Omega_{0}=2/(\pi t_{0})$, is clearly satified: the example in Fig. 1 has
$\Omega_{\rm{sep}}/\Omega_{0}\approx 46$. The soliton-like propagation of this
fundamental two-color soliton pair is demonstrated in Figs. 1(d,e).

Figure 1: Idealized setup yielding a two-color soliton molecule. (a)
Propagation constant, (b) inverse group velocity, and, (c) group-velocity
dispersion. Circle and diamond highlight the two-color solitons subpulse loci.
The domain of normal dispersion is shaded gray. In (c), A1 and A2 label
distinct domains of anomalous dispersion, and N labels the domain of normal
dispersion. (d) Intensity, and, (e) spectrum of a fundamental two-color
soliton molecule for the propagation range $z=0\ldots 15\,z_{0}$ with
$z_{0}=t_{0}^{2}/(2\beta_{2})$. (f) Initial intensity, and, (g) initial
spectrum at $z=0$. In (f), A labels the scaled intensity $|A|^{2}/P_{0}$, and
B labels the squared magnitude of the meta-envelope
$|F|^{2}/P_{0}={\rm{sech}}^{2}(t/t_{0})$.

##  3 Results

To investigate the radiation that emanates from an oscillating two-color
soliton molecule, we need to realize stable $z$-periodic dynamics.
Specifically, in Sect. 3.1 we demonstrate how to induce amplitude and width
oscillations of soliton molecules and we study their characteristics. In Sect.
3.2 we demonstrate that oscillating soliton molecules shed resonant radiation
in the form of multi-peaked spectral bands, and we compare our numerical
results to predictions based on a phase matching analysis. In Sect. 3.3 we
show that possible degeneracies in the resulting radiation can be lifted by
perturbations.

###  3.1 Oscillating two-color soliton molecules

#### Amplitude oscillations.

So as to excite amplitude and width oscillations of two-color soliton
molecules, we consider initial conditions of the form of Eqs. (8,9), and
include an additional dimensionless factor $N$ through $P_{0}\to
P_{0}^{\prime}\equiv N^{2}P_{0}$, allowing to increase thir initial amplitude.
The choice $N=1$ yields the basic soliton molecule of Eq. (8) [Fig. 1(f)], for
which the pulse-enclosing envelope $F$ is coined the _fundamental metasoliton_
in Ref. [24]. In analogy with the nomenclature of usual single-pulse nonlinear
Schrödinger solitons, choosing $N>1$ define higher-order metasolitons [24],
resulting in oscillating soliton molecules. For brevity, we will refer to $N$
simply as the _order_ of the soliton molecule. Figure 2 demonstrates such
oscillating soliton molecules for $N=1.1$ [Fig. 2(a)], $N=1.4$ [Fig. 2(b)],
and $N=1.65$ [Fig. 2(c)]. If the initial perturbation of the two-color soliton
molecule is small ($N=1.1$), the variation of the scaled peak intensity is
largely harmonic [Fig. 2(d)], i.e. a single spatial Fourier mode describes the
observed intensity variation [Fig. 2(g)]. As evident from Figs. 2(e,h) and
2(g,i), increasing the molecule order results in nonlinear peak intensity
oscillations, and, consequently, an increasing number of spatial Fourier
modes. Thereby, the variation of the molecule peak intensity is characterized
by an oscillation period $\Lambda$ that decreases with increasing order $N$.

Figure 2:  Vibration characteristics of two-color soliton molecules. (a-c)
Time-domain propagation dynamics for (a) $N=1.1$, (b) $N=1.4$, and, (c)
$N=1.65$. The horizontal dashed lines indicate $z=(\pi/2)\,z_{0}$, i.e. the
oscillation period of a higher-order standard NSE soliton with dispersion
length $z_{0}$. (d-f) Variation of the scaled peak intensity followed for 8
consecutive oscillation periods for (d) $N=1.1$ with $\Lambda=4.30\,z_{0}$,
(e) $N=1.4$ with $\Lambda=1.93\,z_{0}$, and, (f) $N=1.65$ with
$\Lambda=1.13\,z_{0}$. (g-i) Spatial Fourier modes corresponding to the
intensity variations shown in (d-f).

#### Amplitude oscillation wavelength.

To quantify these dynamic amplitude and width variations we performed a series
of numerical simulations for soliton molecules with different initial duration
and order in the range $[1.05,2]$, see Fig. 3. For $N\to 1$, the observed
oscillation period agrees well with that of standard nonlinear Schrödinger
solitons of order $N$ governed by a standard NSE, indicated by the solid line
(labeled C) in Fig. 3. This finding is reasonable, since, for $N\approx 1$,
the coupling between the molecules subpulses remains effectively incoherent so
that a given subpulse can be described by a standard NSE with modified
coefficients [39]. For $N\to 2$, the oscillation period of standard NSE
solitons approaches the value $z_{\rm{sol}}=(\pi/2)z_{0}$, i.e. the soliton
period of a higher-order (single) soliton with dispersion length $z_{0}$,
indicated by the horizontal line in Fig. 3. We find that the oscillation
period of soliton molecules with large duration $t_{0}=12\,\mathrm{fs}$,
labeled B in Fig. 3, are in excellent agreement with this behavior. This
finding is consistent with results obtained for a $N=2$ metasoliton of
duration $t_{0}=\sqrt{2|\beta_{4}|/(3\beta_{2}\epsilon^{2})}\approx
16.3\,\mathrm{fs}$ ($\epsilon=0.05$), studied in Ref. [24]. In contrast,
soliton molecules with smaller duration $t_{0}=6\,\mathrm{fs}$, labeled A in
Fig. 3, deviate from this behavior for $N\gtrapprox 1.6$. From our numerical
experiments we found that for soliton molecules with small duration, even
small amplitude oscillations suffice to cause a molecule to adiabatically
change its shape, affecting its oscillation period and explaining the observed
trend of the data. Overall, the close correspondence between the oscillatory
behavior of soliton molecules and standard NSE solitons allows to directly
transfer approximate descriptions of the evolution of near-fundamental
standard NSE solitons, obtained in terms of variational approaches [45, 46],
to the present case. A simple qualitative model of the observed subpulse
oscillations can be derived on basis of a coupled system of variational
equations, governing the evolution of the parameters of a hyperbolic-secant
trial pulse [46]. From this, an equation of motion can be obtained that
captures the variation of the pulse width about some preferred, stationary
value (see A). In lowest order approximation, it describes a harmonic
oscillation. The corresponding wavelength

$\Lambda^{\rm{HO}}(N)=\frac{\pi^{2}}{\sqrt{2}N^{2}}z_{0},$ (10)

derived in A, is in qualitative agreement with our numerical results, see the
dashes line (labeled D) in Fig. 3. A less intuitive, but more accurate
estimate of the oscillation period can be obtained by using Gaussian trial
functions to approximate NSE solitons in a variational analysis of the
standard NSE, yielding the approximate oscillation period [45]

$\Lambda^{\rm{var}}(N)=\frac{\pi\sqrt{2}N^{2}}{(\sqrt{2}N^{2}-1)^{3/2}}z_{0},$
(11)

which is in good agreement with our numerical results, as shown by the dash-
dotted line (labeled E) in Fig. 3. A fit of the molecule data to a polynomial
model of the form

$\Lambda^{\rm{fit}}(N)=z_{0}\sum_{n=0}^{3}c_{n}N^{-n},$ (12)

yields the parameters $(c_{0},c_{1},c_{2},c_{3})=(-15.1,64.3,-93.2,56.2)$ for
soliton molecules with initial duration $t_{0}=6\,\mathrm{fs}$, shown as short
dashed line (labeled E) in Fig. 3. For molecules with $t_{0}=12\,\mathrm{fs}$,
a fit yields $(c_{0},c_{1},c_{2},c_{3})=(-3.8,24.7,-47,39)$.

Figure 3:  Change of the oscillation wavelength $\Lambda$ with soliton
molecule order $N$. Data points, labeled A and B, show numerical results for
soliton molecules with initial duration $t_{0}=6\,\mathrm{fs}$ and
$t_{0}=12\,\mathrm{fs}$, respectively. The oscillation period is scaled by
$z_{0}=t_{0}^{2}/(2\beta_{2})$. Solid line, labeled C, shows numerical results
for standard NSE solitons of order $N$. Dashed line, labeled D, shows results
obtained by assuming a harmonic oscillation of the molecules subpulse widths.
Dash-dotted line, labeled E, shows the approximate oscillation period obtained
from a variational approach. Short dashed line, labeled F, shows a fit of Eq.
(12) to the data for molecules with $t_{0}=6\,\mathrm{fs}$.

###  3.2 Degenerated multi-frequency radiation

Below we show that the nonstationary dynamics of soliton molecules, which is
reminiscent of molecular vibrations, results in emission of resonant multi-
frequency radiation.

#### Example for $N=1.65$.

As a first example, we consider a soliton molecule with duration
$t_{0}=6\,\mathrm{fs}$, and order $N=1.65$. Figure 4(a) shows the time-domain
intensity on a logarithmic scale [cf. Fig. 2(c)], indicating that, in response
to the initial perturbation, the localized molecule sheds radiation along
either $t$-direction. The periodic variation of the spectrum, accompanying the
amplitude and width oscillations of the molecule with wavelength
$\Lambda/z_{0}=2.23$, is shown in Fig. 4(b). Distinct, comb-like bands of
frequencies can be seen to develop in the vicinity of the subpulse loci
$\Omega_{1}$ and $\Omega_{2}$. These newly generated spectral bands appear
particularly pronounced when the molecule assumes its maximum time-domain
width during an oscillation period, so that the resonances are not drowned by
the spectrum of the molecule, demonstrated for the specific choice
$z/z_{0}\approx 28.6$ in Figs. 4(b,d). The spectral location of this resonant
radiation can be predicted in terms of wave number matching conditions,
derived for the dynamically evolving, oscillating molecule state. For clarity,
below we summarize the derivation of the resonance conditions detailed in Ref.
[38]. Considering an oscillating soliton molecule, consisting of two subpulses
$U_{1}$ and $U_{2}$, we can account for amplitude and width oscillations by
representing both in terms of a Fourier series [38]

$U_{n}(z,t)=\sum_{m}C_{nm}(t)\,\exp\left(i\kappa_{n}z+im\frac{2\pi}{\Lambda}z\right),\quad
n\in(1,2),~{}m\in\mathbb{Z}$ (13)

wherein $C_{nm}$ are expansions coefficients, $\kappa_{n}$ are the wavenumbers
that govern the $z$-propagation of the subpulses, and $\Lambda$ is the
wavelength of the periodic $z$-variation of the molecule. Below, we abbreviate
the wavenumbers associated with the harmonics of the molecules $z$-oscillation
wavelength $\Lambda$ by $K_{m}=2\pi m/\Lambda$. Under the assumption that
these subpulses satisfy a pair of coupled NSEs [25, 39], and by assessing how
dispersive radiation is excited by their mutual periodic driving [38], the
resonance conditions

$\displaystyle D_{n}(\Omega_{RR})-\kappa_{n}=K_{m},$ $\displaystyle\qquad
n\in(1,2),~{}m\in\mathbb{Z}$ (14) $\displaystyle
D_{n}(\Omega_{RR})-2\kappa_{n}+\kappa_{n^{\prime}}=K_{m},$
$\displaystyle\qquad n,n^{\prime}\in(1,2),~{}n\neq
n^{\prime},~{}m\in\mathbb{Z}$ (15)

can be obtained. Therein,
$D_{n}(\Omega)\equiv\beta(\Omega)-\beta(\Omega_{n})-\beta_{1}(\Omega)(\Omega-\Omega_{n})$,
for $n=1,2$, defines the dispersion profile in the reference frame of each of
the subpulses, and $\Omega_{RR}$ specifies frequencies at which resonant
radiation (RR) is excited. The individual subpulse wavenumbers are related to
the subpulse amplitudes $a_{1}$ and $a_{2}$ through

$\displaystyle\kappa_{1}=\frac{\gamma}{2}(a_{1}^{2}+2a_{2}^{2}),$ (16)
$\displaystyle\kappa_{2}=\frac{\gamma}{2}(a_{2}^{2}+2a_{1}^{2})+\beta(\Omega_{2})-\beta(\Omega_{1}).$
(17)

Equation (14) yields resonance conditions for the generation of Cherenkov
radiation emitted by each of the molecules subpulses, related to the $m$-th
harmonic of the driving period $\Lambda$. Equation (15) defines further
resonance conditions indicative of four-wave mixing (FWM) processes that
involve both subpulses.

Figure 4: Resonant radiation of symmetric two-color soliton molecules. (a-d)
Results for soliton molecules of order $N=1.65$. (a) Evolution of the
intensity on a logarithmic scale. (b) Evolution of the spectrum. Horizontal
short dashed line indicates $z/z_{0}=28.6$. (c) Dispersion profile and
graphical solution of the resonance conditions Eqs. (14,15) for
$z$-oscillation periods of order $m$. (d) Spectrum at $z/z_{0}\approx 28.6$.
Vertical lines in (c,d) indicate loci of resonant frequencies. (e-h) Same for
$N=2.0$. In (h) $z/z_{0}\approx 29.73$. Vertical dashed lines in (d,h)
indicate approximate loci of the resonances for $m=-1$, given by Eq. (19).

In Ref. [38], the subpulse wavenumbers $\kappa_{1}$ and $\kappa_{2}$ had to be
extracted _a posteriori_ from the simulated propagation dynamics. In the
present case, we can establish a reasonable _a priori_ estimate by assuming
the molecule to be close to a two-color soliton pair with subpulse amplitudes
$a_{1}=a_{2}=\sqrt{2\beta_{2}/(3\gamma t_{0}^{2})}$ [39] and
$\beta(\Omega_{1})=\beta(\Omega_{2})=3\beta_{2}^{2}/(2|\beta_{4}|)$, giving

$\kappa_{1}=\kappa_{2}=\frac{\beta_{2}}{t_{0}^{2}}.$ (18)

Consequently, the resonance conditions (14) and (15) are degenerate. The
wavenumber profile in the reference frame of the subpulse at $\Omega_{1}$,
given by $D_{1}(\Omega)-\kappa_{1}$, is shown in Fig. 4(c). In the reference
frame of the subpulse at $\Omega_{2}$, the wavenumber profile
$D_{2}(\Omega)-\kappa_{2}$ looks identical. This is a consequence of the
symmetry of the propagation constant about the point $\Omega=0$, and the
symmetric choice $-\Omega_{1}=\Omega_{2}=\sqrt{6}\,\mathrm{rad/fs}$ of the
supulse loci. A graphical solution of the resonance conditions (14) and (15)
can be obtained by superimposing the wavenumber profile with the wavenumbers
$K_{m}$, associated with the harmonics of the molecules $z$-oscillation
wavelength. This is demonstrated in Fig. 4(c), where $K_{m}$ for $m=-8\ldots
1$ are included as horizontal dashed lines. The crossing points, emphasized by
dots in Fig. 4(c), indicate frequency loci $\Omega_{RR}$ at which the coupling
of energy into dispersive waves is expected to be efficient. As shown by the
vertical lines in Figs. 4(c,d), we find that these resonant frequencies are in
excellent agreement with our numerical results. The large number of resonant
frequencies, excited and continuously fed during propagation of the molecule,
is related to the large number of modes required by a Fourier decomposition of
the nonlinear vibration of the molecule for $N=1.65$ [cf. Figs. 2(f,i)].
Molecules of order $N=1.4$ yield fewer and less pronounced resonances. As
shown by the spectrogram Fig. 6(a), all resonantly excited modes leave the
molecule at a characteristic group-velocity, characteristic for their location
in the vicinity of $\Omega_{1}$ and $\Omega_{2}$: all resonantly excited modes
in $[0,\Omega_{1}]$ and $[\Omega_{2},\infty]$ propagate towards larger
positive times; those in $[0,\Omega_{2}]$ and $[\Omega_{1},\infty]$ propagate
towards increasingly negative times. We can simplify the above resonance
condition to obtain a simple expression that specifies the location of the
first few resonances reasonably well: by making a parabolic approximation of
the propagation constant in the vicinity of $\Omega_{1}$ and $\Omega_{2}$, and
by neglecting the contribution of the subpulse wavenumbers $\kappa_{1}$ and
$\kappa_{2}$, we find

$\Omega_{\rm{RR},m}^{(1,2),\pm}\approx\Omega_{1,2}\pm\sqrt{\frac{|K_{m}|}{\beta_{2}}}.$
(19)

The approximate positions of the resonances for $m=-1$ are indicated by
vertical dashed lines in Figs. 4(d).

#### Example for $N=2$.

According to Eq. (19) it should be possible to control the separation between
consecutive resonances by either decreasing $\beta_{2}$, resulting in a more
“flat” dispersion profile $D_{1}(\Omega)$ in Fig. 4(c), or by increasing the
wavenumbers $K_{m}$, resulting in more widely spaced horizontal dashed lines
in Fig. 4(c). To verify this prediction we next consider a soliton molecule of
order $N=2$. Since such a molecule exhibits $\Lambda(N=2)<\Lambda(N=1.65)$
[Fig. 3], we expect to observe resonant frequencies with larger spacing as in
the previously considered case. As evident from Figs. 4(e,f), the evolution of
the molecule proceeds in two distinct stages: during the initial propagation
stage $z/z_{0}\lessapprox 10$, a large amount of dispersive radiation emanates
from the molecule; then, for $z/z_{0}\gtrapprox 10$, the molecule seems to
have reached a steadily oscillating stage. The initial burst of energy results
in an adiabatic reshaping of the molecule in course of which its temporal
width decreases to about $t_{0}\approx 4\,\mathrm{fs}$, and its spectrum
broadens. Consequently, the molecule also develops a notable spectral density
within the domain of normal dispersion [Fig. 4(f)]. The wavelength of the
$z$-periodic oscillations also decrease slightly from initially
$\Lambda/z_{0}\approx 1.1$ to $\Lambda/z_{0}\approx 0.92$ for
$z/z_{0}\gtrapprox 10$. This causes the wavenumbers $K_{m}$, indicated by
horizontal dashed lines in Fig. 4(g), to shift downwards, and the
corresponding resonant frequencies to shift away from the nearest subulse
location. The associated spectral lines then appear broadened [Fig. 4(h)]. A
similar observation has been made in a study of the resonance mechanism of
higher-order solitons in a NSE with added third order dispersion [8]: adding
the Raman-effect, which causes a self-frequency-shift of solitons, was found
to wash out the resonantly excited radiation bands. The approximate positions
of the resonances for $m=-1$ are again indicated by vertical dashed lines in
Figs. 4(h). Contrasting the dynamics for $N=2$ with that for $N=1.65$ provides
an unexpected observation: since the possible driving wavenumbers $K_{m}$ are
now more widely spaced, see the horizontal dashed lines in Fig. 4(g), the
molecule can support only fewer resonant frequencies; as a result, its
amplitude and width oscillations appear more harmonic. As evident from the
spectrogram Fig. 6(b), due to the aforementioned decrease of $\Lambda$ for
increasing $z$, some resonances, which are present during the initial
“nonlinearly oscillating” propagation stage ($z/z_{0}\lessapprox 10$), die out
and are absent in the later “harmonically oscillating” propagation stage
($z/z_{0}\gtrapprox 10$). Let us note that for both cases considered above,
the particular resonance for $m=0$, indicating the excitation of the usual
Cherenkov radiation in absence of amplitude and width oscillations [47], is
not excited.

In both cases studied above, we find that the amplitude and width oscillations
of the soliton molecules are quasistationary and decay only slowly due to the
emission of resonant radiation.

Figure 5: Resonant radiation of a two-color soliton molecule of order $N=1.65$
for a non-symmetric propagation constant, perturbed by
$\beta_{3}=-0.2\,\mathrm{fs^{3}/\mu m}$. (a) Evolution of the intensity on a
linear scale. (b) Evolution of the spectrum. A filtered view of the spectral
components enclosed by the box (labeled A) in (b) is shown within the box
(labeled A) in (a). (c) Dispersion profile as seen for the subpulse at
$\Omega_{1}$, and graphical solution of the resonance conditions (14,15) for
$z$-oscillation periods of order $m=-5\ldots 5$. (d) Spectrum at
$z/z_{0}\approx 28.5$. Vertical lines in (c,d) indicate loci of resonant
frequencies. B and E label the loci of FWM resonances [Eqs. (15)] expected in
absence of oscillations. Similarly, C and D label those of Cherenkov
resonances [Eqs. (14)] in absence of oscillations.

###  3.3 Lifting the degeneracies by perturbations

In order to lift the degeneracy of the resonance conditions (14,15), which
occurs for $\kappa_{1}=\kappa_{2}$, Eqs. (16,17) suggest two possibilities:
either setting $a_{1}\neq a_{2}$ (see B), i.e. assigning different amplitudes
for both subpulses, or devising a propagation scenario for which
$\beta(\Omega_{1})\neq\beta(\Omega_{2})$. Subsequently, we opt for the latter
and assess the impact of a non-symmetric propagation constant on the dynamics
of a soliton molecule. Specifically, we modify the propagation constant to
$\beta(\Omega)=\frac{\beta_{2}}{2}\Omega^{2}+\frac{\beta_{3}}{6}\Omega^{3}+\frac{\beta_{4}}{24}\Omega^{4}$
using $\beta_{3}=-0.2\,\mathrm{fs^{3}/\mu m}$. This breaks the symmetry
inherent in the propagation scenarios studied above. In Fig. 5 we show how a
soliton molecule with $N=1.65$, and subpulse center frequencies again at
$\Omega_{1}=-\Omega_{2}=\sqrt{6}\,\mathrm{rad/fs}$, responds to such a
perturbation. In this case, the effective values of second order dispersion
felt by both subpulses are
$\beta_{2}(\Omega_{1})\approx-1.51\,\mathrm{fs^{2}/\mu m}$, and
$\beta_{2}(\Omega_{2})\approx-2.49\,\mathrm{fs^{2}/\mu m}$. We then start with
a (fundamental) two-color soliton pair in the form of Eq. (7) with
$a_{1}=\sqrt{(2\alpha-1)|\beta_{2}(\Omega_{1})|/(3\gamma t_{0}^{2})}$ and
$a_{2}=\sqrt{(2\alpha^{-1}-1)|\beta_{2}(\Omega_{2})|/(3\gamma t_{0}^{2})}$,
with $\alpha=\beta_{2}(\Omega_{2})/\beta_{2}(\Omega_{1})\approx 1.65$ [39],
and increase its amplitude by the factor $N$. As a first observation, we note
that the molecule acquires an additional shift in its group-velocity due to
$\beta_{3}\neq 0$. In Fig. 5(a) we compensated for the resulting drift by
transforming to a reference frame with velocity $v_{0}=-1.877\,\mathrm{\mu
m/fs}$, within which the molecule appears stationary. Further, the spectrum in
Figs. 5(b,d) is now very different from that in Figs. 4(b,d). The soliton
molecule exhibits a non-symmetric spectrum, consisting of a strong pulse,
centered at $\Omega_{1}$, and a weaker pulse, centered at $\Omega_{2}$. In
such a scenario, due to a strong entanglement of the subpulses [35], we expect
that the strong pulse drives the weak pulse so that the wavelength of the
intensity and width oscillations are equal for both pulses. From our numerical
simulation we find $\Lambda/z_{0}=3.2$. Below we perform a phase-matching
analysis in the frame of reference of the stronger subpulse (i.e. the one at
$\Omega_{1}$). In this case, we have to resort to measuring the subpulse
wavenumbers $\kappa_{1}$ and $\kappa_{2}$ _a posteriori_ from the simulated
dynamics. A graphical solution of the resonance conditions (14) for
$m=-5\ldots 5$ shows that now there are also multi-peaked spectral bands of
resonances excited around $m=0$, indicated by the red short-dashed line in
Fig. 5(c). A filtered view of the resonant radiation at $\Omega\approx
0.6\,\mathrm{rad/fs}$, enclosed by a box (labeled A) in Fig. 5(b), is shown in
the box (labeled A) in Fig. 5(a). As evident from this filtered view, the
resonant radiation is emitted in a pulse-wise fashion and with decreasing
intensity for increasing propagation distance [48]. The phase-relationship
between theses separated pulses then yields the multi-peaked spectral band in
Fig. 5(d). Additional bands of resonances, associated with Cherenkov
resonances, are visible at $\Omega\approx\Omega_{1}$ and $\Omega\approx
3.5\,\mathrm{rad/fs}$. For clarity, we provide a graphical solution of the
phase-matching conditions (14) for $m=-5\ldots 5$, only. Similar effects for
perturbation by $\beta_{3}$ have previously been reported for higher-order
solitons in the NSE [8], and dissipative solitons in fiber cavities [9].
Additional four-wave mixing (FWM) resonances arising from Eq. (15) can be
identified at $\Omega\approx-3.8\,\mathrm{rad/fs}$. In contrast to an earlier
approach, where such resonances where caused by the nonlinear mixing of
solitons and dispersive waves [37], we here find that these resonances are
caused by a mixing of the subpulses of a soliton molecule.

All of the resonantly generated radiation emanates from the molecule in the
vicinity of the Cherenkov, and FWM resonances specified by Eqs. (14,15). The
pulse-wise emission is clearly evident in the spectrogram Fig. 6(c). In Fig.
6(c) we did not compensate for the non-vanishing group velocity of the
molecule, so that, at $z/z_{0}\approx 28$, it has drifted to
$t\approx-300\,\mathrm{fs}$. In the non-oscillating case, the observed multi-
frequency spectral bands would degenerate into single, sharp spectral lines
[47, 36, 37]. Specifically, the points in the spectrum labeled C and D in Fig.
4(d) correspond to the location of Cherenkov radiation, and the points labeled
B and E correspond to additional FWM processes in absence of oscillations.
Finally, let us note that in the present propagation scenario we do not
observed resonant radiation in the vicinity of the point labeled E in Fig.
5(d). As discussed in Refs. [37, 38], this can be attributed to a low
efficiency of the corresponding resonant process. In comparison to the earlier
two cases studies in Sect. 3.2, we here find that the amount of resonant
radiation emanating from the soliton molecule is larger and that the amplitude
oscillations decay faster. Let us note that lifting the degeneracies between
the resonance conditions (14,15) through weak perturbations first leads to a
split-up of the resonance lines into multiple lines (see B). Only when
increasing the strength of the perturbation, complex spectra as in Fig. 5 are
obtained.

As pointed out above, a viable means to lift existing degeneracies is to
consider pulses at a group velocity matched pair of frequencies for which
$\beta(\Omega_{1})\neq\beta(\Omega_{2})$. Above, this was achieved by
considering a non-symmetric propagation constant. Similarly, this can be
achieved by retaining a symmetric propagation constant but shifting the
subpulses to a group-velocity matched pair of frequencies different from the
specific pair $\Omega_{1,2}$ considered above. In addition, allowing a
frequency dependence of the nonlinear coefficient $\gamma$ in the nonlinear
part of Eq. (1) directly breaks the symmetry of the subpulses that comprise a
molecule. Propagation models of this type were considered in previous studies
in Refs. [23, 25, 35, 38].

Figure 6: Spectrograms of vibrating soliton molecules and their resonantly
shed radiation. (a) Spectrogram of the oscillating molecule for $N=1.65$ at
$z/z_{0}\approx 28.6$ [Figs. 4(a-d)], showing many narrowly spaced resonances.
The spectrogram is computed for $\sigma=24\,\mathrm{fs}$ in Eq. (6). (b)
Spectrogram ($\sigma=24\,\mathrm{fs}$) of a molecule for $N=2$ at
$z_{0}/z=29.73$ [Figs. 4(e-h)], showing few widely spaced resonances. (c)
Spectrogram ($\sigma=16\,\mathrm{fs}$) for a molecules with $N=1.65$ and a
non-symmetric propagation constant [Fig. 5]. The horizontal short-dashed lines
indicate zero-dispersion points. The dashed line indicates
$\beta_{1}(\Omega)\,z$, i.e. the light-cone delimiting temporal position of a
mode at $\Omega$ emitted at $z=0$. Movies of the dynamics are provided as
Supplementary Material.

##  4 Summary and conclusions

We studied two-color soliton molecules, i.e. pulse compounds consisting of two
tightly bound subpulses at vastly distinct frequencies, supported by nonlinear
waveguides with two distinct domains of anomalous dispersion.

First, we analyzed the amplitude and width oscillations of such molecules,
induced by increasing the molecule order in a systematic way. The order of a
molecule is defined similar to the usual order of solitons in the standard
NSE. We contrasted the periodic dynamics exhibited by two-color soliton
molecules with that of single solitons in the standard NSE, finding that, in
the limit of large molecule duration and despite their complex internal
structure, their $z$-oscillation wavelength matches that of a standard NSE
soliton of the same order. By elaborating on a well established variational
analysis of the NSE we were led to expect that near-fundamental soliton
molecules oscillate harmonically, whereas the amplitude oscillations quickly
grow nonlinear with increasing molecule order. This simple and intuitive
picture was found to be in good qualitative agreement with our numerical
results.

Second, we investigated the emission spectrum of soliton molecules, excited by
their periodic amplitude and width oscillations. As soon as a molecule
oscillates, it sheds resonant radiation in the form of multi-frequency
spectral bands, reminiscent of the shape of traditional Japanese _Kushi_
combs. The frequency spacing of resonantly excited modes is determined by the
$z$-oscillation wavelength of the molecule. We demonstrated that the location
of the resonant frequencies can be reliably predicted by means of resonance
conditions resulting from a phase-matching analysis. Specifically, there are
two types of resonance conditions which separately determine the generation of
Cherenkov radiation and further FWM processes. We showed that if the
propagation scenario is entirely symmetric, these resonance conditions are
redundant, so that the generated radiation is degenerate. Breaking the
symmetry of the propagation scenario, which we here accomplished by perturbing
the propagation constant, the degeneracy can be lifted and distinct multi-
frequency spectral bands uniquely associated with Cherenkov radiation, and
further FWM processes, are obtained.

The occurrence of such multi-frequency radiation, especially in the degenerate
form, constitutes a fundamental phenomenon in nonlinear waveguides. The
presented study clarifies the mechanism underlying the generation of resonant
radiation of two-color pulse solitons, observed previously in terms of
theoretical studies [23], and recently also in experiments on multi-wavelength
mode-locked solitons in laser cavities [32]. Our results provide further
insight into the puzzling propagation dynamics and radiative processes in
presence of two-frequency soliton molecules. Due to the manifold of two-
frequency pulse compounds with differing substructure, their emission spectra
manifest in various complex forms. This demonstrates a peculiar nonlinear
optics system for investigating emission of resonant radiation.

OM, SW, IB, UM, and AD acknowledge financial support from Deutsche
Forschungsgemeinschaft (DFG) under Germany’s Excellence Strategy within the
Cluster of Excellence PhoenixD (Photonics, Optics, and Engineering—Innovation
Across Disciplines) (EXC 2122, Project No. 390833453). AY acknowledges
financial support from Priority 2030 Academic Leadership Program and
goszadanie no. 2019-1246. IB also acknowledges support from DFG (Project No.
BA4156/4-2). UM also acknowledges support from DFG (Project No. MO 850-20/1).

##  Appendix A Simplified description of soliton molecule amplitude and width
oscillations

Considering the standard NSE in nondimensional form

$i\partial_{z}A=-\frac{1}{2}\partial_{t}^{2}A-|A^{2}|A,$ (20)

a variational approach based on the trial function

$A(z,t)=\eta\,{\rm{sech}}\left(\frac{t}{w}\right)\exp\left(i\theta+ib\frac{t^{2}}{2w}\right),$
(21)

yields the set of coupled ordinary differential equations

$\displaystyle\frac{d}{dz}(\eta^{2}w)=0,$ (22) $\displaystyle\frac{d}{dz}w=b,$
(23) $\displaystyle\frac{d}{dz}b=\frac{4}{\pi^{2}w^{3}}(1-\eta^{2}w^{2}),$
(24) $\displaystyle\frac{d}{dz}\theta=\frac{5\eta^{2}}{6}-\frac{1}{3w^{3}},$
(25)

describing the $z$-variation of the amplitude ($\eta$), pulse width ($w$),
pulse phase ($\theta$), and, chirp ($b$) of the trial function (21) [46].
Equation (22) implies the constant of motion $\eta^{2}w=\rm{const.}$, linking
intensity and width variations of the trial function (21). For an approximate
analysis of the pulse width variation along $z$, we consider the
$z$-derivative of Eq. (23), which, together with Eq. (24), defines the
nonlinear equation of motion

$\frac{d^{2}}{dz^{2}}w=\frac{4}{\pi^{2}w^{3}}(1-\eta^{2}w^{2})\equiv f(w).$
(26)

Assuming a pulse width $w(z)=w_{0}+\epsilon(z)$ that performs periodic
oscillations with amplitude $\epsilon(z)$ about some preferred value $w_{0}$,
we may expand Eq. (26) up to first order in $\epsilon$ to find
$\frac{d^{2}}{dz^{2}}\epsilon\approx f(w_{0})+f^{\prime}(w_{0})\,\epsilon$
with $f(w_{0})=\alpha w_{0}^{-3}+\beta w_{0}^{-1}$, and
$f^{\prime}(w_{0})\equiv\frac{\partial}{\partial w}f(w)|_{w=w_{0}}=-3\alpha
w_{0}^{-4}-\beta w_{0}^{-2}$, with $\alpha=4\pi^{-2}$ and
$\beta=-\alpha\eta^{2}$. Imposing the condition $f(w_{0})=0$ fixes
$w_{0}=\eta^{-1}$, i.e. a reasonable choice for $w_{0}$ given the trial
function (21), and brings Eq. (26) to the form of a harmonic oscillator

$\displaystyle\frac{d^{2}}{dz^{2}}\epsilon(z)+K^{2}(\eta)\,\epsilon(z)\approx
0,$ (27)

with angular wavenumber $K(\eta)=2\sqrt{2}\eta^{2}/\pi$. The wavelength of
this harmonic oscillation then reads

$\Lambda^{\rm{HO}}(\eta)=\frac{2\pi}{K(\eta)}=\frac{\pi^{2}}{\sqrt{2}\eta^{2}}.$
(28)

We can readily transfer this result to the case of amplitude and width
oscillations of two-color soliton molecules. In the limit where the spectral
extent of a molecules subpulses is small in comparison to the frequency gap
separating them, i.e. under the condition $\pi
t_{0}\sqrt{6\beta_{2}/|\beta_{4}|}\gg 1$, each subpulse can effectively be
modeled in terms of a standard NSE [39]. Then we are led to expect that the
oscillation period of each soliton molecule subpulse, and thus the entire
soliton molecule, is reasonably well described by the approximate result Eq.
(28). We plot $\Lambda^{\rm{HO}}(\eta)/z_{0}$, with $z_{0}=1$ for Eq. (20), as
dashed line, labeled D, in Fig. 3. Let us note that expanding Eq. (26) up to
second order in $\epsilon$ yields a correction term
$\propto[f^{\prime\prime}(w_{0})/2]\epsilon^{2}$ with
$f^{\prime\prime}(w_{0})\equiv\frac{\partial^{2}}{\partial
w^{2}}f(w)|_{w=w_{0}}=40\eta^{5}/\pi^{2}$, indicating that the anharmonicity
of the oscillation increases rapidly with increasing amplitude. This finding
is in qualitative agreement with the numerical results in Sec. 3.1.

Figure 7: Spectrogram ($\sigma=70~{}\mathrm{fs}$) of a vibrating soliton
molecule with slightly different subpulse amplitudes. The main plot shows the
propagation dynamics starting from an initial condition given by Eq. (29) with
$N_{1}=1.59$ and $N_{2}=1.66$. Parts of the main plot enclosed by the box is
detailed in the close-up view labeled A. Due to the unequal subpulse
amplitudes, the degeneracy among the resonances conditions (14,15) is lifted.
For comparison, the close-up view labeled B shows the corresponding plotting
range for the degenerate case with $N_{1}=N_{2}=1.59$.

##  Appendix B Splitting of resonance lines in presence of small amplitude
differences

To complement the strong symmetry breaking perturbations induced by a non-
symmetric propagation constant, discussed in Sec. 3.3, we here demonstrate
degeneracy-lifted multifrequency radiation resulting from a small difference
in the molecules subpulse amplitudes. Specifically, we consider the symmetric
propagation constant
$\beta(\Omega)=\frac{\beta_{2}}{2}\Omega^{2}+\frac{\beta_{4}}{24}\Omega^{4}$
for $\beta_{2}=1\,\mathrm{fs^{2}/\mu m}$ and $\beta_{4}=-1\,\mathrm{fs^{4}/\mu
m}$, and the initial condition

$\displaystyle
A_{0}(t)=a_{1}\,{\rm{sech}}(t/t_{1})\,e^{-i\Omega_{1}t}+a_{2}\,{\rm{sech}}(t/t_{2})\,e^{-i\Omega_{2}t},$
(29)

where $\Omega_{1,2}=\mp\sqrt{6}\,\mathrm{rad/fs}$, $t_{0}=6\,\mathrm{fs}$, and
$a_{1,2}^{2}=N_{1,2}^{2}2\beta_{2}/(\gamma\,t_{0}^{2})$, with $N_{1}=1.57$ and
$N_{2}=1.66$. Correspondingly, $a_{2}\gtrapprox a_{1}$ in Eqs. (16,17). A
spectrogram of the resulting vibrating soliton molecule at $z/z_{0}\approx 80$
[$z_{0}=t_{0}^{2}/(2\beta_{2})$] is shown in Fig. 7. As evident from the
spectrogram range
$[\Omega=-2.3\ldots-1.5~{}\mathrm{rad/fs}]\times[t=-900\ldots-300~{}\mathrm{fs}]$,
highlighted in close-up view A, the degeneracy of the resonances conditions
(14,15) is lifted. The primary resonances, resulting from solutions to Eq.
(14), are now accompanied by pairs of secondary resonances, resulting from
solutions to Eq. (15) for $n=1,2$, which appear on either side of a primary
resonance. To facilitate comparison to the degenerate scenario, close-up view
B shows the corresponding numerical results for $N_{1}=N_{2}=1.59$.

## References

## References

  * [1] Krainov V and Smirnov B M 2019 Atomic and Molecular Radiative Processes (Springer)
  * [2] Mitschke F 2016 Fiber Optics: Physics and Technology (Springer)
  * [3] Agrawal G P 2019 Nonlinear Fiber Optics (Academic Press)
  * [4] Drazin P G and Johnson R S 1989 Solitons: An Introduction (Cambridge University Press)
  * [5] Yang T S and Kath W L 2001 Physica D 149 80–94
  * [6] Pelinovsky D E and Yang J 2004 SIAM J. Appl. Math. 64 1360–1382
  * [7] Conforti M, Trillo S, Mussot A and Kudlinski A 2015 Sci. Rep. 5 9433
  * [8] Driben R, Yulin A V and Efimov A 2015 Opt. Express 23 19112–19117
  * [9] Melchert O, Demircan A and Yulin A 2020 Sci. Rep. 10 8849
  * [10] Stratmann M, Pagel T and Mitschke F 2005 Phys. Rev. Lett. 95 143902
  * [11] Hause A, Hartwig H, Böhm M and Mitschke F 2008 Phys. Rev. A 78 063817
  * [12] Ueda T and Kath W L 1990 Phys. Rev. A 42 563–571
  * [13] Menyuk C R 1988 J. Opt. Soc. Am. B 5 392–402
  * [14] Afanasjev V V, Dianov E M, Prokhorov A M and Serkin V N 1988 Pisma Zh. Eksp. Teor. Fiz. 48 588
  * [15] Trillo S, Wabnitz S, Wright E M and Stegeman G I 1988 Opt. Lett. 13 871–873
  * [16] Afanasjev V, Dianov E and Serkin V 1989 IEEE J. Quantum Electron. 25 2656–2664
  * [17] Afanasyev V V, Kivshar Y S, Konotop V V and Serkin V N 1989 Opt. Lett. 14 805–807
  * [18] Mesentsev V K and Turitsyn S K 1992 Opt. Lett. 17 1497–1499
  * [19] Akhmediev N and Ankiewicz A 2000 Chaos 10 600–612
  * [20] Oreshnikov I, Driben R and Yulin A 2017 Phys. Rev. A 96 013809
  * [21] Krupa K, Nithyanandan K, Andral U, Tchofo-Dinda P and Grelu P 2017 Phys. Rev. Lett. 118 243901
  * [22] Wang Z Q, Nithyanandan K, Coillet A, Tchofo-Dinda P and Grelu P 2019 Nature Communications 10 830
  * [23] Melchert O, Willms S, Bose S, Yulin A, Roth B, Mitschke F, Morgner U, Babushkin I and Demircan A 2019 Phys. Rev. Lett. 123 243905
  * [24] Tam K K K, Alexander T J, Blanco-Redondo A and de Sterke C M 2020 Phys. Rev. A 101 043822
  * [25] Melchert O, Willms S, Morgner U, Babushkin I and Demircan A 2021 Scientific Reports 11 11190
  * [26] Philbin T G, Kuklewicz C, Robertson S, Hill S, König F and Leonhardt U 2008 Science 319 1367–1370
  * [27] Demircan A, Amiranashvili S and Steinmeyer G 2011 Phys. Rev. Lett. 106(16) 163901
  * [28] Faccio D 2012 Cont. Phys. 53 97–112
  * [29] Demircan A, Amiranashvili S, Brée C and Steinmeyer G 2013 Phys. Rev. Lett. 110 233901
  * [30] Plansinis B W, Donaldson W R and Agrawal G P 2015 Phys. Rev. Lett. 115(18) 183901
  * [31] Lourdesamy J P, Runge A F J, Alexander T J, Hudson D D, Blanco-Redondo A and de Sterke C M 2021 Nat. Phys. 18 59–66
  * [32] Mao D, Wang H, Zhang H, Zeng C, Du Y, He Z, Sun Z and Zhao J 2021 Nat. Commun. 12 6712
  * [33] Melchert O, Yulin A and Demircan A 2020 Opt. Lett. 45 2764
  * [34] Moille G, Li Q, Kim S, Westly D and Srinivasan K 2018 Opt. Lett. 43 2772–2775
  * [35] Willms S, Melchert O, Bose S, Yulin A, Oreshnikov I, Morgner U, Babushkin I and Demircan A 2022 Phys. Rev. A 105 053525
  * [36] Yulin A V, Skryabin D V and Russell P S J 2004 Opt. Lett. 29 2411–2413
  * [37] Skryabin D V and Yulin A V 2005 Phys. Rev. E 72 016619
  * [38] Oreshnikov I, Melchert O, Willms S, Bose S, Babushkin I, Demircan A, Morgner U and Yulin A 2022 Cherenkov radiation and scattering of external dispersive waves by two-color solitons arXiv:2207.03541 URL {https://doi.org/10.48550/arXiv.2207.03541}
  * [39] Melchert O and Demircan A 2021 Opt. Lett. 46 5603
  * [40] Heidt A M 2009 J. Lightwave Tech. 27 3984
  * [41] Melchert O and Demircan A 2022 Computer Physics Communications 273 108257
  * [42] Hult J 2007 J. Lightwave Tech. 25 3770
  * [43] Melchert O, Roth B, Morgner U and Demircan A 2019 SoftwareX 10 100275
  * [44] Cohen L 1989 Procedings of the IEEE 77 941–981
  * [45] Anderson D 1983 Phys. Rev. A 27 3135
  * [46] Kath W L and Smyth N F 1995 Phys. Rev. E 51 1484
  * [47] Akhmediev N and Karlsson M 1995 Phys. Rev. A 51 2602–2607
  * [48] Cristiani I, Tediosi R, Tartara L and Degiorgio V 2004 Opt. Express 12 124–135