# A generalized nonlinear Schrödinger Python module implementing different
models of input pulse quantum noise

O. Melchert<EMAIL_ADDRESS>A. Demircan<EMAIL_ADDRESS>hannover.de Leibniz Universität Hannover, Institute of Quantum Optics (IQO),
30167 Hannover, Germany  Cluster of Excellence PhoenixD (Photonics, Optics,
and Engineering - Innovation Across Disciplines), Hannover, Germany

###### Abstract

We provide Python tools enabling numerical simulation and analysis of the
propagation dynamics of ultrashort laser pulses in nonlinear waveguides. The
modeling approach is based on the widely used generalized nonlinear
Schrödinger equation for the pulse envelope. The presented software implements
the effects of linear dispersion, pulse self-steepening, and the Raman effect.
The focus lies on the implementation of input pulse shot noise, i.e. classical
background fields that mimick quantum noise, which are often not thoroughly
presented in the scientific literature. We discuss and implement commonly
adopted quantum noise models based on pure spectral phase noise, as well as
Gaussian noise. Coherence properties of the resulting spectra can be
calculated. We demonstrate the functionality of the software by reproducing
results for a supercontinuum generation process in a photonic crystal fiber,
documented in the scientific literature. The presented Python tools are are
open-source and released under the MIT license in a publicly available
software repository.

###### keywords:

Generalized nonlinear Schrödinger equation , quantum noise , spectral
coherence , Python

###### Contents

  1. 1 Introduction
  2. 2 Scientific problem solved by the software
  3. 3 Software description
    1. 3.1 Software Functionalities
    2. 3.2 Implementation details
      1. 3.2.1 The GNLS data structure
      2. 3.2.2 Implemented noise models
      3. 3.2.3 Functions for calculating coherence properties
  4. 4 Sample results
    1. 4.1 A minimal working example
    2. 4.2 Using GNLS with py-fmas
  5. 5 Impact and conclusions
  6. Acknowledgements

##  1 Introduction

The propagation of laser pulses in nonlinear waveguides supports the
generation of supercontinuum spectra [1, 2, 3]. Starting from a spectrally
narrow input pulse, the interplay of linear and nonlinear effects induces
tremendous spectral broadening, yielding flat spectra that can extend from the
violet to the infrared [4]. Such effects can be achieved, e.g., in photonic
crystal fibers (PCFs) [5, 6], wherein supercontinuum spectra can be produced
using $\sim\\!100\,\mathrm{fs}$-duration pulses, peak powers
$\sim\\!10\,\mathrm{kW}$ and propagation lengths on the order of
$1\,\mathrm{m}$ [4]. The resulting broad, flat spectra with high spectral
density find application, e.g., in optical frequency metrology [7], and
optical technologies [2].

A flexible theoretical framework for studying the complex physical processes
associated with the generation of supercontinuum spectra is provided by the
generalized nonlinear Schrödinger equation (GNLS) [1]. In order to model the
propagation dynamics of laser pulses it combines the effects of linear
dispersion, pulse self-steepening [8, 9], and the Raman effect [10]. This
accounts for various processes that support the generation of widely extended
supercontinuum spectra, such as the modulation istability [11], soliton-
fission [12, 13], and self-frequency shift of Raman solitons [10]. It furthers
forms the basis for modeling optical rouge waves [14, 15, 16, 17], and the
interaction of solitons and dispersive waves in presence of an optical event-
horizon [18, 19, 20, 21]. The GNLS describes the real-valued optical field in
terms of a complex-valued envelope at a fixed reference frequency using a
nonlinear wave equation of first order in the propagation coordinate. It has
proven to yield reliable results even in the single cycle regime [22].
Overall, the GNLS provides a classical description of the propagation dynamics
of input pulses. Basic phenomena described by the GNLS are susceptible to the
presence of noise. When attempting to model noise on the quantum level, a
fundamental issue arises since the rigorous description of quantum noise
sources requires the introduction of quantum mechanical operators [23, 24, 25,
26], which are incompatible with classical fields. When the number of photons
in the input pulse is large, however, semiclassical approaches exist that
account for quantum fluctuations by randomizing physical interactions using a
classical background field of low power, entering the GNLS via the initial
condition. Such input pulse noise is modeled by simply adding to the coherent
input pulse an incoherent, stochastic noise field. We consider commonly
adopted one photon per frequency mode [3, 27, 28, 29, 30], and half a photon
per temporal mode [31, 32, 33] noise models, and further implement a
background field based on a classical analog of the zero-point field of
quantum field theory (QFT) [34]. Including noise by perturbing the initial
conditions is in accord with the Wigner method, an efficient computational
approach for modeling quantum noise in optical fibers [26, 25, 35].
Considering an ensemble of independent simulation runs perturbed by noise,
coherence properties of the simulated spectra can be investigated. Accounting
for shot-to-shot fluctuations in numerical simulations appears essential for
the dynamics and enables a better comparison to actual experiments, employing
multishot measurement techniques [27, 28].

In this article we present software tools allowing to reliably simulate and
analyze basic phenomena described by the GNLS in presence of various types of
input pulse shot noise. We detail the scientific problem addressed by the
provided software in Sec. 2, and discuss its implementation, with emphasis on
input pulse noise models, in Sec. 3. Section 4 reports verification tests,
reproducing results for a supercontinuum generation process in a photonic
crystal fiber, well documented in the scientific literature, where it has been
used for illustrating nonlinear-optics effects [3], and for benchmarking
algorithms [36, 37, 38, 39, 40]. We comment on impact and conclude in Sec. 5.

##  2 Scientific problem solved by the software

The provided Python tools enable simulation and analysis of the dynamics of
ultrashort laser pulses, governed by the generalized nonlinear Schrödinger
equation (GNLS) [1, 3]

$\displaystyle\partial_{z}u=i\sum_{n\geq 2}$
$\displaystyle\frac{\beta_{n}}{n!}(i\partial_{t})^{n}u+i\gamma\left(1+\frac{i\partial_{t}}{\omega_{0}}\right)$
$\displaystyle\times\left[\,u(z,t)\int
R(t^{\prime})|u(z,t-t^{\prime})|^{2}~{}{\rm{d}}t^{\prime}\right],$ (1)

for a complex-valued field envelope $u\equiv u(z,t)$ on a periodic time-domain
of extend $T$ with boundary condition $u(z,-T/2)=u(z,T/2)$ and propagation
distance $z$. An initial condition $u(z=0,t)=u_{0}(t)$ needs to be specified.
In Eq. (1), $\beta_{n}$ (in units of $\mathrm{ps^{n}/\mu m}$) is the $n$-th
order dispersion coefficient, $\gamma$ ($\mathrm{W^{-1}/\mu m}$) is a scalar
nonlinear coefficient, and $\omega_{0}$ ($\mathrm{rad/ps}$) a reference
frequency. $t$ is a retarded time, measured in a reference frame moving with
the group velocity at $\omega_{0}$. The Raman effect is included via the total
response function

$\displaystyle R(t)=(1-f_{\rm{R}})\,\delta(t)+f_{\rm{R}}\,h_{\rm{R}}(t),$ (2)

where the first term results in an instantaneous Kerr-type nonlinear response,
and where the Raman response function $h_{\rm{R}}(t)$ enters with fractional
contribution $f_{\rm{R}}$ [41, 42]. A generic two-parameter model is

$\displaystyle
h_{\rm{R}}(t)=\frac{\tau_{1}^{2}+\tau_{2}^{2}}{\tau_{1}\tau_{2}^{2}}\,e^{-t/\tau_{2}}\,\sin(t/\tau_{1})\,\Theta(t),$
(3)

with the Heaviside step function $\Theta(t)$ ensuring causality [41]. Based on
a single damped harmonic-oscillator approximation of molecular responses, this
simple model reproduces the Raman gain spectrum measured for silica glass
reasonably well. For silica fibers, numerical values for the three relevant
model parameters are $f_{\rm{R}}=0.18$, $\tau_{1}=12.2~{}\mathrm{fs}$, and
$\tau_{2}=32~{}\mathrm{fs}$ [41]. More complex response functions, based on
improved response models and valid for other media, exist [43, 44, 45]. In
this framework, using a discrete sequence of angular frequency detunings
$\Omega\in\frac{2\pi}{T}\mathbb{Z}$ relative to $\omega_{0}$, the expressions

$\displaystyle\mathsf{F}[u(z,t)]\equiv\frac{1}{T}\int_{-T/2}^{T/2}u(z,t)\,e^{i\Omega
t}~{}{\rm{d}}t=u_{\Omega}(z),$ (4a)
$\displaystyle\mathsf{F}^{-1}[u_{\Omega}(z)]\equiv\sum_{\Omega}u_{\Omega}(z)\,e^{-i\Omega
t}=u(z,t),$ (4b)

specify forward [Eq. (4a)], and inverse [Eq. (4b)] Fourier transforms,
relating the field envelope $u(z,t)$ to the spectral envelope $u_{\Omega}(z)$.
Using Parseval’s identity for the transforms (4), the energy in both domains
is

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

with instantaneous power $|u(z,t)|^{2}$ ($\mathrm{W=J/s}$) and power spectrum
$|u_{\Omega}(z)|^{2}$ ($\mathrm{W}$). While Eq. (1) represents the GNLS in its
time-domain formulation, a consistent formulation in the frequency domain is
also possible [46].

Input pulse noise is modeled by simply adding to the coherent input pulse
$u_{0}(t)$ an incoherent, stochastic noise field $\Delta u(t)$ with properties

$\displaystyle\langle u(t)\rangle=0,$ (6a) $\displaystyle\langle\Delta
u(t)\Delta u^{*}(0)\rangle=\sigma^{2}\delta(t),$ (6b)

where $\langle\ldots\rangle$ denotes an ensemble average over independent
instances of noise, and $\sigma^{2}$ is the noise variance. According to Eq.
(6a), at each point $t$, this background field has zero mean. The
$\delta$-function in the autocorrelation (6b) indicates that the noise varies
fast in comparison to any reasonable field $u_{0}(t)$, and that subsequent
actions of the noise are uncorrelated. For instance, considering _half a
photon_ with energy $\hbar\omega_{0}$ per temporal mode [31, 32, 33], the
noise variance is $\sigma^{2}=\hbar\omega_{0}/2$ with reduced Planck constant
$\hbar$. Noise models for treating more specific pump laser spectra exist [47,
48]. Let us note that there also exist formulations of the propagation
dynamics in terms of stochastic nonlinear Schrödinger equations, which are
equivalent to quantum field operator equations and account for vaccuum
fluctuations during propagation [49, 50, 51].

In the frequency domain, shot-to-shot fluctuations of the field can then be
characterized by the spectrally resolved modulus of first order coherence for
zero time-lag

$\displaystyle|g_{12}(\Omega)|=\left|\frac{\langle
u_{\Omega,m}^{\phantom{*}}u_{\Omega,k}^{*}\rangle_{m\neq
k}}{\sqrt{\langle|u_{\Omega,m}|^{2}\rangle\langle|u_{\Omega,k}|^{2}\rangle}}\right|,$
(7)

where the angular brackets specify an average over nonidentical pairs of
fields, labeled $m$ and $k$, obtained from an ensemble of $M$ simulation runs
with independent noise fields [27, 28, 3]. At a given angular frequency
detuning $\Omega$, $0\leq|g_{12}(\Omega)|\leq 1$, where
$|g_{12}(\Omega)|\approx 1$ indicates good stability in amplitude and phase.
Different from this measure of interpulse coherence, the modified intrapulse
coherence

$\displaystyle\tilde{\Gamma}(\Omega_{1},\Omega_{2})=\frac{|\langle
u^{\phantom{*}}_{\Omega_{1},m}u^{*}_{\Omega_{2},m}\rangle|}{\langle|u^{\phantom{*}}_{\Omega_{1},m}u^{*}_{\Omega_{2},m}|\rangle}$
(8)

allows to assess the coherence between different spectral components within a
pulse. With focus on $f$-to-$2f$ ($f=$ frequency) setups, the intrapulse
coherence was shown to play an important role in carrier-envelope phase
measurement and stabilization of ultrashort pulses [52]. Equation (8) modifies
the intrapulse coherence of Ref. [52] by relaxing the $f/2f$ condition.
Instead, two general distinct frequencies $\Omega_{1}$, and $\Omega_{2}$ are
considered.

##  3 Software description

GNLStools is written using the Python programming language [53], and depends
on the functionality of numpy, scipy [54], and matplotlib [55]. It can be
cloned directly from GitHub [56], where it is available under the MIT license.

###  3.1 Software Functionalities

The current version of GNLStools features:

  * 1.

A data structure for the GNLS (1), which can be easily tailored to a specific
nonlinear waveguide.

  * 2.

Functions that implement optional input pulse noise models considering both,
pure spectral phase noise as well as Gaussian noise.

  * 3.

Functions for calculating the coherence properties of simulated spectra.

The software can be used on its own, it includes a basic driver script
implementing a fixed stepsize “fourth-order Runge Kutta in the interaction
picture” (RK4IP) solver [36] (Sec. 4.1), or as extension module for py-fmas
[40] (Sec. 4.2). In the latter case, a user can take advantage of variable
stepsize $z$-propagtion algorithms and more elaborate Raman-response
functions. The presented software is not overly feature-rich, but allows to
reliably simulate the complex physical processes that enable generation of
supercontinuum spectra in nonlinear waveguides [3].

###  3.2 Implementation details

Solving the GNLS in the time-domain formulation Eq. (1) presents several
disadvantages. Higher-oder dispersion is difficult to handle. Even for the
lowest order $t$-derivatives, introduction of finite-difference operations
comes with unwanted truncation errors. It is therefore advantageous to perform
some of the operations in the frequency-domain, allowing to employ spectral
derivatives $\partial_{t}^{n}\,u=\mathsf{F}^{-1}[(-i\Omega)^{n}\,u_{\Omega}]$.
The operating principles of common algorithms for the solution of nonlinear
Schrödinger-type equations, like the split-step Fourier method [57, 58], the
RK4IP method [36], and the variable step-size conservation quantity error
method [37, 38], exploit this. To facilitate implementation, the interval
$-T/2\ldots T/2$ is divided into $N$ equal subintervals with temporal grid
spacing $\Delta t=N/T$, yielding discrete grid points $t_{m}=m\,\Delta t$, and
detuning grid points $\Omega_{m}=m\,\Delta\Omega$, $\Delta\Omega=2\pi/T$, with
$m=-N/2,\ldots,N/2-1$. Subsequently we write $u(z,t)|_{t=t_{m}}=u_{m}(z)$, and
$u_{\Omega}(z)|_{\Omega=\Omega_{m}}=u_{\Omega_{m}}(z)$. Respecting the sign
choice and normalization of the transforms (4), we use the numpy native
discrete Fourier transform (DFT) routine ifft to implement Eq. (4a), and fft
to implement Eq. (4b) [59].

####  3.2.1 The GNLS data structure

Instantiating an instance of the class GNLS requires a user to specify several
input parameters. Below, they are listed as “parameter_name (data type):
description”:

  * 1.

w (array, float): discrete angular frequencies $(\Omega_{m})_{m=-N/2}^{N/2-1}$
in units of $\mathrm{rad/fs}$;

  * 2.

beta_n (array, float): ordered sequence of dispersion parameters
$(\beta_{n})_{n=2}^{N_{\rm{max}}}=(\beta_{2},\beta_{3},\ldots,\beta_{N_{\rm{max}}})$,
with $N_{\rm{max}}\geq 2$ and $\beta_{n}$ in units of $\mathrm{fs^{n}/\mu m}$;

  * 3.

gamma (float): nonlinear parameter $\gamma$ $(\mathrm{W^{-1}/\mu m})$;

  * 4.

w0 (float): reference angular frequency $\omega_{0}$ $(\mathrm{rad/fs})$;

  * 5.

fR (float): Raman contribution $f_{\rm{R}}$ (default: $0.18$);

  * 6.

tau1 (float): time-scale $\tau_{1}$ $(\mathrm{fs})$ (default:
$12.2\,\mathrm{fs}$);

  * 7.

tau2 (float): time-scale $\tau_{2}$ $(\mathrm{fs})$ (default:
$32\,\mathrm{fs}$).

While the angular frequency grid w is passed as positional argument, all other
parameters are passed as keyword arguments. An example is discussed in Sect.
4.1. For using an instance of the GNLS class with a $z$-propagation algorithm,
several instance methods are available. Below, they are listed in the format
“method_name(arg1,arg2,...): description”:

Lw(self):

Using spectral derivatives, the frequency domain representation of the linear
operator $L\equiv i\sum_{n\geq 2}\frac{\beta_{n}}{n!}(i\partial_{t})^{n}$ on
the right-hand-side (rhs) of Eq. (1) can be written as

$\displaystyle L_{\Omega}\equiv
i\sum_{n=2}^{N_{\rm{max}}}\frac{\beta_{n}}{n!}\Omega^{n}.$ (9a)

For practical reasons, the sum has to be truncated at a finite integer number
$N_{\rm{max}}\geq 2$. This class method returns $L_{\Omega}$, evaluated at the
points $(\Omega_{m})_{m=-N/2}^{N/2-1}$.

Input parameters:

  * 1.

The method relies on class attributes only.

Output parameters:

  * 1.

Lw (array): $L_{\Omega}$ [Eq. (9a)] at the points
$(\Omega_{m})_{m=-N/2}^{N/2-1}$.

The method is decorated as `@property`, so it can be conveniently used as
`gnls_instance.Lw`.

Nw(self, uw):

Using the convolution theorem for the transforms (4), the second term on the
rhs of Eq. (1) can be written in a mixed representation as

$\displaystyle N_{\Omega}(u)\equiv
i\gamma(1+\omega_{0}^{-1}\Omega)\,\mathsf{F}\left[(1\\!-\\!f_{\rm{R}})\mathcal{I}_{1}+f_{\rm{R}}\mathcal{I}_{2}\right],$
(9b)

where $\mathcal{I}_{1}=|u|^{2}u$,
$\mathcal{I}_{2}=\mathsf{F}^{-1}[\tilde{h}_{\rm{R}}(\Omega)\mathsf{F}[|u|^{2}]]u$,
and

$\displaystyle\tilde{h}_{\rm{R}}(\Omega)\equiv
T\,\mathsf{F}[h_{\rm{R}}(t)]=\frac{\tau_{1}^{-2}+\tau_{2}^{-2}}{\tau_{1}^{-2}-(\Omega+i\tau_{2}^{-1})^{2}}.$
(9c)

The method takes on input the spectral envelope $u_{\Omega}$, retrieves
$u=\mathsf{F}^{-1}[u_{\Omega}]$, and evaluates and returns Eq. (9b) at
$(\Omega_{m})_{m=-N/2}^{N/2-1}$.

Input parameters:

  * 1.

uw (array): Spectral envelope $u_{\Omega_{m}}$ at the angular frequency grid
points $(\Omega_{m})_{m=-N/2}^{N/2-1}$.

Output parameters:

  * 1.

Nw (array): $N_{\Omega}$ [Eq. (9b)] at the points
$(\Omega_{m})_{m=-N/2}^{N/2-1}$.

claw_Ph(self,i,zi,w,uw):

Class method evaluating a conservation law of the GNLS (1), related to the
classical analog of the photon number [41]. This method considers the energy
(5) in the form $E=\hbar\sum_{\Omega}n_{\Omega}\,(\omega_{0}+\Omega)$, where
the dimensionless quantity $n_{\Omega}\equiv
T|u_{\Omega}|^{2}/[\hbar(\omega_{0}+\Omega)]$ specifies the number of photons
with energy $\hbar(\omega_{0}\\!+\\!\Omega)$. The total number of photons is
then given by

$\displaystyle
C_{\rm{Ph}}(z)\equiv\sum_{\Omega}n_{\Omega}=\frac{2\pi}{\hbar\Delta\Omega}\sum_{\Omega}\frac{|u_{\Omega}(z)|^{2}}{\omega_{0}+\Omega}.$
(9d)

In a rigorous quantum mechanical treatment, the photon number is instead
defined using photon creation and annihilation operators [24]. Let us note
that while the GNLS (1) conserves the total number of photons Eq. (9d), it
does not conserve the pulse energy Eq. (5) [41].

Input parameters:

  * 1.

i (int): Integer label of the current $z$-propagation step.

  * 2.

zi (float): Current propagation distance.

  * 3.

w (array): Angular frequency grid $(\Omega_{m})_{m=-N/2}^{N/2-1}$.

  * 4.

uw (array): Spectral envelope $u_{\Omega_{m}}$ at the angular frequency grid
points $(\Omega_{m})_{m=-N/2}^{N/2-1}$.

Output parameters:

  * 1.

Cph (float): Total number of photons $C_{\rm{Ph}}$ at the current point along
$z$ [Eq. 9d].

The method has the structure of an user-action function for use with py-fmas
[40] (Sec. 4.2).

####  3.2.2 Implemented noise models

We provide optional functions for generating time-domain representations of
input pulse noise, consistent with Eqs. (6). Below, they are listed in the
format “function_name(arg1,arg2,...): description”:

noise_model_01(t,w0,s0):

Function generating an instance of noise by directly sampling in the time-
domain. The underlying noise model assumes normally distributed amplitudes.

Input parameters:

  * 1.

t (array): Time-grid $(t_{m})_{m=-N/2}^{N/2-1}$.

  * 2.

w0 (float): Pulse center frequency $\omega_{0}$.

  * 3.

s0 (int): Integer seed $s_{0}$ for pseudo random number generator.

Output parameters:

  * 1.

du_t (array): Instance $(\Delta u_{m})_{m=-N/2}^{N/2-1}$ of time-domain noise.

An instance of this type of noise is obtained by directly sampling a sequence
$(\Delta u_{m})_{m=-N/2}^{N/2-1}$ of complex-valued noise-amplitudes in the
time domain as

$\displaystyle\Delta u_{m}$
$\displaystyle=\sqrt{\frac{\hbar\omega_{0}}{4\Delta t}}\,(X+iY),$ (10a)

where the independent random variables $X$ and $Y$ yield independent
identically distributed (iid) standard normal random numbers [i.e.
$X,Y\sim\mathcal{N}(0,1)$]. Performing an ensemble average, this type of noise
exhibits the properties

$\displaystyle\langle\Delta u_{m}\rangle=0,$ (10b) $\displaystyle\langle\Delta
u_{m}\Delta u_{0}^{*}\rangle=\frac{\hbar\omega_{0}}{2\Delta t}\,\delta_{m0},$
(10c)

i.e. the average field is zero and the autocorrelation (10c) is equivalent to
Eq. (6b) on a discrete grid with grid spacing $\Delta t$ [31]. An instance of
this noise has average energy $\langle\sum_{m}|u_{m}|^{2}\Delta
t\rangle=N\hbar\omega_{0}/2$, corresponding, on average, to _half a photon_
with energy $\hbar\omega_{0}$ per temporal mode.

noise_model_02(t,w0,s0):

Function generating an instance of time-domain noise by sampling its Fourier
representation. The underlying noise model assumes pure phase noise.

Input parameters:

  * 1.

t (array): Time-grid $(t_{m})_{m=-N/2}^{N/2-1}$.

  * 2.

w0 (float): Pulse center frequency $\omega_{0}$.

  * 3.

s0 (int): Integer seed $s_{0}$ for pseudo random number generator.

Output parameters:

  * 1.

du_t (array): Instance $(\Delta u_{m})_{m=-N/2}^{N/2-1}$ of time-domain noise.

An instance of this type of noise is obtained using a Fourier method by first
sampling a sequence $(\Delta u_{\Omega_{m}})_{m=-N/2}^{N/2-1}$ of random
complex-valued spectral amplitudes

$\displaystyle\Delta
u_{\Omega_{m}}\equiv\sqrt{\frac{\hbar(\omega_{0}\\!+\\!\Omega_{m})}{T}}e^{-i\Phi},$
(11a)

where the random variable $\Phi$ yields iid phase-angles uniformly distributed
in the range $0\ldots 2\pi$ [i.e. $\Phi\sim U(0,2\pi)$]; an inverse Fourier
transform, consistent with Eq. (4b), is then used to obtain the sequence
$(\Delta u_{m})_{m=-N/2}^{N/2-1}$ of time-domain noise amplitudes.

For this type of noise, the magnitude of the spectral amplitudes (11a) is
definite and the energy per mode $\Omega$ is
$T|u_{\Omega}|^{2}=\hbar(\omega_{0}\\!+\\!\Omega)$. With the number of photons
$n_{\Omega}$ as per Eq. (9d), adding noise as $u_{\Omega}\rightarrow
u_{\Omega}+\Delta u_{\Omega}$ results in $n_{\Omega}\rightarrow n_{\Omega}+1$.
The photon occupation number is thus increased by the minimal definite amount
of _one (entire) photon_ per mode. The phase of this photon is, however,
entirely indefinite and the noise in different modes has no particular phase
relationship. An instance of the noise has total energy
$T\sum_{\Omega}|u_{\Omega}|^{2}=N\hbar\omega_{0}$.

noise_model_03(t,w0,s0):

Function generating an instance of time-domain noise by sampling its Fourier
representation. The underlying noise model assumes normally distributed
spectral amplitudes.

Input parameters:

  * 1.

t (array): Time-grid $(t_{m})_{m=-N/2}^{N/2-1}$.

  * 2.

w0 (float): Pulse center frequency $\omega_{0}$.

  * 3.

s0 (int): Integer seed $s_{0}$ for pseudo random number generator.

Output parameters:

  * 1.

du_t (array): Instance $(\Delta u_{m})_{m=-N/2}^{N/2-1}$ of time-domain noise.

An instance of this type of noise is obtained using a Fourier method by first
sampling a sequence $(\Delta u_{\Omega_{m}})_{m=-N/2}^{N/2-1}$ of random
complex-valued spectral amplitudes

$\displaystyle\Delta u_{\Omega_{m}}\equiv$
$\displaystyle\sqrt{\frac{\hbar(\omega_{0}\\!+\\!\Omega_{m})}{T}}\sqrt{I}e^{-i\Phi},$
(12)

where the random variable $I$ obeys an exponential distribution with rate
parameter $2$ [i.e. $I\sim{\mathrm{Exp}}(2)$ with expectation value $\langle
I\rangle=1/2$], and $\Phi$ is uniformly distributed in $0\ldots 2\pi$ [i.e.
$\Phi\sim U(0,2\pi)$]; an inverse Fourier transform, consistent with Eq. (4b),
is then used to obtain the sequence $(\Delta u_{m})_{m=-N/2}^{N/2-1}$ of time-
domain noise amplitudes.

Averaging over the noise in a mode $\Omega$ yields energy $\langle T|\Delta
u_{\Omega}|^{2}\rangle=\hbar(\omega_{0}\\!+\\!\Omega)\langle
I\rangle=\hbar(\omega_{0}\\!+\\!\Omega)/2$, corresponding, on average, to
_half a photon_ per mode. In contrast to the pure phase noise model
implemented by noise_model_02, this type of noise exhibits instance to
instance fluctuations of the energy in each mode. Thereby, the random variable
$I$ takes the role of an occupation number with support $I\in[0,\infty)$. On
average, an instance of this noise has energy $\langle
T\sum_{\Omega}|u_{\Omega}|^{2}\rangle=N\hbar\omega_{0}/2$. Alternatively,
using the Box-Müller transform [60], we may rewrite $\sqrt{I}e^{-i\Phi}$, with
$I\sim{\rm{Exp}}(2)$, $\Phi\sim U(0,2\pi)$, in Eq. (12) as $(X+iY)/2$, with
$X,Y\sim\mathcal{N}(0,1)$ two iid standard normal random variables. This noise
model is based on a classical analog of the zero-point field of QFT, which
exactly reproduces the statistics of the electromagnetic vacuum [34].

Fourier methods that obtain time-domain noise by sampling frequency dependent
spectral noise amplitudes generally result in correlations $\langle\Delta
u_{m}\Delta u^{*}_{0}\rangle\sim|t_{m}|^{-1}$. This is demonstrated in Fig. 1,
where the ensemble averaged field and autocorrelation for the three noise
models are shown. As evident from the inset in Fig. 1(b), both Fourier methods
result in correlations that persist over many grid spacings $\Delta t$.

Let us note that the above noise models differ with respect to details of
implementation and interpretation. Nevertheless, with some degree of
approximation, they can easily be related to each other. For instance, under
the assumption that the required bandwidth of the computational domain is
small compared to the pulse center frequency, i.e.
$N\Delta\Omega\ll\omega_{0}$, we may approximate Eq. (11a) by $\Delta
u_{\Omega_{m}}\approx\sqrt{\hbar\omega_{0}/T}\exp(-\Phi)$. All terms in the
inverse Fourier transform [Eq. (4b)] of these noise spectral amplitudes are
then identically distributed. Approximating these sums using the central limit
theorem yields Eq. (10a) times a factor of $\sqrt{2}$ (this is a consequence
of one model assuming one entire photon per mode, while the other assumes
half-a photon per mode). Let us note that this approximation is accompanied by
difficulties: it corresponds to _one (entire) photon_ per mode if the energy
(5) is written as $E=\hbar\omega_{0}\sum_{\Omega}n_{\Omega}^{\prime}$, with
$n_{\Omega}^{\prime}\equiv T|u_{\Omega}|^{2}/(\hbar\omega_{0})$. Thus,
irrespective of $\Omega$, all photons are assumed to contribute the same
energy $\hbar\omega_{0}$. Put into the context of pulse propagation models,
this is a feature of standard nonlinear Schrödinger-type equations for which
energy conservation and photon number conservation are trivially linked. This,
however, is at odds with the GNLS (1) which conserves the number of photons
but not the energy [41]. Under the same assumption and approximation, Eq. (12)
can easily be related to Eq. (10a): applying the inverse Fourier transform, we
may rewrite the individual terms using the Box-Müller transform [60] to
directly obtain Eq. (10a). Moreover, in terms of a Fourier method and in the
above approximation, practically any iid spectral amplitudes will yield time-
domain noise with properties Eqs. (6) [61]. As a technical detail, let us note
that when using a Fourier method to set up noise, the noise spectral
amplitudes also depend on the normalization of the transform pair (4). For
instance, in Ref. [30], a GNLS on the infinite $t$-domain with a different
transform pair and, consequently, spectral noise amplitudes with different
normalization, was considered.

Figure 1: Ensemble averaged noise moments. (a) Real part of the time-domain
noise amplitudes as function of the grid point index $m=t_{m}/\Delta t$. (b)
Scaled autocorrelation. The inset shows a close up view of the index range
$-6\ldots 6$. Noise is sampled on a grid with temporal extent
$T=4~{}\mathrm{(ps)}$, $N=2^{13}$ grid points, and for
$\omega_{0}=2.2559~{}\mathrm{rad/fs}$. Averages are performed over $10^{4}$
independent instances of noise.

####  3.2.3 Functions for calculating coherence properties

We provide optional functions for calculating the coherence properties of
spectra obtained from pulse propagation simulations in terms of the GNLS (1).
Below, they are listed in the format “function_name(arg1,arg2,...):
description”:

coherence_interpulse(w,uw_list):

Function computing the interpulse coherence Eq. (7).

Input parameters:

  * 1.

w (array): Angular-frequency grid $(\Omega_{m})_{m=-N/2}^{N/2-1}$.

  * 2.

uw_list (array): List `[uw1, uw2, ...]`, comprising independent spectra `uw1,
uw2, ...`, obtained for the same propagation distance $z$ but for different
instances of input pulse noise.

Output parameters:

  * 1.

g12 (array): $g_{12}(\Omega)$ [Eq. (7)] at the points
$(\Omega_{m})_{m=-N/2}^{N/2-1}$.

coherence_intrapulse(w,uw_list,w1,w2):

Function computing the intrapulse coherence Eq. (8).

Input parameters:

  * 1.

w (array): Angular-frequency grid $(\Omega_{m})_{m=-N/2}^{N/2-1}$.

  * 2.

uw_list (array): List `[uw1, uw2, ...]`, comprising independent spectra `uw1,
uw2, ...`, obtained for the same propagation distance $z$ but for different
instances of input pulse noise.

  * 3.

w1 (float): Reference angular frequency $\Omega_{1}$.

  * 4.

w2 (float): Reference angular frequency $\Omega_{2}$.

Output parameters:

  * 1.

G (float): $\tilde{\Gamma}^{\rm{CEP}}(\Omega_{1},\Omega_{2})$ [Eq. 8].

Figure 2: Supercontinuum generation in a photonic crystal fiber. (a)
Instantaneous power, and, (b) spectrum after $14\,\mathrm{cm}$ of propagation.
In (b), result of the presented implementation (labeled A) is compared to
pyNLO [62] (labeled B). (c) Propagation of the instantaneous power, and, (d)
spectrum as function of propagation distance. White boxes in (c) show a close-
up view of the soliton fission process. (d) Normalized deviation of photon
number ($C_{\rm{Ph}}$) and energy ($E$) along the fiber.

##  4 Sample results

As a verification test of the presented Python tools, we consider a
supercontinuum generation process in a photonic crystal fiber, detailed in
Ref. [3]. Specifically, we perform simulations in terms of the GNLS (1), using
the sequence of dispersion coefficients $\beta_{2}=-1.183\times
10^{-2}\,\mathrm{fs^{2}/\mu m}$, $\beta_{3}=8.10383\times
10^{-2}\,\mathrm{fs^{3}/\mu m}$, $\beta_{4}=-9.5205\times
10^{-2}\,\mathrm{fs^{4}/\mu m}$, $\beta_{5}=0.20737\,\mathrm{fs^{5}/\mu m}$,
$\beta_{6}=-0.53943\,\mathrm{fs^{6}/\mu m}$,
$\beta_{7}=1.3486\,\mathrm{fs^{7}/\mu m}$,
$\beta_{8}=-2.5495\,\mathrm{fs^{8}/\mu m}$,
$\beta_{9}=3.0524\,\mathrm{fs^{9}/\mu m}$,
$\beta_{10}=-1.7140\,\mathrm{fs^{10}/\mu m}$, and $\gamma=0.11\times
10^{-6}\,\mathrm{W^{-1}/\mu m}$. For the Raman response we use the standard
values for silica fibers $f_{\rm{R}}=0.18$, $\tau_{1}=12.2\,\mathrm{fs}$, and
$\tau_{2}=32\,\mathrm{fs}$ [41]. As initial condition we take a hyperbolic-
secant pulse $u_{0}(t)=\sqrt{P_{0}}\,{\rm{sech}}(t/t_{0})$, with duration
$t_{0}=28.4\,\mathrm{fs}$, peak power $P_{0}=10\,\mathrm{kW}$, and pump
wavelength $\lambda_{0}=835\,\mathrm{nm}$ corresponding to $\omega_{0}=2\pi
c/\lambda_{0}\approx 2.2559\,\mathrm{rad/fs}$ (with speed of light
$c=0.29979\,\mathrm{\mu m/fs}$). The number of photons [Eq. (9d)] in this
pulse amounts to $C_{\rm{ph}}\approx 2.4\times 10^{9}$. The $z$-propagation
dynamics of the above propagation scenario without noise background, followed
over $14\,\mathrm{cm}$, is shown in Figs. 2(a-d).

###  4.1 A minimal working example

To demonstrate how to initialize and use the GNLS data structure for
performing pulse propagation in terms of Eq. (1), we provide a minimal example
in code-listing 1. The provided code initializes the GNLS data structure in
lines 10–28, sets up the input pulse in line 30, retrieves an instance of
input pulse noise in line 31, and propagates the initial field for
$10\,\mathrm{cm}$ using the fourth-order Runge Kutta in the interaction
picture (RK4IP) method [36]. This type of algorithmic approach is also
referred to as integrating factor method [63], or linearly exact Runge Kutta
method [64]. Pulse propagation is performed in lines 35–42 using a fixed
stepsize of $10\,\mathrm{\mu m}$, and using a temporal domain of extent
$T=7\,\mathrm{ps}$ with $N=2^{13}$ grid points. The script terminates in just
under $15\,\mathrm{sec.}$ (Apple M1 chip @ 3.2 GHz) and reproduces Figs.
2(a,b) [cf. Figs. 18(a,b) of Ref. [3], Fig. 2(a) of Ref. [36], Fig. 1(a) of
Ref. [38], and Fig. 2(b,c) of Ref. [40]]. The script shown in listing 1 is
located in project folder `numExp01` provided with the code [56].

###  4.2 Using GNLS with py-fmas

##### Integration with py-fmas

Figure 2 is produced by using the GNLS data structure in conjunction with py-
fmas [40], a Python package for the numerical simulation of the
$z$-propagation dynamics of ultrashort optical pulses in terms of the analytic
signal of the optical field. GNLStools can be used as an elaboration module,
providing an envelope based model which integrates well with the propagation
algorithms provided by py-fmas. Specifically, pulse propagation in Figs. 2-4
is performed by the variable stepsize conservation quantity error method [37,
38, 40], using Eq. (9d) to guide step-size selection. Figure 2(e) shows the
numerical error accumlated for the photon number [Eq. (9d)] and energy [Eq.
(5)], showing that $C_{\rm{Ph}}$ is conserved up to order $10^{-8}$, while $E$
is not conserved [41].

Figure 3: Shot-to-shot fluctuations of (a) instantaneous power, and, (b)
corresponding spectra, for input pulses with $t_{0}=85.0\,\mathrm{fs}$
duration input pulses and $10\,\mathrm{cm}$ of propagation (cf. Fig. 18 of
Ref. [3]).

Considering $C_{\rm{Ph}}$, the feature at $z\approx 1~{}\mathrm{cm}$ indicates
the onset of soliton-fission (close-up view in Fig. 2c), and the change of the
trend at $z\approx 7~{}\mathrm{cm}$ indicates the onset of the interaction of
a soliton and a dispersive wave which mainly determines the interaction
dynamics and generation of supercontinuum spectra. A small project workflow
with a driver script that performs the simulation, and a postprocessing script
that generates Fig. 2, are located in project folder `numExp02` provided with
the code [56].

Listing 1: Python code for solving the GNLS (1) using the fourth-order Runge
Kutta in the interaction picture method [36].

⬇

1import numpy as np

2import matplotlib.pyplot as plt

3from GNLStools import GNLS, noise_model_01

4

5# \-- SET COMPUTATIONAL GRID

6z, dz = np.linspace(0, 0.1e6, 10000, retstep=True)

7t = np.linspace(-3500, 3500, 2**13, endpoint=False)

8w = np.fft.fftfreq(t.size, d=t[1]-t[0])*2*np.pi

9# \-- INSTANTIATE GENERALIZED NONLINEAR SCHROEDINGER EQUATION

10gnls = GNLS(

11 w, # (rad/fs)

12 beta_n = [

13 -1.1830e-2, # (fs^2/micron) beta_2

14 8.1038e-2, # (fs^3/micron) beta_3

15 -0.95205e-1, # (fs^4/micron) beta_4

16 2.0737e-1, # (fs^5/micron) beta_5

17 -5.3943e-1, # (fs^6/micron) beta_6

18 1.3486, # (fs^7/micron) beta_7

19 -2.5495, # (fs^8/micron) beta_8

20 3.0524, # (fs^9/micron) beta_9

21 -1.7140, # (fs^10/micron) beta_10

22 ],

23 gamma=0.11e-6, # (1/W/micron)

24 w0= 2.2559, # (rad/fs)

25 fR = 0.18, # (-)

26 tau1 = 12.2, # (fs)

27 tau2 = 32.0 # (fs)

28 )

29# \-- SPECIFY INITIAL PULSE

30ut = np.sqrt(1e4)/np.cosh(t/28.4)

31dut = noise_model_01(t, 2.2559, 1)

32uw = np.fft.ifft(ut + dut)

33

34# \-- RK4IP PULSE PROPAGATION

35P = np.exp(gnls.Lw*dz/2)

36for n in range(1,z.size):

37 uw_I = P*uw

38 k1 = P*gnls.Nw(uw)*dz

39 k2 = gnls.Nw(uw_I + k1/2)*dz

40 k3 = gnls.Nw(uw_I + k2/2)*dz

41 k4 = gnls.Nw(P*uw_I + k3)*dz

42 uw = P*(uw_I + k1/6 + k2/3 + k3/3) + k4/6

43

44# \-- PLOT RESULTS

45fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(8, 3))

46I = np.abs(np.fft.fft(uw))**2

47ax1.plot(t, I*1e-3)

48ax1.set_xlim(-200,3200); ax1.set_xlabel(r"Time $t$ (fs)")

49ax1.set_ylim(0,6); ax1.set_ylabel(r"Intensity $|u|^2$ (kW)")

50Iw = np.abs(uw)**2

51ax2.plot(2*np.pi*0.29979/(w+2.2559), 10*np.log10(Iw/np.max(Iw)))

52ax2.set_xlim(0.45,1.4); ax2.set_xlabel(r"Wavelength $\lambda$ (micron)")

53ax2.set_ylim(-60,0); ax2.set_ylabel(r"Spectrum $|u_\lambda|^2$ (dB)")

54fig.tight_layout(); plt.show()

##### Shot-to-shot fluctuations

Figure 3 demonstrates shot-to-shot variations in pulse intensity and spectrum
after $10\,\mathrm{cm}$ of propagation, arising from the inclusion of input
pulse noise for an input pulse with duration $t_{0}=85.1\,\mathrm{fs}$. The
noise is generated via noise model 1, i.e. by direct sampling in the time-
domain. A small project workflow with a driver script that performs the
simulation, and a postprocessing script that generates Fig. 3, are located in
project folder `numExp03_noise_model_01`, provided with the code [56].

##### Coherence properties

The coherence properties of input pulses of different duration, obtained by
performing ensemble averages over $200$ independent simulation runs with
different noise seeds, are shown in Fig. 4. Figure 4 reproduces Fig. 19 of
Ref. [3], wherein a detailed discussion of the effects of input pulse noise
and the interpretation of the coherence can be found. In addition, in Fig. 4
we also show the modified intrapulse coherence [Eq. 8] discussed in Ref. [52].
A postprocessing script that generates Fig. 4 is located in folder
`numExp03_noise_model_01` provided with the code [TBW]. From our experience,
the slight correlations introduced by the Fourier method based noise models,
see Sect. 3.2.2 and Fig. 1, do not affect the coherence properties of the
obtained supercontinuum spectra.

##  5 Impact and conclusions

The presented GNLStools comprise a data structure implementing the generalized
nonlinear Schrödinger equation, two commonly adopted models of quantum noise
and a further noise model based on a classical analog of the zero-point field
of QFT, and functions for assessing the coherence properties of simulated
spectra. It provides all features required for studying basic phenomena
supported by the GNLS in the presence of input pulse shot noise. The provided
software can be used on its own, as demonstrated in Sec. 4.1, or as
elaboration module for use with the propagation algorithms provided by the py-
fmas package [40], as demonstrated in Sec. 4.2.

With the provided software we hope to shed some more light on the
implementation of quantum noise models for use in pulse propagation studies,
an issue that is often not thoroughly presented in the scientific literature.
We hope to directly benefit students and researchers alike, which are new to
the field of nonlinear optics and seek a tutorial-type introduction on how to
perform pulse propagation simulations including quantum noise, and assess the
coherence properties of the resulting spectra. The minimal example in Sec. 4.1
can very well serve as classroom code or as starting point for seminar
projects in computation oriented courses, aimed at demonstrating algorithmic
approaches that go beyond simple split-step Fourier methods commonly used for
solving nonlinear Schrödinger-type equations [1, 57, 58].

Finally, we would like to reference original research applying the presented
GNLS tools. In Ref. [21] we used the presented software to demonstrate an
efficient all-optical switching scheme, based on controlling the features of
soliton fission induced supercontinuum spectra using a time-delayed dispersive
wave.

Figure 4: Coherence properties of supercontinuum spectra after
$10\,\mathrm{cm}$ of propagation (cf. Fig. 19 of Ref. [3]). (a) from top to
bottom: exemplary single-shot output intensity, ensemble averaged spectrum,
and degrees of coherence $|g_{12}|$ [Eq. (7)], and $\tilde{\Gamma}$ [Eq. (8)],
as function of wavelength $\lambda=2\pi c/(\omega_{0}+\Omega)$. For the
calculation of the modified intrapulse coherence
$\tilde{\Gamma}(\Omega_{1},\Omega_{2})$, $\Omega_{1}=\Omega$, and
$\Omega_{2}=-0.96\,\mathrm{rad/fs}$, corresponding to $\lambda_{1}=\lambda$,
and $\lambda_{2}=1.45\,\mathrm{nm}$. Duration of the input pulses is
$t_{0}=28.4\,\mathrm{fs}$. (b) same for $t_{0}=56.7\,\mathrm{fs}$, and (c)
same for $t_{0}=85.0\,\mathrm{fs}$.

## Declaration of competing interests

The authors confirm that there are no known conflicts of interest associated
with this publication.

## Acknowledgements

We acknowledge support from the Deutsche Forschungsgemeinschaft (DFG) under
Germany’s Excellence Strategy within the Cluster of Excellence PhoenixD
(Photonics, Optics, and Engineering – Innovation Across Disciplines) (EXC
2122, projectID 390833453).

## References

  * [1] G. P. Agrawal, Nonlinear Fiber Optics, Academic Press, 2019.
  * [2] F. Mitschke, Fiber Optics: Physics and Technology, Springer, 2016.
  * [3] J. M. Dudley, G. Genty, S. Coen, Supercontinuum generation in photonic crystal fiber, Rev. Mod. Phys. 78 (2009) 1135.
  * [4] J. K. Ranka, R. S. Windeler, A. J. Stentz, Visible continuum generation in air–silica microstructure optical fibers with anomalous dispersion at 800 nm, Opt. Lett. 25 (2000) 25.
  * [5] J. C. Knight, Photonic crystal fibres, Nature 424 (2003) 847.
  * [6] T. A. Birks, J. C. Knight, P. S. J. Russell, Endlessly single-mode photonic crystal fiber, Opt. Lett. 22 (1997) 961.
  * [7] D. J. Jones, S. A. Diddams, J. K. Ranka, A. Stentz, R. S. Windeler, J. L. Hall, S. T. Cundiff, Carrier-Envelope Phase Control of Femtosecond Mode-Locked Lasers and Direct Optical Frequency Synthesis, Science 288 (2000) 635.
  * [8] F. DeMartini, C. H. Townes, T. K. Gustafson, P. L. Kelley, Self-Steepening of Light Pulses, Phys. Rev. 164 (1967) 312.
  * [9] J. R. de Oliveira, M. A. de Moura, J. M. Hickmann, A. S. L. Gomes, Self-steepening of optical pulses in dispersive media, J. Opt. Soc. Am. B 9 (1992) 2025.
  * [10] J. P. Gordon, Theory of the soliton self-frequency shift, Opt. Lett. 11 (1986) 662.
  * [11] A. Demircan, U. Bandelow, Supercontinuum generation by the modulation instability, Opt. Commun. 244 (2005) 181–185.
  * [12] A. V. Husakou, J. Herrmann, Supercontinuum Generation of Higher-Order Solitons by Fission in Photonic Crystal Fibers, Phys. Rev. Lett. 87 (2001) 203901.
  * [13] A. Demircan, U. Bandelow, Analysis of the interplay between soliton fission and modulation instability in supercontinuum generation, Appl. Phys. B 86 (2007) 31.
  * [14] D. R. Solli, C. Ropers, P. Koonath, B. Jalali, Optical rouge waves, Nature 450 (2007) 1054–1057.
  * [15] D. R. Solli, C. Ropers, B. Jalali, Active Control of Rogue Waves for Stimulated Supercontinuum Generation, Phys. Rev. Lett. 101 (2008) 233902.
  * [16] A. Demircan, S. Amiranashvili, C. Brée, C. Mahnke, F. Mitschke, G. Steinmeyer, Rogue events in the group velocity horizon, Sci. Rep. 2 (2012) 850.
  * [17] J. M. Dudley, F. Dias, M. Erkintalo, G. Genty, Instabilities, breathers and rogue waves in optics, Nat. Photonics 8 (2014) 755–764.
  * [18] A. Demircan, S. Amiranashvili, G. Steinmeyer, Controlling Light by Light with an Optical Event Horizon, Phys. Rev. Lett. 106 (2011) 163901.
  * [19] A. Demircan, S. Amiranashvili, C. Brée, U. Morgner, G. Steinmeyer, Supercontinuum generation by multiple scatterings at a group velocity horizon, Opt. Express 22 (2014) 3866–3879.
  * [20] D. V. Skryabin, A. V. Gorbach, Colloquium: Looking at a soliton through the prism of optical supercontinuum, Rev. Mod. Phys. 82 (2010) 1287–1299.
  * [21] O. Melchert, C. Brée, A. Tajalli, A. Pape, R. Arkhipov, S. Willms, I. Babushkin, D. Skryabin, G. Steinmeyer, U. Morgner, A. Demircan, All-optical supercontinuum switching, Communications Physics 3 (2020) 146.
  * [22] T. Brabec, F. Krausz, Nonlinear Optical Pulse Propagation in the Single-Cycle Regime, Phys. Rev. Lett. 78 (1997) 3282.
  * [23] H. A. Haus, Y. Lai, Quantum theory of soliton squeezing: a linearized approach, J. Opt. Soc. Am. B 7 (1990) 386.
  * [24] C. H. Henry, R. F. Kazarinov, Quantum noise in photonics, Rev. Mod. Phys. 68 (1996) 801.
  * [25] P. D. Drummond, J. F. Corney, Quantum noise in optical fibers. I. Stochastic equations, J. Opt. Soc. Am. B 18 (2) (2001) 139.
  * [26] J. F. Corney, P. D. Drummond, Quantum noise in optical fibers. II. Raman jitter in soliton communications, J. Opt. Soc. Am. B 18 (2001) 153.
  * [27] J. M. Dudley, S. Coen, Numerical Simulations and Coherence Properties of Supercontinuum Generation in Photonic Crystal and Tapered Optical Fibers, IEEE J. Sel. Top. Quantum Electron. 8 (2002) 651.
  * [28] J. M. Dudley, S. Coen, Coherence properties of supercontinuum spectra generated in photonic crystal and tapered optical fibers, Opt. Lett. 27 (2002) 1180.
  * [29] S. M. Kobtsev, S. V. Smirnov, Modelling of high-power supercontinuum generation in highly nonlinear, dispersion shifted fibers at CW pump, Opt. Express 13 (2005) 6912–6918.
  * [30] E. Genier, P. Bowen, T. Sylvestre, J. M. Dudley, P. Moselund, O. Bang, Amplitude noise and coherence degradation of femtosecond supercontinuum generation in all-normal-dispersion fibers, J. Opt. Soc. Am. B 36 (2019) A161–A167.
  * [31] R. Paschotta, Noise of mode-locked lasers (Part I): numerical model, Appl. Phys. B 79 (2004) 153.
  * [32] A. Ruehl, M. J. Martin, K. C. Cossel, L. Chen, H. McKay, B. Thomas, C. Benko, L. Dong, J. M. Dudley, M. E. Fermann, I. Hartl, J. Ye, Ultrabroadband coherent supercontinuum frequency comb, Phys. Rev. A 84 (2011) 011806.
  * [33] B. Zhou, M. Bache, Invited Article: Multiple-octave spanning high-energy mid-IR supercontinuum generation in bulk quadratic nonlinear crystals, APL Photonics 1 (2016) 050802.
  * [34] M. Ibison, B. Haisch, Quantum and classical statistics of the electromagnetic zero-point field, Phys. Rev. A 54 (1996) 2737–2744.
  * [35] M. J. Werner, M. G. Raymer, M. Beck, P. D. Drummond, Ultrashort pulsed squeezing by optical parametric amplification, Phys. Rev. A 52 (1995) 4202.
  * [36] J. Hult, A Fourth-Order Runge–Kutta in the Interaction Picture Method for Simulating Supercontinuum Generation in Optical Fibers, J. Lightwave Tech. 25 (12) (2007) 3770.
  * [37] A. M. Heidt, Efficient Adaptive Step Size Method for the Simulation of Supercontinuum Generation in Optical Fibers, J. Lightwave Tech. 27 (18) (2009) 3984.
  * [38] A. A. Rieznik, A. M. Heidt, P. G. Konig, V. A. Bettachini, D. F. Grosz, Optimum integration procedures for supercontinuum simulation, IEEE Photonics Journal 4 (2012) 552–560.
  * [39] S. Balac, A. Fernandez, SPIP: A computer program implementing the Interaction Picture method for simulation of light-wave propagation in optical fibre, Comp. Phys. Commun. 199 (2016) 139.
  * [40] O. Melchert, A. Demircan, py-fmas: A python package for ultrashort optical pulse propagation in terms of forward models for the analytic signal, Computer Physics Communications 273 (2022) 108257.
  * [41] K. J. Blow, D. Wood, Theoretical Description of Transient Stimulated Raman Scattering in Optical Fibers, IEEE J. Quantum Electron. 25 (1989) 2665.
  * [42] R. H. Stolen, J. P. Gordon, W. J. Tomlinson, Raman response function of silica-core fibers, J. Opt. Soc. Am. B 6 (1989) 1159.
  * [43] Q. Lin, G. P. Agrawal, Raman response function for silica fibers, Opt. Lett. 31 (2006) 3086.
  * [44] D. Hollenbeck, C. D. Cantrell, Multiple-vibrational-mode model for fiber-optic Raman gain spectrum and response function, J. Opt. Soc. Am. B 19 (2002) 2886\.
  * [45] C. Agger, C. Petersen, S. Dupont, H. Steffensen, J. K. Lyngsø, C. L. Thomsen, J. Thøgersen, S. R. Keiding, O. Bang, Supercontinuum generation in zblan fibers—detailed comparison between measurement and simulation, J. Opt. Soc. Am. B 29 (2012) 635.
  * [46] P. L. François, Nonlinear propagation of ultrashort pulses in optical fibers: total field formulation in the frequency domain, J. Opt. Soc. Am. B 8 (1991) 276–293.
  * [47] M. H. Frosz, O. Bang, A. Bjarklev, Soliton collision and Raman gain regimes in continuous-wave pumped supercontinuum generation, Opt. Express 14 (2006) 9391\.
  * [48] M. H. Frosz, Validation of input-noise model for simulations of supercontinuum generation and rogue waves, Opt. Express 18 (2010) 14778.
  * [49] S. J. Carter, P. D. Drummond, M. D. Reid, R. M. Shelby, Squeezing of quantum solitons, Phys. Rev. Lett. 58 (1987) 1841–1844.
  * [50] P. D. Drummond, S. J. Carter, Quantum-field theory of squeezing in solitons, J. Opt. Soc. Am. B 4 (10) (1987) 1565–1573.
  * [51] E. Brainis, D. Amans, S. Massar, Scalar and vector modulation instabilities induced by vacuum fluctuations in fibers: Numerical study, Phys. Rev. A 71 (2005) 023808.
  * [52] N. Raabe, T. Feng, T. Witting, A. Demircan, C. Brée, G. Steinmeyer, Role of Intrapulse Coherence in Carrier-Envelope Phase Stabilization, Phys. Rev. Lett. 119 (2017) 123901.
  * [53] G. Rossum, Python Reference Manual, Tech. rep., Amsterdam, The Netherlands, The Netherlands (1995).
  * [54] P. Virtanen, R. Gommers, T. E. e. a. Oliphant, SciPy 1.0: fundamental algorithms for scientific computing in Python, Nature Methods 17 (2020) 261.
  * [55] J. D. Hunter, Matplotlib: A 2D graphics environment, Computing in Science & Engineering 9 (2007) 90.
  * [56] O. Melchert, GNLStools – Data structures and functions for simulation and analysis of the propagation dynamics of laser pulses in nonlinear waveguides, https://github.com/omelchert/GNLStools.git, [Online; accessed 2022-06-07] (2022).
  * [57] T. Taha, M. Ablowitz, Analytical and Numerical Aspects of Certain Nonlinear Evolution Equations. II. Numerical, Nonlinear Schrödinger Equation, J. Comp. Phys. 55 (1984) 203.
  * [58] P. DeVries, Application of the Split Operator Fourier Transform method to the solution of the nonlinear Schrödinger equation, AIP Conf. Proc. 160 (1987) 269\.
  * [59] On the discrete computational grid we identify $\mathsf{F}$ with ifft, and $\mathsf{F}^{-1}$ with fft. Consequently, the spectral envelope in the provided implementation differs from the one defined by Eq. (4a) by an unimportant phase factor.
  * [60] G. E. P. Box, M. E. Muller, A Note on the Generation of Random Normal Deviates, The Annals of Mathematical Statistics 29 (2) (1958) 610 – 611.
  * [61] G. H. Goedecke, Stochastic Electrodynamics. I. On the Stochastic Zero-Point Field, Foundations of Physics 13 (1983) 1101.
  * [62] G. Ycas, pynlo: Nonlinear optics modeling for python, https://github.com/pyNLO/PyNLO (2013).
  * [63] L. Trefethen, Spectral Methods in MATLAB, SIAM, Philadelphia, 2000.
  * [64] B. García-Archilla, Some Practical Experience with the Time Integration of Dissipative Equations, J. Comp. Phys. 25 (1995) 122.


# NAS-FM: Neural Architecture Search for Tunable and Interpretable Sound
Synthesis based on Frequency Modulation

Zhen Ye  Hong Kong Baptist University  Wei Xue∗ Hong Kong University of
Science and Technology  Xu Tan  Microsoft Research Asia  Qifeng Liu  Hong Kong
Institute of Science & Innovation, Chinese Academy of Sciences  Hong Kong
University of Science and Technology  Yike Guo Corresponding authors: Wei Xue
<EMAIL_ADDRESS>Yike<EMAIL_ADDRESS>Hong Kong University of Science and
Technology

###### Abstract

Developing digital sound synthesizers is crucial to the music industry as it
provides a low-cost way to produce high-quality sounds with rich timbres.
Existing traditional synthesizers often require substantial expertise to
determine the overall framework of a synthesizer and the parameters of
submodules. Since expert knowledge is hard to acquire, it hinders the
flexibility to quickly design and tune digital synthesizers for diverse
sounds. In this paper, we propose “NAS-FM”, which adopts neural architecture
search (NAS) to build a differentiable frequency modulation (FM) synthesizer.
Tunable synthesizers with interpretable controls can be developed
automatically from sounds without any prior expert knowledge and manual
operating costs. In detail, we train a supernet with a specifically designed
search space, including predicting the envelopes of carriers and modulators
with different frequency ratios. An evolutionary search algorithm with
adaptive oscillator size is then developed to find the optimal relationship
between oscillators and the frequency ratio of FM. Extensive experiments on
recordings of different instrument sounds show that our algorithm can build a
synthesizer fully automatically, achieving better results than handcrafted
synthesizers. Audio samples are available at https://nas-fm.github.io/.

##  1 Introduction

Creating and rendering music has become increasingly convenient with the help
of digital sound synthesizers which simulate the timbre of real instruments
that are potentially expensive and rare. Lots of sound synthesizers are
designed whose commercial values are widely recognized, while the complexity
of the synthesizers also increases dramatically.

Synthesizer | Tunable | Interpretable |  | Automatic
---
Design
| Digital
---
Synthesizer
✓ | ✓ | ✗
| Neural Network
---
Synthesizer
✗ | ✗ | ✓
Our NAS-FM | ✓ | ✓ | ✓

Table 1: Comparison between the proposed NAS-FM with other strategies for
digital audio synthesis.

Instead of carefully designing delicate structures of highly diversified
instruments, our target is to design a general framework which can flexibly
construct digital synthesizers simply from recordings and allow further tuning
of the timbres in a controllable and interpretable way.

Early approaches to sound synthesis are parametric, designing a set of
components (e.g., oscillators and modulators) based on digital signal
processing (DSP), and ultimately creating complex timbres with specific
structures combining these components. Typical works include subtractive
synthesis Huovilainen and Välimäki (2005), additive synthesis Serra and Smith
(1990), frequency modulation (FM) Chowning (1973) and wavetable synthesis
Bristow-Johnson (1996), among which FM-based methods are widely used in sound
synthesis due to their flexibility in tuning timbre with only a few
parameters.

While satisfying sounds can be produced, designing the parametric structure of
the synthesizer as well as parameter values requires considerable expertise.
This is due to the non-linear interactions between synthesizer parameters, as
well as the wide range of possible values for each parameter. Achieving a
desired sound often involves a process of iterative adjustments and fine-
tuning, which can be time-consuming and require a deep understanding of
digital signal processing techniques. Although some efforts have been made to
determine the parameters of FM synthesizer based on estimation theory Justice
(1979), genetic algorithm Horner (1998), LSTM-based Yee-King et al. (2018),
VAE-based Le Vaillant et al. (2021) and dilated CNN based Chen et al. (2022)
methods, these methods are limited by static spectra or an audio clip
conditioned on an entire ADSR envelope under a specific pitch value and a
fixed duration length. Thus, these methods cannot be applied to the audio with
continuously varying pitch and loudness on dynamic spectra.

Neural synthesis methods have been developed recently to learn a deep neural
network to produce audio in the data-driven scheme and achieved impressive
results in terms of audio fidelity. However, for fully-neural generative
models such as WaveGlow Prenger et al. (2019), and HiFi-GAN Kong et al.
(2020a), sufficient data is usually required to train a complex model with
non-interpretable parameters, which greatly limits the controllability of the
synthesizer as required by the music industry, specifically when the generated
non-perfect sound needs tuning. Although the harmonic-plus-noise model-based
differential DSP (DDSP) is integrated into the network to improve
controllability Engel et al. (2020), it only allows the transfer of timbres
between different instruments rather than explicit adjustment. FM-based neural
synthesizer DDX7 Caspe et al. (2022) is then developed. However, this method
can only learn the time-varying parameters of the submodules under the
assumption that the overall structure has been manually defined.

The large reliance on expert knowledge and extensive time cost greatly limit
the feasibility of building digital synthesizers for general instruments.
Moreover, the handcrafted framework may make the resulting synthesizer
suboptimal in modelling the real instrument. This paper proposes NAS-FM, a
neural architecture search (NAS) based FM synthesizer. It can be seen as an
effort to build interpretable neural generative models. A key of the proposed
method is NAS. With a carefully designed search procedure, different
structures and oscillator sizes are included in a universal search space of
the supernet, and an evolutionary search algorithm is developed to find the
optimal structure between oscillators and the frequency ratio of FM. The
advantages of the proposed NAS-FM are described below:

  * •

The audio synthesizer can be built based on recordings without any expert
knowledge, largely simplifying the pipeline and reducing the cost of audio
synthesizer construction. It is also possible to quickly build new variants of
the target sound with flexible adjustments;

  * •

It returns a tunable and interpretable conventional FM-based interface with a
few parameters, making it easily embedded into existing audio workstations;

  * •

Extensive experiments on recordings of different instruments demonstrate that
synthesizers fully automatically built by the proposed NAS-FM can achieve
better results to carefully handcrafted synthesizers.

##  2 Related Work

###  2.1 Sound Synthesis

Sound synthesis contains digital signal processing (DSP) methods and neural
network methods. DSP methods have been integrated into the digital audio
workstation used by musicians. More specifically, DSP methods start from
several simple waves such as sine, sawtooth, square and triangle generated by
an oscillator. The additive synthesizer Serra and Smith (1990) generates new
sounds by adding various simple waves. The Subtractive synthesizer Huovilainen
and Välimäki (2005) filters a simple wave from white noises. FM synthesizers
Chowning (1973) rely on simple waves to modulate frequency to create complex
timbre. Wavetable synthesis manipulate a collection of short samples to design
new sounds. These traditional methods need users to determine the
configuration manually for a given sound.

Neural network synthesis models adopt the deep neural network to learn the
mapping function between audio and given input, for instance, pitch and
loudness. The early exploration begins with auto-regressive models such as
WaveRNN Kalchbrenner et al. (2018) and SampleRNN Mehri et al. (2016). The
following works Engel et al. (2019) Kong et al. (2020b) are based on various
generative models to further improve the quality of synthesis sound. However,
the above methods may lead to glitch problems because of the lack of phase
continuity. Therefore, DDSP Engel et al. (2020) rely on the harmonic-plus-
noise model to keep the phase continuous and also make the sound can be
directly controlled by pitch and loudness. While these methods can be
optimized automatically with the help of gradient descent to obtain the model,
there are few control factors to help users manipulate the synthesized result
directly. Therefore, our approach aims to introduce the FM synthesizer with
controllable factors to help the user interact with the synthesized audio.

###  2.2 FM Parameter Estimation

FM parameter estimation also called FM matching is adopted to determine the
configuration of an FM synthesizer. The early approach considers this problem
as a searching problem that uses the genetic algorithm(GA) Jh (1975) or its
variants to find a best-fit configuration in a specific searching space.
Horner employs GA solving a sound matching problem for different FM algorithms
such as Formant FM Horner et al. (1993), Double FM Horner (1996), Nested FM
and Feedback FM Horner (1998). These methods can achieve very close re-
synthesizing results with a static target spectrum when selecting an
appropriate FM algorithm as prior. The following methods leverage an open-
source FM synthesizer Dexed synthesizer Gauthier (2013) to construct a large
number of pair data between presets111Preset means a full FM configuration for
specific sound designed by userscollected on the Internet and synthesized
audio clip generated by Dexed. By reversing this process, these works employ
LSTM Yee-King et al. (2018), VAE Le Vaillant et al. (2021) and dilated CNN
Chen et al. (2022) to estimate the preset. Since the synthesized sound is
generated by pressing a note with specific velocity, duration and pitch using
Dexed synthesizer, the result of the model for a realistic sound audio clip
without any prior is unpredictable. DDX7 Caspe et al. (2022) predict the
envelopes of oscillators using the widely-used TCN Oord et al. (2016) decoder
with algorithm and frequency ratios as prior. Therefore, our method is
designed to construct an FM synthesizer without prior knowledge.

###  2.3 Neural Architecture Search

Neural Architecture Search (NAS) can automatically find the best neural
architecture for a specific task. Many works focus on computer vision Liu et
al. (2018)or natural language processing Xu et al. (2021) tasks. More
recently, one-shot NAS Guo et al. (2020) Bender et al. (2018) train a shared
supernet once which includes all candidate architectures. Then, the supernet
is used as an estimator to evaluate every possible architecture in the search
space. This method has been widely used on various applications to determine a
good structure, for example, image classification Guo et al. (2020), object
detection Liang et al. (2021), 3d scene understanding Tang et al. (2020), BERT
compression Xu et al. (2022). These methods introduce neural architecture
search to each specific task which gets a better model architecture than
manual design. Although NAS has been widely used in lots of areas, the
applications are mainly focused on neural networks. Actually, we are the first
to bring NAS to fully automatically build the sound synthesizer.

##  3 Review of FM Synthesizer

###  3.1 Basics of FM

FM was originally proposed in Chowning (1973) for sound synthesis and
different timbres are produced by controlling a set of parameters. With two
sound sources, the modulator oscillator $\sin(2\pi{f_{m}}t)$ and the carrier
oscillator $\sin(2\pi{f_{c}}t)$, FM basically generates a time-domain signal
$y(t)$:

$y(t)=a(t)\sin(2\pi{f_{c}}t+I\sin(2\pi{f_{m}}t)),$ (1)

where $f_{c}$ and $f_{m}$ are the carrier frequency and modulation frequency
respectively, $I$ is the modulation index, and $a(t)$ is the amplitude
envelope of the carrier. $y(t)$ can be further decomposed by using Bessel
functions of the first kind as

$y(t)=a(t)\sum\limits_{n=-\infty}^{n=+\infty}{{J_{n}}(I)\sin(2\pi({f_{c}}+n{f_{m}})\cdot
t)}$ (2)

which shows that the sidebands of $y(t)$ distribute evenly around $f_{c}$ with
spacing as $f_{m}$, and the spectra is harmonic when the frequency ratio
$r\in\mathbb{Q}$ where $r=f_{c}/f_{m}$. ${J_{n}}(I)$ is a Bessel function of
the modulation index $I$, and dynamic spectra can be generated if the $I$
becomes a time-variant function $I(t)$.

Figure 1: Different FM algorithms: (a) Single FM; (b) Nested FM; (c) Formant
FM; (d) Double FM. The green box refers to the carrier and the blue box refers
to the modulator Figure 2: Digital synthesizers in the commercial product
YAMAHA DX7. The lower part is the user interface of a synthesizer. The upper
part shows the FM algorithms for audio synthesis, with green and blue boxes
denoting carriers and modulators, respectively.

###  3.2 FM algorithms

The (1) explains the basic module of the FM synthesizer to generate sounds.
More diversified and expressive FM synthesizers can be further developed by
designing complicated topologies, which are called “FM algorithms”, on
connecting the carrier and modulator oscillators.

Typical FM algorithms are shown in Fig. 1. Fig. 1(a) denotes the single FM
expressed by (1). The Nested FM Justice (1979), formant FM Horner et al.
(1993), and double FM Schottstaedt (1977) are illustrated in Fig. 1(b)-(d),
whose outputs are calculated by

$\displaystyle y(t)=$ $\displaystyle
a(t)\sin(2\pi{f_{c}}t+{I_{m1}(t)}\sin[2\pi{f_{m1}}t$
$\displaystyle+{I_{m2}(t)}\sin(2\pi{f_{m2}}t)]),$ (3) $\displaystyle y(t)=$
$\displaystyle{a_{1}}(t)\sin[2\pi{f_{c1}}t+{I_{m}(t)}\sin(2\pi{f_{m}}t)]$
$\displaystyle+{a_{2}}(t)\sin[2\pi{f_{c2}}t+{I_{m}(t)}\sin(2\pi{f_{m}}t)],$
(4)

and

$\displaystyle y(t)=$ $\displaystyle
a(t)\sin[2\pi{f_{c}}t+{I_{m1}(t)}\sin(2\pi{f_{m1}}t)$
$\displaystyle+{I_{m2}(t)}\sin(2\pi{f_{m2}}t)],$ (5)

respectively. These FM algorithms produce different timbres, and as shown in
Fig. 2222We do not consider the feedback FM in our work as the same reason in
Caspe et al. (2022), by selecting different FM algorithms.

To briefly summarize, a set of parameters controls the FM output, which can be
time-variant or fixed. In our work, time-variant parameters include the
oscillators’ envelopes, i.e., $a(t)$ and $I(t)$, which are determined by the
time-variant input, such as $f_{0}$ and loudness. Fixed parameters are the
chosen FM algorithm and the frequency ratio of each oscillator, which need to
be determined by the expected timbre.

Figure 3:  The overall architecture of NAS-FM. The learnable Reverb module is
optional depending on whether the real room effect should be simulated for
real-world audios.

##  4 Proposed NAS-FM synthesizer

In this section, we aim to fully automatize the designing of an FM synthesizer
in a data-driven manner, thus eliminating the reliance on expertise and labour
to tune the timbres.

A NAS-FM synthesizer is proposed, with an overall framework shown in Fig. 3 on
the next page. For input audio, the pitch and loudness are extracted, and
these two features are fed into an oscillator envelope prediction network to
estimate the envelopes of the oscillators. We develop architecture search
methods to determine the optimal FM algorithm, and given the FM algorithm and
oscillator envelopes, sounds with a specific timbre can be produced. Depending
on whether the sound is from the real environment, a learnable reverb module
Engel et al. (2020) can be optionally used to simulate the room effects. The
framework is optimized in an auto-encoder setting, i.e., seeking to recover
the original signal in the final output.

Besides the FM algorithm search, the framework is similar to the framework in
DDX7 Caspe et al. (2022), which estimates pitch and loudness by CREPE Kim et
al. (2018) and A-weighting loudness Moore et al. (1997), designs the
oscillator envelope prediction network as a temporal convolutional network
(TCN) Bai et al. (2018). However, DDX7 requires a prior FM configuration
designed by experts. By adopting NAS, the proposed framework designs the FM
synthesizer in the fully-automated data-driven pipeline and also returns a
tunable and interpretable interface used by musicians. In the following,
details of the proposed NAS-FM will be introduced, which include a) converting
the FM algorithm to a graph, b) designing the search space, c) training the
supernet, and d) selecting the FM algorithm.

###  4.1 Directed Acyclic Graph of FM Algorithm

To facilitate discussion in this section, we convert an FM synthesizer, with
examples shown in Fig. 1, to a directed acyclic graph (DAG) Pham et al. (2018)
with an ordered sequence of $N$ nodes, where $N$ represents the number of
oscillators. Each node $x^{(i)}$ refers to the oscillator’s output. Each
directed edge $(i,j)$ indicates whether node $x^{(j)}$ is the modulating node
of $x^{(i)}$. The topology of the graph is associated with the FM algorithm.
Thus, an intermediate node $x^{(i)}$ can be expressed as

${x^{(i)}(t)}={a_{i}}(t)\sin(2\pi{f_{i}}t+\sum\limits_{j\in\mathbb{M}}{{x^{(j)}(t)}}),$
(6)

where $\mathbb{M}$ is a set of oscillators modulating node $x^{(i)}$. When $M$
is empty, $x^{(i)}$ outputs standard sine wave. The output of FM is calculated
through the sum of carrier nodes

$y(t)=\sum\limits_{i\in\mathbb{C}}{{x^{(i)}}(t)},$ (7)

where $\mathbb{C}$ is the set of carriers.

###  4.2 Search Space Design

There are a huge number of possible configurations for the FM synthesizer. By
converting the FM algorithm to the DAG, the principles of NAS can be applied
to design the FM algorithm. In NAS for neural networks, all possible
architectures of the network, which span the “search space”, can be
represented by a general DAG, with each candidate architecture as a sub-graph.
Similarly, we design the search space for FM algorithm here.

We visualize the operation in (6) in Fig. 4, where two issues should be solved
to define an oscillator: a) what is the frequency ratio? The frequency ratio
is defined originally in (2), and is the ratio between $f_{i}$ in (6) and
fundamental frequency $F_{0}$ here; b) which modulators will be connected to
the current oscillator? As depicted in Fig. 4, we define the frequency ratio
set $\mathbb{FR}=\\{1,2,3,..,K\\}$ which consists of $K$ integers, and assumes
that there are $N$ candidates in the modulator set $\mathbb{M}$. Actually, the
interval in the frequency ratio set can be reduced to 0.1 or even smaller, to
0.01 instead of 1, for a more refined exploration.

Traversing the original search space is impossible considering the exponential
number of possibilities for oscillator connections. Now let us analyze the
design of the search space. Previous works adopted evolutionary search on
frequency ratio set under a fixed FM algorithm such as double FM Horner (1996)
and nested FM Horner (1998). Therefore, the challenge is how to design an
appropriate space including various FM algorithms, instead of enumerating all
possible FM algorithms. In Guo et al. (2020) and Xu et al. (2022), the weight
sharing is proposed to reduce the search space, which forces sub-graph
candidates to have the same weights in the graph nodes commonly shared by
different candidates. This also makes it possible to directly train a supernet
including all possible candidates. Here we propose a novel envelope-sharing
strategy for FM synthesizers to make all FM configuration candidates trained
in the supernet.

Specifically, we construct a search space as shown in Fig. 3. The search space
can be divided into a carrier layer and several modulator layers. We set the
oscillator at the same layer with the same frequency ratios sharing the same
envelope. Due to the envelope-sharing strategy, There are following rules in
our search space. (1) A oscillator in a certain layer can only be modulated by
the upper layer; (2) The sum of the output of selected oscillators at the
carrier layer forms the final signal;(3) A oscillator is discarded when there
is no connection with other oscillators or final output. In addition,
according to our experience, we find two modulator layers are enough. The
number of candidate oscillators in a layer depends on the oscillator number in
the expected FM synthesizer.

Figure 4:  The oscillator in NAS-FM. The output of an oscillator conditioning
on $F0$ depends on the choice of frequency ratio, envelope and modulator. The
choices are determined by NAS.

###  4.3 Supernet Training with Proxy Oscillator

Although the proposed envelope-sharing strategy makes all possible FM
configurations can be trained using a supernet that includes all candidates,
we still encounter the challenge of a large search space. The huge search
space makes it hard to evaluate each FM configuration in the supernet during
training. To further improve the training efficiency, we propose to use the
“proxy oscillator” as the proxy for all oscillators in each layer to determine
the envelopes of the oscillators.

Specifically, during supernet training, a fixed Nested FM in (3.2) with a
carrier and two modulators are chosen, which has three layers in accordance
with the settings in the Sec. 4.2. For each oscillator, uniform sampling Guo
et al. (2020) is adopted to determine the frequency ratio. With these two
configurations fixed, the envelopes of the oscillators are learned, and after
training, the learned envelopes are shared within the same layer for a
complicated search space. Utilizing the technique, we can expand the width of
the search space flexibly.

###  4.4 FM Algorithm Selection

After the supernet training and search space design, we use the evolutionary
algorithm to conduct FM algorithm selection. Specifically, we put all $N$
candidate oscillators $ol$ in order. a certain FM configuration is encoded as
an individual can be formulated as

$\\{{f_{o{l_{1}}}},{f_{o{l_{2}}}},...,{f_{o{l_{N}}}},{l_{o{l_{1}}}}{l_{o{l_{2}}}},...,{l_{o{l_{N}}}}\\}$
(8)

where $ol_{i}$ is the $i_{th}$ oscillator, corresponding ${f_{o{l_{i}}}}$ and
${l_{o{l_{i}}}}$ indicate the selected frequency ratio and connection
relationship. The connection relationship is the relation between the lower
layer. If the oscillator belongs to the carrier, the lower layer is the output
signal. If the connection relationship of this oscillator is none means it is
discarded. Firstly, a random population is initialized within the initial
space. Secondly, we evaluate the fitness score of the generated individuals
and select the top individuals. The fitness function will be introduced in the
experiment section. Thirdly, crossover and mutation are used to generate new
individuals to update the population, until meeting the stopping criterion.

##  5 Experiments

To evaluate the proposed NAS-FM approach which aims to automatically learn the
FM synthesizer that is applicable to the music industry, we aim to answer the
following questions: a) Given sound recordings, can the FM synthesizers
learned by the proposed method be comparable to the manually designed
synthesizers? b) We fuse different FM algorithms into one universal search
space, is this strategy better than separately searching the best
configuration for each FM algorithm? c) Can we controllably tune the timbre of
the produced sounds or create new instrument timbres by modifying the
parameters of the learned FM synthesizer?

###  5.1 Dataset

We conduct experiments on the benchmark URMP dataset Li et al. (2018). Sound
recordings of three real instruments, which are violin, flute, and trumpet,
are chosen and the “optimal” manually designed FM algorithms of these
instruments are given as in Caspe et al. (2022). The three instruments are
also the most typical instrument of the strings, woodwinds, and brass. Each
instrument recording is divided into 4-second segments, and silent segments
are discarded. The loudness and pitch are estimated by the A-weighting
loudness Moore et al. (1997) and CREPE Kim et al. (2018) methods. Each audio
clip is resampled to the 16kHz sampling rate and analyzed with a frame size of
2048 and a hop size of 64, yielding 1000 frames. We split the dataset into
train, validation, and test sets with proportions of 0.75, 0.125, and 0.125,
respectively.

###  5.2 Experimental Setup

For training, a supernet contains $c*m*m$ paths where $c$ is the number of
frequency ratios of the carrier, and $m$ is the number of frequency ratios in
modulators. We set $c$ as 15 and $m$ as 5. The model adopts the stack TCN
architecture Bai et al. (2018) as a sequence-to-sequence model to predict
oscillator envelope through pitch and loudness. We first stack 4 TCN
architecture to extract a hidden temporal feature. Then, we construct
$(c+m+m)$ TCN architectures with different weights to predict the envelope of
the corresponding oscillator from the hidden feature. In fact, other sequence-
to-sequence models can also be employed to model the temporal relationship in
our pipeline. In addition, uniform sampling of oscillators on each layer with
different ratios is adopted. We use Adam optimizer with an initial learning
rate of 3e-4. For regularize, we set the maximum value of the oscillator for
the carrier and modulator as 1 and 2, respectively. The exponential decay
strategy with a decreasing factor of 0.98 every 10k steps is adopted for the
learning rate .The whole training steps are 500k, and the batch size is 16.

For searching, we use the evolutionary algorithm with a population size of
$P$, crossover size of $P/2$, and mutation size of $P/2$ with mutation
probability and max iterations of $T$ depending on the search space size. In
addition, since we use the fixed number of the candidate oscillator, the
discarded candidate oscillator with different frequency ratios will cause the
same fitness score. Therefore, we force the generated candidates with a unique
fitness score to be legal.

After searching, we directly extract the weight of the sub-model with the best
fitness score from the supernet as our final model. Actually, fine-tuning the
model or training from scratch with the searched FM configuration may further
improve the performance. However, our aim is not to focus on the resynthesis
performance but to get a controllable synthesizer close to the target sound so
that musicians can further utilize it to tweak the sound or do other more
interesting things such as sound morphing and sound interpolation.

###  5.3 Evaluation Metric

We use Fréchet Audio Distance (FAD) Kilgour et al. (2018) as the evaluation
metric to measure the distance between real and generated sound. This method
extracts the embedding of the audio using a pre-trained VGG-like model. The
distance is calculated as follows:

$FAD={\left\|{{\mu_{r}}-{\mu_{g}}}\right\|^{2}}+tr({\Sigma_{r}}+{\Sigma_{g}}-2\sqrt{{\Sigma_{r}}{\Sigma_{g}}})$
(9)

where $r$ and $g$ are the real audio and generated audio respectively,
${\mu_{r}}$ and ${\mu_{g}}$ are the mean vector of embedding and
${\Sigma_{r}}$ and ${\Sigma_{g}}$ are the covariances of embedding. A smaller
Fréchet Audio Distance indicates a higher level of similarity between the
distributions of real and generated data. In our experiment, we calculate the
Fréchet Audio Distance between the real data and synthesize validation data as
the fitness score during searching. And calculating the FAD between the
synthesized test data as the final evaluation result.

| Fréchet Audio Distance (↓)
---|---
Model | Flute | Violin | Trumpet
Test Data | 1.180 | 0.308 | 0.554
DDX7 | 7.841 | 3.497 | 4.442
NAS-FM | 7.077 | 3.255 | 3.384
Table 2: FAD of resynthesis results using 6 oscillators

###  5.4 Comparison with Manually-Designed FM Synthesizers

In this part, we are interested to know if our method could achieve comparable
results to manually designed FM synthesizers. The baseline is DDX7 Caspe et
al. (2022) which the author retrieves on the web manually to find the patch
with the most similar sound to the target sound. The patches they find have
six oscillators, and we adopt the same oscillator number with them. Since the
authors did not open source their evaluation code, we train the model
following their method ten times with a random seed and took the best one as
their result.

In our method, we use a $3*3$ candidate oscillator search space forcing three
oscillators discarded to ensure six oscillators. During the search, we follow
the above procedure and set the population size as 1000 and max iterations as
50. Results are shown in Table 2. The Test Data line means the Fréchet Audio
Distance between the entire real data from the real test data. We can see that
our NAS-FM outperforms the hand-designed DDX7 baseline across all musical
instrument recordings. The results show that our NAS-FM can search for more
comparative FM configurations than manually designed FM synthesizers.

| Fréchet Audio Distance (↓)
---|---
Model | Flute | Violin | Trumpet
Nested FM | 12.75 | 6.02 | 8.24
Formant FM |  14.48 | 7.56 | 8.68
Double FM | 11.16 | 7.27 | 9.91
Single+ FM | 11.20 | 12.07 | 9.28
NAS-FM | 11.16 | 6.02 | 8.24
Table 3: Ablated of search space using 3 oscillators

###  5.5 Ablation Study of Search Space

The previous methods of FM parameter estimation Chen et al. (2022) Mitchell
and Sullivan (2005) Horner et al. (1993) were usually constrained to a
specific FM algorithm. However, our approach puts the different FM algorithms
in one search space. To demonstrate the significance of the method, we first
want to answer whether different timbres have their own FM algorithms that are
more suitable for them when operator size is fixed. If the answer is yes, can
our method find the algorithm that best fits this timbre? We conduct our
experiments conditioned on three oscillators and enumerate all possible FM
algorithms:

  * •

Nested FM: a carrier with two nested modulators as shown in Fig. 1(b);

  * •

Formant FM: two carriers sharing a modulator as shown in Fig. 1(c);

  * •

Double FM: a carrier with two modulators in the same row as shown in Fig.
1(d);

  * •

Single FM+: a single oscillator adds a single FM. Single FM is shown in Fig.
1(a);

  * •

NAS-FM (Ours): Constructing a $3*2$ candidate oscillator search space forcing
three oscillators discarded to search both FM algorithms and frequency ratios.

During the search process, the search space of a fixed FM algorithm equals
only contains the frequency ratios of each oscillator. In addition, we set the
population size as 30 and max iterations as 20 for each method for a fair
comparison. The results are shown in Table 3.

To answer the first question, we compare the result across all FM algorithms.
The Nested FM achieves the best performance on violin and trumpet. However,
this algorithm yields the second-worst performance on the Flute. The double FM
performs the best on the Flute but worst on the trumpet. Therefore, we find
that no fixed FM algorithm achieves the best results for every musical
instrument recording. However, our method can always find the best FM
configuration for different instruments under the same searching setting. This
proves the merit of our well-designed search space.

###  5.6 Tuning of Learned FM synthesizers

In this section, we check the tuning ability of NAS-FM. Two important tasks in
sound synthesis are considered: sound morphing and timbre interpolation.

Figure 5: Examples of sound morphing. The top figure is the synthesized sound
of the trumpet. In the middle and bottom figures, we morph the sound by tuning
a certain parameter in our NAS-FM.

####  5.6.1 Sound Morphing

Due to the timbre of an instrument being controlled by a set of frequency
ratio parameters, we want to show the ability of our NAS-FM tweaks the timbre
from a trained synthesizer by adjusting a few parameters. Begin with a
synthesizer designed for the trumpet, Fig. 5 (a) is the synthesized recording
of the trumpet using our NAS-FM. More specifically, the searched FM algorithm
is the same as the 7-th algorithm shown in Fig. 2 without the feedback module.
It consists of two parts. The left part is a single FM in which the carrier
frequency ratio is three and the modulator is one. The right part includes
four oscillators which consist of a double FM with one more nested modulator.
The frequency ratios among a carrier, double modulators and a nested modulator
are 7,1,2 and 1, respectively. In Fig. 5(b), we modify the second formant
position from the 7th harmonic to the 10th harmonic by modifying the frequency
ratio of the carrier in the right part from seven to ten. In addition, we
decrease the second and fourth harmonic by tuning the frequency ratio of the
modulator in the left part from one to two, as shown in Fig. 5(c). We can find
that due to the adjustment of one parameter in NAS-FM, the entire distribution
and details of frequency can be easily changed.

####  5.6.2 Timbre Interpolation

Another meaningful application is timbre interpolation. Since our approach
determines the FM configuration by searching for the Fréchet Audio Distance to
the target audio clip. We find that the target audio clips can produce new
timbre. Therefore, we can change our fitness score to a variant form. For
instance, we set the FAD from synthesized audio to violin as $d_{v}$ and the
FAD from synthesized audio to trumpet as $d_{t}$. Next, we define the novel
fitness function as ${d_{v}}+{d_{a}}+\left|{{d_{v}}-{d_{a}}}\right|$ aiming to
find a synthesizer to generate sound regarded as an intermediate timbre
between violin and trumpet. Surprisingly, we use the supernet trained from
violin recordings to conduct the evolutionary algorithm and get an interesting
result. As shown in Fig 6, we obtain a new instrument recording similar to
both violin and trumpet.

Figure 6: An example of timbre interpolation. The top and bottom figures show
the linear spectrogram of the synthesized sound of the trumpet and violin,
respectively. The middle figure shows the interpolation result between the
violin and the trumpet.

##  6 Conclusion

We present NAS-FM, a tunable and interpretable sound synthesizer based on
frequency modulation. Given a target sound, we prove that our method can
automatically design an FM synthesizer instead of spending a huge time
designing it manually. Meanwhile, our auto-designed synthesizer can achieve
comparable results to the handcrafted one. Furthermore, Our NAS-FM leverages
the widely-used FM synthesizer as the main component in our framework. This
makes our method can be directly understood by musicians that can be used to
create new sounds without extra knowledge of the neural networks.

## Acknowledgments

The research was supported by the Theme-based Research Scheme (T45-205/21-N)
and Early Career Scheme (ECS-HKUST22201322), Research Grants Council of Hong
Kong.

## References

  * Bai et al. [2018] Shaojie Bai, J Zico Kolter, and Vladlen Koltun.  An empirical evaluation of generic convolutional and recurrent networks for sequence modeling.  arXiv preprint arXiv:1803.01271, 2018.
  * Bender et al. [2018] Gabriel Bender, Pieter-Jan Kindermans, Barret Zoph, Vijay Vasudevan, and Quoc Le.  Understanding and simplifying one-shot architecture search.  In Proc. Intl. Conf. Machine Learning (ICML), pages 550–559. PMLR, 2018.
  * Bristow-Johnson [1996] Robert Bristow-Johnson.  Wavetable synthesis 101, a fundamental perspective.  In Audio engineering society convention 101. Audio Engineering Society, 1996.
  * Caspe et al. [2022] Franco Caspe, Andrew McPherson, and Mark Sandler.  DDX7: Differentiable FM Synthesis of Musical Instrument Sounds.  Proc. Intl. Soc. for Music Information Retrieval Conf. (ISMIR), 2022.
  * Chen et al. [2022] Zui Chen, Yansen Jing, Shengcheng Yuan, Yifei Xu, Jian Wu, and Hang Zhao.  Sound2synth: Interpreting sound via fm synthesizer parameters estimation.  Proc. Intl. Joint Conf. on Artificial Intelligence (IJCAI), 2022\.
  * Chowning [1973] John M Chowning.  The synthesis of complex audio spectra by means of frequency modulation.  J. Audio Eng. Soc., 21(7):526–534, 1973.
  * Engel et al. [2019] Jesse Engel, Kumar Krishna Agrawal, Shuo Chen, Ishaan Gulrajani, Chris Donahue, and Adam Roberts.  GANsynth: Adversarial neural audio synthesis.  Proc. Intl. Conf. on Learning Representations (ICLR), 2019.
  * Engel et al. [2020] Jesse Engel, Lamtharn Hantrakul, Chenjie Gu, and Adam Roberts.  DDSP: Differentiable digital signal processing.  Proc. Intl. Soc. for Music Information Retrieval Conf. (ISMIR), 2020.
  * Gauthier [2013] Pascal Gauthier.  dexed.  https://github.com/ asb2m10/dexed, 2013.
  * Guo et al. [2020] Zichao Guo, Xiangyu Zhang, Haoyuan Mu, Wen Heng, Zechun Liu, Yichen Wei, and Jian Sun.  Single path one-shot neural architecture search with uniform sampling.  In Proc. European Conf. on Computer Vision (ECCV), pages 544–560. Springer, 2020.
  * Horner et al. [1993] Andrew Horner, James Beauchamp, and Lippold Haken.  Machine tongues XVI: Genetic algorithms and their application to FM matching synthesis.  Computer Music Journal, 17(4):17–29, 1993.
  * Horner [1996] Andrew Horner.  Double-modulator fm matching of instrument tones.  Computer Music Journal, 20(2):57–71, 1996.
  * Horner [1998] Andrew Horner.  Nested modulator and feedback fm matching of instrument tones.  IEEE Trans. Speech Audio Process., 6(4):398–409, 1998.
  * Huovilainen and Välimäki [2005] Antti Huovilainen and Vesa Välimäki.  New approaches to digital subtractive synthesis.  In Proc. Intl. Computer Music Conf. (ICMC), 2005.
  * Jh [1975] Holland Jh.  Adaptation in natural and artificial systems.  Ann Arbor, 1975.
  * Justice [1979] James Justice.  Analytic signal processing in music computation.  IEEE Trans. Acoust., Speech, Signal Process., 27(6):670–684, 1979\.
  * Kalchbrenner et al. [2018] Nal Kalchbrenner, Erich Elsen, Karen Simonyan, Seb Noury, Norman Casagrande, Edward Lockhart, Florian Stimberg, Aaron Oord, Sander Dieleman, and Koray Kavukcuoglu.  Efficient neural audio synthesis.  In Proc. Intl. Conf. Machine Learning (ICML), pages 2410–2419. PMLR, 2018.
  * Kilgour et al. [2018] Kevin Kilgour, Mauricio Zuluaga, Dominik Roblek, and Matthew Sharifi.  Fréchet audio distance: A metric for evaluating music enhancement algorithms.  arXiv preprint arXiv:1812.08466, 2018.
  * Kim et al. [2018] Jong Wook Kim, Justin Salamon, Peter Li, and Juan Pablo Bello.  Crepe: A convolutional representation for pitch estimation.  In Proc. Intl. Conf. on Acoustics, Speech, and Signal Processing (ICASSP), pages 161–165. IEEE, 2018.
  * Kong et al. [2020a] Jungil Kong, Jaehyeon Kim, and Jaekyoung Bae.  Hifi-GAN: Generative adversarial networks for efficient and high fidelity speech synthesis.  Proc. Advances in Neural Information Processing Systems (NIPS), 33:17022–17033, 2020.
  * Kong et al. [2020b] Zhifeng Kong, Wei Ping, Jiaji Huang, Kexin Zhao, and Bryan Catanzaro.  Diffwave: A versatile diffusion model for audio synthesis.  Proc. Intl. Conf. on Learning Representations (ICLR), 2020.
  * Le Vaillant et al. [2021] Gwendal Le Vaillant, Thierry Dutoit, and Sébastien Dekeyser.  Improving synthesizer programming from variational autoencoders latent space.  In Proc. Conf. on Digital Audio Effects, pages 276–283. IEEE, 2021\.
  * Li et al. [2018] Bochen Li, Xinzhao Liu, Karthik Dinesh, Zhiyao Duan, and Gaurav Sharma.  Creating a multitrack classical music performance dataset for multimodal music analysis: Challenges, insights, and applications.  IEEE Trans. Multimedia, 21(2):522–535, 2018.
  * Liang et al. [2021] Tingting Liang, Yongtao Wang, Zhi Tang, Guosheng Hu, and Haibin Ling.  Opanas: One-shot path aggregation network architecture search for object detection.  In Proc. IEEE Conf. on Computer Vision and Pattern Recognition (CVPR), pages 10195–10203, 2021.
  * Liu et al. [2018] Hanxiao Liu, Karen Simonyan, and Yiming Yang.  Darts: Differentiable architecture search.  Proc. Intl. Conf. on Learning Representations (ICLR), 2018.
  * Mehri et al. [2016] Soroush Mehri, Kundan Kumar, Ishaan Gulrajani, Rithesh Kumar, Shubham Jain, Jose Sotelo, Aaron Courville, and Yoshua Bengio.  Samplernn: An unconditional end-to-end neural audio generation model.  In Proc. Intl. Conf. on Learning Representations (ICLR), 2016\.
  * Mitchell and Sullivan [2005] Thomas Mitchell and Charlie Sullivan.  Frequency modulation tone matching using a fuzzy clustering evolution strategy.  In Audio Engineering Society Convention 118. Audio Engineering Society, 2005.
  * Moore et al. [1997] Brian CJ Moore, Brian R Glasberg, and Thomas Baer.  A model for the prediction of thresholds, loudness, and partial loudness.  J. Audio Eng. Soc., 45(4):224–240, 1997.
  * Oord et al. [2016] Aaron van den Oord, Sander Dieleman, Heiga Zen, Karen Simonyan, Oriol Vinyals, Alex Graves, Nal Kalchbrenner, Andrew Senior, and Koray Kavukcuoglu.  Wavenet: A generative model for raw audio.  arXiv preprint arXiv:1609.03499, 2016.
  * Pham et al. [2018] Hieu Pham, Melody Guan, Barret Zoph, Quoc Le, and Jeff Dean.  Efficient neural architecture search via parameters sharing.  In Proc. Intl. Conf. Machine Learning (ICML), pages 4095–4104. PMLR, 2018.
  * Prenger et al. [2019] Ryan Prenger, Rafael Valle, and Bryan Catanzaro.  Waveglow: A flow-based generative network for speech synthesis.  In Proc. Intl. Conf. on Acoustics, Speech, and Signal Processing (ICASSP), pages 3617–3621. IEEE, 2019.
  * Schottstaedt [1977] Bill Schottstaedt.  The simulation of natural instrument tones using frequency modulation with a complex modulating wave.  Computer Music Journal, pages 46–50, 1977.
  * Serra and Smith [1990] Xavier Serra and Julius Smith.  Spectral modeling synthesis: A sound analysis/synthesis system based on a deterministic plus stochastic decomposition.  Computer Music Journal, 14(4):12–24, 1990.
  * Tang et al. [2020] Haotian Tang, Zhijian Liu, Shengyu Zhao, Yujun Lin, Ji Lin, Hanrui Wang, and Song Han.  Searching efficient 3d architectures with sparse point-voxel convolution.  In Proc. European Conf. on Computer Vision (ECCV), pages 685–702. Springer, 2020.
  * Xu et al. [2021] Jin Xu, Xu Tan, Renqian Luo, Kaitao Song, Jian Li, Tao Qin, and Tie-Yan Liu.  Nas-bert: task-agnostic and adaptive-size bert compression with neural architecture search.  In Proc. ACM SIGKDD Intl. Conf. on Knowledge Discovery and Data Mining (SIGKDD), pages 1933–1943, 2021.
  * Xu et al. [2022] Jin Xu, Xu Tan, Kaitao Song, Renqian Luo, Yichong Leng, Tao Qin, Tie-Yan Liu, and Jian Li.  Analyzing and mitigating interference in neural architecture search.  In Proc. Intl. Conf. Machine Learning (ICML), pages 24646–24662. PMLR, 2022.
  * Yee-King et al. [2018] Matthew John Yee-King, Leon Fedden, and Mark d’Inverno.  Automatic programming of vst sound synthesizers using deep networks and other techniques.  IEEE Trans. Emerging. Topics in Comput. Intel., 2(2):150–159, 2018.