# Early-time resonances in the three-dimensional wall-bounded axisymmetric
Euler and related equations

Sai Swetha Venkata Kolluru  Rahul Pandit  Department of Physics, Indian
Institute of Science, Bengaluru, India - 560012

###### Abstract

We investigate the complex-time analytic structure of solutions of the
3D-axisymmetric, wall-bounded, incompressible Euler equations, by starting
with the initial data proposed in Ref. 1, to study a possible finite-time
singularity. We use our pseudospectral Fourier-Chebyshev method [2], with
quadruple-precision arithmetic, to compute the time-Taylor series coefficients
of the flow fields, up to a high order. We show that the resulting
approximations display early-time resonances; the initial spatial location of
these structures is different from that for the tygers, which we have obtained
in Ref. 2. We then perform asymptotic analysis of the Taylor-series
coefficients, by using generalised ratio methods, to extract the location and
nature of the convergence-limiting singularities and demonstrate that these
singularities are distributed around the origin, in the complex-$t^{2}$ plane,
along two curves that resemble the shape of an eye. We obtain similar results
for the 1D wall-approximation (of the full 3D-axisymmetric Euler equation)
called the 1D HL model, for which we use Fourier-pseudospectral methods to
compute the time-Taylor series coefficients of the flow fields. Our work
examines the link between tygers, in Galerkin-truncated pseudospectral
studies, and early-time resonances, in truncated time-Taylor expansions of
solutions of PDEs, such as those we consider.

##  I Introduction

The regularity problem for hydrodynamical partial differential equations
(PDEs), such as the three-dimensional (3D) Euler and the 3D Navier-Stokes
equations [3, 4, 5], has been the subject of several theoretical and numerical
investigations over the last few decades. However, it is still not known
whether the solutions of these PDEs, starting from smooth initial data, retain
smoothness, for all time, or whether they lose regularity, at a finite time,
leading to a finite-time singularity (often called blowup in common parlance)
[6]. The breakdown of the regularity of solutions of the 3D Euler equations
has been conjectured to play a key role in understanding of the structure of
Navier–Stokes turbulence [see, e.g., Refs. 7, 8, 9, 10, 11].

The singularity problem for hydrodynamical PDEs has often been investigated by
examining the analytic structure of solutions in the complex-time domain. The
solution of the governing PDE is represented as a Taylor series in integer
powers $n$ of the time $t$. The asymptotic behavior of these Taylor
coefficients is determined by the dominant complex-time singularities that lie
closest to the point about which the Taylor expansion is made. Early
investigations of the finite-time-singularity problem for the 3D Euler
equation, starting from Taylor–Green initial data [12], used methods of
asymptotic analysis that were developed originally for the study of
singularities in critical phenomena [13], e.g., the ratio test [14] and Padé-
approximant methods [15]. These methods were used in conjunction with high-
precision arithmetic to study the analytic structure of Euler flows [12, 16,
17, 18], high-Reynolds-number turbulent flows [19], and related hydrodynamical
problems [20, 21, 22].

Over the years, direct-numerical-simulation (DNS) searches for possible
singular solutions of hydrodynamical PDEs have grown apace with advances in
high-performance computing and state-of-the-art numerical schemes have been
developed to capture the growth of singular solutions. Complex-space
singularities in solutions of such PDEs can be tracked by studying the
temporal evolution of spectra by using the analyticity-strip technique [23].
However, as the solution proceeds towards finite-time blowup, the
approximations of singular solutions, obtained by using finite-resolution and
finite-precision numerical methods, often suffer from loss of accuracy or
stability or both. This is because the flow fields develop increasingly small-
scale structures for times close to the time of blowup [6].

The fields are approximated by employing spectral methods in space; and the
governing PDEs are solved via finite-difference time-marching schemes [24,
25]. The resulting numerical approximation is the exact solution of the
spectrally truncated version of the governing PDE; but it deviates from the
true solution of the PDE as the solution approaches the time of singularity
(or potential singularity). These deviations first appear as localised,
oscillatory structures in real space; they have been called tygers. Tygers
were first reported in Fourier pseudospectral approximations of the one-
dimensional (1D) inviscid Burgers quation [26], where they were attributed to
nonlinear wave-particle resonances in regions of flow with positive strain.
The spontaneous emergence of tygers in real space is accompanied by complex-
space singularities crossing into the analyticity strip of width $\Delta_{x}$,
where $\Delta_{x}$ is the smallest grid spacing in the domain. Tygers have
also been reported in spectral studies of solutions of the 3D-axisymmetric,
wall-bounded, incompressible Euler PDE [2] and of the 3D incompressible Euler
PDE [27].

Recently, the complex-time analytic structure of solutions of the 1D inviscid
Burgers equation has been reinvestigated by using high-order time-Taylor-
series methods [28]. Unlike earlier work, which concentrates on the
singularity on the real time axis, Ref. 28 charts the landscape of
singularities in the complex-time plane, by investigating the convergence of
time-Taylor series of the 1D Burgers velocity field at each spatial point. For
single- and multi-mode sine-wave initial data, the authors of Ref. 28 find
several singularities, which are arranged in the form of an eye that is
centered at the origin. Furthermore, they find that the truncated time-Taylor-
series approximation of the velocity field shows spatially localised
oscillatory structures, called early-time resonances, that are reminiscent of
the tygers [26] mentioned above. These early-time resonances emerge at the
spatial points at which the convergence of the Taylor-series expansion is
first lost; and the corresponding complex-time singularities, situated off of
the real time axis, are the convergence-limiting singularities [28].

We investigate the complex-time analytic structure of solutions of the
3D-axisymmetric, wall-bounded, incompressible Euler equations, starting from
the initial data proposed in Ref. 1, to study a possible finite-time
singularity. We use our pseudospectral Fourier-Chebyshev method [2], with
quadruple-precision arithmetic to compute the time-Taylor series coefficients
of the flow fields, up until a high order. The resulting approximations
display early-time resonances; the initial spatial location of these
structures is different from that for the tygers, which we have obtained in
Ref. 2. We then perform asymptotic analysis of the Taylor-series coefficients,
by using the Darboux-type ratio method due to Mercer and Roberts [29], to
extract the location and nature of the convergence-limiting singularities. We
find that these singularities are distributed around the origin, in the
complex-$t^{2}$ plane, along two curves that resemble the shape of an eye. We
obtain similar results for the 1D wall-approximation, of the full
3D-axisymmetric Euler equation, called the 1D HL model [1, 30]. We use a
Fourier pseudospectral method [2], with quadruple-precision arithmetic, to
compute the time-Taylor-series coefficients of the flow fields. Our work
examines early-time resonances, in truncated time-Taylor expansions of
solutions of these PDEs, and compares them with tygers, in Galerkin-truncated
pseudospectral studies [26, 2].

The remaining part of this paper is organized as follows. For simplicity, we
first discuss the 1D HL model and obtain and analyse time-Taylor series for
its solutions in Sec. II.1. We then discuss and analyse the time-Taylor series
for solutions of the 3D-axisymmetric wall-bounded incompressible Euler
equations in Sec. II.2. We present our results for the 1D HL model in Sec.
III.1, including the complex-singularity landscape. We present similar results
for the 3D-axisymmetric wall-bounded incompressible Euler equations in Sec.
III.2. Finally, in Sec. IV, we conclude with a discussion of our results. In
Appendices A and B we give, respectively, a Mathematica code for the symbolic
computation of time-Taylor coefficients for the 1D HL model and some details
of the Mercer-Roberts method.

##  II Models and Methods

A 1D wall-approximation for the full 3D-axisymmetric incompressible Euler
(3DAE) equation was developed by Luo and Hou [1]. It is now called the 1D HL
model [30]. We discuss this 1D HL model first, in Sec. II.1.1, because it is
much simpler than the full 3DAE (see below). We construct the time-Taylor
series representation for solutions of the 1D HL model and derive recursion
relations for these time-Taylor coefficients in Sec. II.1.2. We describe the
Fourier pseudospectral methods and the numerical computation of the
coefficients for this model in Sec. II.1.3. In Section II.2.1, we describe the
3DAE; and we derive the recursion relations for the time-Taylor coefficients
of the series for the fields in Sec. II.2.2. Fourier-Chebyshev pseudospectral
methods, which we use for the computation of these coefficients, are presented
in Sec. II.2.3.

###  II.1 1D HL model

####  II.1.1 Model

The 1D HL model approximates the singular dynamics of the 3DAE [see Eq. (7)
below], at the wall at $r=1$, as discussed in Refs. 1, 30, 2. This model is

$\displaystyle\partial_{t}u+v\ \partial_{z}u$ $\displaystyle=0,$ (1a)
$\displaystyle\partial_{t}\omega+v\ \partial_{z}\omega$
$\displaystyle=\partial_{z}u,$ (1b) $\displaystyle\partial_{z}v$
$\displaystyle=\mathcal{H}[\omega]\,,$ (1c)

where $\mathcal{H}[.]$ is the Hilbert transform operator; $u$, $\omega$, and
$v$ are related, respectively, to the angular components of the velocity,
vorticity, and stream function in the 3DAE (see below) at $r=1$; Eq. (1c) is
the Biot-Savart-type law for this model. The domain is periodic in $z$ with
the periodicity length $L=2\pi$.

####  II.1.2 Time-Taylor series expansion and recursion relations

We write the time-Taylor expansions

$\displaystyle u(z,t)$ $\displaystyle=$
$\displaystyle\sum_{n=0}^{\infty}u_{n}(z)t^{n}\,,$ (2a)
$\displaystyle\omega(z,t)$ $\displaystyle=$
$\displaystyle\sum_{n=0}^{\infty}\omega_{n}(z)t^{n}\,,$ (2b) $\displaystyle
v(z,t)$ $\displaystyle=$ $\displaystyle\sum_{n=0}^{\infty}v_{n}(z)t^{n}\,,$
(2c)

where $u_{n}(z),\,\omega_{n}(z),\,$ and $v_{n}(z)$ are time-Taylor coefficient
of order $n$, which are functions of $z$. Hereafter, we suppress the
dependence on $z$ for notational simplicity. To construct recursion relations
for the time-Taylor coefficients, we substitute the series expansions
(2a)-(2c) in Eqs. (1a)-(1c). By comparing coefficients of terms of degree
$t^{n+1}$, we obtain

$\displaystyle u_{n}$ $\displaystyle=\frac{-1}{n}\sum_{l+m=n-1}v_{l}\
\partial_{z}u_{m},\;\;n>0\,,$ (3a) $\displaystyle\omega_{n}$
$\displaystyle=\frac{-1}{n}\sum_{l+m=n-1}v_{l}\
\partial_{z}\omega_{m}+\frac{1}{n}\partial_{z}u_{n-1},\;\;n>0\,,\;\;{\rm{and}}$
(3b) $\displaystyle v_{n}$
$\displaystyle=\int\mathcal{H}(\omega_{n}(z^{\prime}))dz^{\prime},\;\;n>0\,.$
(3c)

The initial conditions are $u_{0}(z)=u(z,t=0)$, $\omega_{0}(z)=\omega(z,t=0)$,
and $v_{0}(z)=\int\mathcal{H}(\omega_{0}(z^{\prime}))dz^{\prime}$.

Once the initial data are specified, we can use the recursion relations in
Eqs. (3) to derive closed-form expressions for $u_{n}(z),\ \omega_{n}(z)$ and
$v_{n}(z)$, analytically. However, the derivation of high-order coefficients
quickly becomes cumbersome, even if we use symbolic-computation software [see
Appendix A for details]. Therefore, we compute the time-Taylor coefficients on
a finite grid by using quadruple-precision Fourier pseudospectral numerical
methods.

####  II.1.3 Fourier pseudospectral methods

In the $L$-periodic domain, we define a uniform collocation grid of $N$ points
given by $\mathbf{X}_{N}=\\{z_{j}=j\Delta z\ :j=0,1,...,N-1\\}$ where $\Delta
z=L/N$. The Fourier pseudospectral projection operation $P_{N}$ is:

$\displaystyle P_{N}u(z,t)$ $\displaystyle=\sum_{|k|\leq
N/2}\hat{u}(k,t)e^{i\frac{2\pi}{L}kz}\,;$ (4a) $\displaystyle\hat{u}(k,t)$
$\displaystyle=\frac{1}{N}\sum_{j=0}^{N-1}u(z_{j},t)e^{-i\frac{2\pi}{L}kz_{j}}\,.$
(4b)

Given the initial data $u_{0}(z),\ \omega_{0}(z)$ and $v_{0}(z)$, we first
perform a pseudospectral projection to obtain $P_{N}u_{0}$, $P_{N}\omega_{0}$,
and $P_{N}v_{0}$ by using the FFTW3 library [31]. To compute
$u_{n}(\mathbf{X}_{N})\coloneqq\\{u_{n}(z):z\in\mathbf{X}_{N}\\}$, we set
$n=1$ in Eq. (3a). We evaluate the derivative $\partial_{z}u_{n}$ on the
right-hand side (RHS) of Eq. (3a) in Fourier spectral space. The nonlinear
terms are computed in real space. Similarly, we evaluate
$\omega_{1}(\mathbf{X}_{N})$ from Eq. (3b). We use
$\omega_{1}(\mathbf{X}_{N})$ to compute $v_{1}(\mathbf{X}_{N})$ from Eq. (3c);
the integration is performed in Fourier spectral space. The values of
$u_{1}(\mathbf{X}_{N})$, $\omega_{1}(\mathbf{X}_{N})$, and
$v_{1}(\mathbf{X}_{N})$ can then be used to compute the values of
$u_{2}(\mathbf{X}_{N})$, $\omega_{2}(\mathbf{X}_{N})$, and
$v_{2}(\mathbf{X}_{N})$. We repeat this procedure for $n\geq 2$.

###  II.2 3D-axisymmetric wall-bounded incompressible Euler equation

####  II.2.1 Model

We now consider the 3D Euler equation in the vorticity–stream function
formulation:

$\displaystyle\mbox{\boldmath$\omega$}_{t}+\mbox{\boldmath$u$}\cdot\nabla\mbox{\boldmath$\omega$}=\mbox{\boldmath$\omega$}\cdot\nabla\mbox{\boldmath$u$}\,,$
(5)

where $\mbox{\boldmath$\omega$}=\nabla\times\mbox{\boldmath$u$}$ is the
vorticity, $\mbox{\boldmath$u$}=\nabla\times\mbox{\boldmath$\psi$}$ is the
velocity field, and $\psi$ is the vector-valued stream function. The Poisson
equation $\mbox{\boldmath$\omega$}=-\nabla^{2}\mbox{\boldmath$\psi$}$ gives
the relation between the vorticity and the stream function.

To represent axisymmetric fields, we use $\mbox{\boldmath$u$}(r,z)=u^{r}(r,z)\
\mathbi{\hat{e}_{r}}+u^{\theta}(r,z)\ \mathbi{\hat{e}_{\theta}}+u^{z}(r,z)\
\mathbi{\hat{e}_{z}}$, where $\ \mathbi{\hat{e}_{r}},\
\mathbi{\hat{e}_{\theta}}$, and $\ \mathbi{\hat{e}_{z}}$ are unit vectors in
the cylindrical coordinate system. Furthermore, we define new variables in
terms of the angular components $u^{\theta}$, $\omega^{\theta}$, and
$\psi^{\theta}$:

$u^{1}=u^{\theta}/r;\qquad\omega^{1}=\omega^{\theta}/r;\qquad\psi^{1}=\psi^{\theta}/r\,.$
(6)

The 3D-axisymmetric wall-bounded incompressible Euler (3DAE) equations are
then given in terms of the new variables (6) by

$\displaystyle u^{1}_{t}+u^{r}u^{1}_{r}+u^{z}u^{1}_{z}$
$\displaystyle=2u^{1}\psi_{z}^{1},$ (7a)
$\displaystyle\omega_{t}^{1}+u^{r}\omega_{r}^{1}+u^{z}\omega_{z}^{1}$
$\displaystyle=((u^{1})^{2})_{z},$ (7b)
$\displaystyle-\Big{(}\partial_{r}^{2}+\frac{3}{r}\partial_{r}+\partial_{z}^{2}\Big{)}\psi^{1}$
$\displaystyle=\omega^{1},$ (7c)

with $u^{r}=-r\psi^{1}_{z}$ and $u^{z}=2\psi^{1}+r\psi^{1}_{r}$. For a
detailed derivation, see Refs. 32, 1. So long as the solutions to Eq. (7) are
smooth [$C^{\infty}(\mathbf{R}\times\bar{\mathbf{R}}^{+})$],
$u^{\theta},\omega^{\theta}$, and $\psi^{\theta}$ all vanish at $r=0$; this
ensures that the coordinate singularity at $r=0$ does not enter the evolution
equations (7) directly [33]. Here, $\mathbf{R}$ represents the set of real
numbers and $\bar{\mathbf{R}}^{+}$ the set of affinely extended positive real
numbers.

We solve Eq. (7) in the domain $D(1,L)=\\{(r,z):0\leq r\leq 1,0\leq z\leq
L\\}$, with $L$-periodic boundary conditions along $\hat{e}_{z}$, the no-flow
conditions at $r=1$ (8a), and the following pole conditions at $r=0$ (8b):

$\displaystyle\psi^{1}(r=1,z,t)$ $\displaystyle=$ $\displaystyle 0\,;$ (8a)
$\displaystyle u^{1}_{r}(r=0,z,t)=\omega^{1}_{r}(r=0,z,t)$ $\displaystyle=$
$\displaystyle\psi^{1}_{r}(r=0,z,t)=0\,.$ (8b)

The 3D-axisymmetric field variables in Eq. (7) at $r=1$ can be approximated by
the 1D-HL-model (1) fields in Eq. (1) as discussed in Ref. 1:

$u(z)\rightarrow(u^{1})^{2}(r=1,z),\qquad\omega(z)\rightarrow\omega^{1}(r=1,z),\qquad
v(z)\rightarrow\partial_{r}\psi^{1}(r=1,z).$ (9)

####  II.2.2 Time-Taylor series expansion and recursion relations

We write the time-Taylor series expansions for the fields $u^{1}(r,z,t),\
\omega^{1}(r,z,t)$, and $u^{z}(r,z,t)$ as

$\displaystyle u^{1}(r,z,t)=\sum_{n=0}^{\infty}u^{1}_{n}(r,z)\
t^{n}\,,\qquad\omega^{1}(r,z,t)=\sum_{n=0}^{\infty}\omega^{1}_{n}(r,z)\
t^{n}\,,\;\;{\rm{and}}\qquad\psi^{1}(r,z,t)=\sum_{n=0}^{\infty}\psi^{1}_{n}(r,z)\
t^{n}\,.$ (10)

Similarly, we write $u^{r}(r,z,t)=\sum_{n=0}^{\infty}u^{r}_{n}(r,z,n)t^{n}$
and $u^{z}(r,z,t)=\sum_{n=0}^{\infty}u^{z}_{n}(r,z,n)t^{n}$ for the radial and
axial velocities, respectively. By substituting these series expansions for
the fields in Eqs. (7a)-(7c) and equating the coefficients of $t^{n+1}$ on
both sides, we obtain:

$\displaystyle
u^{1}_{n}(r,z)=-\frac{1}{n}\sum_{l+m=n-1}\Big{(}u^{r}_{l}\partial_{r}u^{1}_{m}+u^{z}_{l}\partial_{z}u^{1}_{m}-2u^{1}_{l}\partial_{z}\psi^{1}_{m}\Big{)}\,;$
(11a)
$\displaystyle\omega^{1}_{n}(r,z)=-\frac{1}{n}\sum_{l+m=n-1}\Big{(}u^{r}_{l}\partial_{r}\omega^{1}_{m}+u^{z}_{l}\partial_{z}\omega^{1}_{m}-\partial_{z}(u^{1}_{l}u^{1}_{m})\Big{)}\,.$
(11b)

Once we determine $\omega^{1}_{n}(r,z)$, we use the Poisson equation (7c) to
obtain $\psi^{1}_{n}(r,z)$ as follows:

$\displaystyle-\left(\partial_{r}^{2}+\frac{3}{r}\partial_{r}+\partial_{z}^{2}\right)\psi^{1}_{n}$
$\displaystyle=\omega^{1}_{n}\,;$ (11c)
$\displaystyle\partial_{r}\psi^{1}_{n}(r=0,z)=0\,;\qquad$
$\displaystyle\qquad\psi^{1}_{n}(r=1,z)=0\,.$ (11d)

Here, the boundary conditions in Eqs. (8a)-(8b) are applied to
$\psi^{1}_{n}(r,z)$ in Eq. (11d). We can then obtain the time-Taylor
coefficients for the radial and axial velocity fields:

$\displaystyle u^{r}_{n}$ $\displaystyle=(-r\partial_{z})\psi^{1}_{n}\,;$
(11e) $\displaystyle u^{z}_{n}$
$\displaystyle=(2+r\partial_{r})\psi^{1}_{n}\,.$ (11f)

The initial data are used to define $u^{1}_{0},\,\omega^{1}_{0}$ and
$\psi^{1}_{0}$ [cf. Subsec. II.1.2 for the 1D HL model]. For $n=1$,
$u^{1}_{1}(r,z)$ and $\omega^{1}_{1}(r,z)$ can be obtained analytically, as
closed-form expressions, by using Eqs. (11a)-(11b). However, unlike the
Biot–Savart law for the 1D HL model in Eq. (1c), the Poisson equation in Eqs.
(11c)-(11d) cannot be solved analytically for $\psi^{1}_{1}(r,z)$ via symbolic
computation. Thus, we use our Fourier-Chebyshev pseudospectral scheme [2] with
quadruple-precision arithmetic. We describe our methods briefly in the next
Subsection II.2.3 [for details see Ref. 2].

####  II.2.3 Fourier-Chebyshev pseudospectral method for the computation of
time-Taylor coefficients

We discretize the domain by using the collocation grid
$\mathbf{X}_{N,M}\coloneqq\\{(r_{i},z_{j}):i=1,...,N\text{ and
}j=0,...,M-1\\}$, where the axial nodes
$\\{z_{j}=\tfrac{Lj}{M}\\}_{j=0,...M-1}$ are uniformly spaced. The radial
nodes $\\{r_{i}=\frac{1}{2}(1+\cos(\frac{\pi(i-0.5)}{N})\\}_{i=1,...N}$, the
roots of the highest-order Chebyshev polynomial $T_{N}(2r-1)$ in the spectral
basis, are non-uniformly spaced and cluster near $r=0$ and $r=1$. The Fourier-
Chebyshev pseudospectral approximation for the field $u^{1}(r,z,t)$ is then
given by

$\displaystyle
P_{N,M}u^{1}(r,z,t)=\sum_{|k|<M/2}\sum_{l=0}^{N}\hat{u}^{1}(k,l,t)\ e^{ikz}\
T_{l}(2r-1)\,,$ (12)

where $T_{l}(x)$ is the order-$l$ Chebyshev polynomial of the first kind.

We compute the derivative terms in spectral space. The Poisson problem in Eqs.
(11c)-(11d) is solved by using the Fourier-Chebyshev Tau Poisson solver
described in Refs. 2, 34. The nonlinear terms are computed in real space.
Given the initial data $u^{1}_{0}(\mathbf{X}_{N,M}),\
\omega^{1}_{0}(\mathbf{X}_{N,M})$ and $\psi^{1}_{0}(\mathbf{X}_{N,M})$, we
compute $u^{1}_{1}(\mathbf{X}_{N,M})$ and $\omega^{1}_{1}(\mathbf{X}_{N,M})$
from Eq. (11a)-(11b). We then solve the Poisson problem (11c)-(11d) for
$\psi^{1}_{1}(\mathbf{X}_{N,M})$. We repeat this process to obtain the time-
Taylor coefficients for $n=2,3,...,N_{t}$, where $N_{t}$ is the order at which
we truncate the time-Taylor expansion in Eq. (10). The choice of $N_{t}$ plays
a crucial role in ensuring that the results of our asymptotic analysis are
accurate [see Sec. III.2 for details]: time-Taylor coefficients, computed for
$n>N_{t}$, suffer from errors that arise because of (a) fixed-precision
arithmetic and (b) the finite number of modes retained in any implementation
of the pseudospectral method.

##  III Results

We compute and examine the time-Taylor coefficients for special choices of
initial data that lead to (possible) finite-time singularities [2]. We then
perform asymptotic analysis [35, 29, 14] of these coefficients to estimate the
position and the nature of the convergence-limiting singularities. In Section
III.1, we chart out the pattern of these singularities for a singular solution
of the 1D HL model (1). In Section III.2, we repeat the analysis for the
potentially singular solution of the 3DAE equation (7) with the initial
condition used in Refs. 1, 2, 36, 37.

###  III.1 1D HL-model

We start with the following initial data:

$\displaystyle u_{0}(z)=\sin^{2}(z)\,;\qquad\omega_{0}(z)=0.$ (13)

It has been proved [30] that the flow, which evolves from Eq. (13), develops a
finite-time singularity. For a scaled version of the initial data used in Ref.
2, the time of singularity lies near $t_{*}\simeq 0.0035$.

$\begin{array}[]{|c|c|c|c|}\hline\cr{n}&u_{n}(z)&{\omega_{n}}(z)&{v}_{n}(z)\\\
\hline\cr&&&\\\ 0&\sin^{2}(z)&0&0\\\ 1&0&\sin(2z)&\sin(z)(-\cos(z))\\\
2&\sin^{2}(z)\cos^{2}(z)&0&0\\\
3&0&\frac{1}{3}\sin(4z)&-\frac{1}{12}\sin(4z)\\\
4&\frac{1}{12}\sin(2z)\sin(4z)&0&0\\\
5&0&\frac{1}{15}(2\sin(6z)-\sin(2z))&\frac{1}{90}(3\sin(2z)-2\sin(6z))\\\
6&\frac{1}{270}\sin^{2}(2z)(17\cos(4z)+7)&0&0\\\
7&0&\frac{109\sin(8z)-86\sin(4z)}{1890}&\frac{172\sin(4z)-109\sin(8z)}{15120}\\\
.&.&.&.\\\ .&.&.&.\\\ \hline\cr\end{array}$

Table 1: Closed-form symbolic expressions of the time-Taylor coefficients
$u_{n}(z),\omega_{n}(z)$, and $v_{n}(z)$ for low-order $n\leq 7$, obtained by
using the recursion relations in Eq. (3) and starting from the initial data in
Eq. (13) [see Appendix A for details].

We now apply the recursion relations in Eqs. (3a)- (3c), starting from the
initial data in Eq. (13). We use symbolic computation to compute the closed-
form expressions, for low-order ($n\leq 7$) coefficients $u_{n}(z),\
\omega_{n}(z)$ and $v_{n}(z)$, that we list in TABLE 1. Appendix A contains
our Mathematica code for the determination of these coefficients. The odd-
order coefficients for $u(z)$ and the even-order coefficients for $\omega(z)$
and $v(z)$ vanish identically because of the symmetry of the initial data
(13). The difficulty of symbolic computation and the storage of high-order
coefficients increases rapidly with the order $n$. Therefore, we use the
Fourier pseudospectral methods outlined in Sec. II.1.3 to compute the time-
Taylor coefficients up to order $n\leq N_{t}=100$, on a uniform grid with
$N=256$ points. For accuracy, it is important to use quadruple-precision
arithmetic. In FIG. 1, we present a sign-coded heat map of the absolute values
of the coefficients $u_{n}(z),\ \omega_{n}(z),$ and $\psi_{n}(z)$ for $n\leq
100$ evaluated at the collocation points $z\in\mathbf{X}_{N=256}$; the
coefficients with positive (negative) signs are shown in red (blue); the color
bars use a log scale (base 10). White regions appear in these plots where the
coefficients vanish identically [e.g., for odd orders in $u_{n}(z)$] or fall
below $10^{-15}$.

(0,0) node[inner sep=0]; (-0.4,1.74) node (a); (0,0) node[inner sep=0];
(-0.4,1.74) node (b); (0,0) node[inner sep=0]; (-1.35,1.74) node (c);

Figure 1:  Sign-coded heat maps of the absolute values of the time-Taylor
coefficients (a) $u_{n}(z)$, (b) $\omega_{n}(z)$, and (c) $v_{n}(z)$, for
$n\leq 100$ evaluated at $z\in\mathbf{X}_{N=256}$. The coefficients with
positive (negative) signs are shown in red (blue); $u_{n}(z)$ vanishes
identically for odd $n$ (white bands); $\omega_{n}(z)$ and $v_{n}(z)$ vanish
identically for even $n$, because of the symmetry in the initial condition. We
use a log scale (base 10) for the color bar.

In FIG. 2, we plot the truncated series approximations
$u_{N_{t}}\coloneqq\sum_{n=0}^{N_{t}}{u}_{n}t^{n}$,
$\omega_{N_{t}}\coloneqq\sum_{n=0}^{N_{t}}{\omega}_{n}t^{n}$ and
$v_{N_{t}}\coloneqq\sum_{n=0}^{N_{t}}{v}_{n}t^{n}$ at time $t=1.13$, for
$N_{t}=10\text{ (yellow)},30$ (blue), and $100$ (red). We observe the
emergence of localised oscillatory structures at
$z=\\{\tfrac{\pi}{4},\tfrac{3\pi}{4},\tfrac{5\pi}{4},\tfrac{7\pi}{4}\\}$, for
all the fields; we refer to these structures as early-time resonances, a term
introduced in Ref. 28 for the 1D inviscid Burgers equation in a periodic
domain. Similar oscillatory structures were also reported for the Fourier-
pseudospectral study of the 1D HL model in Ref. 2, near the time of
singularity; they were identified as nonlinear wave-particle resonances called
tygers. However, the spatial locations of the tygers, observed in Ref. 2 for
$u(z,t)$ and $\omega(z,t)$, are different from what we see in FIG. 2: In
particular, for the initial condition (13) the Fourier-pseudospectral study in
Ref. 2 yields tygers at $z=\\{\tfrac{\pi}{2},\tfrac{3\pi}{2}\\}$ for $u(z,t)$
and $\omega(z,t)$. This dissimilarity between early-time resonances and
tygers, which we observe here, generalises significantly its first observation
in Ref. 28 for the 1D inviscid Burgers equation.

(0,0) node[inner sep=0]; (-3.1,-0.8) node (a); (2.1,-0.8) node (b); (7.5,-0.8)
node (c);

Figure 2: Plots vs. $z$ of the truncated time-Taylor series approximations for
the fields (a) $u_{N_{t}}(z,t)$, (b) $\omega_{N_{t}}(z,t)$, and (c)
$v_{N_{t}}(z,t)$ at time $t=1.13$, summed for $N_{t}=10\text{ (yellow)},\ 30$
(blue) and $100$ (red). Early-time resonances, which occur in all three
fields, are oscillatory structures localised around
$z=\\{\tfrac{\pi}{4},\tfrac{3\pi}{4}\\}$; their amplitude is largest for
$\omega_{N_{t}}$.

In FIG. 2, the emergence of early-time resonances signals the loss of
convergence of the Taylor-series approximation. For an analytic function
$u(t)$, the convergence of the Taylor series about a point $t=t_{0}$ is
limited by the nearest singularity in the complex-$t$ plane. If $R$ is the
distance of the nearest singularity from the point of expansion $t_{0}$, the
Taylor series loses convergence for $|t-t_{0}|>R$. The nature and position of
this singularity (or singularities) governs the large-$n$ behaviour of the
Taylor coefficients. Conversely, we can use the large-$n$ behaviour of
$u_{n}(z),\ \omega_{n}(z)$, and $v_{n}(z)$ in FIG. 1 to chart the positions of
complex-time singularities.

Given our initial data (13), the time-Taylor series expansion for $u(z,t)$
[TABLE 1] contains only even-order terms, so we use $q\equiv t^{2}$ as the
expansion variable. In contrast, the expansions for $\omega(z,t)$ and $v(z,t)$
do not have even-order terms in $t$; this allows us to rewrite their series in
terms of $q\equiv t^{2}$, with a prefactor of $t=q^{1/2}$ as follows:

$\displaystyle u_{N_{t}}(z,q)$ $\displaystyle=$
$\displaystyle\sum_{n=0}^{N_{t}/2}u_{2n}(z){q}^{n}\,;\;\;$ (14a)
$\displaystyle\omega_{N_{t}}(z,q)$ $\displaystyle=$ $\displaystyle q^{1/2}\
\sum_{n=0}^{N_{t}/2}\omega_{2n+1}(z){q}^{n}\,;\;\;$ (14b) $\displaystyle
v_{N_{t}}(z,q)$ $\displaystyle=$ $\displaystyle q^{1/2}\
\sum_{n=0}^{N_{t}/2}v_{2n+1}(z){q}^{n}\,.$ (14c)

Furthermore, in FIG. 1 we observe that $u_{2n}(z=0)$ and $u_{2n}(z=\pi)$
display alternating signs with a long period; the period over which the sign
of the coefficient alternates is reduced as $z\rightarrow\\{\pi/2\\}$. Given
these alternating signs, we cannot use the standard Domb–Sykes ratio method
[14, 35] for determining the radius of convergence of the series
$u_{N_{t}}(z,t)$. This method can only be applied when the coefficients are
all of the same sign or have strictly alternating signs. When the signs
alternate with a period greater than $1$, the sign of the ratio
$u_{2(n+1)}(z)/u_{2n}(z)$ changes with $n$, so $u_{2(n+1)}(z)/u_{2n}(z)$ does
not have a well-defined limit as $n\to\infty$.

Mercer and Roberts [29] (MR) generalised the Domb–Sykes method for series with
coefficients whose signs change with order $n$ and have a period of
alternation greater than $1$. In the MR method, the pattern of alternating
signs is attributed to a pair of complex-conjugate singularities that lie
close to the real domain [see Appendix B]. For a given value of $z$, the
Taylor series $u_{N_{t}}(z,q)$ is approximated by the MR-model function
$\mathfrak{u}(z,\mathbf{q})$, in complex-$\mathbf{q}$ space, where
$\mathbf{q}$ is the complexified version of the variable $q=t^{2}$:

$\displaystyle\mathfrak{u}(z,\mathbf{q})=\left(1-\frac{\mathbf{q}}{\mathbf{q_{*}}}\right)^{\nu(z)}+\left(1-\frac{\mathbf{q}}{\overline{\mathbf{q}}_{*}}\right)^{\nu(z)}\,;\qquad\mathbf{q_{*}}(z)\coloneqq
R(z)e^{i\theta(z)}.$ (15)

Here, $\mathbf{q}_{*}$ denotes the position of the convergence-limiting
complex singularity and $\overline{\mathbf{q}}_{*}$ its complex conjugate. The
unknown functions $R(z)>0,\ \theta(z)\in[0,\pi]$, with $\nu(z)\in\mathbb{R}$
neither zero nor a positive integer. For each value of $z\in\mathbf{X}_{N}$,
we construct the MR coefficients $B_{k}^{2}(z)$ given in Eq. (16a). We then
perform a nonlinear fit for $B^{2}_{k}$ with the functional forms given in Eq.
(16b) below:

$\displaystyle B^{2}_{k}(z)$ $\displaystyle=\frac{u_{k+1}(z)\cdot
u_{k-1}(z)-u^{2}_{k}(z)}{u_{k}(z)\cdot u_{k-2}(z)-u^{2}_{k-1}(z)}\,;\qquad
2\leq k\leq(N_{t}-1)\,;$ (16a) $\displaystyle u_{n}(z)$
$\displaystyle=2(-1)^{n}\binom{\nu(z)}{n}R(z)^{-n}\cos(n\theta(z))\,;\qquad
0\leq n\leq N_{t}\,.$ (16b)

From these we obtain estimates for the functions $R(z),\ \theta(z),$ and
$\nu(z)$ [FIGs. 3 (a)-(c), respectively], for $z\in\mathbf{X}_{256}$, by using
the LMFIT package [38] in Python. In Appendix B, we discuss the $\tfrac{1}{k}$
linear fit for $B_{k}(z)$; this linear fit works for
$z=\\{0,\,\tfrac{\pi}{8}\,\tfrac{\pi}{4}\\}$; for other values of $z$ we must,
perforce, use the nonlinear fit provided by LMFIT. The estimates from Eqs.
(14a), (14b), and (14c) are shown, respectively, in yellow, pink, and blue in
FIGs. 3 (a)-(c).

(0,0) node[inner sep=0]; (-1.1,1.45) node (a); (0,0) node[inner sep=0];
(-1.2,1.45) node (b); (0,0) node[inner sep=0]; (-1.1,1.45) node (c);

Figure 3: Plots vs. $z$ of the estimates of (a) $R$, (b) $\theta$, and (c)
$\nu$, obtained from $u_{n}\text{(orange)},\ \omega_{n}$(pink) and
$v_{n}$(blue), for $z\in\mathbf{X}_{256}$ and $n\leq N_{t}=100$, by using the
MR method. A full nonlinear fit of $B_{k}^{2}$ in Eq. (16a) is performed by
using the LMFIT package in Python.

The convergence-limiting singularities for $z=\\{0,\pi\\}$ occur at
$\theta(0)=\theta(\pi)=0$; and for $z=\tfrac{\pi}{2}$ they occur at
$\theta(\tfrac{\pi}{2})=\pi$ [panel (b) of FIG. 3]; i.e., these singularities
are on the positive and negative real axes in the complex-$\mathbf{q}$ plane,
respectively. They lie farthest from the origin because the radius of
convergence $R(z)$ assumes its maximum value at these points [panel (a) of
FIG. 3]; $R(z)$ assumes its minimum value $R_{m}$ at
$z=\\{\tfrac{\pi}{4},\tfrac{3\pi}{4}\\}$; here, $\theta(z)=\tfrac{\pi}{2}$.
Thus, singularities for $z=\\{\tfrac{\pi}{4},\tfrac{3\pi}{4}\\}$ are closest
to the origin and are positioned on the imaginary axis $\Re(\mathbf{q})=0$. In
particular, when $t^{2}>R_{m}$, where
$R_{m}=R(z=\\{\tfrac{\pi}{4},\tfrac{3\pi}{4}\\})$, the fields in FIG. 2
develop early-time resonances. The minimum radius of convergence $R_{m}\simeq
1.17$ or $t_{m}\simeq 1.08$ for all fields [see FIG. 3 and FIG. 4]; however,
early-time resonances appear more prominently in $\omega_{N_{t}}$ than in
$u_{N_{t}}$ and $v_{N_{t}}$ [see FIG. 2].

In FIG. 4, we plot $(R(z),\theta(z))$, for different values of $z$ in the
complex $\mathbf{q}$ plane, by using polar coordinates. These singularities
form of an approximately elliptical eye, which is centered at the origin and
compressed along the $\Im(\mathbf{q})$ axis. The arrangement of such complex
singularities in the form of an eye was first reported for time-Taylor
expansions of the 1D inviscid Burgers equation [28], where the eye occurs in
the complex-$t$ plane. Given the symmetries of initial conditions that we use,
these eyes occur in the complex-$\mathbf{q}$ plane [see FIG. 4]. In panel (c)
of FIG. 3, we plot the exponent $\nu(z)$ of the convergence-limiting
singularities for all the fields; here, we observe that $\nu$ takes non-
integer values that depend on $z$. This $z$ dependence is significant only
near $z=\\{0,\tfrac{\pi}{2},\pi\\}$; otherwise $\nu\simeq\tfrac{1}{2}$ for
$u$, and $\nu=-\tfrac{1}{2}$ for $\omega$ and $v$. We cannot obtain accurate
estimates for singularities that lie arbitrarily close to the real axis in
FIG. 4. This is because the MR method is effective only when the series has
settled down to its final regular behaviour. Given finite-resolution and
finite-precision errors, the truncated Taylor series that we use in Eq. (14)
does not contain enough terms to capture this behaviour. [See Appendix B for
details.]

Figure 4: Polar plot, in the complex-$\mathbf{q}$ plane, of $R(z)$ vs.
$\theta(z)$, depicting the positions of convergence-limiting complex
singularities $\mathbf{q}_{*}=Re^{i\theta}$, as well as their complex-
conjugates, for $u$ (yellow), $\omega$ (pink), and $v$ (blue) at
$z\in\mathbf{X}_{N=256}$. The singularities are arranged in the shape of an
eye, which is squashed along the $\Im(\mathbf{q})$ axis and is centered at the
origin.

###  III.2 3D-axisymmetric wall-bounded incompressible Euler equation

We use the singular initial condition proposed by Ref. 1 for this model to
obtain

$\displaystyle
u^{1}_{0}(r,z)=e^{-30(1-r^{2})^{4}}\sin{\left(\tfrac{2\pi}{L}z\right)}\,,$
(17a)
$\displaystyle\omega^{1}_{0}(r,z)=0\,,\qquad{\rm{and}}\qquad\psi^{1}_{0}(r,z)=0\,,$
(17b)

where $L=2\pi$ 111Note that Refs. 1, 2 use $L=1/6$.. For every point
$(r,z)\in\mathbf{X}_{N,M}$, we use our Fourier-Chebyshev pseudospectral method
[2] to compute $u^{1}_{n}(\mathbf{X}_{N,M})$ and
$\omega^{1}_{n}(\mathbf{X}_{N,M})$ via the recursion relations (11a) and
(11b), respectively. We then solve the Poisson problem (11c)-(11d) for
$\psi^{1}_{n}(\mathbf{X}_{N,M})$ by using the Tau Poisson solver in Fourier-
Chebyshev spectral space [2]. The time-Taylor coefficients for the radial and
axial velocities, $u^{r}_{n}(\mathbf{X}_{N,M})$ and
$u^{z}_{n}(\mathbf{X}_{N,M})$, are then computed from
$\psi^{1}_{n}(\mathbf{X}_{N,M})$ [see Sec. II.2.2].

(0,0) node[inner sep=0]; (-0.8,1.66) node (a); (0,0) node[inner sep=0];
(-0.8,1.66) node (b); (0,0) node[inner sep=0]; (-1.65,1.66) node (c);

Figure 5:  Sign-coded heat maps of the absolute values of the time-Taylor
coefficients (a) $u^{1}_{n}(r_{N},z)$, (b) $\omega^{1}_{n}(r_{N},z)$, and (c)
$u^{z}_{n}(r_{N},z)$, for $n\leq N_{t}=86$ and evaluated at $r_{N},\
z\in\mathbf{X}_{N,M}$, where $(N=128,M=256)$. The coefficients with positive
(negative) signs are shown in red (blue); $u^{1}_{n}(r_{N},z)$ vanish
identically for odd $n$ (white bands); and $\omega^{1}_{n}(r_{N},z)$ and
$u^{z}(r_{N},z)$ vanish identically for even $n$, because of the symmetry of
the initial data. We use a log scale (base 10) for the color bar.

In FIG. 5, we present the sign-coded heat map of the absolute values of the
coefficients $u^{1}_{n},\,\omega^{1}_{n}$ and $\psi^{1}_{n}$ for $n\leq 86$
evaluated at the collocation points $(r_{N},z)\in\mathbf{X}_{N,M}$ where
$(N=128,M=256)$; $r_{N}$ is the Chebyshev node that lies closest to the
boundary at $r=1$. As in FIG. 1, the coefficients with positive (negative)
signs are shown in red (blue), with log scale (base 10) color bars; the
coefficients vanish identically [e.g., for odd orders in $u^{1}_{n}$] or fall
below $10^{-40}$ in white regions. The odd-order coefficients for $u^{1}_{n}$
and the even-order coefficients for $\omega^{1}_{n}$ and $u^{z}_{n}$ vanish
identically because of the symmetry of the initial condition (17).

In FIG. 6, we plot, as a function of $(r,z)$, the time-Taylor series truncated
at $N_{t}=86$ and evaluated at time $t=2.6$ for the fields (a)
$u^{1}_{N_{t}}(\mathbf{X}_{N,M})$, (b) $\omega^{1}_{N_{t}}(\mathbf{X}_{N,M})$,
and (c) $u^{z}_{N_{t}}(\mathbf{X}_{N,M})$. Early-time resonances emerge at the
wall ($r=1$) in all fields. In FIG. 7, we plot these fields versus $z$, at
$r=r_{N}$. These plots show clearly the development of oscillations at
$z=\\{\tfrac{\pi}{4},\tfrac{3\pi}{4}\\}$ (in FIG. 6, also at
$\\{\tfrac{5pi}{4},\tfrac{7\pi}{4}\\}$). In FIG. 4. of Ref. 2 we have reported
the emergence of tygers, $z=\\{\tfrac{\pi}{2},\tfrac{3\pi}{2}\\}$ in all
fields in this model, while using traditional Fourier-Chebyshev pseudospectral
methods, for these initial data. We again emphasize that (as in the 1D HL
model) these tygers are different from the early-time resonances in FIG. 7.

(0,0) node[inner sep=0]; (-3.7,1.4) node (a); (1.4,1.4) node (b); (6.8,1.4)
node (c);

Figure 6:  Surface plots vs. $(r,z)$ of the time-Taylor series approximations
of the fields (a) $u^{1}$, (b) $\omega^{1}$, and (c) $u^{z}$, truncated at
$N_{t}=86$ and evaluated at time $t=2.6$ on a Fourier-Chebyshev collocation
grid $\mathbf{X}_{N,M}$, where $(N=128,M=256)$. We see the development of
early-time resonances in all fields near the wall at $r=1$. These resonant
oscillations are localised around
$z=\\{\tfrac{\pi}{4},\tfrac{3\pi}{4},\tfrac{5\pi}{4},\tfrac{7\pi}{4}\\}$; they
grow in amplitude and spread outwards in time.

(0,0) node[inner sep=0]; (-6.2,-0.8) node (a); (-1,-0.8) node (b); (4.3,-0.8)
node (c);

Figure 7: Plots vs. $z$ of the truncated time-Taylor series approximations for
the fields (a) $u^{1}_{N_{t}}(r_{N},z,t)$, (b)
$\omega^{1}_{N_{t}}(r_{N},z,t)$, and (c) $u^{z}_{N_{t}}(r_{N},z,t)$, evaluated
at time $t=2.6$ and summed for three orders of truncation, namely,
$N_{t}=10\text{ (yellow)},\ 30$ (blue) and $100$ (red). Early-time resonances
occur in all three fields as oscillatory structures localised around
$z=\\{\tfrac{\pi}{4},\tfrac{3\pi}{4}\\}$. The oscillations have the largest
amplitude for $\omega^{1}_{N_{t}}$.

Given that the even-order coefficients for $u^{1}_{n}(r_{N},z)$ and the odd-
order coefficients for $\omega^{1}_{n}(r_{N},z),u^{z}_{n}(r_{N},z)$ are
identically zero [FIG. 5], the time-Taylor series can be rewritten for the
variable $q=t^{2}$ [cf. Eq. (14) for the 1D HL model]:

$\displaystyle u^{1}_{N_{t}}(r_{N},z,q$
$\displaystyle)=\sum_{n=0}^{N_{t}/2}u^{1}_{2n}(r_{N},z){q}^{n}\,;$ (18a)
$\displaystyle\omega^{1}_{N_{t}}(r_{N},z,q)=q^{1/2}\
\sum_{n=0}^{N_{t}/2}\omega^{1}_{2n+1}(r_{N},z){q}^{n}\,;$ $\displaystyle\qquad
u^{z}_{N_{t}}(r_{N},z,q)=q^{1/2}\
\sum_{n=0}^{N_{t}/2}u^{z}_{2n+1}(r_{N},z){q}^{n}\,.$ (18b)

For a given field, we also note that the non-zero time-Taylor coefficients,
computed at a point $z$ near the wall $r_{N}$ by using the initial data
((17)), have a pattern of alternating signs [cf. Sec. III.1 for the 1D HL
model]. Thus, we can use the MR methods as we did for the 1D HL model in Eq.
(15). For each value of $z\in\mathbf{X}_{N,M}$ near the wall at $r_{N}$, we
construct $B_{k}^{2}(z)$ as in Eq. (16a). Our estimates for $R,\ \theta$ and
$\nu$, which we obtain by using the LMFIT package [38] in Python as for the 1D
HL model, are shown in panels (a)-(c) of FIG. 8, where we superpose the
estimates for $u^{1}_{N_{t}}(r_{N},z_{j})$ (yellow),
$\omega^{1}_{N_{t}}(r_{N},z_{j})$ (pink) and $u^{z}_{N_{t}}(r_{N},z_{j})$
(blue).

(0,0) node[inner sep=0]; (-1.2,1.45) node (a); (0,0) node[inner sep=0];
(-1.3,1.45) node (b); (0,0) node[inner sep=0]; (-1.1,1.45) node (c);

Figure 8: Plots vs. $z$ of the estimates of (a) $R$, (b) $\theta$, and (c)
$\nu$, obtained from $u^{1}_{n}(r_{N},z)\text{(orange)},\
\omega^{1}_{n}(r_{N},z)$(pink) and $u^{z}_{n}(r_{N},z)$(blue), for
$z\in\mathbf{X}_{N,M}$, where $(N=128,M=256)$ and $n\leq N_{t}=86$, using the
MR method. A full nonlinear fit of $B_{k}^{2}$ in Eq. (16a) is performed by
using the LMFIT package in Python. For more details on the MR method, see
Appendix B.

The convergence-limiting singularities for the series, evaluated at the wall
for $z=\\{0,\pi\\}$, are located on the positive real axis because
$\theta(z)=0$ [cf. FIG. 4 for the 1D HL model]; these singularities are
situated farthest from the origin in the complex-$\mathbf{q}$ plane. The
estimates for $R(z)$ in FIG. 8(a) are not arranged symmetrically about
$z=\\{\tfrac{\pi}{4},\tfrac{3\pi}{4}\\}$ [unlike in FIG. 3(a) for the 1D HL
model]. Here, the convergence-limiting singularities with the smallest radius
of convergence [$R_{m}\simeq 5.95$] occur at $z\simeq 0.88>\tfrac{\pi}{4}$ and
$z\simeq 2.26<\tfrac{3\pi}{4}$ in FIG. 8 (a). In the complex-$\mathbf{q}$
plane, the arrangement of the convergence-limiting singularities is in the
form of an eye [FIG. 9], whose center lies on the real axis but is displaced
to the right of the origin, for all three fields: $u^{1}\text{ (orange)},\
\omega^{1}$ (pink) and $u^{z}$ (blue), because of the asymmetry of $R(z)$
mentioned above. In panel (c) of FIG. 8, we plot the order $\nu$ of the
convergence-limiting singularities obtained for the fields $u^{1}$ (yellow),
$\omega^{1}$ (pink) and $u^{z}$ (blue). For $u^{1}$, $\nu$ takes non-integer
values which depend on $z$; for $z$ near $\tfrac{\pi}{2}$, it is negative but
otherwise remains positive. For $\omega^{1}$, $\nu\simeq-1/3$; and for
$u^{z}$, $\nu\in(0.5,1)$. [See Appendix B for estimates for $R$ obtained at
different distances $r$. ]

Figure 9: Polar plot in complex-$\mathbf{q}$ plane of $R(z)$ vs. $\theta(z)$
depicting the position of convergence-limiting complex singularities
$\mathbf{q}_{*}=Re^{i\theta}$ as well as their complex-conjugates, for $u^{1}$
(yellow), $\omega^{1}$ (pink) and $u^{z}$ (blue) at $r=r_{N}$ and
$z\in\mathbf{X}_{N,M}$, where $(N=128,M=256)$. The singularities are arranged
in the shape of an eye, which is squashed along the $\Im(\mathbf{q})$ axis;
the center of the eye is displaced to the right of the origin.

##  IV Conclusions

We have shown how to adapt the methods introduced in Ref. 28 to investigate
early-time resonances in two PDEs that have attracted considerable attention
recently [1, 39, 30, 2, 36, 37] in the context of possible finite-time
singularities in ideal hydrodynamical systems. The first is the 1D HL model
and the second is the 3DAE; the former is related approximately to the latter
[see, e.g., Refs. 1, 30, 2]. We have used initial data from earlier studies of
singularities for both these PDEs and employed our pseudospectral Fourier-
Chebyshev method [2], with quadruple-precision arithmetic, to compute the
time-Taylor series coefficients of the flow fields, up to a high order
[$N_{t}=100$ (1DHL) and $N_{t}=86$ (3DAE)]. The resulting approximations
display early-time resonances, which we have studied in detail. In particular,
we have demonstrated that the initial spatial location of these resonances is
different from that for the tygers, which we have obtained in Ref. 2. We have
analysed the time Taylor-series coefficients, by using the MR method, to
extract $R(z)$, $\theta(z)$, and $\nu(z)$, and thence the location and nature
of the convergence-limiting singularities. We have found that these
singularities are distributed around the origin, in the complex-$t^{2}$ plane,
along two curves that resemble an approximately elliptical eye.

The connection between complex-space and complex-time singularity landscapes
is a challenging problem that is now being explored for ideal hydrodynamical
PDEs such as the inviscid, unforced 1D Burgers equation [28]. Our work has
enlarged the scope of such studies by investigating such connections for the
1DHL and 3DAE systems. It is interesting to explore the relationship between
the singularity detection for complex-space singularities (e.g, via the
analyticity-strip method [23]) and time-Taylor-series methods that we use
here. It is important to understand the link (if any) between tygers, in
Galerkin-truncated pseudospectral studies, and the early-time resonances, in
truncated time-Taylor expansions of solutions of PDEs. The development of
numerical schemes for the mitigation of tygers or early-time resonances also
poses interesting questions that are being addressed in several recent studies
[40, 28, 41, 42].

###### Acknowledgements.

We thank SERB, CSIR, and NSM (India) for their support and SERC (IISc) for
computational resources. We thank J.K. Alageshan, U. Frisch, K.V. Kiran, T.
Matsumoto, C. Rampf, and S.S. Ray for useful discussions.

## References

  * [1] G. Luo and T. Y. Hou, “Potentially singular solutions of the 3D axisymmetric Euler equations,” Proceedings of the National Academy of Sciences, vol. 111, no. 36, pp. 12968–12973, 2014.
  * [2] S. S. V. Kolluru, P. Sharma, and R. Pandit, “Insights from a pseudospectral study of a potentially singular solution of the three-dimensional axisymmetric incompressible euler equation,” Physical Review E, vol. 105, no. 6, p. 065107, 2022.
  * [3] P. Constantin and C. Foiaş, Navier-stokes equations.  University of Chicago Press, 1988.
  * [4] C. R. Doering and J. D. Gibbon, Applied analysis of the Navier-Stokes equations.  Cambridge university press, 1995.
  * [5] P. Constantin, “Euler equations, navier-stokes equations and turbulence,” in Mathematical Foundation of Turbulent Viscous Flows: Lectures given at the CIME Summer School held in Martina Franca, Italy, SEptember 1-5, 2003, pp. 1–43, Springer, 2005.
  * [6] J. Gibbon, “The three-dimensional Euler equations: how much do we know?,” Physica D: Nonlinear Phenomena, vol. 237, no. 14-17, pp. 1894–1904, 2008\.
  * [7] L. Onsager, “Statistical hydrodynamics,” Nuovo Cimento, vol. 6 (Suppl.), p. 279–287, 1949.
  * [8] P. Constantin, E. Weinan, and E. S. Titi, “Onsager’s conjecture on the energy conservation for solutions of euler’s equation,” Communications in Mathematical Physics, vol. 165, no. 1, pp. 207–209, 1994.
  * [9] G. L. Eyink and K. R. Sreenivasan, “Onsager and the theory of hydrodynamic turbulence,” Reviews of modern physics, vol. 78, no. 1, p. 87, 2006.
  * [10] G. L. Eyink, “Dissipative anomalies in singular Euler flows,” Physica D: Nonlinear Phenomena, vol. 237, no. 14, pp. 1956 – 1968, 2008.  Euler Equations: 250 Years On.
  * [11] C. De Lellis and L. Székelyhidi Jr, “Dissipative euler flows and onsager’s conjecture,” Journal of the European Mathematical Society, vol. 16, no. 7, pp. 1467–1505, 2014.
  * [12] R. H. Morf, S. A. Orszag, and U. Frisch, “Spontaneous singularity in three-dimensional inviscid, incompressible flow,” Physical Review Letters, vol. 44, no. 9, p. 572, 1980.
  * [13] D. Gaunt and A. Guttmann, “Asymptotic analysis of coefficients,” Phase transitions and critical phenomena, vol. 3, pp. 181–243, 1974.
  * [14] C. Domb and M. F. Sykes, “On the susceptibility of a ferromagnetic above the curie point,” Proceedings of the Royal Society of London. Series A. Mathematical and Physical Sciences, vol. 240, no. 1221, pp. 214–228, 1957.
  * [15] C. M. Bender and S. A. Orszag, Advanced mathematical methods for scientists and engineers I: Asymptotic methods and perturbation theory.  Springer Science & Business Media, 2013.
  * [16] M. E. Brachet, D. I. Meiron, S. A. Orszag, B. Nickel, R. H. Morf, and U. Frisch, “Small-scale structure of the taylor–green vortex,” Journal of Fluid Mechanics, vol. 130, pp. 411–452, 1983.
  * [17] M. E. Brachet, D. Meiron, S. Orszag, B. Nickel, R. Morf, and U. Frisch, “The taylor-green vortex and fully developed turbulence,” Journal of Statistical Physics, vol. 34, no. 5, pp. 1049–1063, 1984.
  * [18] R. Pelz and Y. Gulak, “Evidence for a real-time singularity in hydrodynamics from time series analysis,” Physical review letters, vol. 79, no. 25, p. 4998, 1997.
  * [19] R. Morf, S. Orszag, D. Meiron, M. Meneguzzi, and U. Frisch, “Analytic structure of high reynolds number flows,” in Seventh International Conference on Numerical Methods in Fluid Dynamics, pp. 292–298, Springer, 1981\.
  * [20] U. Frisch and R. Morf, “Intermittency in nonlinear dynamics and singularities at complex times,” Physical review A, vol. 23, no. 5, p. 2673, 1981.
  * [21] D. I. Meiron, G. R. Baker, and S. A. Orszag, “Analytic structure of vortex sheet dynamics. part 1. kelvin–helmholtz instability,” Journal of Fluid Mechanics, vol. 114, pp. 283–298, 1982.
  * [22] K. Ohkitani and J. D. Gibbon, “Numerical study of singularity formation in a class of euler and navier–stokes flows,” Physics of Fluids, vol. 12, no. 12, pp. 3181–3194, 2000.
  * [23] C. Sulem, P.-L. Sulem, and H. Frisch, “Tracing complex singularities with spectral methods,” Journal of Computational Physics, vol. 50, no. 1, pp. 138–161, 1983.
  * [24] R. M. Kerr, “Evidence for a singularity of the three-dimensional, incompressible Euler equations,” Physics of Fluids A: Fluid Dynamics, vol. 5, no. 7, pp. 1725–1746, 1993.
  * [25] C. Cichowlas and M. Brachet, “Evolution of complex singularities in Kida–Pelz and Taylor–Green inviscid flows,” Fluid Dynamics Research, vol. 36, no. 4-6, p. 239, 2005.
  * [26] S. S. Ray, U. Frisch, S. Nazarenko, and T. Matsumoto, “Resonance phenomenon for the Galerkin-truncated Burgers and Euler equations,” Phys. Rev. E, vol. 84, p. 016301, 2011.
  * [27] S. D. Murugan and S. S. Ray, “Genesis of thermalization in the three-dimensional, incompressible, galerkin-truncated euler equation,” Physical Review Fluids, vol. 8, no. 8, p. 084605, 2023.
  * [28] C. Rampf, U. Frisch, and O. Hahn, “Eye of the tyger: Early-time resonances and singularities in the inviscid burgers equation,” Physical Review Fluids, vol. 7, no. 10, p. 104610, 2022.
  * [29] G. Mercer and A. Roberts, “A centre manifold description of contaminant dispersion in channels with varying flow properties,” SIAM Journal on Applied Mathematics, vol. 50, no. 6, pp. 1547–1565, 1990.
  * [30] K. Choi, T. Y. Hou, A. Kiselev, G. Luo, V. Sverak, and Y. Yao, “On the finite-time blowup of a one-dimensional model for the three-dimensional axisymmetric euler equations,” Communications on Pure and Applied Mathematics, vol. 70, no. 11, pp. 2218–2243, 2017.
  * [31] M. Frigo and S. G. Johnson, “The design and implementation of fftw3,” Proceedings of the IEEE, vol. 93, no. 2, pp. 216–231, 2005.
  * [32] A. J. Majda and A. L. Bertozzi, Vorticity and incompressible flow, vol. 27.  Cambridge University Press, 2002.
  * [33] J.-G. Liu and W.-C. Wang, “Characterization and regularity for axisymmetric solenoidal vector fields with application to Navier–Stokes equation,” SIAM Journal on Mathematical Analysis, vol. 41, no. 5, pp. 1825–1850, 2009\.
  * [34] R. Peyret, Spectral methods for incompressible viscous flow, vol. 148.  Springer Science & Business Media, 2013.
  * [35] M. Van Dyke, “Analysis and improvement of perturbation series,” The Quarterly Journal of Mechanics and Applied Mathematics, vol. 27, no. 4, pp. 423–450, 1974.
  * [36] D. Barkley, “A fluid mechanic’s analysis of the teacup singularity,” Proceedings of the royal society of london. Series A, Containing papers of a mathematical and physical character, vol. 476, 2020.
  * [37] T. Hertel, N. Besse, and U. Frisch, “The cauchy-lagrange method for 3d-axisymmetric wall-bounded and potentially singular incompressible euler flows,” Journal of Computational Physics, p. 110758, 2021.
  * [38] M. Newville, T. Stensitzki, D. B. Allen, M. Rawlik, A. Ingargiola, and A. Nelson, “Lmfit: Non-linear least-square minimization and curve-fitting for python,” Astrophysics Source Code Library, pp. ascl–1606, 2016.
  * [39] G. Luo and T. Y. Hou, “Toward the finite-time blowup of the 3D axisymmetric Euler equations: a numerical investigation,” Multiscale Modeling & Simulation, vol. 12, no. 4, pp. 1722–1776, 2014.
  * [40] S. D. Murugan, U. Frisch, S. Nazarenko, N. Besse, and S. S. Ray, “Suppressing thermalization and constructing weak solutions in truncated inviscid equations of hydrodynamics: Lessons from the Burgers equation,” Physical Review Research, vol. 2, no. 3, p. 033202, 2020.
  * [41] V. S. S. Kolluru, N. Besse, and R. Pandit, “Novel spectral methods for shock capturing and the removal of tygers in computational fluid dynamics,” arXiv preprint arXiv:2402.17688, 2024.
  * [42] N. Besse, “Analysis of the spectral relaxation method and the spectral purging method for nonlinear conservation laws.” Submitted.

##  Appendix A Code

The following MATHEMATICA 12 code computes analytical closed form expressions
for the time-Taylor coefficients; $u_{n}(z)$ is stored in u[n,z],
$\omega_{n}(z)$ in w[n,z] and $v_{n}(z)$ in v[n,z]. Here, Ntmax represents the
order of truncation $N_{t}$. The initial conditions are given in Lines 4-6.
Recursion relations given in Eqs. (3) are given in Lines 8-10. The
coefficients are then appended to listc which can ultimately be stored
symbolically and later retrieved.

Listing 1: Mathematica 12 code

⬇

1listc={}; Ntmax=70; f[n_] = Floor[(n - 1)/2] ;

2HilbertTransform[f_, x_, X_] := Module[ {fp = FourierParameters -> {1, -1},
k}, InverseFourierTransform[ -I (2 HeavisideTheta[k] - 1) *
FourierTransform[f, x, k, fp], k, X, fp]];

3

4u[0, z_] := (Sin[z])^2;

5w[0, z_] := 0;

6v[0, z_] := Integrate[HilbertTransform[w[0, x], x, z], z];

7

8w[n_, z_] := w[n, z] = If[EvenQ[n], 0, (- Sum[ If[i != (n - 1 - i), v[i, z]
D[w[n - 1 - i, z], z] + v[n - 1 - i, z] D[w[i, z], z], v[i, z] D[w[i, z], z]],
{i, 0, f[n], 1}] + D[u[n - 1, z], z])/n ];

9u[n_, z_] := u[n, z] = If[OddQ[n], 0,

10 (-Sum[ If[i != (n - 1 - i), v[i, z] D[u[n - 1 - i, z], z] + v[n - 1 - i, z]
D[u[i, z], z], v[i, z] D[u[i, z], z]], {i, 0, f[n], 1}])/n];

11v[n_, z_] := v[n, z] = If[EvenQ[n], 0, Integrate[HilbertTransform[w[n, x],
x, z], z]];

12

13Do[nn = n; AppendTo[listc, {n, u[n, z], w[n, z], v[n, z]}], {n, 0, Ntmax}];

##  Appendix B Mercer-Roberts method

For perturbation series where the coefficients’ signs follow a non-trivial
pattern, Mercer and Roberts [29] generalized the Domb-Sykes method to allow
for a pair of complex conjugate singularities.

Consider a model function that has complex conjugate singularities at
$\mathfrak{t}_{*}$ and $\overline{\mathfrak{t}_{*}}$

$\mathfrak{u}(\mathfrak{t})=\left(1-\frac{\mathfrak{t}}{\mathfrak{t_{*}}}\right)^{\nu}+\left(1-\frac{\mathfrak{t}}{\overline{\mathfrak{t_{*}}}}\right)^{\nu};\quad\text{where
}\mathfrak{t_{*}}\coloneqq Re^{i\theta}$ (19)

For $|\mathfrak{t}|<R$, the model function $\mathfrak{u}$ has the following
Taylor expansion around $t=0$:

$\mathfrak{u}(\mathfrak{t})=\sum_{n=0}^{\infty}2(-1)^{n}\binom{\nu}{n}R^{-n}\cos(n\theta)\mathfrak{t}^{n}$
(20)

Mercer and Roberts showed that $R,\theta$ and $\nu$ can be determined by
relating the model function in Eq. (20) to the original series
$u(t)=\sum_{n=0}^{\infty}u_{n}t^{n}$. For each $4$-tuple of $u_{n}$, we
construct:

$\displaystyle B^{2}_{k}$
$\displaystyle=\frac{u_{k+1}u_{k-1}-u^{2}_{k}}{u_{k}u_{k-2}-u^{2}_{k-1}}$
(21a) $\displaystyle cos\theta_{k}$
$\displaystyle=\frac{1}{2}\left(\frac{u_{k-1}B_{k}}{u_{k}}+\frac{u_{k+1}}{u_{k}B_{k}}\right)\,,$
(21b)

where $k=2,3,\ldots$. The leading-order behaviour is then obtained by
substituting Eq. (20) into Eq. (21):

$\displaystyle B_{k}=$
$\displaystyle\frac{1}{R}-\left(\frac{\nu+1}{R}\right)\frac{1}{k}+\left(\frac{\nu+1}{2R}\frac{\sin(2k-1)\theta}{\sin\theta}\right)\frac{1}{k^{2}}+\mathcal{O}(\frac{1}{k^{3}})$
(22a) $\displaystyle\cos\theta_{k}$
$\displaystyle=\cos\theta+\left[\cos\theta\ (\nu+1)\
\left(1-\frac{\cos(2k-1)\theta}{\cos\theta}\right)\right]\frac{1}{k^{2}}+\mathcal{O}(\frac{1}{k^{3}})\,.$
(22b)

For large $k$, we can extract $1/R$ from the intercept on the vertical axis,in
the plot of $B_{k}$ vs. $1/k$, and $\cos\theta$ from the plot of
$\cos\theta_{k}$ vs. $\frac{1}{k^{2}}$. In FIG. 10(a), we show the plots
versus $1/k$ of $B_{k}$, obtained from $u_{n}(z)$ in Eqs. (3a) for the 1D HL
model, at the spatial locations $z=\\{0,\tfrac{\pi}{8},\tfrac{\pi}{4}\\}$ in
yellow, red, and blue, respectively. For these values of $z$, the leading-
order linear behaviour dominates asymptotically. For $\theta$ close to
$\\{0,\pm\pi\\}$, the estimators in Eq. (22) become unreliable as high-orders
terms are not negligible.

In FIG. 10(b), we show the eye, including the MR estimates obtained by using
the linear forms (FIG. 10(a)), exactly at the singular $z$ points.

In FIG. 11 (a), we show the estimates for $R$, obtained as above, at different
distances $r$ from the wall. In particular, we compute the time-Taylor
approximations at the wall $r=1$, at $r=r_{N}$ (as before), and at an interior
point $r=0.96$. All these values of $r$ yield eyes with shifted centers of
$R(z)$ about $z=\\{\tfrac{\pi}{4},\tfrac{3\pi}{4}\\}$ persists. In FIG. 11
(b), we plot the associated eyes in the complex-$\mathbf{q}$ plane; the center
of the eye is shifted from the origin for all the choices of $r$ above. To the
best of our knowledge, the eye associated with the 3DAE is the first with a
shifted center; this shift remains to be understood.

(0,0) node[inner sep=0]; (-2.2,2.5) node (a); (0,0) node[inner sep=0];
(-2.5,2) node (b);

Figure 10: For the 1D HL model: (a) Plots vs. $1/k$ of $B_{k}(z)$ computed
from $u_{n}(z)$ using Eq. (16a), for $z=0\text{ (yellow)},\pi/8\text{
(red)},\pi/4\text{ (blue)}$. (b) Extended polar plot in complex-$\mathbf{q}$
plane of $R(z)$ vs. $\theta(z)$ depicting the position of convergence-limiting
complex singularities $\mathbf{q}_{*}=Re^{i\theta}$ as well as their complex-
conjugates, for $u$ (yellow), $\omega$ (pink) and $v$ (blue). Here, we note
the singularities on the real-axis obtained using MR method for this model.
The time-Taylor coefficients are estimated on a uniform grid of
$\mathbf{X}_{N=256}$ using quadruple-precision Fourier pseudospectral methods.

(0,0) node[inner sep=0]; (-2.1,3) node (a); (0,0) node[inner sep=0] ;
(-2.5,2.6) node (b);

Figure 11: For 3DAE: (a) Plots vs. $z$ of the estimates of $R$ and, (b) polar
plot, in the complex-$\mathbf{q}$ plane, of $R(z)$ vs. $\theta(z)$ depicting
the position of convergence-limiting complex singularities
$\mathbf{q}_{*}=Re^{i\theta}$ as well as their complex-conjugates, obtained
from $u^{1}_{n}(r_{N},z)$ for $z\in\mathbf{X}_{N,M}$, where $(N=128,M=256)$
and $n\leq N_{t}=86$, using the MR method; the estimates are shown for
different distances from the wall at $r\in\\{0.96,\,r_{N},\,1\\}$ (see
legend). Here, we see that the eye is shifted for all $r$ and the MR estimates
worsen far from the wall (where the singularity precipitates).


# Music Recommendation Based on Audio Fingerprint

D.S.U.Diego Saldaña-Ulloa  , , ,

###### Abstract

This work combined different audio features to obtain a more robust
fingerprint to be used in a music recommendation process. The combination of
these methods resulted in a high-dimensional vector. To reduce the number of
values, PCA was applied to the set of resulting fingerprints, selecting the
number of principal components that corresponded to an explained variance of
$95\%$. Finally, with these PCA-fingerprints, the similarity matrix of each
fingerprint with the entire data set was calculated. The process was applied
to 200 songs from a personal music library; the songs were tagged with the
artists’ corresponding genres. The recommendations (fingerprints of songs with
the closest similarity) were rated successful if the recommended songs’ genre
matched the target songs’ genre. With this procedure, it was possible to
obtain an accuracy of $89\%$ (successful recommendations out of total
recommendation requests).

audio fingerprint,

music recommendation,

spectral audio features,

MFCC,

chroma features,

tempo features,

###### keywords:

## Introduction

The use of digital systems that operate with music tracks has increased
considerably in recent years. Different tools exist for analyzing audio files,
and applications focus on correctly identifying a musical record to provide
additional services [1]. Audio identification has become an essential part of
different tasks on a day-to-day basis. This process is carried out employing
audio fingerprints. In this scope, a fingerprint is defined as a unique
representation of an audio track formed by different descriptive audio
features.

Audio fingerprinting techniques are sometimes confused with audio watermarking
techniques. Audio watermarking aims to embed information in an audio track,
for example, general information about the track for its correct
identification. In this way, audio fingerprints can have many immediate
applications, such as identification of copyrighted tracks, music
recommendations, or voice recognition, to name a few [2]. An audio
identification process through watermarking would involve directly extracting
the embedded message. However, audio watermarking could introduce noise since
the audio track information must be altered [1].

An audio identification system using fingerprints must have specific
characteristics that make the system robust and capable of performing correct
matches [1], among which we find: the number of accurate identifications, the
security of the system against noisy signals, the granularity related to
fragments audio files as short as possible, the independence of the audio
formats, the scalability over the number of fingerprints in the system and the
complexity related to the computational cost.

The audio identification system generally comprises an extraction block and a
block that stores the information for later identification. A system like this
can develop numerous applications; one of the most striking at present is the
one related to music recommendation systems [3]. The central idea relies on
obtaining the fundamental characteristics of a musical piece (fingerprint) and
recommending similar pieces.

In this way, the musical recommendation problem depends to a large extent on
the construction of the fingerprint, that is, the set of techniques used to
define the fingerprint of a song. The process generally involves extracting
time-frequency features based on modifications to the Fourier transform [4],
and there are different methods to achieve this objective [5] [6] [7].

In this work, we combine different methods that extract information from an
audio track based on its spectral, chromatic, and tempo features to extract
the most significant amount of structural information and provide robustness
to the process. Subsequently, we apply this procedure to a set of predefined
songs to fully characterize them; the intention is to develop a recommendation
system based on an individual’s musical tastes.

The structure of this work is as follows: Section 1 deals with previous
related works, section 2 introduces the theory for the extraction of spectral,
chromatic, and tempo features given an audio signal; Section 3 presents the
design proposal of the feature extraction system, and the recommendation
process, Section 4 focuses on the experimental results taking into account the
predefined data set; finally, the conclusions are detailed.

##  1 Related Works

The first version of audio fingerprinting was used in the 90s to detect
advertisements in broadcast streams [8]. This first version used the signal’s
energy information to match the signal waves’ envelope. The objective focused
on reusing transmissions and the automatic placement of other advertisements.

Currently, the processes to obtain the fingerprint of an audio signal focus on
the extraction of structural features based on spectral analysis. For example,
it is common to use variations of the Fourier transform, such as the Short
Time Fourier Transform STFT [4][9][10], which aims to extract the
transformation coefficients associated with the frequencies of the signal.
Unlike the standard Fourier transform, STFT operates on time windows,
resulting in the final coefficients being associated with interval frequencies
of the signal and not the complete signal. Another spectral technique used for
fingerprinting is the Mel Frequency Cepstral Coefficients (MFCC) [11][12][13].
The MFCC is obtained by mapping the STFT coefficients to the Mel scale,
defined as a perceptual scale based on pitches and frequencies perceptible by
the human ear.

There are also features related to the pitch classes of Western music called
chroma features. These variables capture the melodic characteristics of a
piece of music. Different works use this type of variable [7][14], and, like
the MFCC, they work by calculating the STFT coefficients, combining them with
binning techniques.

On the other hand, it is common to use features related to the tempo or beats
of the audio signal. This feature type has been used as a complementary
process in audio analysis to obtain fingerprints related to musical aspects
such as tempo or rhythm [15][16]. Generally, this approach works using the
spectral parameters obtained by STFT to identify the onset of the signal
corresponding to the beats’ fundamental characteristics.

The types of features that we have mentioned correspond to spectral,
chromatic, tempo, and Mel Scale-based. Some works have focused on combining
some [10][11] [13] [15] to obtain a more robust fingerprint capable of
capturing more information about the audio signal and achieving correct
identification. Despite the above, we did not find a work that is precisely
dedicated to combining all the previous methods to obtain a fingerprint.

Regarding musical recommendation systems, some works focus exclusively on the
similarities of songs based on spectral features [17][18]. Other works are in
charge of combining this type of feature with more contextual information
about the song, for example, the name, the duration, or the genre [19][20].
The purpose is to create vector arrays that fully identify a piece of music,
making it as distinctive as possible. With these arrays, songs can be compared
to recommend similar items. Currently, music recommendation systems are used
by large streaming music companies and use a combination of methods for
constructing the features as well as different algorithms for comparing each
element [3].

##  2 Audio Features

Next, we will present the theory related to different types of audio features
that will be used to develop a robust fingerprint. These features represent
the audio signal’s spectral, chromatic, Mel scale-based, and tempo
characteristics.

###  2.1 Short Time Fourier Transform

The Fourier transform is useful for analyzing different types of signals. This
transformation aims to decompose the signal into frequency components
$w\in\mathbb{R}$ associated with coefficients $d_{w}\in\mathbb{R}$. The
coefficients indicate the degree to which a sinusoidal component with
frequency $w$ adjusts the signal. In this way, the Fourier transform
represents the signal in the frequency domain compared to its original form,
which represents the time domain. For example, the Fourier transform of a
musical extract indicates which notes (frequencies) have been played but not
when they are played [5].

In order to make use of the temporal information associated with a signal, the
so-called short-time Fourier transform (STFT) is used. This transform type
operates on small portions of the signal (time windows) compared to the
original transform. Let $x\in\mathbb{R}$ be a discrete signal,
$r:[0:N-1]\rightarrow\mathbb{R}$ a time window function of length
$N\in\mathbb{N}$, and a parameter $H\in\mathbb{N}$ called hop size that
indicates the intervals in which the function $r$ moves along the signal; The
STFT of the signal is given by

$\xi(m,k)=\sum^{N-1}_{n=0}x(n+mH)r(n)exp(-2\pi ikn/N),$ (1)

where $m\in\mathbb{Z}$ and $k\in[0:N/2]$. $\xi(m,k)$ is the $k^{\text{th}}$
fourier coefficient for the $m^{\text{th}}$ time frame.

Since the results of STFT are complex numbers, we can transform those values
into a real two-dimensional representation by squaring their magnitude

$\Theta(m,k)=|\xi(m,k)|^{2}.$ (2)

This is called the spectrogram and is helpful for visualizing the
transformation. The horizontal axis represents the time, and the vertical axis
represents frequency.

###  2.2 Mel Frequency Cepstral Coefficients

This type of coefficient is based on the structural analysis of timbre.
Initially, the coefficients were used for voice recognition tasks, and
obtaining them involves mapping the frequencies to the Mel scale [5]. The Mel
scale is a perceptual scale that tries to mimic the hearing tones of the human
ear. It divides the frequencies into equidistant intervals on the final scale.
The following function gives the frequency transformation to the Mel scale

$F_{\text{Mel}}(x)=1127\ \text{ln}(1+\frac{x}{700}).$ (3)

To compute the MFCC, a series of steps are used [13] [21] [22] that can be
summarized as 1) calculate the STFT coefficients (as in the previous section);
2) make a mapping between the calculated coefficients and the Mel scale. A
series of triangular filter banks along the Mel scale usually do this
transformation. The objective is a binarization along the scale, depending on
the number of filter banks used. The Mel coefficients have the following form

$S(n)=\sum^{N-1}_{k=0}|\xi(k)|^{2}J_{m}(k),$ (4)

where $0\leq m\leq M-1$, $m$ is the number of triangular filters, and $J$ is
the weight applied to the $k^{\text{th}}$ energy spectrum bin that contributes
to the $m^{\text{th}}$ output band. This result is the energy spectrum of the
Mel frequencies. 3) Finally, the Discrete Cosine Transform (DCT) is calculated
on the previous coefficients to decorrelate them and obtain only a determined
number of coefficients

$\Lambda_{m}(n)=\sum^{M-1}_{m=0}S(n)\ \text{cos}(\frac{\pi n(m-0.5)}{M}),$ (5)

where $n\in\mathbb{N}$ is the number of MFCC desired, in this work 13
coefficients were used because is a typical value according to [13] [21].

###  2.3 Chroma Coefficients

This type of feature considers properties related to the harmony and melody of
a signal. In order to obtain these features, tempered scales are considered,
mainly the one related to Western music that divides the scale into 12
different tones $\\{C,C\\#,D,D\\#,E,F,F\\#,G,G\\#,A,A\\#,B\\}$. To compute the
chroma coefficients, we must calculate the STFT coefficients and obtain a log
frequency spectrogram

$Y(n,p)=\sum_{k\in P(p)}|\xi(n,k)|^{2},$ (6)

where $p\in[0:127]$ is the MIDI note number and represents the note’s pitch,
and $P(p)$ is a function that maps $k$ to the MIDI numbers. In short, the MIDI
note number encodes the musical pitches along different octaves. Once this
spectrogram is obtained, the chroma coefficients or chromagram are calculated
by adding all the pitch coefficients that belong to the same chroma

$\Psi(n,c)=\sum_{p\in[0:127]:p\text{mod}12=c}Y(n,p),$ (7)

where $c\in[0:11]$.

###  2.4 Tempo and Beat Features

The beat corresponds to the pulse in an audio signal and is an additional
characteristic that can be appreciated when listening to a song. The tempo is
the rate at which the beats occur, corresponding to a fundamental aspect of
music. Like chromatic coefficients, these features are related to the melodic
aspects of an audio signal.

In order to identify the beats of a piece of music, one of the approaches used
is beat tracking by dynamic programming. The goal is to generate a sequence of
beat times that correspond to perceived onsets in the audio signal and that
have a regular rhythm [23]. In order to do this, an objective function is used
that combines a part related to the locations of the beats and another
function that penalizes the deviations of irregular intervals

$\Xi({t_{i}})=\sum^{N}_{i=1}O(t_{i})+\alpha\sum^{N}_{i=2}F(t_{i}-t_{i-1},\tau_{p})$
(8)

where ${t_{i}}$ is the sequence of $N$ beats, $O(t)$ is an onset strength
envelope of the audio signal, $\alpha$ is a parameter to measure the
importance of the penalty function $F$, and $\tau_{p}$ is an ideal beat
spacing.

To calculate the best value of $\Xi$, we take the recursive relation

$\Xi^{*}(t)=O(t)+\max_{\tau=0,...,t}\\{\alpha
F(t-\tau,\tau_{p})+\Xi^{*}(\tau)\\}$ (9)

In other words, the best value $\Xi$ for time t is the sum of the local onset
strength plus the best value of $\Xi$ at time $\tau$. Since this function
always operates with the preceding $\tau$ value (to get the instants of time
$t$ that mark the beats), the iteration of this procedure is from the end
towards the start of the signal.

Table 1: Set of Songs Genre | Number of songs
---|---
Metal Variants | 28
Rock, Alternative, Indie, New Wave | 27
Ballad, Pop | 15
House, Electronic, Dance, Trance | 14
Synthpop, Electropop, Technopop | 12
Pop, Ballad, Rock | 12
Pop, Rock, Folk, Indie | 11
Jazz, Blues, Ballad | 10
Pop, Punk, Rock | 10
Electropop, Electronic, Hip hop | 9
Pop, Dance, Electropop | 9
Nortech, Electronic | 8
Rock, Pop, Urbano | 7
Classical, Soundtrack | 6
Pop, Indie, Rock | 6
Rock, Metal | 5
Pop, Rap, Dance | 4
Rock, Blues, Jazz | 4
Hip Hop, Rap, Rock | 3

##  3 Audio Features and Music Recomendation

In practical terms, the methods described in the previous section produce
matrix coefficients that describe essential aspects of a musical composition.
For example, the coefficients resulting from STFT, also called a spectrogram,
give an overview of an audio signal since we obtain information about the
frequencies that compose it. MFCCs provide information about timbre related to
speech. Through the chroma coefficients, we can obtain properties of the
harmony and melody related to the Western musical scale (chromagram). Finally,
we can identify the exact moments corresponding to a song’s beats (pulses)
related to its tempo.

Although most of the mentioned methods start initially from the STFT results,
each differs in the sequence of operations to obtain more information
regarding an audio signal or song. This work focuses on obtaining the
previously mentioned coefficients of a set of songs and combining them to
obtain a more robust fingerprint that considers particular and essential
aspects of a song.

To carry out this objective, we define the vector $\Omega$ as the vector that
includes the coefficients of each of the methods described above, such that,

$\Omega=\bar{\Theta}\oplus\bar{\Lambda}\oplus\bar{\Psi}\oplus\bar{\Xi},$ (10)

where $\bar{\Theta}$, $\bar{\Lambda}$, $\bar{\Psi}$ denotes the matrix of the
mean values of the rows, $\bar{\Xi}$ is a vector of size of $\Omega$ with
tempo values, and $\oplus$ indicates the concatenation of these matrices. With
the mean values of the rows of each matrix, we reduce the dimension
encapsulating the variations of each frequency along all the time.

Defining this vector or fingerprint focuses on condensing as much information
as possible from an audio signal. However, due to the concatenation of the
matrices, the resulting vector can have high dimensionality. Since the
objective is to obtain musical recommendations, reducing the dimensionality of
the vector is essential to reduce the complexity of the process.

One widely used technique to reduce a data set’s dimensionality is principal
component analysis (PCA) [24]. This method creates linear combinations of the
original variables called principal components. In this way, the principal
components retain part of the original information (variance) of the original
data. This reduction method aims to choose the number of principal components
that retain a considerable percentage of the total variance of the original
data [24].

Let $\Omega=\\{\Omega^{\\{1\\}},\Omega^{\\{2\\}},...,\Omega^{\\{n\\}}\\}$ be
the matrix formed by the union of different musical fingerprints and
$\Omega^{\prime}$ the resulting matrix of the principal components of
$\Omega$, where
$\Omega^{\prime}$={$\Omega^{\prime}_{1},\Omega^{\prime}_{2},...,\Omega^{\prime}_{n}$}
correspond to the principal components of each of the musical fingerprints; we
can compute the distance between fingerprints $i,j$ by some distance metric,
like euclidean distance

$S_{e}=||\Omega_{i}^{\prime}-\Omega_{j}^{\prime}||$ (11)

For this case, the $S_{e}$ matrix corresponds to the distance matrix between
each fingerprint. This matrix expresses the similarities between music tracks
and is widely used in music recommendation systems [3]. A music recommendation
system aims to find the top K of the closest fingerprints (smallest distance)
corresponding to the musical recommendations given a set of target songs. The
target songs correspond to songs presumably added to a list of favorite songs.
This type of recommendation is called content-based and is used to give the
most similar items the user has liked [25].

Figure 1: The number of components needed to explain the $95\%$ of the
variance

##  4 Experimental Results

For the experimental part, 200 songs from a personal music library composed of
different genres and artists were used. Each song was tagged with the genres
corresponding to the artist, obtained from Wikipedia. Table 1 shows the number
of songs corresponding to a set of genres. Additionally, the bitrate of each
audio file was kept constant at 128 kbps (CD quality). In order to compute the
coefficient described in the previous sections, 60 seconds of audio were used
starting at second 60, i.e., from the second 60 to the 120. Each audio file
was processed in a mono signal. For the methods that use STFT, the sampling
rate was set to 22050 Hz (samples per second), a window size of 2048, and a
cosine window (hann).

The first part of the process focused on fingerprinting the sound files, as
described in the previous section. In our case, each resulting vector
(fingerprint) comprised 1062 coefficients. PCA was applied to each vector to
reduce its dimensionality, and the number of components selected corresponded
to a retained variance of $95\%$ Figure 1. This reduced the size of the vector
from 1062 to 29, that is, a reduction of $97\%$ on the number of dimensions.

Figure 2 shows a graph of the two principal components of each vector song. It
can be seen that the genres of the songs tend to cluster: in the center, we
find the genres of rock, pop, alternative, and indie; in the lower right part,
we find the metal variants; in the upper right part, the genres related to
electronic music and in the central left part classical, ballads, jazz and
indie. This shows us that chromatic, tempo, and spectral coefficients are
useful in describing a music track.

The similarity between each pair of songs was computed using the 29 principal
components for each song and the Euclidean distance, as described in the
previous section. Then, the top-K of the most similar elements for each track
was obtained using the similarity matrix. For this case, the top 3 most
similar songs were retrieved, and a successful evaluation was set for any
recommendation that corresponded to the genre of the target song. This type of
evaluation system is commonly used within recommendation systems [3]. With
this consideration, it was possible to obtain an accuracy of $89\%$
(successful recommendations out of total recommendation requests), reaffirming
the ability of the combination of chroma, tempo, MFCC, and STFT coefficients
to obtain a robust fingerprint and to be able to categorize a song.

Figure 2: The two principal components and the distribution of songs in this
two dimensional space

##  5 Conclusions

This work combined different audio features to obtain a robust musical
fingerprint of a set of songs from a personal music library of 200 songs. Each
song was tagged according to the musical genres of the artist. This
fingerprint was later used in the process of recommending similar songs. The
features used to define the fingerprint correspond to spectral, chromatic,
tempo, and Mel scale-based features. Since each of the methods described above
generally results in a matrix of values, the average values of the rows were
obtained to capture the variation of each frequency throughout the signal.
This process enabled obtaining a simplified vector (fingerprint) that combines
the audio information described above. Because the size of the resulting
vector contained 1062 elements, PCA was applied to the set of fingerprints of
each song, and the 29 main components corresponding to a retained variance of
$95\%$ were selected. With this information (PCA-fingerprints), the similarity
matrix was calculated, defined as the matrix of values that reflect the
similarity of each song with the rest of the data set. For each song, the top
3 most similar songs were obtained, and the recommendation was qualified as
successful if the genre of any of the recommendations matched the genre of the
target song. With this process, it was possible to obtain $89\%$ accuracy in
the recommendation process, demonstrating that the union of different audio
features allows for a more robust fingerprint that captures relevant song
information.

## References

  * [1] Cano, P., Batlle, E., Kalker, T., Haitsma, J.: A review of audiofingerprinting. VLSI Signal Processing 41, 271–284 (2005)
  * [2] Lerch, A.: Audio fingerprinting. In: An Introduction to Audio Content Analysis Applications in Signal Processing and Music Informatics, 1st edn., p. 163. John Wiley & Sons, United States of America (2012)
  * [3] Schedl, M., Zamani, H., Chen, C., Deldjoo, Y., Elahi, M.: Current challenges and visions in music recommender systems research. Int J Multimed Info Retr 7, 95–116 (2018)
  * [4] Bellittini, C., Mazzini, G.: A framework for robust audio fingerprinting. Journal of Communications 5, 409–424 (2010)
  * [5] Muller, A.: Content-based audio retrieval .g. In: Fundamentals of Music Processing, 1st edn., p. 163. Springer, Switzerland (355)
  * [6] Haitsma, J., Kalker, T.: A highly robust audio fingerprinting system with an efficient search strategy. Journal of New Music Research 32, 211–221 (2003)
  * [7] Kim, S., Unal, E., Narayanan, S.: Music fingerprint extraction for classical music cover song identification. In: IEEE International Conference on Multimedia and Expo, pp. 1261–1264 (2008)
  * [8] Lourens, J.: Detection and logging advertisements using its sound. Journal of New Music Research 36, 231–236 (1990)
  * [9] Ramona, M., Peeters, G.: Audio identification based on spectral modeling of bark-bands energy and synchronization through onset detection. In: ICASSP, pp. 477–480 (2011)
  * [10] Ramalingam, A., Krishnan, S.: Gaussian mixture modeling of short-time fourier transform features for audio fingerprinting. In: IEEE Transactions on Information Forensics and Security, pp. 457–463 (2006)
  * [11] Thiruvengatanadhan, R., Dhanalakshmi, P.: Indexing and retrieval of music using gaussian mixture model techniques. In: International Journal of Computer Applications, pp. 35–41 (2016)
  * [12] Abu-El-Quran, A., Goubran, R., Chan, A.: Music fingerprint extraction for classical music cover song identification. In: IEEE International Conference on Multimedia and Expo, pp. 1261–1264 (2008)
  * [13] Khasana, Z.: Feature extraction techniques in music information retrieval. PhD thesis, TALLINN UNIVERSITY OF TECHNOLOGY (2020)
  * [14] Miotto, R., Orio, M.: A music identification system based on chroma indexing and statistical modeling. In: International Society for Music Information Retrieval Conference (2008)
  * [15] Ellis, D., Poliner, G.: Identifying ‘cover songs’ with chroma features and dynamic programming beat tracking. In: IEEE International Conference on Acoustics, Speech and Signal Processing - ICASSP (2007)
  * [16] Kurth, F., Gerhmann, T.: The cyclic beat spectrum: Tempo-related audio features for time-scale invariant audio identification. In: International Conference on Music Information Retrieval (2006)
  * [17] Badu, D., Suprita, P.: Music emotion recommender system using spectral features-a malayalam cine music deployment. In: International Conference on Electronics and Sustainable Communication Systems (ICESC), pp. 1306–1310 (2021)
  * [18] Han, H., Luo, X., Yang, T., Shi, Y.: Music recommendation based on feature similarity. In: IEEE International Conference of Safety Produce Informatization, pp. 650–654 (2018)
  * [19] Han, B., Rho, S., Jun, S., Hwang, E.: Music emotion classification and context-based music recommendation. Multimedia Tools and Applications 47, 433–460 (2010)
  * [20] Deldjoo, Y., Schedl, M., Knees, P.: Content-driven music recommendation: Evolution, state of the art, and challenges. Arxiv (2021)
  * [21] Logan, B.: Mel frequency cepstral coefficients for music modeling. In: International Society for Music Information Retrieval Conference (2000)
  * [22] Wibawa, I., Darmawan, D.: Implementation of audio recognition using mel frequency cepstrum coefficient and dynamic time warping in wirama praharsini. In: Tenth International Conference and Workshop on High Dimensional Data Analysis (2021)
  * [23] Ellis, D.: Beat tracking by dynamic programming. Journal of New Music Research 36(1), 51–60 (2007)
  * [24] Ringnér, M.: What is principal component analysis? Nature Biotechnology 26(3), 303–304 (2008)
  * [25] Lops, M..S.G. P & De Geminis: Content-based recommender systems: State of the art and trends. In: Recommender Systems Handbook. Springer, . (2011)