# Instantaneous normal modes in liquids: a heterogeneous-elastic-medium
approach

Stefano Mossa  Univ. Grenoble Alpes, CEA, IRIG-MEM-LSim, 38054 Grenoble,
France  Taras Bryk  Institute for Condensed Matter Physics,National Academy of
Sciences of Ukraine, UA-79011 Lviv, Ukraine  Institute of Applied Mathematics
and Fundamental Sciences, Lviv National Polytechnic University, UA-79013 Lviv,
Ukraine  Giancarlo Ruocco  Center for Life Nano Science @Sapienza, Istituto
Italiano di Tecnologia, 295 Viale Regina Elena, I-00161, Roma, Italy
Dipartimento di Fisica, Universita’ di Roma “La Sapienza”, P’le Aldo Moro 5,
I-00185, Roma, Italy  Walter Schirmacher  Center for Life Nano Science
@Sapienza, Istituto Italiano di Tecnologia, 295 Viale Regina Elena, I-00161,
Roma, Italy  Institut für Physik, Staudinger Weg 7, Universität Mainz, D-55099
Mainz, Germany

###### Abstract

Since decades, the concept of vibrational density of states in glasses has
been mirrored in liquids by the instantaneous-normal-mode spectrum. While in
glasses instantaneous configurations correspond to minima of the potential-
energy hypersurface and all eigenvalues of the associated Hessian matrix are
therefore positive, in liquids this is no longer true, and modes corresponding
to both positive and negative eigenvalues exist. The instantaneous-normal-mode
spectrum has been numerically investigated in the past, and it has been
demonstrated to bring important information on the liquid dynamics and
transport properties. A systematic deeper theoretical understanding is now
needed. Heterogeneous-elasticity theory has proven to be particularly
successful in explaining many details of the low-frequency excitations in
glasses, ranging from the thoroughly studied boson peak, down to the more
elusive non-phononic excitations observed in numerical simulations at the
lowest frequencies. Here we present an extension of heterogeneous-elasticity
theory to the liquid state, and show that the outcome of the theory agrees
well to the results of extensive molecular-dynamics simulations of a model
liquid at different temperatures. We find that the spectrum of eigenvalues
$\rho(\lambda)$ has a sharp maximum close to (but not at) $\lambda=0$, and
decreases monotonically with $|\lambda|$ on both its stable and unstable side.
We show that the spectral shape strongly depends on temperature, being
symmetric at high temperatures and becoming rather asymmetric at low
temperatures, close to the dynamical critical temperature. Most importantly,
we demonstrate that the theory naturally reproduces a surprising phenomenon, a
zero-energy spectral singularity with a cusp-like character developing in the
vibrational spectra upon cooling. This feature, known from a few previous
numerical studies, has been generally overlooked in the past due to a
misleading representation of the data. We provide a thorough analysis of this
issue, based on both very accurate predictions of our theory, and
computational studies of model liquid systems with extended size.

##  I Introduction

Clarifying the properties of the potential-energy hypersurface (or potential-
energy landscape, PEL) and its connection with the (slow) dynamics and
transport properties in supercooled liquids is a hot topic since decades (see,
among many others, Seeley and Keyes (1989); Stratt (1995); Bembenek and Laird
(1995); Keyes (1997); Sciortino and Tartaglia (1997); Sastry  _et al._ (2001);
Taraskin and Elliott (2002); Debenedetti and Stillinger (2001); Cavagna
(2001); Broderix  _et al._ (2000); Angelani  _et al._ (2000a, b, c); Grigera
_et al._ (2002); Fabricius and Sariolo (2002); Ciliberti and Grigeranz (2003);
Sampoli  _et al._ (2003); Angelani  _et al._ (2003); Berthier and Garrahan
(2003); Doye and Wales (2002, 2003); Wales and Doye (2003); Doliwa and Heuer
(2003); Grigera (2006); Clapa  _et al._ (2012); Coslovich  _et al._ (2019);
Zhang  _et al._ (2019)). The idea traces back to Goldstein Goldstein (1969)
who described the PEL as an ensemble of basins of attraction, each one
associated with a local minimum (the inherent structure, IS). The dynamics was
thus described as the interplay of (quasi-)harmonic motion confined around
minima, and the jumps (more or less rare, depending on temperature, $T$)
between the basins of attraction associated with different minima.

This picture was further developed by Stillinger and Weber Stillinger and
Weber (1982); Stillinger (1995). They defined the basin of attraction of a
given IS (i. e., the set of points connected to the minimum through a
steepest-descent path) in such a way that the basins of attraction generate a
complete description of the configurational space. Next, they wrote the
partition function as $Z=\sum_{\\{b\\}}Z_{b}$, where $\\{b\\}$ indicates the
set of all the basins, and $Z_{b}$ is the partition function associated to the
dynamics within the basin $b$. As a consequence, the free energy of a liquid
could be expressed Stillinger and Weber (1982); Stillinger (1995) as the sum
of the free energy of a ”typical” basin, ${\cal F}_{b}(T)$, and a term
accounting for the multiplicity of basins accessible by the system, ${\cal
N}(T)=\exp(\Sigma(T)/k_{B})$, with $\Sigma(T)$ the configurational entropy.
According to these ideas, the thermodynamics of the liquids could therefore be
described in terms of the minima of the PEL.

A further step forward was proposed by Keyes Seeley and Keyes (1989); Keyes
(1997), who used the instantaneous normal modes (INM) to predict not only the
very short time dynamics within a basin, but also the long time diffusion
processes of crossing the basins’ borders. In the INM approach, the local
curvature of the PEL is calculated along equilibrated molecular-dynamics (MD)
trajectories. The idea is to recover the global shape of the landscape from a
collection of local curvatures, specifically those that are sampled during the
equilibrium dynamics. For systems evolving close to a local minimum, the local
curvatures will be all positive. In contrast, points near the borders
separating basins will have at least one negative curvature, most likely more
than one. Intuitively, one can expect that the larger the amount of negative
curvatures present, the easier the inter-basins jumps, thus inducing a larger
diffusion coefficient, $D$. Indeed, although the exact functional form is
system dependent, the INM approach suggests that $D$ is a monotonic function
of the fraction of unstable modes.

The connection between the properties of the INM and the diffusion coefficient
in supercooled liquids immediately points towards the physics of the
structural arrest, and the theories for the glass transition(see, among other
reviews, Berthier and Biroli (2011)). Among them, we mention the mode-coupling
theory (MCT) Götze (2009), a mean-field approach which predicts the existence
of a finite dynamical critical temperatures, $T_{MCT}$, where $D=0$. The study
of the INM reveals Angelani  _et al._ (2000b); Broderix  _et al._ (2000) that
near $T_{MCT}$ the fraction of unstable modes vanishes. One can therefore
conclude that at $T_{MCT}$ the system ceases to evolve near the borders of the
basins of attraction, starting to spend most of the time close to minima.

An analysis made by Kurchan and Laloux Kurchan and Laloux (1996) (KL) may help
to elucidate the above process. KL focused on the properties of the borders of
the basins of attraction in high-dimensional configuration spaces. They
concluded that in an $N$-dimensional space a basin has a volume at least of
the order $\ell^{N}$, where $\ell$ is the typical distance between
configurations pertaining to the same basin. This has as a consequence that in
the thermodynamic limit, $N\rightarrow\infty$, most of the volume associated
to the basins is located near their borders, away from the minima. KL then
hypothesized that, at high $T$, the system populates regions of borders
belonging to many basins. Near $T_{MCT}$, in contrast, it is mainly situated
in regions of borders of a few basins, which obviously encompass smaller
volumes compared to the former case.

Interestingly, the KL analysis points towards a hierarchical organization of
the configuration space. Indeed, let us assume that we were able to study the
potential energy $U$ restricted to the borders, $\partial b$, of the basins at
low $T$. This restricted PEL is a $(N-1)$-dimensional manifold where all the
stationary points have their order decreased by one. As a consequence, saddles
of order one of the original PEL become local minima of the restricted PEL.
Similarly, at higher $T$ the trajectories will be localized close to the
$(N-2)$-dimensional borders of $(N-1)$-dimensional basins of attractions. One
can therefore iterate the procedure $k$ times until reaching the relevant
$(N-k)$-dimensional manifold in the vicinity of which equilibrium trajectories
are “almost always” within some small distance. This analysis overall
demonstrates the importance of saddle points (of increasing order on
increasing $T$) for the description of the liquid dynamics, and, specifically,
the importance of saddles of order one for the dynamics close to the arrest
temperature.

Following this line of thought, Angelani et al. Angelani  _et al._ (2000b) and
Broderix et al. Broderix  _et al._ (2000) studied the role played by saddles
in the supercooled liquid dynamics using a trick to find saddles of order one,
i.e., by searching for the stationary points of the function
$\Phi=|\vec{\nabla}{\cal V}|^{2}$, with ${\cal V}$ the system potential
energy. They found that the order of the saddles is a decreasing function of
the temperature and it vanishes at $T_{MCT}$. This numerical result pinpoints
the physical meaning of $T_{MCT}$: it is the temperature marking the cross-
over from a dynamics taking place near the basins’ borders, to a dynamics
mainly localized within the basins.

Despite a large body of literature focusing on the INMs, and on the pivotal
role they play in structural arrest, a robust theory describing the main
features of the density of the instantaneous normal modes (including a
surprising zero-energy spectral singularity at $T_{MCT}$, see below) is still
lacking. Here, following our recent work Schirmacher  _et al._ (2022), we
present a theoretical description of the INM spectra in liquids, accompanied
by a extensive numerical simulation of a model system. In our theory, a
generalization of the heterogeneous-elastic theory (HET) solved in the self-
Consistent Born approximation (SCBA) Schirmacher (2006); Schirmacher  _et al._
(2007); Marruzzo  _et al._ (2013); Schirmacher  _et al._ (2014), the
instantaneous-normal-mode spectrum of the liquid is described as that of an
elastic medium, where the local shear moduli exhibit strong spatial
fluctuations, including a large number of negative values. Realizing that
these fluctuations are merely snapshots of thermal fluctuations Mizuno  _et
al._ (2020); Mizuno and Mossa (2019), we assume that they obey Gaussian
statistics, with a variance proportional to the temperature, thus introducing
the temperature dependence of the INM in a very natural way. The outcome of
the theory is shown to favorably compare with the results of the extensive MD
simulations.

##  II Theory

We consider a system formed by $N$ particles with a potential energy ${\cal
V}(\\{{\bf r}_{i}(t)\\}_{i=1..N})$, which depends on the instantaneous
positions ${\bf r}_{i}(t)$ of the particles at time $t$. We assume that it can
be expressed in terms of pairwise potentials $\phi(r_{ij}(t))$ as,

${\cal V}(\\{{\bf r}_{i}(t)\\}_{i=1..N})=\frac{1}{2}\sum_{i\neq
j}\phi(r_{ij}(t)),$ (1)

with $r_{ij}(t)=|{\bf r}_{i}(t)-{\bf r}_{j}(t)|$. In the INM procedure one
considers virtual small displacements, ${\bf u}_{i}(t)$, around the
instantaneous particle positions, ${\bf r}_{i}(t)$, obtained in a MD
simulation. The second-order coefficients of a Taylor expansion of ${\cal V}$
form the Hessian matrix,

$H_{ij}^{\alpha\beta}(t)=\frac{\partial^{2}}{\partial r_{i}^{\alpha}\partial
r_{j}^{\beta}}{\cal V}(\\{{\bf r}_{1}(t)\dots\bf r_{N}(t)\\}),$ (2)

($\alpha,\beta=1,\ldots,3$), which can be diagonalized to obtain the
eigenvalues $\lambda_{p}$, with $p=1,\dots,3(N-1)$ 1113 eigenvalues are
trivially zero, due to the conditions $\sum_{i}H_{ij}^{\alpha\beta}=0$.,
together with the associated spectrum, $\rho(\lambda_{p})$.

If the system is close enough to a minimum, all the curvatures are positive
and so are all the $\lambda_{p}$. These are related to the square of the
vibrational frequencies, $\lambda_{p}=\omega_{p}^{2}$, and one can consider
the Hessian as the counterpart of the dynamical matrix of a solid. In this
regime one can therefore identify the spectrum with the density of vibrational
states (DOS) of the system,

$g(\omega)=2\,\omega\rho(\lambda(\omega))=\frac{\textstyle 1}{\textstyle
3N-3}\sum_{p=1}^{3N-3}\delta(\omega-\omega_{p})\,.$ (3)

If, in contrast, the system is away from a minimum, one or more negative
curvatures exist, and the non-trivial eigenvalues of the Hessian are both
positive and negative. In this case the identification of $\lambda_{p}$ with
the square of a ”frequency” is obviously no longer possible. For the unstable
part of the spectrum ($\lambda_{p}\leq 0$) it has become common practice to
define the positive quantity $\tilde{\omega}=-i\sqrt{\lambda}$, together with
the corresponding “DOS” as,

$g(\tilde{\omega})\equiv 2\,\tilde{\omega}\rho(\lambda(\tilde{\omega})).$ (4)

This function is usually plotted on the negative $\omega$-axis, and the stable
$g(\omega)$, according to Eq. (3), on the positive $\omega$-axis.

This representation, although commonly used as an alternative to directly
looking at the $\rho(\lambda)$, is misleading as unfortunately hides many
crucial details of the INM spectrum. Indeed, the former is obtained from the
latter by multiplication by $|\omega|$, a procedure, which strongly suppresses
all low-frequency details of the spectrum $\rho(\lambda)$ Schirmacher  _et
al._ (2022). Based on this observation, Sastry et al. Sastry  _et al._ (2001),
and Taraskin and Elliott Taraskin and Elliott (2002) showed that the INM
spectrum of liquids, when plotted directly as $\rho(\lambda)$ exhibits a
characteristic cusp-like maximum at $\lambda=0$, which cannot be seen in an
analysis based on Eqs. (3) and (4).

In order to obtain a theoretical description of the spectrum $\rho(\lambda)$,
we start from the equations of motion of linear elasticity with a spatially
fluctuating shear modulus, $G({\bf r})$, similar to HET in glasses. HET has
been rather successful in explaining several low-frequeny anomalies in the
vibrational spectrum of glasses Schirmacher (2006); Schirmacher  _et al._
(2007); Marruzzo  _et al._ (2013); Schirmacher  _et al._ (2014); Pan  _et al._
(2021), including the so-called Boson peak, an enhancement of the DOS at
finite frequencies over the Debye prediction $g(\omega)\propto\omega^{2}$. We
now aim to explore the application of this theory in the unstable regime, and
compare the results with the INM spectrum of a simulated soft-sphere liquid
(see below).

We consider the equations of motion for the continuum version of the virtual
displacements, ${\bf u}_{i}(t)\rightarrow{\bf u}({\bf r},t)$. The time
derivative of this quantity, ${\bf v}({\bf r},t)=\dot{\bf u}({\bf r},t)$, is
the Eulerian instantaneous local velocity field. The HET equation of motion
is,

$\rho\ddot{\bf u}({\bf r},t)=\nabla M({\bf r})\nabla\cdot{\bf u}({\bf
r})-\nabla\times G({\bf r})\nabla\times{\bf u}({\bf r}),$ (5)

where $M({\bf r})=K_{0}+\frac{4}{3}G({\bf r})$ is the longitudinal modulus.
Here, $K_{0}$ is the bulk modulus, which is considered not to exhibit relevant
spatial fluctuations. In contrast, the fluctuating shear modulus is written as
$G({\bf r})=G_{0}-\Delta({\bf r})$, where $G_{0}$ is the spatial average of
$G({\bf r})$, and $\Delta({\bf r})$ is a local deviation from the average,
which we assume to follow a normal distribution with a standard deviation
$\sigma^{2}$ given by

$\sigma^{2}=\langle\Delta({\bf r})^{2}\rangle$ (6)

In the HET theory for a glass Schirmacher (2006); Schirmacher  _et al._
(2007); Marruzzo  _et al._ (2013); Schirmacher  _et al._ (2014); Pan  _et al._
(2021), the standard deviation $\sigma$ is, of course, temperature
independent. In the present treatment for a liquid, in contrast, one expects a
temperature dependent $\sigma$, as the features of the visited regions of the
PEL depend on $T$. By comparing with the simulation data we shall see that
$\sigma^{2}$ is proportional to the temperature, a feature already observed
earlier in a simulation of a glass-forming liquid at higher temperatures
Mizuno and Mossa (2019).

To proceed further, we decompose the displacements into the longitudinal
($\nabla\times{\bf u}_{L}=0$) and transverse ($\nabla\cdot{\bf u}_{T}=0$)
components, and transform to frequency space, leading to the stochastic
Helmholtz equations,

$\omega^{2}{\bf u}_{L,T}({\bf r},\omega)=\nabla v_{L,T}^{2}({\bf
r})\nabla^{2}{\bf u}_{L,T}({\bf r},\omega).$ (7)

Here, $v_{L}^{2}({\bf r})=M({\bf r})/{\rho}$ and $v_{T}({\bf r})^{2}=G({\bf
r})/{\rho}$ are the locally fluctuating squared sound velocities, longitudinal
(L) and transverse (T), respectively, and $\rho$ is the mass density.
Considering Gaussian fluctuations for $\Delta({\bf r})$, we can use the
replica trick to calculate the average spectrum McKane and Stone (1981); John
_et al._ (1983); Schirmacher (2006).

###  II.0.1 The Self-Consistent Born Approximation (SCBA)

From a saddle-point approximation of the resulting effective action, we obtain
an effective-medium theory where the fluctuating quantity $\Delta({\bf r})$ is
replaced by a complex frequency-dependent self-energy, $\Sigma(z)$, with
$z=\omega^{2}+i\epsilon=\lambda+i\epsilon$, for $\epsilon\rightarrow 0_{+}$.
$\Sigma(z)$ obeys the self-consistent SCBA Eq. (10) Schirmacher (2006);
Schirmacher  _et al._ (2007); Marruzzo  _et al._ (2013); Schirmacher  _et al._
(2014), see below.

In three dimensions, the longitudinal and transverse Green’s functions for
Eqs. (7) are,

${\cal G}_{L,T}(k,z)=\frac{1}{-z+k^{2}v_{L,T}(z)^{2}},$ (8)

where the $v_{L,T}(z)$ are the effective complex frequency-dependent sound
velocities,

$\displaystyle v^{2}_{L}(z)$ $\displaystyle=$
$\displaystyle\frac{1}{\rho}\big{[}K_{0}+\frac{4}{3}\big{(}G_{0}-\Sigma(z)\big{)}\big{]}$
(9) $\displaystyle v^{2}_{T}(z)$ $\displaystyle=$
$\displaystyle\frac{1}{\rho}\big{[}G_{0}-\Sigma(z)\big{)}\big{]}.$

We choose $G_{o}$ as the reference shear modulus $G_{0}$, which is the average
shear modulus in the stable glass. $G_{0}$ might be approximately correspond
to $G_{\infty}$, the liquid high-frequency shear modulus, which appears in the
Maxwell relation $\tau=\eta/G_{\infty}$ between the relaxation time $\tau$ and
the shear viscosity $\eta$ Hertzfeld and Litovitz (1959). The self-consistent
HET-SCBA equation for the self energy $\Sigma(z)$ is,

$\displaystyle\Sigma(z)$ $\displaystyle=$
$\displaystyle\sigma^{2}V_{c}\,{\textstyle\left(\frac{1}{2\pi}\right)^{3}}\;\frac{1}{\rho}$
(10) $\displaystyle\times$ $\displaystyle\int_{|{\bf k}|\leq k_{\xi}}d^{3}{\bf
k}\;k^{2}\bigg{(}\frac{2}{3}{\cal G}_{L}(k,z)+{\cal G}_{T}(k,z)\bigg{)},$

where $V_{c}$ is the coarse-graining volume, and the upper cutoff, $k_{\xi}$,
is inversely proportional to the correlation length, $\xi$, of the
fluctuations. We now introduce the dimensionless variables, $q$=$k/k_{\xi}$,
$\widetilde{\Sigma}(z)$=$\Sigma(z)/G_{0}$, $\tilde{K}_{0}$=$K_{0}/G_{0}$,
$\tilde{M}(\lambda)$=$\tilde{K}+\frac{4}{3}[1-\widetilde{\Sigma}(\lambda)]$,
and the (also dimensionless) disorder parameter,

$\gamma=\frac{1}{G_{0}^{2}}\sigma^{2}V_{c}k_{\xi}^{3}/6\pi^{2}\,.$ (11)

In terms of these quantities, Eq. (10) takes the form,

$\widetilde{\Sigma}(z)=\;\gamma\int_{0}^{1}dq\;q^{4}\bigg{(}2\;\tilde{\cal
G}_{L}(q,z)+3\;\tilde{\cal G}_{T}(q,z)\bigg{)},$ (12)

where,

$\tilde{\cal
G}_{L,T}(q,z)=\frac{1}{-\frac{z\rho}{G_{0}k^{2}_{\xi}}+\frac{q^{2}\rho}{G_{0}}v_{L,T}(z)^{2}}.$
(13)

By further defining the dimensionless squared frequency,
$\tilde{z}$=$z\rho/G_{0}k_{\xi}^{2}$, and, correspondingly,
$\tilde{\lambda}$=$\lambda\rho/G_{0}k_{\xi}^{2}$, we can write,

$\widetilde{\Sigma}(\tilde{z})=\;\gamma\int_{0}^{1}dq\;q^{4}\bigg{(}2\;\tilde{\cal
G}_{L}(q,\tilde{z})+3\;\tilde{\cal G}_{T}(q,\tilde{z})\bigg{)},$ (14)

with

$\displaystyle\tilde{\cal G}_{L}(q,\tilde{z})$ $\displaystyle=$
$\displaystyle\frac{1}{-\tilde{z}+q^{2}\big{[}\tilde{K}_{0}+\frac{4}{3}\big{(}1-\widetilde{\Sigma}(\tilde{z})\big{)}\big{]}}$
(15) $\displaystyle\tilde{\cal G}_{T}(q,\tilde{z})$ $\displaystyle=$
$\displaystyle\frac{1}{-\tilde{z}+q^{2}\big{[}1-\widetilde{\Sigma}(\tilde{z})\big{]}}.$

The self-consistent Eqs. (14) and (15), must be solved numerically for the
self energy, $\widetilde{\Sigma}(\tilde{z})$, and the Green’s functions,
$\tilde{\cal G}_{L,T}(q,\tilde{z})$. The latter can be used in turn to
calculate the dimensionless density of eigenvalues as,

$\tilde{\rho}(\tilde{\lambda})=\frac{1}{3\pi}\Im{\Big{\\{}{\cal\tilde{G}}_{L}(\tilde{z})+2{\cal\tilde{G}}_{T}(\tilde{z})\Big{\\}}},$
(16)

with the local Green’s functions,

${\cal\tilde{G}}_{L,T}(\tilde{z})=3\frac{k_{\xi}^{3}}{k_{D}^{3}}\int_{0}^{k_{D}/k_{\xi}}dqq^{2}{\cal\tilde{G}}_{L,T}(q,\tilde{z}),$
(17)

and the Debye cutoff, $k_{D}=\sqrt[3]{6\pi^{2}N/V}$. In what follows, for the
sake of simplicity, we set $k_{\xi}=k_{D}$, obtaining,

${\cal\tilde{G}}_{L,T}(\tilde{z})=3\int_{0}^{1}dq\;q^{2}{\cal\tilde{G}}_{L,T}(q,\tilde{z}).$
(18)

The numerical solution of the self-consistent equations provides the density
of eigenvalues shown in Fig. 1, for the indicated values of the parameter
$\gamma$ (Eq. (11)). We chose these values to correspond to the chosen
temperatures in the simulation (Fig. 4), according to $\gamma 3.5T$, where $T$
is given in Lennard-Jones (LJ) units.

Figure 1: INM spectra, $\rho(\lambda)$, of our unstable-elasticity model,
calculated in the self-consistent Born approximation (SCBA) as described in
the main text. We have considered different values of the disorder parameter
$\gamma$, all of them in the liquid side ($\gamma>\gamma_{c}$=1/4), as
indicated in the legend. For better comparison with the numerical data of Fig.
4 we converted the eigenvalues to LJ units according to
$\lambda=40\,\tilde{\lambda}$.

###  II.0.2 The transverse contribution

A simple but robust approximation, at least in the low-frequency region,
consists in only considering the transverse contribution to the DOS. In this
case, the HET-SCBA equation, Eq. (12) becomes,

$\displaystyle\widetilde{\Sigma}(\tilde{z})$ $\displaystyle=$
$\displaystyle\;\gamma\int_{0}^{1}dq\;q^{4}\frac{3}{-\tilde{z}+q^{2}\big{[}1-\widetilde{\Sigma}(\tilde{z})\big{]}}=$
$\displaystyle=$
$\displaystyle\;\gamma\frac{1}{[1-\widetilde{\Sigma}(\tilde{z})]}\int_{0}^{1}dq\;3q^{2}\frac{q^{2}}{q^{2}-\frac{\tilde{z}}{[1-\widetilde{\Sigma}(\tilde{z})]}},$

which, defining the auxiliary function
$f(\tilde{z})\doteq\tilde{z}/[1-\widetilde{\Sigma}(\tilde{z})]$, can be also
expressed as,

$\displaystyle\widetilde{\Sigma}(\tilde{z})(1-\widetilde{\Sigma}(\tilde{z}))=\gamma\int_{0}^{1}dq\;3q^{2}\frac{q^{2}}{q^{2}-f(\tilde{z})}=$
$\displaystyle=$
$\displaystyle\;\gamma\left\\{1+3f(\tilde{z})\left[1+f(\tilde{z})^{1/2}\frac{1}{2}\ln\Big{(}\frac{f(\tilde{z})^{1/2}-1}{f(\tilde{z})^{1/2}+1}\Big{)}\right]\right\\}.$

It is now convenient to write Eq. (II.0.2) as,

$\widetilde{\Sigma}(\tilde{z})(1-\widetilde{\Sigma}(\tilde{z}))=\gamma+\tilde{z}{\cal{F}}(\tilde{z}),$
(21)

where we have defined,

$\displaystyle{\cal{F}}(\tilde{z})$ $\displaystyle\doteq$
$\displaystyle\frac{3\gamma}{1-\widetilde{\Sigma}(\tilde{z})}\left[1+f(\tilde{z})^{1/2}\frac{1}{2}\ln\Big{(}\frac{f(\tilde{z})^{1/2}-1}{f(\tilde{z})^{1/2}+1}\Big{)}\right]=$
(22) $\displaystyle=$ $\displaystyle\gamma\;{\cal\tilde{G}}_{T}(\tilde{z}).$

As the function ${\cal\tilde{G}}_{T}(\tilde{z})$ reaches a finite value for
$|\tilde{z}|\rightarrow 0$, Eq. (21) makes clear that the low-frequency limit
of the self-energy is also finite. We may now solve the quadratic Eq. (21),
and obtain 222We take the solution with the plus sign which, if one goes into
the glass with $\gamma$ very small, recovers the Born approximation result
Ganter and Schirmacher (2010).,

$\widetilde{\Sigma}(\tilde{z})=\frac{1}{2}+\sqrt{\gamma_{c}-\gamma-\tilde{z}{\cal{F}}(\tilde{z})},$
(23)

with $\gamma_{c}=\frac{1}{4}$. The unstable (i.e., liquid) situation is
characterized by $\gamma>\gamma_{c}$. In this case it is convenient re-write
Eq. (23) as,

$\widetilde{\Sigma}(\tilde{z})=\frac{1}{2}+i\sqrt{\gamma-\gamma_{c}+\tilde{z}{\cal{F}}(\tilde{z})}.$
(24)

The self consistent equations for the transverse contribution to the level
density are, therefore, Eqs. (22) and (24).

###  II.0.3 The low-frequency limit

It is worth noting that the self energy $\widetilde{\Sigma}(\tilde{z})$
satisfies the following properties for vanishing argument,

$\displaystyle\Sigma_{0}$ $\displaystyle=$
$\displaystyle\frac{1}{2}+i\sqrt{\gamma-\gamma_{c}}$
$\displaystyle\Sigma_{0}^{*}$ $\displaystyle=$ $\displaystyle 1-\Sigma_{0}$
(25) $\displaystyle|\Sigma_{0}|^{2}$ $\displaystyle=$
$\displaystyle\Sigma_{0}\big{(}1-\Sigma_{0}\big{)}=\gamma,$

where we have denoted, $\Sigma_{0}\doteq\widetilde{\Sigma}(0)$). By defining
$\xi=\xi^{\prime}+i\xi^{\prime\prime}=\Sigma_{0}^{1/2}$, we next obtain,

$\displaystyle\xi^{\prime}$ $\displaystyle=$
$\displaystyle\frac{1}{2}\sqrt{\sqrt{1+4(\gamma-\gamma_{c})}+1}$ (26)
$\displaystyle\xi^{\prime\prime}$ $\displaystyle=$
$\displaystyle\frac{1}{2}\sqrt{\sqrt{1+4(\gamma-\gamma_{c})}-1},$

and by combining both equations,

$\xi^{\prime\prime}=\frac{1}{\xi^{\prime}}\Sigma_{0}^{\prime\prime}=\frac{1}{\xi^{\prime}}\sqrt{\gamma-\gamma_{c}}\,.$
(27)

To finalize the list of properties of the self energy at small frequency, we
observe that, besides the constant $\Sigma_{0}$, the leading term is linear,
as one can conclude from the relation,

$\widetilde{\Sigma}(\tilde{z})=\Sigma_{0}\Big{(}1+i\frac{3}{2}\frac{1}{\sqrt{\gamma-\gamma_{c}}}\tilde{z}\Big{)}+\mathcal{O}(\tilde{z}^{3/2}).$
(28)

The density of eigenvalues $\lambda$, considering the transverse contributions
only, is therefore given by,

$\tilde{\rho}(\tilde{\lambda})=\frac{2}{3\pi}\Im{\Big{\\{}{\cal
G}_{T}(\tilde{z})\Big{\\}}}=\frac{2}{3\pi\gamma}\Im{\Big{\\{}{\cal
F}_{T}(\tilde{z})\Big{\\}}}.$ (29)

We finally analyze the imaginary part of ${\cal{F}}(\lambda)$ in Eq. (22) at
low frequency. For small $\lambda$, the argument of the logarithm becomes -1,
therefore its real part vanishes while the imaginary part become constant and
equal to $\pi$,

$\displaystyle{\cal{F}}^{\prime\prime}(\lambda)$ $\displaystyle\approx$
$\displaystyle\Im{\Big{\\{}\frac{3\gamma}{[1-\widetilde{\Sigma}(\lambda)]}\Big{(}1+i\frac{\pi}{2}f(\lambda)^{1/2}\Big{)}\Big{\\}}}=$
$\displaystyle=$
$\displaystyle\Im{\Big{\\{}\frac{3\gamma}{[1-\widetilde{\Sigma}(\lambda)]}\Big{(}1+i\frac{\pi}{2}\frac{\lambda^{1/2}}{[1-\widetilde{\Sigma}(\lambda)]^{1/2}}\Big{)}\Big{\\}}}.$

###  II.0.4 The zero-frequency value and the cusp

It is now straightforward to obtain the density of eigenvalues at $\lambda=0$
as,

$\tilde{\rho}_{o}\doteq\tilde{\rho}(\tilde{\lambda}\\!=\\!0)=\frac{2}{\pi}\Im\bigg{\\{}\frac{1}{1-\Sigma_{o}}\bigg{\\}}=\frac{1}{\pi}\frac{\sqrt{4\gamma-1}}{\gamma}.$
(31)

Interestingly, this value is not approached linearly by our theory, but rather
in a non-analytic way, producing a characteristic cusp already observed in
Taraskin and Elliott (2002).

Figure 2: Prefactor of the $\lambda^{1/2}$ term in the density of level
$\rho(\lambda)$ (Eq. (35)) for positive (blue line) and negative (red line)
values of $\lambda$. Only the contributions of the transverse modes to the
level distribution are considered here. The prefactor is always negative (thus
the cusp is upward) in the negative $\lambda$ region. For $\lambda>0$ the
direction of the cusp depends on the value of $\gamma$. At $\gamma$ close to
$\gamma_{c}$(=1/4), corresponding to low temperatures, the cusp is downward,
on increasing $\gamma$ (in our approximation, at $\gamma$=1) the cusp changes
direction.

Indeed, beside the constant value reached at $\tilde{\lambda}$=0, the leading
term of the function ${\cal{F}}(\tilde{\lambda})$ is proportional to
$\tilde{\lambda}^{1/2}$. In the low frequency expansion of Eq. (II.0.3) we can
therefore safely substitute $\widetilde{\Sigma}(\tilde{\lambda})$ with
$\Sigma_{0}$, obtaining,

$\displaystyle{\cal{F}}^{\prime\prime}(\lambda)$ $\displaystyle\approx$
$\displaystyle\Im{\Big{\\{}\frac{3\gamma}{[1-\Sigma_{0}]}\Big{(}1+i\frac{\pi}{2}\frac{\tilde{\lambda}^{1/2}}{[1-\Sigma_{0}]^{1/2}}\Big{)}\Big{\\}}}$
$\displaystyle=$
$\displaystyle\Im{\Big{\\{}3\Sigma_{0}\Big{(}1+i\frac{\pi}{2}{\Sigma_{0}^{1/2}}\frac{1}{\sqrt{\gamma}}{\tilde{\lambda}^{1/2}}\Big{)}\Big{\\}}}$
$\displaystyle=$
$\displaystyle\Im{\Big{\\{}3(\Sigma^{\prime}_{0}+i\Sigma^{\prime\prime}_{0})\Big{(}1+i\frac{\pi}{2}(\xi^{\prime}+i\xi^{\prime\prime})\frac{1}{\sqrt{\gamma}}{\tilde{\lambda}^{1/2}}\Big{)}\Big{\\}}}.$

We see that, for $\tilde{\lambda}$ small but $>0$,

${\cal{F}}^{\prime\prime}(\tilde{\lambda})=3\Sigma^{\prime\prime}_{0}+\tilde{\lambda}^{1/2}\frac{3\pi}{2\sqrt{\gamma}}\Big{\\{}\xi^{\prime}\Sigma_{0}^{\prime}-\xi^{\prime\prime}\Sigma_{0}^{\prime\prime}\Big{\\}},$
(33)

while, for negative $\tilde{\lambda}$ we find,

${\cal{F}}^{\prime\prime}(\tilde{\lambda})=3\Sigma^{\prime\prime}_{0}+|\tilde{\lambda}|^{1/2}\frac{3\pi}{2\sqrt{\gamma}}\Big{\\{}-\xi^{\prime}\Sigma_{0}^{\prime\prime}-\xi^{\prime\prime}\Sigma_{0}^{\prime}\Big{\\}}.$
(34)

Summing up, for $\gamma_{c}$=1/4 we can write,

$\tilde{\rho}(\tilde{\lambda})=\frac{2}{3\pi\gamma}{\cal{F}}^{\prime\prime}(\tilde{\lambda})=\tilde{\rho}_{o}+h_{\pm}(\gamma)\;|\tilde{\lambda}|^{1/2},$
(35)

where the function $h_{+}(\gamma)$ ($h_{-}(\gamma)$) must be used for the
positive (negative) $\lambda$-regions, and are defined as,

$h_{\pm}(\gamma)=\frac{\gamma^{-3/2}}{4}\left[-\sqrt{(4\gamma-1)(2\gamma^{1/2}\mp
1)}\pm\sqrt{(2\gamma^{1/2}\pm 1)}\right].$ (36)

A non-analytical behavior of the density of level is therefore clear in Eq.
(35), with an asymmetric square-root cusp centered at $\lambda=0$. The
functions $h_{\pm}(\gamma)$, determining both direction and sharpness of the
cusp, are shown in Fig. 2. We see that for negative values of $\lambda$, we
always obtain a upward $\lambda^{1/2}$ cusp. For $\lambda>0$, in contrast, we
find a positive rise proportional to $\lambda^{1/2}$ in the region where
$\gamma$ is small (i.e., at low temperatures), whereas for large $\gamma$ an
upwards $\lambda^{1/2}$ cusp is obtained. The transition between these two
regimes occurs at $\gamma=1$.

The validity of the low-frequency approximation of Eqs. (35) and (36) is
demonstrated in Fig. 3, where the ”exact” solution of the self-consistent Eqs.
(24) and (22)(full line) is compared with the approximated one (dashed line),
for the two indicated values of $\gamma$. We conclude that the description of
the low frequency level density with a cusp is valid and robust.

Figure 3: Test of the low frequency approximation for the function
${\cal{F}}^{\prime\prime}(\lambda)$. The dashed lines are approximate solution
from Eqs. (35) and (36), while the full lines are obtained by a numerical
solution of the self-consistent Eqs. (24) and (22). The red lines refer to
$\gamma$=0.5, the green to $\gamma$=2.

The origin of the non-analytic singular behavior discussed above can be traced
back to the sum rule $\sum_{j}H_{ij}^{\alpha\beta}=0$ for the Hessian Taraskin
and Elliott (2002), due to the global translation invariance. In our elastic
model the sum rule is transformed to the double spatial derivative in the
elastic wave Eq. (5), or in the Helmholtz Eqs. (7), again reflecting global
translation invariance. Technically, the $\lambda^{1/2}$ singularity stems
from vanishing values of the wave-number $k$ in the integrand of the SCBA Eq.
(10), corresponding to density and stress fluctuations of very large extent.
We convinced ourselves that the singularity at $\lambda=0$ would disappear if
an infrared cutoff $k_{0}$ is introduced at the lower limit of the integral in
Eq. (10), which mimics the effect of a finite sample size. As a consequence,
in numerical simulations at finite system sizes, one expects a finite slope at
$\lambda=0$, at variance with the $\lambda=0$ singularity of the SCBA. We will
address this point below.

Concluding the presentation of the theoretical results, and specifically of
the existence of a cusp in the vibrational levels spectrum, we emphasize that
this result is an additional evidence of the intrinsic weakness of exceedingly
crude approximations sometimes employed in the literature, including the
prescription $\rho(\lambda)\sim\textit{const}$ (or, equivalently,
$g(\omega)\sim|\omega|$) Zaccone and Baggioli (2021); Stamper  _et al._
(2022). Indeed, in considering this approximation many details of the liquid
dynamics are plainly ignored.

##  III Numerical methods

In the following Sections we will compare the above theoretical predictions to
the outcome of extensive computer simulations. Here, we give details about the
atomic liquid model used, the molecular-dynamics (MD) simulations performed,
and the numerical tools employed in the INM analysis. Below, all observables
are expressed in Lennard Jones units by setting the fundamental quantities,
$m$, $\sigma$, $\epsilon$ (mass, length, and energy, respectively), and the
Boltzmann constant $k_{B}$ to 1, while $\sigma\sqrt{m/\epsilon}$ is the unit
of frequency.

We have considered systems formed by equimolar binary mixtures of $N$ soft-
spheres of type $A$ (small) and $B$ (large), respectively, with
$N_{A}=N_{B}=N/2$. The spheres interact via the pair-wise potential,

$V(r)/\epsilon=\sigma^{12}\left(\frac{1}{r^{12}}-\frac{1}{r_{c}^{12}}\right)\times\left[1-\left(\frac{r}{r_{c}}\right)^{2}\right]^{2},$
(37)

which is the product of the usual soft sphere potential ($\propto r^{-12}$),
cutoff and shifted at $r_{c}=1.5$, and a tapering function that ensures
continuity of both potential and forces at the cutoff. Here $r=r_{ij}$ is the
distance between the particles $i$ and $j$,
$\epsilon_{AA}=\epsilon_{BB}=\epsilon_{AB}=1$, $\sigma_{AA}=0.90$,
$\sigma_{BB}=1.09\simeq 1.2\,\sigma_{AA}$, and
$\sigma_{AB}=(\sigma_{AA}+\sigma_{BB})/2\simeq 1$. Masses are
$m_{A}=m_{B}=m=1$. All simulations have been conducted by using the high-
performance-computing simulation code LAMMPS Thompson  _et al._ (2022).

We have considered 7 system sizes, $N=10^{3}\times 2^{n}$ with $n=0,\ldots,7$
($n=6$ excluded), ranging from the situation where the entire set of the
Hessian eigenvalues can be calculated comfortably (small $N$), to that where
we are only able to obtain the vibrational density of state in a limited
frequency range around $\lambda\simeq 0$ (large $N$), although with high
statistical accuracy. We have fixed the simulation box sizes, $L$, by matching
the number density $\hat{\rho}=N/L^{3}=1$, for all systems. We considered 24
values of temperature, $T$, in the range $[0.08:1]$. Starting at $T=1$, the
system was slowly annealed to $T=0.08$ in steps of $\delta T=0.04$. At each
temperature, a total of $12\times 10^{6}$ time steps were performed,
consisting of a $2\times 10^{6}$ thermalization period, followed by the
production run with the dumping of the system configurations used in the
subsequent analysis. The runs were carried out in the $(NVE)$ ensemble with a
numerical integration time step $\delta t=0.002$.

We have implemented the calculation of the Hessian matrix, Eq. (2), in JAX
Bradbury  _et al._ (2018), which can automatically differentiate native Python
and NumPy functions. $H^{\alpha\beta}_{ij}$ has been evaluated at each value
of $T$ and $N$ on a number of independent istantaneous configurations
sufficient to reach the desired statistical accuracy, and subsequently
diagonalized to extract the eigenvalues. We have employed the standard
scipy.linalg Virtanen  _et al._ (2020) linear algebra functions at small
values of $N$, while for $N\geq 16000$ we used RALEIGH Ovtchinnikov (2008), a
Python implementation of the block Jacobi-conjugated gradients algorithm for
computing eigenvalues and eigenvectors of large-scale real symmetric and
Hermitian problems in selected eigenvalue ranges. We have eventually estimated
the INM spectra, $\rho(\lambda)$, by constructing histograms of the obtained
eigenvalues.

Figure 4: Simulated INM spectra, $\rho(\lambda)$ , for $N=4\times 10^{3}$ at
the indicated values of $T$. For the investigated soft sphere system,
$T_{MCT}\simeq 0.05$. All quantities are expressed in LJ units. These data are
the numerical counterpart of those shown in Fig. 1. The number of unstable
modes decreases significantly upon cooling and, as a consequence,
distributions which are quite asymmetric around $\lambda=0$ at high $T$,
becomes consistently increasingly asymmetric by decreasing $T$, as expected.
These data are discussed at length in the main text. Figure 5: Effect of the
finite system size on the low-$\lambda$ INM spectral shape. Simulated INM
spectra, $\rho(\lambda)$, at the lowest investigated temperature, $T=0.08$,
for the indicated values of $N=10^{3}\times 2^{n}$, with $n=0,\ldots,7$ ($n=6$
excluded). At the smaller values of $N$ the $\rho(\lambda)$ intensities
increase continuously on the unstable side, crossing $\lambda=0$ and
eventually reaching maxima at positive values of $\lambda$. The positions of
the maxima, however, continuously shift toward $\lambda=0$ on increasing $N$.
At the highest investigated values of $N$ the cusp-like non-analytical
behavior starts to develop, as discussed into details in the main text.

##  IV Numerical Results

We now present a comparison of the simulated spectra of the soft-sphere liquid
with the predictions obtained by the unstable version of HET-SCBA theory,
together with an in-depth discussion of the data especially in light of the
discussion of Sect. II.0.4. In particular, we focus on the finite-size effects
on the distribution of the eigenvalues of the Hessian matrix, and on the
details of the shape of $\rho(\lambda)$ in different conditions, both in the
limit $\lambda\rightarrow 0$ and for intermediate values of $\lambda$.

###  IV.0.1 Overall INM spectrum

We show in Fig. 4 the calculated INM spectra of the investigated soft-spheres
liquids, for the indicated values of $T$, for $N=$ 4000\. The lowest
considered temperatures are just above the arrest temperature (we estimate,
$T_{MCT}\simeq$ 0.05), while the highest values of $T$ encompass the high-
fluidity regime, up to $T\simeq 20\;T_{MCT}$. We observe that the spectrum
changes from an almost stable situation with a very small fraction of negative
eigenvalues at low $T$ (black symbols), to a distribution which is almost
symmetric with respect to positive and negative $\lambda$’s. These
$T$-dependent data are very similar to those predicted by the HET-SCBA shown
in Fig. 1 for the indicated values of the disorder parameter, $\gamma$. As
said above, we chose these values according to $\gamma=3.5T$. While a more
detailed comparison of the numerical and analytical $\rho(\lambda)$ curves
would give a slight deviation from this linear law, the overall trend is given
correctly. We already mentioned that in our earlier simulation Mizuno and
Mossa (2019), in which the temperature dependence of the variance
$\sigma^{3}\propto\gamma$ was evaluated, an increase with temperature was
observed at the highest considered temperatures.

We turn now to a comparison between the simulated $\rho(\lambda)$ curves (Fig.
4), and those predicted by HET-SCBA (Fig. 1). The only notable difference
between the two sets of data is that in the spectra of Fig. 1 the
negative-$\lambda$ tails are truncated, a typical effect of the SCBA
approximation, which provides spectra extending on a support of finite width.
No cutoff is visible in the simulation data, as expected.

HET-SCBA therefore seems to satisfactorily reproduce the numerical INM spectra
modifications with $T$, including the transition from a stable situation at
$T\simeq T_{MCT}$ to a quite symmetric distribution at high $T$, while, as
expected due to the normalization, the peaks heights decrease with $T$. We
also note that the maximum of the distribution is always found at
$\lambda_{M}>0$, in contrast to the intuitive expectation $\lambda_{M}=0$. We
will better discuss this point below.

Figure 6: Finite-size scaling of values and derivatives of the INM spectra at
$\lambda=0$. The data have been determined by fitting linearly the data of
Fig. 5 around $\lambda=0$. Top: Values of the INM spectra $\rho(\lambda=0)$ at
$T=0.08$, as a function of $N^{-1/2}$. Data are indicated by the symbols, the
solid line is a linear guide for the eyes. Bottom: Inverse of the derivative
of the INM spectra at $\lambda=0$, represented as a function of $N^{-1/3}$.
The details of the observed scaling with the system size are discussed in the
main text.

###  IV.0.2 Finite-size effects

We now quantify the effect due to the finite size of the simulated systems on
the $\rho(\lambda)$. In Fig. 5 we show the details of $\rho(\lambda)$ in the
very low-$|\lambda|$ region, for the indicated system sizes $N$, at the lowest
investigated temperature, $T=0.08$ (corresponding to the black symbols in Fig.
4). We observe a strong $N$-dependence of the spectral shape, with a slope for
$\lambda<0$ which consistently increases by increasing $N$. Importantly, the
theory predicts a derivative which diverges at $\lambda$=0 (see, for instance,
Eq. (35)), which cannot be directly observed in our simulation. Our data,
however, point in the correct direction, allowing us to qualitatively conclude
that our simulation is consistent with the theory in the thermodynamic limit,
$N\rightarrow\infty$.

To better emphasize this important point, we have extracted the main features
of $\rho(\lambda)$ close to $\lambda=0$ by a linear fit of the data of Fig. 5
restricted to the negative (unstable) $\lambda$-branch, and up to the position
of the distribution maximum, always localized at positive values of $\lambda$,
as already noted above. We show in Fig. 6 (top) the obtained values for
$\rho(\lambda=0)$, represented as a function of $N^{-1/2}$, and for the
inverse derivative $1/\rho^{\prime}(0)$ shown as a function of $N^{-1/3}$
(bottom).

In the $N\rightarrow\infty$ limit, we find that $\rho(0)$ extrapolates to a
finite value $\simeq 1.686\times 10^{-2}$ (see Eq. (31)), while
$\rho^{\prime}(0)$ diverges, as expected from the presence of the cusp, Eq.
(35), for $\gamma$ close to criticality. HET-SCBA therefore truly succeeds in
identifying quantitatively the zero-energy singularity in the simulated energy
spectrum. We now demonstrate that it also provides quantitative information
about the cusp-like character of the singularity. Here, the possibility of
analyzing quite large system sizes is crucial.

Figure 7: Details of the low frequency region of the INM spectra. We plot
$\rho(\lambda)$ for N=128000, at the four lowest values of $T$. All details of
the data, both in the unstable and stable $\lambda$-regions, are discussed in
the main text.

###  IV.0.3 The cusp at $\lambda$=0

We are now in the position to analyze quantitatively the shape of the
simulated $\rho(\lambda)$ curves at small values of $|\lambda|$ on both the
stable and unstable regions, and its dependence on temperature or,
equivalently, on the disorder parameter, $\gamma$. We show in Fig. 7 the
details of the $\rho(\lambda)$ at small $|\lambda|$ values for the largest
investigated system size $N$=128000, at the indicated low temperatures.

With reference to Eqs. (35) and (36), at the lowest $T=0.08$ (red symbols,
corresponding to $\gamma\approx\gamma_{c}$), we find that $h_{+}(\gamma)>0$,
while $h_{-}(\gamma)$, in contrast, is negative in all cases, as predicted by
Eq. (36). Also, on increasing $T$ ($\gamma$), the sign of $h_{+}(\gamma)$
reverses, as expected, indicating that already at $T$=0.12 we have $\gamma>1$.
As illustrated in Fig. 8, by plotting $\rho(\lambda)$ against
$\mathrm{sign}(\lambda)|\lambda|^{1/2}$, the curves of Fig. 7 are rectified
(according to Eq. (35)), and the slopes for $\lambda$ positive and negative
are different. This difference tends to vanish on increasing $T$ ($\gamma$),
where the distribution becomes smoother (see blue points in Fig. 8).

We conclude that the simulations performed on a soft sphere model confirm both
qualitatively and quantitatively the prediction of the INM made on the basis
of the HET extended to the liquid case. In particular, we find that the
critical behavior of the distributions $\rho(\lambda)$ in the vicinity of the
zero-energy spectral singularity is very satisfactorily described by the
theory. HET-SCBA seems indeed to adequately grasp the behavior of a realistic
model of an atomic fluid.

Figure 8: Same data of Fig. 7, now plotted as a function of
$\text{sgn}(\lambda)|\lambda|^{1/2}$. The consequences of this representation
are discussed in details in the main text. Figure 9: Local slopes of
$\rho(\lambda)$ close to $\lambda=0$, both in the stable and unstable regions.
Left, Main panel: Local slope for $N=4000$ as a function of temperature, $T$.
The slopes on the stable (squares) and unstable (circles) regions around
$\lambda=0$ have been determined by a simultaneous non-linear fit of the two
regions, as discussed in the text. In the unstable region the slope is always
positive, and consistently increases upon cooling. In the stable region, in
contrast, the slope is always negative, increasing (in absolute value) on
lowering $T$, going through a maximum at $T\simeq 0.2$ and eventually
decreasing.) Left, Inset: Same that in the main panel, for $N=128000$. The
behavior is similar but the decreased impact of finite-size effects indicates
that the slope of the stable region becomes positive at $T\simeq$ 0.1. Right:
Slopes differences $\Delta S=S_{\text{s}}-S_{\text{u}}$ for both values of
$N$, together with the $\propto T^{-2/3}$ dependence (dashed lines). All data
are discussed at length in the main text.

###  IV.0.4 The slope of $\rho(\lambda$) at intermediate frequency

An interesting detail of the spectra obtained in numerical simulations, which
cannot be easily extracted from analytical calculations, are the strength and
sign of the slope of the $\rho(\lambda)$ curves at small but finite
$|\lambda|$ (thus away from the cusp region), in both the stable and unstable
regions. We have estimated both slopes by a single non-linear fit to an
empirical function of $\lambda$, around $\lambda=$ 0 333We have used the model
$\rho(\lambda)=[e^{-\alpha\lambda}(c_{0}+c_{1}\lambda)+e^{\alpha\lambda}(c_{2}+c_{3}\lambda)]/(e^{-\alpha\lambda}+e^{\alpha\lambda})$.
The values of the slopes are directly provided by $c_{1}$ and $c_{3}$. We have
checked that fitting procedures based on different model functions, including
linear models, provide very similar results..

We plot in the left panel of Fig. 9 the estimated slopes as a function of $T$
for $N=$ 4000 (main panel) and 128000 (inset). In the unstable region
(circles) the slope is always positive, and consistently decreases upon
heating. In the stable region, in contrast, the slope is negative.
Interestingly, it increases (in absolute value) on cooling, goes through a
maximum at $T\simeq$ 0.2 and eventually decreases approaching $T_{MCT}$. (Very
similar results have been reported in Schirmacher  _et al._ (2022)). The data
for large $N$ in the inset provide substantially the same picture, with the
additional information that the slope of the stable region becomes positive at
$T\simeq$ 0.1, as already noted above.

Unfortunately we do not have an intuition of the mechanism behind the observed
behavior of the stable region slope, which we postpone to future work. As more
work is needed to extract from the theory a clear explanation of the data
shown in the right panel of Fig. 9. Here, we plot the difference $\Delta
S=S_{s}-S_{u}$ of the stable and unstable slopes, on a double-logarithmic
scale as a function of $T$ for both sets of data. As already observed by some
of us Schirmacher  _et al._ (2022), the data follow a $T^{-2/3}$ law, in the
entire investigated temperature range. A satisfactorily explanation of this
observation is lacking.

##  V Conclusions

Melting of a crystal is accompanied by the suppression of sharp features of
the vibrational spectrum. For instance, the Van Hove singularities in the
density of states, $g(\omega)$, (associated to the critical points of the
dispersion, $\omega(q)$) disappear, and are substituted by the appearance of a
broad distribution of vibrational modes, including both stable and unstable
excitations. Interestingly, the reverse process can be very different. This is
the case for systems that can be supercooled, i.e., where it is possible to
avoid crystallization by a succession of equilibrium liquid states extending
below the melting temperature. In this process, the high fraction of unstable
modes contained in the spectrum of the Hessian matrix of the liquid vanishes
at $T_{MCT}$, which is associated to a dynamical phase transition, and marks
the boundary where a saddle-dominated PEL liquid-like kinetics gives way to a
minima-dominated flow, eventually leading to glass formation. The limiting
shape of the vibrational spectrum approaching $T_{MCT}$ turns out to be
surprising, with the appearance of a zero-energy spectral singularity, now
unrelated to the $\omega(q)$, with a cusp-like character, and an origin which
can be traced back to the most obvious of the symmetries, translational
invariance.

Despite the generality of this phenomenon also discovered, for instance, in
Anderson hamiltonians with off-diagonal disorder Eilmes  _et al._ (1998),
where a similar cusp in the electronic DOS at $E=0$ is observed, an analytic
description of this phenomenon until now has not been provided. Here, we have
addressed in depth the characterization of the zero-energy singularity, by
putting together insight coming from theoretical advances and extensive
computer simulations. We have presented a generalization to the liquid state
of the heterogeneous elastic theory, by integrating in the elastic medium
description of the instantaneous normal modes spectrum an explicit temperature
dependence of the spatial shear modulus fluctuations. Specifically we have
been able not only to demonstrate the existence of a disorder-induced zero-
energy spectral singularity, but also fully characterize the disorder
dependent spectral shape, in both stable and unstable energy regions. We have
next tested the theoretical predictions against extensive Molecular Dynamics
simulations of an atomic glass-forming liquid, confirming the relevance of the
developed mean-field theory for realistic systems.

In a nutshell, we have demonstrated that focusing to the original INM spectrum
$\rho(\lambda)$ reveals much more interesting information than by transforming
it to the $\omega=\sqrt{\lambda}$ spectrum, which completely suppress most
relevant zero-energy features of the spectra. Our findings also pose serious
constraints on the allowed shape of $\rho(\lambda)$ employed in theoretical
developments, sometimes at variance with superficial recent proposals. They
finally constitute a very general step forward for all condensed matter
systems described by non-negligible off-diagonal disorder.

##  VI acknowledgments

We are grateful to Jean-Louis Barrat, Matteo Paoluzzi and Richard M. Stratt
for helpful discussions and suggestions. S. M. acknowledges support by the
project Heatflow (Grant No. ANR-18-CE30-0019-01) funded by the French national
funding agency “Agence Nationale de la Recherche”.

## References

  * Seeley and Keyes (1989) G. Seeley and T. Keyes, “Normal-mode analysis of liquid-state dynamics,” J. Chem. Phys. 91, 5581 (1989).
  * Stratt (1995) R. Stratt, “The instantaneous normal modes of liquids,” Acc. Chem. Res. 28, 1 (1995).
  * Bembenek and Laird (1995) S. D. Bembenek and B. B. Laird, “Instantaneous normal modes and the glass transition,” Phys. Rev. Lett. 74, 936 (1995).
  * Keyes (1997) T. Keyes, “Instantaneous normal mode approach to liquid state,” J. Phys. Chem. 101, 2921 (1997).
  * Sciortino and Tartaglia (1997) F. Sciortino and P. Tartaglia, “Harmonic dynamics in supercooled liquids: The case of water,” Phys. Rev. Lett. 78, 2385 (1997).
  * Sastry  _et al._ (2001) S. Sastry, N. Deo, and S. Franz, “Spectral statistics of instantaneous normal modes in liquids and random matrices,” Phys. Rev. E 64, 016305 (2001).
  * Taraskin and Elliott (2002) S. N. Taraskin and S. R. Elliott, “Disorder-induced zero-energy spectral singularity for random matrices with correlations,” Phys. Rev. B 65, 052201 (2002).
  * Debenedetti and Stillinger (2001) Pablo G Debenedetti and Frank H Stillinger, “Supercooled liquids and the glass transition,” Nature 410, 259–267 (2001).
  * Cavagna (2001) A. Cavagna, “Fragile vs. strong liquids: A saddle-ruled scenario,” Europhys. Lett. 53, 490 (2001).
  * Broderix  _et al._ (2000) K. Broderix, K. K. Bhatacharya, A. Cavagna, A. Zippelius, and I. Giardina, “Energy landscape of a Lennard-Jones liquid: Stastics of stationary points,” Phys. Rev. Lett. 85, 5360 (2000).
  * Angelani  _et al._ (2000a) L. Angelani, G. Parisi, G. Ruocco, and G. Viliani, “Potential energy landscape and long-time dynamics in a simple model glass,” Phys. Rev. E 61, 1681–1691 (2000a).
  * Angelani  _et al._ (2000b) L. Angelani, R. Di Leonardo, G. Ruocco, A. Scala, and F. Sciortino, “Saddles in the energy landscape probed by supercooled liquids,” Phys. Rev. Lett. 85, 5356 (2000b).
  * Angelani  _et al._ (2000c) L. Angelani, R. Di Leonardo, G. Ruocco, A. Scala, and F. Sciortino, “Quasisaddles as relevant points of the potential energy surface in the dynamics of supercooled liquids,” J. Chem. Phys. 116, 10297 (2000c).
  * Grigera  _et al._ (2002) T. S. Grigera, A. Cavagna, I. Giardina, and G. Parisi, “Geometric approach to the dynamic glass transition,” Phys. Rev. Lett. 88, 055502 (2002).
  * Fabricius and Sariolo (2002) G. Fabricius and D. A. Sariolo, “Dictance between inherent structures and the influence of saddles on approaching the mode couling transition in a simple glass former,” Phys. Rev. E 66, 031501 (2002).
  * Ciliberti and Grigeranz (2003) S. Ciliberti and T. S. Grigeranz, “Localization threshold of instantaneous normal modes from level-spacing statistics,” Phys. Rev. E 70, 061502 (2003).
  * Sampoli  _et al._ (2003) M. Sampoli, P. Benassi, R. Eramo, L. Angelani, and G. Ruocco, “The potential energy landscape in the Lennard-Jones binary mixture model,” J. Phys.: Condens. Matter 15, S1227 (2003).
  * Angelani  _et al._ (2003) L. Angelani, G. Ruocco, M. Sampoli, and E. Sciortino, “General features of the energylandscape in Lennard-Jones-like model liquids,” J. Chem. Phys. 119, 2120 (2003).
  * Berthier and Garrahan (2003) L. Berthier and J. P. Garrahan, “Real space origin of temperature crossover in supercooled liquids,” Phys. Rev. E 68, 041201 (2003).
  * Doye and Wales (2002) J. P. K. Doye and D. J. Wales, “Saddle points and dynamics of Lennard-Jones clusters, solids and supercooled liquids,” J. Chem. Phys. 116, 3777 (2002).
  * Doye and Wales (2003) J. P. K. Doye and D. J. Wales, “Comment on “quasisaddles as relevant points of the potential energy surface in the dynamics of supercooled liquids”,” J. Chem. Phys. 118, 5263 (2003).
  * Wales and Doye (2003) D. J. Wales and J. P. K. Doye, “Stationary points and dynamics in high-dimensional systems,” J. Chem. Phys. 119, 12409 (2003).
  * Doliwa and Heuer (2003) B. Doliwa and A. Heuer, “Energy barriers and activated dynamics in a supercooled Lennard-Jones liquid,” Phys. Rev. E 67, 031506 (2003).
  * Grigera (2006) T. S. Grigera, “Geometric properties of the potential energy of the soft-sphere binary mixture,” J. Chem. Phys. 124, 064502 (2006).
  * Clapa  _et al._ (2012) V. I. Clapa, T. Kottos, and E. W. Starr, “Localization transition of instantaneous normal modes and liquid diffusion,” J. Chem. Phys. 136, 144504 (2012).
  * Coslovich  _et al._ (2019) D. Coslovich, A. Ninarello, and L. Berthier, “A localization transition underlies the mode-coupling crossiver of glasses,” SciPost Phys. 7, 077 (2019).
  * Zhang  _et al._ (2019) W. Zhang, J. F. Douglas, and F. W. Starrman, “What does the instantaneous normal mode spectrum tell us about dynamical heterogeneity in glass-forming fluids?” J. Chem. Phys. 151, 184904 (2019).
  * Goldstein (1969) M. Goldstein, “Viscous liquids and the glass transition: A potential energy barrier picture,” J. Chem. Phys. 51, 3728 (1969).
  * Stillinger and Weber (1982) Frank H. Stillinger and Thomas A. Weber, “Hidden structure in liquids,” Phys. Rev. A 25, 978 (1982).
  * Stillinger (1995) F. H. Stillinger, “A topographic view of supercooled liquids and glass formation,” Science 267, 1935 (1995).
  * Berthier and Biroli (2011) Ludovic Berthier and Giulio Biroli, “Theoretical perspective on the glass transition and amorphous materials,” Reviews of modern physics 83, 587 (2011).
  * Götze (2009) Wolfgang Götze,  _Complex dynamics of glass-forming liquids. A mode-coupling theory_ , Vol. 143 (Oxford University Press, 2009).
  * Kurchan and Laloux (1996) J. Kurchan and L. Laloux, J. Phys. A: Math. Gen. 29, 1929 (1996).
  * Schirmacher  _et al._ (2022) W. Schirmacher, T. Bryk, and G. Ruocco, “Modeling the instantaneous normal mode spectra of liquids as that of unstable elastic media,” Proc. Nat. Acad. Sci. 119, e2119288199 (2022).
  * Schirmacher (2006) Walter Schirmacher, “Thermal conductivity of glassy materials and the “boson peak”,” Europhys. Letters 73, 892 (2006).
  * Schirmacher  _et al._ (2007) Walter Schirmacher, Giancarlo Ruocco, and Tullio Scopigno, “Acoustic attenuation in glasses and its relation with the Boson Peak,” Physical Review Letters 98, 025501 (2007).
  * Marruzzo  _et al._ (2013) Alessia Marruzzo, Walter Schirmacher, Andrea Fratalocchi, and Giancarlo Ruocco, “Heterogeneous shear elasticity of glasses: the origin of the boson peak,” Scientific reports 3, 1 (2013).
  * Schirmacher  _et al._ (2014) W. Schirmacher, T. Scopigno, and G. Ruocco, “Theory of vibrational anomalies in glasses,” J. Noncryst. Sol. 407, 133 (2014).
  * Mizuno  _et al._ (2020) Hideyuki Mizuno, Giancarlo Ruocco, and Stefano Mossa, “Sound damping in glasses: Interplay between anharmonicities and elastic heterogeneities,” Physical Review B 101, 174206 (2020).
  * Mizuno and Mossa (2019) Hideyuki Mizuno and Stefano Mossa, “Impact of elastic heterogeneity on the propagation of vibrations at finite temperatures in glasses,” Condensed Matter Physics 22, 43604 (2019).
  * Pan  _et al._ (2021) Z. Pan, O. Benzine ans S. Sawamura, R. Limbach, A. Koike, T. D. Bennett, G. Wilde, W. Schirmacher, and L. Wondraczek, “Disorder classification of the vibrational spectra of modern glasses,” Phys. Rev. B (2021), accepted.
  * McKane and Stone (1981) A. J. McKane and M. Stone, “Localization as an alternative to Goldstone’s theorem,” Ann. Phys. (N. Y.) 131, 36 (1981).
  * John  _et al._ (1983) S. John, H. Sompolinky, and M. J. Stephen, “Localization in a disordered elastic medium near two dimensions,” Phys. Rev. B 28, 5592 (1983).
  * Hertzfeld and Litovitz (1959) K. F. Hertzfeld and T. A. Litovitz,  _Absorption and Dispersion of Ultrasonic Waves_ (Academic Press, New York, 1959).
  * Ganter and Schirmacher (2010) C. Ganter and W. Schirmacher, “Rayleigh scattering, long-time tails and the harmonic spectrum of topologically disordered systems,” Phys. Rev. B 82, 094205 (2010).
  * Zaccone and Baggioli (2021) A. Zaccone and M. Baggioli, Proc. Nat. Acad. Sci. 118, e2022303118 (2021).
  * Stamper  _et al._ (2022) C. Stamper, D. Cortie, Z. Yue, X. Wang, and D. Yu, J. Chem. Phys. Lett. 13, 3105 (2022).
  * Thompson  _et al._ (2022) Aidan P. Thompson, H. Metin Aktulga, Richard Berger, Dan S. Bolintineanu, W. Michael Brown, Paul S. Crozier, Pieter J. in T Veld, Axel Kohlmeyer, Stan G. Moore, Trung Dac Nguyen, Ray Shan, Mark J. Stevens, Julien Tranchida, Christian Trott, and Steven J. Plimpton, “LAMMPS - a flexible simulation tool for particle-based materials modeling at the atomic, meso, and continuum scales,” Computer Physics Communications 271, 108171 (2022).
  * Bradbury  _et al._ (2018) James Bradbury, Roy Frostig, Peter Hawkins, Matthew James Johnson, Chris Leary, Dougal Maclaurin, George Necula, Adam Paszke, Jake VanderPlas, Skye Wanderman-Milne, and Qiao Zhang, “JAX: composable transformations of Python+NumPy programs,”  (2018).
  * Virtanen  _et al._ (2020) Pauli Virtanen, Ralf Gommers, Travis E. Oliphant, Matt Haberland, Tyler Reddy, David Cournapeau, Evgeni Burovski, Pearu Peterson, Warren Weckesser, Jonathan Bright, Stéfan J. van der Walt, Matthew Brett, Joshua Wilson, K. Jarrod Millman, Nikolay Mayorov, Andrew R. J. Nelson, Eric Jones, Robert Kern, Eric Larson, C J Carey, İlhan Polat, Yu Feng, Eric W. Moore, Jake VanderPlas, Denis Laxalde, Josef Perktold, Robert Cimrman, Ian Henriksen, E. A. Quintero, Charles R. Harris, Anne M. Archibald, Antônio H. Ribeiro, Fabian Pedregosa, Paul van Mulbregt, and SciPy 1.0 Contributors, “SciPy 1.0: Fundamental Algorithms for Scientific Computing in Python,” Nature Methods 17, 261–272 (2020).
  * Ovtchinnikov (2008) Evgueni E Ovtchinnikov, “Computing several eigenpairs of hermitian problems by conjugate gradient iterations,” Journal of Computational Physics 227, 9477–9497 (2008).
  * Eilmes  _et al._ (1998) A. Eilmes, R. A. Römer, and M. Schreiber, “The two-dimensional Anderson model of localization with random hopping,” Eur. Phys. J. B 1, 2938 (1998).


# Phonetic Posteriorgrams based End-to-End Adversarial Many-to-Many Singing
Voice Conversion with Non-Parallel Data

Haohan Guo1, Heng Lu2, Lei Xie1
1 School of Computer Science and Engineering, Northwestern Polytechnical
University, Xi’an, China
2 Tecent AI Lab, Beijing, China
<EMAIL_ADDRESS><EMAIL_ADDRESS>Work performed as an intern at
Tencent AI Lab. Lei Xie is the corresponding author.

###### Abstract

Many-to-many singing voice conversion aims to convert the timbre of a song to
anyone’s voice. In this paper, we propose an adversarial training based
embedding-to-audio many-to-many SVC framework (EA-SVC) using a non-parallel
singing corpus to achieve high conversion quality. In this framework, the
audio is decomposed to content, pitch, and timbre represented by phonetic
posteriorgram extracted from the automatic speech recognition (ASR) model, F0,
and speaker embedding extracted from the speaker recognition (SR) model,
respectively. Inspired by MelGAN, we construct a model with the encoder-
upsampler struture to directly generate audio samples from these embeddings,
and train it with a two-stage training strategy, STFT training and adversarial
training. In our experiments, the MOS test shows that our proposed apporach
has obvious improvement over the conventional framework and the WaveNet based
approach in singing quality and singer similarity. In ablation studies, our
proposed timbre tokens help achieve better singer similarity by the preference
ratio 50%, and the conditional multi-scale discriminator can improve both
singing quality and singer similarity by the preference ratio 50% and 50%
respectively. We also compare our two-stage training algorithm with other to
show its important in EA-SVC model training. In addition, we also make an
analysis on pitch control to investigate its controllability in pitch.

## Introduction

Many-to-Many singing voice conversion (SVC) is a technology that converts the
timbre of a song to any target singer. It has been highly concerned due to its
potential applied in various interaction scenarios and products for better
interactive experience and entertainment. Compared with voice conversion (VC,
mainly for speech), the high variety of pitch in SVC brings a higher
requirement for the model capability to ensure stable and high-quality
conversion. To achieve it, many algorithms are proposed based on the non-
parallel data with the large data size and high diversity.

The most many-to-many SVC approaches are based on a generalized autoencoder
framework, which firstly extracts features representing content (lyrics),
pitch and timbre respectively from the audio, then inverses them back in
training. The conversion can be completed by replacing the timbre features
with the target one. So, to generate the expected audio, it is essential to
keep these features accurate and disentangled, i.e. no overlap of information
between any two features. For example, unsupervised learning based methods
aims to extract timbre and pitch independent features in the encoder via
adversarial training Nachmani2019UnsupervisedSV; 9054199 or VAE
luo2020singing; kameoka2019acvae, then combine it with singer ID/embedding and
pitch to reconstruct the audio in the decoder. It keeps sufficient information
from the audio in the encoder to ensure accurate reconstruction, but the high-
purity decomposition still cannot be achieved due to the strong entanglement
of different information in audio, hence degrades the controllability.

Different from it, the PPG-based approach does not extract the features in a
self-learning manner, but use the knowledge learned by models in other tasks.
Phonetic posteriorgrams (PPG), the embedding extracted from automatic speech
recognition (ASR) model trained with large-scale speech corpus, is a more
robust phonetic representation independent of acoustic features. It has been
shown as a reliable and effective approach in both voice conversion
sun2016phonetic; tian2018average; Polyak2019TTSSS and singing voice conversion
Chen2019SingingVC; gao2020personalized. However, compressing the waveform to
low-resolution PPG features also brings the loss of information. So the
generative model is necessarily used to compensate for it to keep high
fidelity and naturalness. As the state-of-the-art generative model in image
generation karras2020analyzing, generative adversarial network (GAN) has also
shown its strong generative performance in audio generation as a neural
vocoder, MelGAN kumar2019melgan. It indicates the great potential applied in
PPG-based SVC. In this paper, we purpose to improve many-to-many SVC with GAN,
and demonstrate that the model trained adversarially in a proper way can
achieve high-quality conversion.

Our main contributions are as follows:

  * •

We propose a PPG-based end-to-end adversarial many-to-many SVC framework. The
generator uses CNN and BLSTM to encode PPG and F0, and a singer-token based
timbre encoder to better extract more generalized timbre representation from
the speaker embedding, then upsamples them to the waveform via a stack of
transposed convolutional layers. Besides the multi-scale discriminator, we
also propose the conditional multi-scale discriminator to make the generated
audio better match the expected pitch and timbre.

  * •

We optimize the training process with a two-stage training strategy, which
firstly trains the generator using multi-resolution STFT loss to achieve a
stable mapping relationship, then use adversarial training to improve its
generative performance.

  * •

We conduct subjective evalution for our approach. It shows the best
performance in singing quality and singer similarity with the MOS of 4.0 & 4.2
over the conventional framework 3.0 & 3.0 that reconstructs wavefrom from the
predicted Mel spectrogram, and the WaveNet based apporach 2.0 & 2.0. Ablation
studies shows the effectiveness of our proposed modules.

Figure 1: The architecture of the conventional SVC framework

## PPG-based Singing Voice Conversion

In this section, we will firstly introduce the framework of PPG-based many-to-
many SVC, including two modules, features extraction, and audio generation.

### Features Extraction

The singing audio usually includes three main components, content, pitch, and
timbre. So the audio can be reconstructed approximately if providing their
corresponding features, and the conversion can be achieved by giving the
expected timbre representation.

Similar with voice conversion, the content in singing voice can be represented
by PPG extracted from the ASR model. Although the singing voice has higher
variety, the ASR model trained with the large-scale speech corpus can still
deal with it, and provide robust, pure content representation. Timbre is the
feature strongly related to the speaker or singer. In the same way, we extract
the timbre representation from the speaker recognition (SR) model trained with
the corpus containing a variety of speakers, recording environments. Pitch
(F0) contains the note of a song and the singing style of a singer, like
trill. In our corpus without serious noise and reverberation, it can be
accurately extracted using signal processing based tools from the waveform.

### Audio Generation

For PPG-based SVC, the generation of the target audio can be composed of two
steps, mapping PPG features to the acoustic space, and upsampling the low-
resolution inputs to the waveform with a high sample rate. For example, in the
conventional framework, the generation is achieved by two models. An acoustic
model is used to predict low-dimensional spectral features, e.g. Mel
spectrogram, from input features. Then the waveform is generated by a
generative model based neural vocoder van2016wavenet;
kalchbrenner2018efficient; prenger2019waveglow. This approach can ensure both
stable space mapping relationship, and high-quality audio generation. So it
has been widely used in various related tasks, such as singing synthesis
kim2018korean and SVC gao2020personalized.

However, there are also two limitation in the framework. Firstly, as the
intermediate audio representation, the low-dimensional spectral feature can
help model a stable mapping relationship, but also remove some timbre and
fine-grained pitch related information. It leads to the limitation of accurate
timbre and pitch reconstruction in the neural vocoder. Secondly, the vocoder
is trained with the recording data, but generates audio according to the
spectral features predicted from the acoustic model. The gap between the
predicted features and the real one may be amplified to degrade sound quality.
Even if the problem can be alleviated by training the vocoder with predicted
features, the acoustic model cannot be trained to adapt the vocoder.

Fortunately, these problems can be avoided if we directly map PPG features to
the waveform in one model. To achieve it, it is essential to propose a better
generative model with good generative capability and stability to keep both
stable space conversion and high-quality audio generation. MelGAN has shown
great potential in audio generation, which can generate high-quality audio
samples in parallel. Besides, rather than the auto-regressive model that has
to be trained in teacher-forcing mode, its training process is more
controllable. To further exploit its effect, an end-to-end adversarial SVC
approach is proposed.

Figure 2: The model structure of the generator in EA-SVC. “x3” means the three
blocks in cascade. The dimension written beside the layer donates its output
channel size.

## End-to-End Adversarial Singing Voice Conversion

In this section, the details of our proposed framwork, EA-SVC, will be
introduced, including the model architecture, training objectives, and the
two-stage training strategy.

### Model Architecture

The framework contains three modules. Feature extration provides 218-dim PPG,
1-dim F0 and 128-dim speaker embedding (SE), as the input features. The
generator converts them to the corresponding audio samples through its encoder
and upsampler. The discriminator is used in the adversarial training process
to improve the generative performance.

#### Encoder

In the generator, before upsampling the input features, the encoder is firstly
used to map these frame-level features to the acoustic space. As shown in the
left part of figure 2, it encodes the three features separately to help better
extract sequential information of the feature itself. For PPG, it firstly
embeds them into a 128-dim space using a convolutional layer with the kernel
size 1. Three residual stacks containing two convolutional layers with the
kernel size 3 are used to extract more high-level and contextual information.
After it, the BLSTM layer is used to extract the 256-dim vectors with more
sequential information. Finally, the convolutional layer with the kernel size
1 compresses them to 64-dim bottleneck features as stable and effective
content representation. The 1-dim F0 sequence is also processed in the same
way.

The timbre representation is extracted in a different way. To perform well on
unseen singers, it is necessary to model a more robust and controllable timbre
space. GST-Tacotron wang2018style can effectively control the speaking style
by controlling the weights on different trainable style tokens. Inspired by
it, we propose timbre tokens during encoding speaker embeddings, which is
shown as follows:

$\displaystyle H_{S}=W_{2}(\sigma_{L}((W_{1}(S)+b_{1}))+b_{2}$ (1)
$\displaystyle Y=Softmax(H_{S})*E_{1:N}.$ (2)

where $W_{1}$, $W_{2}$, $b_{1}$, $b_{2}$ donate the trainable weight matrixs
and biases, $\sigma_{L}$ donates “LeakyReLU”, and we set $N=64$ in our work.
The speaker embedding $S$ is non-linearly transformed to a 128-dim vector
firstly, and then converted to N-dim weight distribution through a Softmax
layer. According to the weight vector, the N trainable 128-dim embeddings
$E_{1:N}$ are weighted added together, then mapped to a 64-dim vector $Y$ as
the timbre representation. Training these embeddings can be seen as a process
of clustering or disentangle. The N embeddings can represent the expected
timbre by combining in the specific proportions. Notably, we adopt Synthesizer
(Dense version) tay2020synthesizer to calculate weight distribution, which is
a simplified version of GST.

To concatenate these features together, SE is replicated to the same length
with other features. We keep the dimension of these encoded features the same
to prevent the model from overfitting to one kind of feature during the
training. Finally, a convolutional layer with the kernel size 7 is used to
fuse these features to a 512-dim contextual vector as the encoder output.

#### Upsampler

The structure of the upsampler is similar with MelGAN, which is composed of a
stack of transposed convolutional layers followed by residual blocks. The
right part of figure 2 shows its detailed architecture. The encoder outputs
with the frameshift 10ms are at 240x lower temporal resolution than the
waveform with the sample rate 24,000. So three transposed convolutional layers
with stride 8, 6, 5 respectively are used to upsample them to the same sample
rate. Their kernel sizes are double of the stride to alleviate Checkerboard
artifacts odena2016deconvolution, which may result in audible high frequency
hissing noise. The deep residual block containing four residual layers with
dilation $3^{0},3^{1},3^{2},3^{3}$ respectively follows each transposed
convolutional layer to enhance the model capability. The deeper model
structure can help the generator capture more high-level information, which
makes the details of the audio more realistic. And the kernel size and wider
diliation size can efficient increase the receptive field, enhance contextual
correlation, thus improve the naturalness and smoothness of the generative
audio. The final convolutional layer with the kernel size 7 and activation
function “Tanh” generates the waveform normalized between -1 and 1.

There is no noise vector as input in the generator, because diversity is not
necessary in our task. We only purpose to improve the naturalness and fidelity
via adversarial training. In addition, to stablize the training process, we
apply LeakyReLU with $\alpha=0.2$ as the default activation function, and
weight normalization for all convolutional layers.

#### Discriminator

The discriminator consists of multiple discriminator blocks, which can
adversarially train the generator from different aspects.

Similar with MelGAN, three multi-scale discriminators (MSDs) are used to deal
with the audio at different sample rate, 24000, 12000, and 6000, to capture
different features at different frequency range. The audio is downsampled via
an average pooling layer with kernal size 4 and stride 2. As shown in the left
of figure 3, each MSD is a window-based discriminator, which compresses the
input audio to small chunks via four strided convolutional layers. Due to the
large kernel size in the strided convolutional layer, grouped convolutions are
used to keep number of parameters small.

The unconditional discriminator can reduce the gap of distribution between the
generated audio and the real one, but cannot keep sufficient correlation
between the audio and input features. In our experiments, the generated audio
can represent content information well via adversarial training, but is still
lack of stable performance in pitch and timbre. So we propose a conditional
multi-scale discriminator (CMSD) inspired by GAN-TTSbinkowski2019high, which
classifies according to F0 and speaker embedding, to make outputs better match
them. The structure is is shown in the right part of figure 3, which has the
same structure except for the CondBlock. It is inserted after two strided
convolutional layers, whose inputs have a similar frame rate with conditions.
In the block, the output from the preceding convolutional layer is given to
the strided convolutional layer without activation function, the condition is
processed by a convolutional layer with LeakyReLU activation function and a
strided convolutional layer. Then these two features are added together and
processed by Leaky ReLU as the output.

There are two CMSD blocks in the discriminator to deal with F0 and SE
respectively. Because the pitch information mainly affects the low-frequency
part, we firstly use the average pooling to downsample the audio to 6k sample
rate. Then we upsample F0 sequence with a frame rate of 100 fps to 375 fps via
nearest neighbor interpolation for length alignment. For SE, the audio is not
downsampled, because the timbre information exists in all bands. But we use
larger strides with $4,4,8,16$ to capture this global feature. In addition, we
also apply weight normalization for all layers in the discriminator for stable
training process.

Figure 3: The structure of the discriminator block. $K,S,G$ donates the kernel
size, stride, and groups, respectively. The dimension written at the right of
the convolutional layer donates the output channels of the layer.

### Training Objective

To achieve more stable and efficient training process, we combine adversarial
loss and multi-resolution STFT loss as the training objective, which has shown
improvement in yang2020multi.

#### Adversarial loss

The least-squares GAN version loss is used in our work. We also try the hinge
version loss, both works, and no obvious difference in the results. For each
discriminator block $D_{k}$, its task is to minimize the loss defined as
follows:

$\displaystyle L_{adv}(D_{k})=(1-D_{k}(x))^{2}+D_{k}(\hat{x})^{2}$ (3)

where $x$ and $\hat{x}$ denote the real waveform and the generated waveform,
respectively.

For generator, except for the least-squares loss, we also apply feature
matching proposed in MelGAN for more effective adversarial training. It
minimizes the L1 distance between the feature maps of the real and generated
audio in each hidden layer of the discriminator.

$\displaystyle
L_{FM}(G,D_{k})=\sum^{N_{k}}_{i=0}\frac{1}{N_{k}}\left\|D^{(i)}_{k}(x)-D^{(i)}_{k}(G(c))\right\|_{1}$
(4) $\displaystyle L_{FM}(G,D)=\frac{1}{K}\sum^{K}_{k=1}L_{FM}(G,D_{k})$ (5)

where $c$ denote the input features, $K$ and $N_{k}$ denote the number of
discriminator blocks and the number of the hidden layers (except for the
output layer) of the k-th discriminator block, respectively. So the
adversarial loss for the generator is defined as:

$\displaystyle L_{adv}(G,D)=\frac{1}{K}(1-G(c))^{2}+\lambda_{FM}L_{FM}(G,D)$
(6)

where $\lambda_{FM}$ is set to 10 in our experiments.

#### Multi-resolution STFT (MR-STFT) loss

It computes STFT loss between the real samples $x$ and the generated samples
$\hat{x}$ on different resolutions. STFT loss is composed of two loss
functions, spectral convergence $L_{sc}$, which can emphasize highly on large
spectral components to help early stage of training, and log STFT magnitude
loss $L_{mag}$, which can help fit small spectral components to help the late
stage of training arik2018fast. The functions are defined as below:

$\displaystyle L_{stft}(x,\hat{x})=L_{sc}(x,\hat{x})+L_{mag}(x,\hat{x})$ (7)
$\displaystyle
L_{sc}(x,\hat{x})=\frac{\left\||STFT(x)|-|STFT(\hat{x})|\right\|_{F}}{\left\||STFT(x)|\right\|_{F}}$
(8) $\displaystyle
L_{mag}(x,\hat{x})=\frac{1}{N}\left\|log|STFT(x)|-log|STFT(\hat{x})|\right\|_{1}$
(9)

where $|STFT(\cdot)|$ and $N$ donate the STFT magnitudes and its dimension
respectively, $\left\|\cdot\right\|_{F}$ and $\left\|\cdot\right\|_{1}$ donate
the Frobenius and $L_{1}$ norms respectively.

Multi-resolution loss is achieved by averaging multiple STFT losses with
different analysis parameters, i.e. frameshift, frame length, and FFT size. It
can help the model generate audio signals with more stable performance in
frequency domain. In our work, three STFT losses are used, which are shown in
Table 1. Finally, the loss function for the generator is:

$\displaystyle L_{G}(G,D)=L_{stft}(x,G(c))+\lambda_{adv}L_{adv}(G,D)$ (10)

where $\lambda_{adv}$ is set to 4 in our experiments.

Table 1: The analysis parameters of the STFT losses in our work. Hanning window is used in the signal analysis. Frame shift | Frame length | FFT size
---|---|---
50 | 240 | 512
120 | 600 | 1024
240 | 1200 | 2048

### Two-stage Training Algorithm

In our experiments, MR-STFT loss shows its great importance in improving
stability, making the generated audio better match input features. To make
full use of it, and further stablize the adversarial training process, we
propose a two-stage training algorithm based on yang2020multi. In the first
stage, only $L_{stft}(G,D)$ is used to train the generator until convergence.
In the second stage, the generator and the discriminator are trained
alternately with $L_{G}(G,D)$ as the training objective.

We purpose to achieve an accurate and stable mapping relationship between
input features and the generated audio via stage 1, which can generate smooth,
intelligible audio with accurate tone and distinguishable timbre. It can also
be seen as a process of modeling the mean value of the target distribution.
The adversarial training stage is responsible to model the variance of the
target distribution to replace the averaged value with the waveform with the
highest probability in the distribution. So the audio generated by the model
trained with adversarial loss is more natural and realistic. In this way, we
can effectively improve output quality more stably.

Figure 4: The T-SNE visualization of distribution of speaker embeddings. The
yellow and purple dots donate the speaker embeddings from NUS-48E and Jingming
dataset, respectively.

## Experiments

### Experimental Configurations

#### Corpus

A Mandarin mutli-singer singing corpus, Jingming, is used in our experiments.
Compared with the open-source standard English singing corpus, NUS-48E
duan2013nus, which contains 169-min audio recordings with 12 singers and four
songs for each singer, it is a larger-scale corpus. It contains nearly 3800
clean singing segments of Chinese popular songs with total 25 hours and about
20 seconds per segment. Figure 4 visualizes the distribution of speaker
embeddings extracted from each singing segments in these two datasets.
Jingming (purple part) has wider distribution and more clusters, which shows
more singers and higher diversity in timbre than NUS-48E. However, it brings a
higher requirement for modelling capability, because it is scratched from the
user data in a sing App, which contains inaccurate pronunciation, uneven
singing capability, and various recording environments.

#### Feature extraction

In our experiments, PPG is the bottleneck feature extracted from a ASR model
trained with a 20,000-hour multi-speaker speech corpus. Its frameshift is
10ms. WORLD morise2016world is used to extract the F0 with frameshift 10ms.
Each singing segment has a speaker embedding, which is a X-Vector extracted
from the TDNN model snyder2018x trained with a large-scale corpus containing
8800 speakers. We downsample all audios to the sample rate 24,000 for model
training.

Figure 5: The training loss curves in stage1.

#### Training configuration

In our experiments, all models are trained on a single P40 GPU with batch size
32. Each sequence in a batch has 24,000 samples (1 second), which is randomly
selected from the audio set. We combine RAdam and LookAhead as the optimizer
to speed up and stabilize the whole training process. As shown in Figure 5,
RAdam liu2019variance helps speed up training convergence in the early stage
via rectifying the variance of the adaptive learning rate. LookAhead
zhang2019lookahead can further help stablize the late training process by
updating parameters with the slow and fast weights trajectory. In our
experiments, the beta of RAdam is set to (0.5, 0.9), the K in LookAhead is set
to 5. The learning rate is exponentially decayed based on the training step
$t$ using this function:

$l=\left\\{\begin{matrix}l_{i}&t<t_{w}\\\ l_{i}*0.5^{(t/t_{w})}&t>t_{w}\\\
\end{matrix}\right.$ (11)

where $l_{i}$, $t_{w}$ donate the initial learning rate and the warmup
training steps respectively. In stage 1, we train the generator for 200,000
steps with $l_{i}=10^{-3},t_{w}=50,000$. In stage 2, we train the generator
and the discrimiantor for 800,000 steps in adversarial way. We set
$l_{i}=10^{-4},t_{w}=200,000$ for the generator, and a higher learning rate
$l_{i}=5*10^{-4}$ for the discriminator to encourage it to learn faster to
capture more distinguishable information. In addition, to stabilize the
parameter update of the generator, we also apply the gradient clip with the
max 2-norm value of 0.5.

### Mean Opinion Score (MOS) Test

We conduct MOS test to evaluate our proposed EA-SVC approach. There are 4
template songs (2 males, 2 females) in the test set. Each template song is
converted to 8 unseen singers (4 males, 4 females), so 32 test cases are
included in the set totally. 20 Chinese native speakers are involved in this
test. They need to give a score for each sample in a test case according to
the criterion: 1 = Bad; 2 = Poor; 3 = Fair; 4 = Good; 5 = Excellent. The final
score for each model is calculated by averaging the collected results.

Except for EA-SVC, two other models involves in this test, the conventional
framework, and a WaveNet based approach. They are trained with Jingming
dataset.

#### Conventional framework

The conventional framework has shown stable conversion performance in SVC
Chen2019SingingVC; gao2020personalized. In this work, we construct a similar
system with a MelGAN-based vocoder to show the improvement of the end-to-end
framework.

The acoustic model has a similar structure with the encoder in EA-SVC. It
inserts two 256-dim BLSTM layers before the final convolutional layer. The
output is 80-dim Mel-spectrogram, which is extracted from the 1025-dim log-
magnitude spectrogram with frameshift 10ms and frame length 50ms. The training
configuration is similar to stage 1, which trains the model for 200,000 steps
with a batch size of 32, each sequence in a batch has the length of 4 seconds.

We construct a MelGAN vocoder with a convolutional layer followed by an
upsampler in the generator. The discriminator and training configuration is
the same as the EA-SVC without CMD blocks. The model is trained for 1,000,000
steps, which is well converged.

#### WaveNet based approach

The idea directly mapping PPG features to waveform has shown improvement in
voice conversion by training a WaveNet model tian2019vocoder. In our work, we
also implement a silimiar WaveNet-based SVC model to compare the performance
of different models in SVC. It replaces the upsampler in EA-SVC with WaveNet
Model. The structure and configuration of WaveNet follows the Nvidia
version111The code links: https://github.com/NVIDIA/nv-wavenet. The model is
trained for 800,000 steps on four P40 GPUs with a batch size of $16*24000$
samples totally to achieve good convergence.

#### Results

Table2 shows the MOS results of these three SVC models. EA-SVC achieves the
best performance in both singer quality and singer similarity by the MOS of
4.0 and 4.0. The WaveNet shows the worst performance in both singer quality
and singer similarity by the MOS of 3.0 and 2.0. In addition to the serious
generated noise and poor singer similarity, the intelligible problems also
exist. We infer that WaveNet is seriously disturbed by exposure bias when
trained with the high-diversity dataset. It also shows EA-SVC is much more
stable model. In conventional framework, the audio generated from the
predicted spectral features has audible buzzy sound, which degrades sound
quality. Singer similarity is also limited by the low-dimensional Mel
spectrogram. 222The audio sample are availble at
https://hhguo.github.io/demo/publications/ConvTTS/index.html.

Table 2: The results of the MOS test models | singing quality | singer similarity
---|---|---
Conventional | 4.0 | 4.0
WaveNet | 4.0 | 4.0
EA-SVC | 4.0 | 4.0

### Ablation Study

#### Timbre encoder

The timbre token is proposed tp improve its generalization on unseen singers.
To show its improvement in singer similarity, we conduct a preference test to
compare the model with or without singer tokens, which replaces it with a
dense layer. The listeners need to choose the better one in each pair of
samples generated by these two models respectively according to the overall
quality and singer similarity. These models are trained without CMSD to remove
its affection.

Table 3: Preference test: w/ or w/o timbre tokens | w/o | neutral | w/ | p-value
---|---|---|---|---
Quality | 20% | 30% | 50% | 0.01
Similarity | 20% | 30% | 50% | 0.01
Table 4: Preference test: w/ or w/o CMSD | w/o | neutral | w/ | p-value
---|---|---|---|---
Quality | 20% | 30% | 50% | 0.01
Similarity | 20% | 30% | 50% | 0.01

As shown in table 3, the timbre encoder with singer tokens can improve the
performance in these two aspects. Especially for the similarity, the
improvement is more significant by the preference ratio of 60%.

#### Conditional multi-scale discriminator (CMSD)

We also make a preference test for the models with or without CMSD. The
results show that CMSD further improves the performance in both quality and
similarity. In our analysis, we find that the problem of unsmooth pitch is
effectively alleviated in the model with CMSD, which is shown in figure 6. It
shows its effectiveness in enhancing the relation between the generated audio
and the input features.

Figure 6: The spectrogram of the generated sustains. (a) and (b) are generated
by the models w/ and w/o CMSD respectively. The blue box of (a) shows the
phenomena of the unsmooth pitch track.

#### Training strategy

To show the important of the two-stage training algorithm and MR-STFT loss, we
compare the models trained with different training strategies:

  * •

$M_{0}$: Trained in stage2 with $L_{adv}$

  * •

$M_{1}$: Trained in stage2 with $L_{stft}+\lambda_{adv}L_{adv}$

  * •

$M_{2}$: Trained in stage1 + stage2 with $L_{stft}+\lambda_{adv}L_{adv}$

As shown in Figure7, both $M_{0}$ and $M_{1}$ easily generate fuzzy formant,
causing serious problems in intelligibility. $M_{2}$ produces more smooth and
clear formant, hence ensures the fidelity of the generative audio.

Figure 7: The Mel spectrograms of the audios generated by the models: (a)
$M_{0}$, (b) $M_{1}$, (c) $M_{2}$

#### Pitch control

We input the pitch estimated by the signal processing based method
morise2016world to represent music note. To show its controllability on the
generated audio, we investigate its performance on pitch control.

The key of a song can be changed higher or lower by globally multiplying F0
with a factor, $\gamma$. For example, the pitch should be an octave (12 keys)
higher when $\gamma=2.0$. To verify the effect of our model in global pitch
control, we rebuild a song with different $\gamma$, which is shown in Figure8.
The results show that our model performs well with $\gamma$ ranging from 0.8
to 1.2, but generated unstable jitters when the pitch is adjusted too much,
such as $\gamma=0.5$ and $\gamma=1.5$. We infer that this situation can be
alleviated by data augmentation, which covers more pitch contours in the train
set.

Figure 8: The pitch contours of the same singing segment converted with
different scales.

## Conclusion

This paper proposes a new framework for PPG-based SVC, EA-SVC. The model can
directly map PPG features to the waveform in parallel via adversarial
training. Based on MelGAN, we construct a generator with an encoder, which
uses two CNN-BLSTM modules to encode PPG and F0 respectively, and a timbre
encoder to extract timbre information from the speaker embedding. To enhance
its generalization for unseen singers, we propose timbre tokens to disentangle
timbre into $N$ components to better represent it. In the adversarial
training, the conditional multi-scale discriminator is proposed to help
generate the audio better matching the input features. And we adopt a two-
stage training strategy to ensure a more stable and effective training
process. We conduct a MOS test to compare EA-SVC with other approaches. The
score of 4.5 shows that our approach achieves better performance in both
singing quality and singer similarity than the conventional framework and the
WaveNet based approach. We also make ablation studies to show the correction
of our choice of the model architecture.