KCL-PH-TH/2024-14

# Multimessenger signals from compact axion star mergers

Liina Chung-Jukkoa<EMAIL_ADDRESS>Eugene A. Lima David J. E. Marsha
aTheoretical Particle Physics and Cosmology Group, Physics Department, King’s
College London, Strand, London WC2R 2LS, United Kingdom

###### Abstract

Axion dark matter can form stable, self-gravitating, and coherent
configurations known as axion stars, which are rendered unstable above a
critical mass by the Chern-Simons coupling to electromagnetism. We study,
using numerical relativity, the merger and subsequent decay of compact axion
stars. We show that two sub-critical stars can merge, and form a more massive,
excited and critical star, which survives for a finite period before rapidly
decaying via electromagnetic radiation. We find a rich multimessenger signal,
composed of gravitational waves, electromagnetic radiation, and axion
radiation. The gravitational wave signal is broken into two parts: a weak and
broad signal from the merger, followed by a much stronger signal of almost
fixed frequency from the decay. The electromagnetic radiation follows only the
gravitational waves from the decay, while the axion signal is continuous
throughout the process. We briefly discuss the detectability of such a signal.

##  I Introduction

Axion stars are self-gravitating compact objects formed of the QCD axion [1,
2, 3, 4, 5, 6, 7, 8, 9, 10, 11] or axion-like particles [12, 13, 14, 15, 16,
17]. If the axion forms part or all of dark matter (DM), these stars would be
expected to be abundant inside DM halos (for initial formation mechanisms, see
Refs. [18, 19, 20, 21, 22]). As a time-periodic solution of the Einstein-
Klein-Gordon equations, axion stars belong to the solitonic class of
oscillatons [23, 24]. In the regime of strong gravity, these stars are
compact, relativistic objects, where the mass of the compact star goes as
$M_{s}\sim m_{\mathrm{Pl}}^{2}/m$. Throughout this work, units of the axion
mass, $m$, are absorbed by choice of length and time units, and we furthermore
use geometrized units with $G=c=1$ 111To convert from our simulation units, as
displayed on the figures, to physical units, the unit of $m$ has to be
multiplied by $m_{\mathrm{Pl}}^{-2}$ before doing the usual conversion from
geometrized units with factors of $G$ and $c$..

Axions are coupled to electromagnetism (EM) through a Chern-Simons term with
coupling constant $g_{a\gamma}$. In Ref. [25] we studied the instability of
compact stars induced by $g_{a\gamma}$. Compact axion stars decay into a less
massive and less compact remnant for masses $M_{s}$ larger than a critical
value, which is set by the coupling constant,
$g_{a\gamma}^{\mathrm{crit}}\propto M_{s}^{-1.35}$. We demonstrated how this
decay process is due to parametric resonance between the axion star and the EM
field, which causes the star to lose energy by emitting EM radiation at its
characteristic frequency. More specifically, when the characteristic frequency
of the oscillating axion star lies within the parametric resonance band of the
underlying scalar field, $\phi$, the oscillations of the star drive photons
out of the electromagnetic field, causing a strong burst of EM radiation. We
discuss the shutting down of the parametric resonance process in more detail
in Appendix A. The instability of compact axion stars in the strong gravity
regime studied in Ref. [25] is in agreement with previous studies of the
instability of non-compact axion stars in the weak field regime [26, 27], and
other related work (e.g. Refs. [8, 9, 28, 29, 30]).

Hence, the existence of compact axion stars is restricted by their coupling to
electromagnetism, which naturally raises the possibility of mergers as a
formation mechanism for compact axion stars, as any “pristine” super-critical
stars created in the early universe e.g. following Ref. [20], would rapidly
decay. Therefore in the following we investigate the EM interaction of an
axion star created by a head-on collision of two less compact but
gravitationally stable stars. The original stars belong to the sub-critical
compactness 222In this paper, we use the same definition for the compactness
parameter $C$ as in [31]. For an axion star mass $M_{s}$, $C=GM_{s}/R$, where
$R$ is the effective radius of the star (includes $95\%$ of its mass). region
found in Ref. [31] where the merger leads to the formation of another excited
but stable star (for more compact stars the merger leads to the prompt
formation of a black hole). Compact axion stars are gravitational objects, and
thus one expects their merger [32] and subsequent decay to produce
gravitational waves, in addition to the EM and axion (scalar) radiation seen
already in Ref. [25]. We extract, for the first time, the full multimessenger
signal in GW, EM and axions caused by axion star collisions and decays. An
example signal extracted from the clean decay of an isolated compact axion
star [25] is shown in Fig. 1 (the more complex signal from mergers will be
shown later in the paper).

Figure 1: Multimessenger signal of axion star decay by electromagnetic
instability, for an isolated super-critical star: gravitational power
(luminosity) (top), total flux of electromagnetic energy (middle), and
corresponding total flux of the scalar energy (bottom). Note that power is
dimensionless in geometrized units. Signals are extracted over a sphere of
radius $r=80m^{-1}$, and the power in electromagnetic and scalar radiation is
obtained by integrating the corresponding flux. We note that the
electromagnetic and GW signal reach the observer at the same time, as both
travel at the speed of light, whereas the scalar emission is delayed as it is
massive. The different frequency wave modes in the scalar emission are
responsible for the long tail in the total scalar flux, as lower frequency
modes travel more slowly away from the decaying star. The EM and scalar fluxes
are calculated according to the method from Eq. 11 in Ref. [33], and the
gravitational power with the usual gravitational luminosity formula (see
Appendix 5).

Our calculations are performed in the strong gravity regime, using the 3+1
numerical relativity code GRChombo [34, 35, 36], hence capturing all
backreactions and other gravitational effects. We show that two stars
originally in the stable region of the plane given by the two dimensionful
parameters $(M_{s},g_{a\gamma})$ (i.e. below the critical line
$g_{a\gamma}^{\mathrm{crit}}(M_{s})$) can merge and form a highly excited,
more compact star that correspondingly has a lower critical coupling, and is
thus unstable electromagnetically. This merged star then decays through
parametric resonance just as the single stars studied in [25], emitting a
strong EM burst. 333As in Ref. [25], an EM seed is required to trigger the
resonance: in our simulations such a seed is given by a weak propagating EM
wave, modelling ambient radiation in the galaxy. The GW signal of the decay in
Fig. 1 does not exhibit any ringdown hierarchy as no black hole is formed.
Instead, its frequency matches the axion mass, implying that it is primarily
sourced by the quadrupole moment during the decay process. The GW signal from
the merger is largely unaffected by the EM coupling, and agrees with the
results of Ref. [31] (which included axion self-interactions but no EM
coupling). Thus, the multimessenger signal of compact axion star mergers has a
distinct portrait: an initial weak GW signal from the merger, followed by a
louder GW burst and accompanying EM signal, and finally a delayed and
prolonged axion signal.

This paper is organised as follows: Section II introduces our oscillaton model
and the initial set up of our simulations. Section III outlines the story of
two axion star colliding head-on, merging, and the merged star decaying
through the electromagnetic instability demonstrated in [25]. Section IV shows
the gravitational wave signal from the electromagnetic decay of axion stars,
and the combined signal from the merger and decay process, followed by the
full multimessenger signal in Section V. Finally, we discuss the observational
consequences of our findings in Section VI.

##  II Theory

Our work in this paper uses the same oscillaton model as in our previous study
[25]. The total action is

$\displaystyle S=\int d^{4}x\sqrt{-g}$
$\displaystyle\left[\frac{m_{\mathrm{Pl}}^{2}}{16\pi}R-\frac{1}{2}\partial_{\mu}\phi\partial^{\mu}\phi-\frac{1}{2}m^{2}\phi^{2}\right.$
(1)
$\displaystyle\qquad\left.-\frac{1}{4}F_{\mu\nu}F^{\mu\nu}-\frac{g_{a\gamma}}{4}\phi
F_{\mu\nu}\tilde{F}^{\mu\nu}\right],$

where $\phi$ denotes the axion field, and $R$ the Ricci scalar. The
electromagnetic field strength tensor and its dual are defined as

$F_{\mu\nu}=\partial_{\mu}A_{\nu}-\partial_{\nu}A_{\mu},\quad\widetilde{F}^{\mu\nu}=\frac{1}{2\sqrt{-g}}\varepsilon^{\mu\nu\rho\sigma}F_{\rho\sigma},$
(2)

with the totally antisymmetric Levi-Civita symbol
$\varepsilon^{\mu\nu\rho\sigma}$, defined through $\varepsilon^{0123}=+1$.

The axion pseudoscalar field is coupled to the field strength tensor through
the Chern-Simons term (last term in Eq. 1), where the coupling constant
$g_{a\gamma}$ sets the strength of the interaction. Moreover, the stress-
energy tensor can be calculated from the action in Eq. (1) to solve Einstein’s
equations $G_{\mu\nu}=8\pi m_{\mathrm{Pl}}^{-2}T_{\mu\nu}$ (see e.g. Ref.
[37]), where we note that the Chern-Simons topological term does not enter the
stress-energy tensor.

The parametric resonance process, as described in Section I, is driven by the
sourced Maxwell’s equation for the evolution of the field strength tensor, Eq.
(4),

$\displaystyle\nabla^{\mu}\nabla_{\mu}\phi-m^{2}\phi$
$\displaystyle=\frac{g_{a\gamma}}{4}F_{\mu\nu}\tilde{F}^{\mu\nu},$ (3)
$\displaystyle\nabla_{\mu}F^{\mu\nu}$ $\displaystyle=-g_{a\gamma}J^{\nu},$ (4)

where the current $J^{\nu}$ is given by
$J^{\nu}=\partial_{\mu}\phi\tilde{F}^{\mu\nu}$. In general, resonance can
happen provided that the photon frequency lies within the parametric resonance
band. As we discussed in [25], the axion mass implies that the oscillations of
the star go as $\omega\sim m\sim R_{s}^{-1}$, hence setting a bandwidth for
the photon frequencies that can drive the parametric resonance process.

Our system of coupled axion star, electromagnetism and the metric is solved
using GRChombo [35, 36, 34]. We follow the same procedure to find initial
conditions for the axion star and the electromagnetic field as in our previous
work [25]. The axion star initial set up is built following Refs. [23, 31, 38,
39, 40, 23, 41, 42], while the implementation of the electromagnetic field
follows the methodology introduced in Refs. [43, 44, 45, 46, 47].
Specifically, for the superposition of two axion stars as initial conditions,
we use the profiles introduced in Ref. [31, 48, 49] where the spatial metric
is renormalized to reduced gauge artifacts. The EM initial conditions were set
up as outlined in Ref. [25], with the initial EM amplitude kept the same for
all stars ($C_{x}=0.001m_{\mathrm{Pl}}$), but a frequency of $k^{(x)}\equiv
2\pi/\lambda\sim 0.10m$ for the single, more massive star and $k^{(x)}\sim
0.12m$ for the colliding stars.

As we solve the fully non-linear evolution equations numerically, our
simulations include all backreactions. We implement periodic boundary
conditions in the propagation direction of our electromagnetic plane wave seed
to keep the radiation in the simulation box, and Sommerfeld boundary
conditions in the remaining two directions. The Sommerfeld boundary conditions
absorb outgoing radiation at the boundary of our simulation box which improves
the accuracy of gravitational wave signal extraction from the simulation. The
evolution of the scalar field is not affected by the choice of boundary
conditions, as the main effect of implementing Sommerfeld boundary conditions
is to absorb outgoing gravitational and electromagnetic radiation. The
convergence of our simulation, along with the convergence of the extracted
gravitational waves, is discussed in Appendix B.

##  III Mergers and electromagnetic instability

Our story of two axion stars colliding head-on and merging, is captured in
Fig. 2, which demonstrates the evolution of the energy density in the scalar
and EM fields over time in a slice through our simulation box. We start with
two gravitationally and electromagnetically stable stars that collide and form
a highly excited merged star that is stable gravitationally, but which after
some time decays electromagnetically through parametric resonance. From now
on, we will call the star resulting from the collision and merger of the two
original axion stars “merged star”, as opposed to the “single star” which is
directly set up by the initial conditions.

Figure 2: Evolution of the energy density in the scalar (blue panels) and EM
field (black panels) as a slice through the center of our simulation box. We
start with two stable stars of mass $M_{s}=0.36m_{\mathrm{Pl}}^{2}/m$ at
$t=0m^{-1}$, which merge to form a highly excited, but stable star (around
$t\sim 180m^{-1}$). Our coupling value is
$g_{a\gamma}=20m_{\mathrm{Pl}}^{-1}$, implying that the final merger star is
super-critical according to Ref. [25], resulting in the decay of the merged
star through parametric resonance (see the lower set of panels). This happens
because, due to the higher mass, the critical coupling of the merged star is
lower than the critical coupling of the original stars, allowing the
parametric resonance process after the merger. The initial seed
electromagnetic field (not visible on this scale) is polarised in the
$x$-direction and propagating from the right to the left. A movie of our
simulations can be found through this link.

More specifically, we start with compact axion stars in the sub-critical
region demonstrated in [31], where the head-on merger of these stars leads to
a highly excited, more compact star. Our stars are initially separated by a
distance of $40m^{-1}$, and have a mass of $M_{s}=0.36m_{\mathrm{Pl}}^{2}/m$.
Using the scaling relation between the critical coupling of a compact star and
its mass $g_{a\gamma}^{\mathrm{crit}}\propto M_{s}^{-1.35}$, which we
demonstrated in Ref. [25], this gives our initial stars a critical coupling of
$g_{a\gamma}^{\mathrm{crit}}\sim 24.6m_{\mathrm{Pl}}^{-1}$.

Hence, we collide these stars with a coupling value of
$g_{a\gamma}=20m_{\mathrm{Pl}}^{-1}$, which ensures that the interaction
between the stars and the electromagnetic field is negligible before the
merger, since the coupling value lies in the stable region of the
$(M_{s},g_{a\gamma})$ plane (under the critical line
$g_{a\gamma}^{\mathrm{crit}}(M_{s})$), as demonstrated in Ref. [25]. In other
words, the decay of the original stars through the parametric resonance
process, as outlined in Section I, is forbidden.

However, the stars collide and merge to form a single highly excited star, as
shown in [31]. This merged star potentially possesses a critical coupling
below the coupling set in the simulation
($g_{a\gamma}=20m_{\mathrm{Pl}}^{-1}$). We note that evaluating the critical
coupling of this merged star analytically is not straightforward, as the star
is highly excited with some “breathing mode” oscillation patterns and
gravitational cooling by axion emission [24], and hence does not map directly
to the simple compact stars with critical line found in Ref. [25].
Nonetheless, approximating the merged star’s effective radius by the location
where the scalar energy density has dropped to $5\%$ of its central value, we
compute the ADM mass of the merged star, which gives us an estimate of its
critical coupling from the critical line in [25]. We find
$g_{a\gamma}^{\mathrm{crit}}\sim 16m_{\mathrm{Pl}}^{-1}$, which is below the
simulated value. Hence, we expect parametric resonance to kick in after the
merger and cause the merged star to decay.

Figure 3: Total simulation box energy in the scalar (solid line) and
electromagnetic field (dashed line) for the merger and electromagnetic decay
process. The merger timescale is shown in green, and the decay through
parametric resonance of the resulting highly excited merged star in grey. The
energy is in units of the original colliding stars’ mass. The electromagnetic
energy $\mathrm{EM}_{0}$ immediately before the decay has been subtracted from
the total; the decrease in the EM energy is due to the absorbing boundary
conditions. Strong amplification of the EM field during the decay of the
excited merged star is clearly visible ($\Delta E_{\gamma}$). The unphysical
bump in the total scalar energy corresponds to the lapse being driven to its
minimum value due to the large curvatures present at merger.

This is indeed what we observe in our simulations, with the merging of the
original stars shown in the upper set of panels of Fig. 2, corresponding to
$t=0m^{-1}$ to $t=340m^{-1}$, and the parametric decay of the merged star at
later times $t=630m^{-1}$ to $t=700m^{-1}$ in the lower set of panels. The
electromagnetic decay leaves behind a less compact remnant, and we discuss the
corresponding shutting down of the resonance process in Appendix A.

The same merger and decay process is shown in terms of the total energy in the
scalar and electromagnetic field in Fig. 3. The green area corresponds to the
merging process, while the grey area shows the time span of the parametric
resonance decay of the merged star. We have scaled the energies to the mass of
the original star ($M_{s}$), and hence the solid line representing the scalar
energy starts at $2M_{s}$. The total scalar energy in the simulation box stays
roughly constant before and after the merger of the stars, but decays through
the parametric resonance process after $t\sim 600m^{-1}$, as expected.
Accordingly, the electromagnetic energy in the box is strongly amplified
during parametric resonance. The decay process pictured in Fig. 3 is in
perfect agreement with the decay of single compact stars that we studied in
[25].

We use Sommerfeld boundary conditions for our simulation box and hence the
electromagnetic radiation is constantly being absorbed at the boundaries of
our box, which explains its decrease until the parametric resonance kicks in.
Furthermore, note the _total_ EM energy of the box is larger than that of the
scalar – this is due to the large box size we used in our simulations since
the EM seed wave fills the entire domain. The key physical effect is the EM
energy “bump” occuring at $t=600m^{-1}$, i.e. $\Delta E_{\gamma}\sim\Delta
E_{\mathrm{\phi}}$ where scalar energy is converted into EM energy via
parameteric resonance.

Figure 4: The power spectrum of the electromagnetic field at the start of the
simulation ($t=0m^{-1}$) and after the decay of the merged star
($t=700m^{-1}$). A clear excitation in all the EM components around the
characteristic frequency of the star, $k\sim 0.5m$, is visible. The $E_{y}$
component of the emission does not acquire energy at the seed frequency due to
our Sommerfeld boundary conditions.

The power spectrum of the electromagnetic field at the start of the simulation
and right after the decay of the merged star has happened is shown in Fig. 4.
The dashed black line demonstrates the original seed (polarised in the
x-direction), while the coloured lines corresponding to $t=700m^{-1}$ show all
the components of the EM field ($E_{x}$, $E_{y}$, $E_{z}$), with a clear
excitation around $k\sim 0.5m$, which corresponds to the characteristic
frequency of the star, and demonstrates parametric resonance. Similar to the
power of a single star decay discussed in [25], the EM power is approximately
isotropic. Interestingly, there is a small double peak around $k\sim 0.5m$,
which is likely due to a slight difference in excitation frequency of the
oscillations around the collision axis and that perpendicular to it.

##  IV Gravitational waves

In Fig. 5 we show the gravitational wave signal corresponding to the
electromagnetic decay of a single compact star with mass
$M_{s}=0.60m_{\mathrm{Pl}}^{2}/m$, using $g_{a\gamma}=16m_{\mathrm{Pl}}^{-1}$.
The two polarisation modes of GWs are captured through the real and imaginary
parts of the Newton-Penrose Weyl scalars $\Psi_{4}$ (see, e.g. Ref.[50]). Fig.
5 shows the real and imaginary part of the $l=2,m=0$ mode, although the
$l=2,m=\pm 2$ modes were observed to be of similar amplitudes. The signal is
extracted at $r=80m^{-1}$.

Figure 5: Gravitational wave waveform of the $l=2,m=0$ mode for the
electromagnetic decay of a single star. The waveform is extracted at
$r=80m^{-1}$ at refinement level 2. This gravitational wave, corresponding to
a burst of electromagnetic energy, does not exhibit a ringdown. We note that
the $l=2,m=\pm 2$ modes were of similar order of magnitude. Figure 6:
Gravitational wave waveforms for the merger and subsequent decay process,
showing the dominant quadrupole mode $l=2,m=0$. We note that the $l=2,m=\pm 2$
modes were of the same order of magnitude than the mode shown. The wavefronts
were extracted at radius $r=80m^{-1}$, refinement level 2. The first
excitation around $t=350m^{-1}$ in the real part corresponds to the merger
signal, and the second excitation around $t=750m^{-1}$ is the signal from the
electromagnetic decay of the star. In the case of the decay, both the real and
imaginary parts are excited, meaning that both gravitational polarisation
modes are present. A zoomed in picture of the signal from the merger of the
stars is depicted as a green wavefront, corresponding to the shaded green
region in the main plot.

As the decay of the star is not symmetric (for example, our electromagnetic
seed is polarised in one direction and propagating in another), both
gravitational polarisation modes are excited, corresponding to the real and
imaginary parts of the Weyl scalar. In contrast, the merger signal is
axisymmetric, and hence only a single polarisation mode is excited. We also
notice that this decay signal does not have a ringdown hierarchy like black
hole mergers, creating a distinct signature.

Moreover, in Fig. 6 we map the whole gravitational wave waveform for the
merger of two stars and the subsequent parametric decay process, outlined in
Section III. The gravitational wave signal from the merger at around $t\sim
350m^{-1}$ has already been studied in [31]. We note that the GW signal from
the decay of the merged star into photons around $t\sim 750m^{-1}$ closely
matches the signal from the decay of a single star shown in Fig. 5.

The frequency of the oscillations in the GW signal were computed by taking a
Fourier transform, and found to be around $\omega\sim m$ for the decay, which
matches the frequency of the EM flux through the same surface (see Section V),
and moreover the mass of the scalar field. The frequency match $\omega\sim m$
can be understood from linearising the axion-GW interaction. The frequency of
the merger signal, on the other hand, is around one order of magnitude lower,
$\omega\sim 0.1m$. For an axion mass of $m=10^{-12}\text{ eV}$, the decay
frequency corresponds to $f\sim 10^{3}\mathrm{Hz}$, implying that both the
decay and the merger frequencies lie close to the current peak sensitivity of
the LIGO band.

Figure 7: Cumulative gravitational energy in each $l$ mode for the decay of
the merged star (corresponding to the signal at $t\sim 750m^{-1}$ in Fig. 6).
The dashed lines indicate the maximum energy of the corresponding mode for the
single star decay. We note that the $l=2$ mode dominates, as expected, and
that the spectrum of the merged star is more varied than the single one, where
most of the energy is in the dominant $l=2$ mode.

Finally, we study the gravitational wave energy resulting from the
electromagnetic decay of the star. Fig. 7 shows the cumulative energy in the
gravitational wave signal from the decay of the merged star for the dominant
$l$ modes, extracted from the signal between $t=600m^{-1}$ to $t=900m^{-1}$.
This spans the decay process of the merged star. As expected from other
compact objects, the dominant mode is the quadrupole, $l=2$. Comparing these
modes to the modes excited by the decay of the single star (maximum energy in
each mode shown as dashed lines), we note that the spectrum of the merged star
has more power in higher multipoles as the merger excites asymmetric spatially
oscillating modes in the star (the single star, on the other hand, is
spherically symmetric before the decay process). Concretely, less than 1% of
the GW energy from the decay of a single star is found in modes $l>2$, whereas
the merged star has around 10% of the energy in the $l=3$ and $l=4$ modes.

The total gravitational energy ($E_{GW}$) from the decay is $\sim 0.03\%$ of
the total energy of the original stars ($2M_{s}$), comparable to the ratio in
black hole head-on mergers. We expect this efficiency to be similar for non-
head on mergers as it is driven by the decay process. On the other hand, the
electromagnetic decay of the merged axion star carries around ten times more
energy in gravitational waves than the merger of axion stars.

Our simulations also demonstrate that the total gravitational energy emitted
by the single star decay scales roughly linearly with compactness, with an
increase of over $90\%$ in the energy emitted from $C=0.10$ to $C=0.14$,
keeping the coupling constant fixed at $g_{a\gamma}=16m_{\mathrm{Pl}}^{-1}$.
However, as the value of the axion-photon coupling affects the decay process
significantly, a multi-parameter study would be required for a full
understanding of the energy scaling.

##  V Multimessenger signal

Figure 8: Multimessenger signal from the merger of two axion stars and
subsequent decay of the merged star, caused by an electromagnetic instability:
the gravitational power (top), total flux of electromagnetic energy extracted
over a sphere at the same radius as the GW signal, $r=80m^{-1}$ (middle), and
corresponding total flux of the scalar energy (bottom), where the power in
electromagnetic and scalar radiation is obtained by integrating the
corresponding flux. The different frequency wave modes in the scalar emission
are responsible for the long tail in the scalar flux, as lower frequency modes
travel slower away from the decaying star. The system starts shedding scalar
waves already due to the merger which happens around $t\sim 300m^{-1}$,
resulting in a complex pattern of scalar radiation after the electromagnetic
decay (around $t\sim 750m^{-1}$).

As the previous sections have explored in detail, the electromagnetic
instability of the axion star results in a multimessenger signal, with both an
electromagnetic and a gravitational component. The gravitational wavefront and
the electromagnetic energy ejected are expected to reach an observer away from
the source at the same time, travelling at the speed of light. This is
represented in Fig. 1, where the upper two panels represent the gravitational
power (in black) from the decay of a single star (corresponding to the
wavefront in Fig. 5), and the corresponding total flux of electromagnetic
energy (in orange), integrated over a sphere at the same radius
($r=80m^{-1}$).

The bottom panel (in blue) represents the scalar (axion) waves ejected by the
star, captured by the total scalar flux over the same sphere. These total
energy fluxes are calculated by integrating over the flux in Eq. 11 of Ref.
[33]. Firstly, we observe that the scalar energy is around two orders of
magnitude lower than its electromagnetic counterpart (the power units are the
same in both panels), hence being a subdominant effect. Secondly, the total
scalar flux is delayed compared to the EM and GW signals, as massive modes
travel slower than light 444To confirm the delay of the massive modes compared
to the EM and GW signal, we computed the expected time delay for the single
star decay using the relativistic dispersion relation and $k\sim 0.5m$, which
corresponds to the characteristic frequency of the star, assuming it matches
the maximum group velocity. This gives us $\Delta t\sim 100m^{-1}$, which
matches the time delay in Fig.1.. In addition, in the decay process, different
wave modes of the scalar field are excited as the star transitions from one
state to another. These scalar waves of different wavelengths travel at
different speeds, resulting in a long tail of total scalar flux as shown in
Fig. 1, as longer wavelength modes travel slower (recall that
$\omega^{2}=k^{2}+m^{2}$ for a free scalar wave) resulting in the
characteristic “dispersion” of the wave packet.

In Fig. 8 we show the corresponding GW, EM and scalar signals for the merger
simulations. The GW and EM signals from the decay phase are strikingly similar
to the single star case, with a relatively small GW bump from the merger. On
the other hand, the merged star emits scalar waves immediately following the
merger at $t\sim 300m^{-1}$, due to being in an excited state [24] ending in a
burst during the electromagnetic decay phase.

Finally, we comment on the energy scales of the system: the largest source of
emitted energy in the parametric, electromagnetic decay of an axion star is
the electromagnetic emission from the star, which is around two orders of
magnitude larger than the scalar emission. The gravitational energy emitted is
a further two orders of magnitude smaller than the scalar emission, hence
creating a clear hierarchy in the magnitude of the signals emitted. For a
scalar field mass of $m=10^{-12}\text{ eV}$, the total electromagnetic energy
ejected by the parametric decay of the merged star is around
$10^{49}\mathrm{J}$, which compares to the energy scale of a supernova.

##  VI Discussion and Conclusions

In this work we have studied, in fully non-linear simulations, the
multimessenger signal produced by the collision of two compact axion stars
close to the critical mass defined by the axion-photon coupling. The initial
merger produces a GW signal consistent with that expected for real scalar
boson stars [32, 31]. This signal is relatively small in amplitude, due to the
low compactness of $C<0.14$. The merger forms an excited, more massive axion
star that crosses the threshold of electromagnetic stability. This star
subsequently decays into electromagnetic radiation via parametric resonance,
triggered by ambient photons. The excited super-critical star first cools by
emission of axions, and then explodes, emitting axions, GWs, and photons. The
GWs from the decay are almost monochromatic, at a frequency one order of
magnitude larger than the GWs from collision, and do not follow a ringdown
hierarchy. GWs from the decay carry over one order of magnitude more energy
than those from the merger.

The length of time of the delay between merger and decay is fixed by the
coupling $g_{a\gamma}$, the mass of the merged star (which defines its
critical coupling), as well as the amplitude of the EM seed through a
logarithmic dependence, as we confirmed in [25]. The dominant energy in the
decay is carried by photons, since criticality is governed by parametric
resonance of the Chern-Simons interaction between axions and photons. The
energy in photons is two orders of magnitude larger than the energy in axions,
which itself is two orders of magnitude more than the energy in GWs.

This paper has been concerned with the theory and simulation of axion star
mergers and explosions, computing for the first time the decay channels. We
discuss briefly here the observational consequences, but leave full study to a
future work. By colliding two axion stars of mass
$M_{s}=0.36m_{\mathrm{Pl}}^{2}/m$, $C=0.03$ each, with axion parameters
$m=10^{-12}\text{ eV},g_{a\gamma}=20m_{\mathrm{Pl}}^{-1}$, selected to give
rise to GWs from merger and decay in the LIGO frequency band, we determined
that the GWs would be visible up to distances of $\sim 110\text{ Mpc}$ and
$270\text{ Mpc}$, for merger and decay respectively. The calculations follow
Eq. 13 in [44], assuming a detected strain $h=10^{-22}$, which is the current
LIGO peak sensitivity for the decay frequency range. If GW emission and photon
emission from the decay could both be observed, then with the right waveform
model and timing data this may allow the weaker merger GW signal to be
extracted from noise in post-processing.

An axion star explosion from higher mass axions produces high frequency GWs,
for example $m\approx 10^{-6}\text{ eV}$ gives GWs in the GHz band, although
the luminosity would be smaller as the axion star mass scales inversely to
$m$. Such high frequency GWs can be detected using microwave cavities [51]
(such as the ones commonly used in axion dark matter haloscopes) or phononic
materials [52] . We have not considered this possibility in detail for GWs
from axion star explosions, as the signal is likely too weak to give
sensitivity to large astronomical distances, but this may warrant further
investigation.

For GWs in the LIGO band, the EM radiation is emitted in the kHz band with
total energy comparable to a supernova. Some consequences of the EM radiation
emitted by an axion star explosion in this band are discussed in Refs. [53,
54]. However, for a multimessenger signal including GWs, there are several
important differences to consider for the photon emission compared to Refs.
[53, 54]. Compact stars giving a GW signal require a different formation
mechanism: the core halo-mass relation [18] for axion stars in dark matter
halos predicts that compact axion stars are extremely rare, if present at all.
One possible mechanism to form compact axion stars is outlined in Ref. [20],
which involves a spike in the axion primordial power spectrum, which could
emerge due to phase transitions and string network decay after inflation (e.g.
as axiton formation in Ref. [55]), or blue inflationary isocurvature [56]. In
such mechanisms, compact axion stars can form at rare density peaks, much like
primordial black holes. Another possibility involves an extended matter
dominated epoch [57]. The physics of an individual explosion will be the same
as studied in Refs. [53, 54], but the population level effects will be absent.
Furthermore, a multimessenger signal involving GWs will probe only the local
Universe at $D\lesssim 300\text{Mpc}$, rather than the early universe $z>10$
considered in Refs. [53, 54].

The EM emission from the explosion of an axion star in the kHz band is
efficiently absorbed by neutral hydrogen, causing it to ionize [58]. Thus, we
would observe the EM emission indirectly from the heating of the surrounding
gas. Following Ref. [53, 54] we find that in the Milky Way, taking values for
the free electron density and baryon temperature in virial equilibrium, the
absorption length for photons in the kHz band is short, on the order of
parsecs. Thus the photons from an exploding axion star in the local Universe
heat the surrounding gas and rapidly ionize and shock it; the resulting
emission in optical frequencies from the hot expanding gas may resemble a
supernova with no associated production of heavy elements.

It may also be possible to detect the semi-relativistic axions produced. Axion
dark matter haloscope data is normally analysed assuming a search for the cold
galactic axions of the standard halo model [59]. However, the physics of
axion-photon conversion can happen also for relativistic axions, such as those
of the “cosmic axion background” [60]. In particular, the broadband DMRadio
proposal [61] covers the correct frequency range for axions that are part of a
multimessenger signal of an axion star explosion with accompanying GWs in the
LIGO band.

## Acknowledgements

We would like to thank Lewis Croney, Miguel Escudero, Tamara Evstafyeva and
Thomas Helfer for useful discussions. We also thank members of the GRTL
Collaboration for technical support and help. LCJ is supported by a
studentship funded by the Science and Technologies Facilities Council (UK).
DJEM is supported by an Ernest Rutherford Fellowship from the Science and
Technologies Facilities Council (ST/T004037/1). This work used the
DiRAC@Durham Cosma facility managed by the Institute for Computational
Cosmology on behalf of the STFC DiRAC HPC Facility (www.dirac.ac.uk), under
DiRAC grant ACTP238. The equipment was funded by BEIS capital funding via STFC
capital grants ST/P002293/1, ST/R002371/1 and ST/S002502/1, Durham University
and STFC operations grant ST/R000832/1.

## References

  * Peccei and Quinn [1977] R. D. Peccei and H. R. Quinn, CP Conservation in the Presence of Instantons, Phys. Rev. Lett. 38, 1440 (1977).
  * Weinberg [1978] S. Weinberg, A New Light Boson?, Phys. Rev. Lett. 40, 223 (1978).
  * Wilczek [1978] F. Wilczek, Problem of Strong p and t Invariance in the Presence of Instantons, Phys. Rev. Lett. 40, 279 (1978).
  * Kim [1979] J. E. Kim, Weak Interaction Singlet and Strong CP Invariance, Phys. Rev. Lett. 43, 103 (1979).
  * Shifman  _et al._ [1980] M. A. Shifman, A. I. Vainshtein, and V. I. Zakharov, Can Confinement Ensure Natural CP Invariance of Strong Interactions?, Nucl. Phys. B 166, 493 (1980).
  * Dine  _et al._ [1981] M. Dine, W. Fischler, and M. Srednicki, A Simple Solution to the Strong CP Problem with a Harmless Axion, Phys. Lett. B 104, 199 (1981).
  * Zhitnitsky [1980] A. R. Zhitnitsky, On Possible Suppression of the Axion Hadron Interactions. (In Russian), Sov. J. Nucl. Phys. 31, 260 (1980).
  * Preskill  _et al._ [1983] J. Preskill, M. B. Wise, and F. Wilczek, Cosmology of the Invisible Axion, Phys. Lett. B 120, 127 (1983).
  * Abbott and Sikivie [1983] L. F. Abbott and P. Sikivie, A Cosmological Bound on the Invisible Axion, Phys. Lett. B 120, 133 (1983).
  * Dine and Fischler [1983] M. Dine and W. Fischler, The Not So Harmless Axion, Phys. Lett. B 120, 137 (1983).
  * Di Luzio  _et al._ [2020] L. Di Luzio, M. Giannotti, E. Nardi, and L. Visinelli, The landscape of QCD axion models, Phys. Rept. 870, 1 (2020), arXiv:2003.01100 [hep-ph] .
  * Turner [1983] M. S. Turner, Coherent Scalar Field Oscillations in an Expanding Universe, Phys. Rev. D 28, 1243 (1983).
  * Witten [1984] E. Witten, Some Properties of O(32) Superstrings, Phys. Lett. B 149, 351 (1984).
  * Khlopov  _et al._ [1985] M. Khlopov, B. A. Malomed, and I. B. Zeldovich, Gravitational instability of scalar fields and formation of primordial black holes, Mon. Not. Roy. Astron. Soc. 215, 575 (1985).
  * Svrcek and Witten [2006] P. Svrcek and E. Witten, Axions In String Theory, JHEP 06, 051, arXiv:hep-th/0605206 .
  * Arias  _et al._ [2012] P. Arias, D. Cadamuro, M. Goodsell, J. Jaeckel, J. Redondo, and A. Ringwald, WISPy Cold Dark Matter, JCAP 06, 013, arXiv:1201.5902 [hep-ph] .
  * Marsh [2016] D. J. E. Marsh, Axion Cosmology, Phys. Rept. 643, 1 (2016), arXiv:1510.07633 [astro-ph.CO] .
  * Schive  _et al._ [2014] H.-Y. Schive, T. Chiueh, and T. Broadhurst, Cosmic Structure as the Quantum Interference of a Coherent Dark Wave, Nature Phys. 10, 496 (2014), arXiv:1406.6586 [astro-ph.GA] .
  * Levkov  _et al._ [2018] D. G. Levkov, A. G. Panin, and I. I. Tkachev, Gravitational Bose-Einstein condensation in the kinetic regime, Phys. Rev. Lett. 121, 151301 (2018), arXiv:1804.05857 [astro-ph.CO] .
  * Widdicombe  _et al._ [2018] J. Y. Widdicombe, T. Helfer, D. J. E. Marsh, and E. A. Lim, Formation of Relativistic Axion Stars, JCAP 10, 005, arXiv:1806.09367 [astro-ph.CO] .
  * Eggemeier and Niemeyer [2019] B. Eggemeier and J. C. Niemeyer, Formation and mass growth of axion stars in axion miniclusters, Phys. Rev. D 100, 063528 (2019), arXiv:1906.01348 [astro-ph.CO] .
  * Chen  _et al._ [2021] J. Chen, X. Du, E. W. Lentz, D. J. E. Marsh, and J. C. Niemeyer, New insights into the formation and growth of boson stars in dark matter halos, Phys. Rev. D 104, 083022 (2021), arXiv:2011.01333 [astro-ph.CO] .
  * Seidel and Suen [1991] E. Seidel and W. M. Suen, Oscillating soliton stars, Phys. Rev. Lett. 66, 1659 (1991).
  * Seidel and Suen [1994] E. Seidel and W.-M. Suen, Formation of solitonic stars through gravitational cooling, Phys. Rev. Lett. 72, 2516 (1994), arXiv:gr-qc/9309015 .
  * Chung-Jukko  _et al._ [2023] L. M. Chung-Jukko, E. A. Lim, D. J. E. Marsh, J. C. Aurrekoetxea, E. de Jong, and B.-X. Ge, Electromagnetic instability of compact axion stars, Phys. Rev. D 108, L061302 (2023), arXiv:2302.10100 [hep-ph] .
  * Hertzberg and Schiappacasse [2018] M. P. Hertzberg and E. D. Schiappacasse, Dark Matter Axion Clump Resonance of Photons, JCAP 11, 004, arXiv:1805.00430 [hep-ph] .
  * Levkov  _et al._ [2020] D. G. Levkov, A. G. Panin, and I. I. Tkachev, Radio-emission of axion stars, Phys. Rev. D 102, 023501 (2020), arXiv:2004.05179 [astro-ph.CO] .
  * Kephart and Weiler [1987] T. W. Kephart and T. J. Weiler, LUMINOUS AXION CLUSTERS, Phys. Rev. Lett. 58, 171 (1987).
  * Boskovic  _et al._ [2019] M. Boskovic, R. Brito, V. Cardoso, T. Ikeda, and H. Witek, Axionic instabilities and new black hole solutions, Phys. Rev. D 99, 035006 (2019), arXiv:1811.04945 [gr-qc] .
  * Ikeda  _et al._ [2019] T. Ikeda, R. Brito, and V. Cardoso, Blasts of Light from Axions, Phys. Rev. Lett. 122, 081101 (2019), arXiv:1811.04950 [gr-qc] .
  * Helfer  _et al._ [2019a] T. Helfer, E. A. Lim, M. A. G. Garcia, and M. A. Amin, Gravitational Wave Emission from Collisions of Compact Scalar Solitons, Phys. Rev. D 99, 044046 (2019a), arXiv:1802.06733 [gr-qc] .
  * Giudice  _et al._ [2016] G. F. Giudice, M. McCullough, and A. Urbano, Hunting for Dark Particles with Gravitational Waves, JCAP 10, 001, arXiv:1605.01209 [hep-ph] .
  * Clough [2021] K. Clough, Continuity equations for general matter: applications in numerical relativity, Class. Quant. Grav. 38, 167001 (2021), arXiv:2104.13420 [gr-qc] .
  * Clough  _et al._ [2015] K. Clough, P. Figueras, H. Finkel, M. Kunesch, E. A. Lim, and S. Tunyasuvunakool, GRChombo : Numerical Relativity with Adaptive Mesh Refinement, Class. Quant. Grav. 32, 245011 (2015), arXiv:1503.03436 [gr-qc] .
  * Andrade  _et al._ [2021] T. Andrade _et al._ , GRChombo: An adaptable numerical relativity code for fundamental physics, J. Open Source Softw. 6, 3703 (2021), arXiv:2201.03458 [gr-qc] .
  * Radia  _et al._ [2022] M. Radia, U. Sperhake, A. Drew, K. Clough, P. Figueras, E. A. Lim, J. L. Ripley, J. C. Aurrekoetxea, T. França, and T. Helfer, Lessons for adaptive mesh refinement in numerical relativity, Class. Quant. Grav. 39, 135006 (2022), arXiv:2112.10567 [gr-qc] .
  * Gorbar  _et al._ [2021] E. V. Gorbar, K. Schmitz, O. O. Sobol, and S. I. Vilchinskii, Gauge-field production during axion inflation in the gradient expansion formalism, Phys. Rev. D 104, 123504 (2021), arXiv:2109.01651 [hep-ph] .
  * Helfer  _et al._ [2017] T. Helfer, D. J. E. Marsh, K. Clough, M. Fairbairn, E. A. Lim, and R. Becerril, Black hole formation from axion stars, JCAP 03, 055, arXiv:1609.04724 [astro-ph.CO] .
  * Alcubierre  _et al._ [2003] M. Alcubierre, R. Becerril, S. F. Guzman, T. Matos, D. Nunez, and L. A. Urena-Lopez, Numerical studies of Phi**2 oscillatons, Class. Quant. Grav. 20, 2883 (2003), arXiv:gr-qc/0301105 .
  * Michel and Moss [2018] F. Michel and I. G. Moss, Relativistic collapse of axion stars, Phys. Lett. B 785, 9 (2018), arXiv:1802.10085 [gr-qc] .
  * Urena-Lopez  _et al._ [2002] L. A. Urena-Lopez, T. Matos, and R. Becerril, Inside oscillatons, Class. Quant. Grav. 19, 6259 (2002).
  * Urena-Lopez [2002] L. A. Urena-Lopez, Oscillatons revisited, Class. Quant. Grav. 19, 2617 (2002), arXiv:gr-qc/0104093 .
  * Zilhão  _et al._ [2015] M. Zilhão, H. Witek, and V. Cardoso, Nonlinear interactions between black holes and Proca fields, Class. Quant. Grav. 32, 234003 (2015), arXiv:1505.00797 [gr-qc] .
  * Helfer  _et al._ [2019b] T. Helfer, J. C. Aurrekoetxea, and E. A. Lim, Cosmic String Loop Collapse in Full General Relativity, Phys. Rev. D 99, 104028 (2019b), arXiv:1808.06678 [gr-qc] .
  * Gundlach  _et al._ [2005] C. Gundlach, J. M. Martin-Garcia, G. Calabrese, and I. Hinder, Constraint damping in the Z4 formulation and harmonic gauge, Class. Quant. Grav. 22, 3767 (2005), arXiv:gr-qc/0504114 .
  * Palenzuela  _et al._ [2010] C. Palenzuela, L. Lehner, and S. Yoshida, Understanding possible electromagnetic counterparts to loud gravitational wave events: Binary black hole effects on electromagnetic fields, Phys. Rev. D 81, 084007 (2010), arXiv:0911.3889 [gr-qc] .
  * Hilditch [2013] D. Hilditch, An Introduction to Well-posedness and Free-evolution, Int. J. Mod. Phys. A 28, 1340015 (2013), arXiv:1309.2012 [gr-qc] .
  * Helfer  _et al._ [2022] T. Helfer, U. Sperhake, R. Croft, M. Radia, B.-X. Ge, and E. A. Lim, Malaise and remedy of binary boson-star initial data, Class. Quant. Grav. 39, 074001 (2022), arXiv:2108.11995 [gr-qc] .
  * Evstafyeva  _et al._ [2023] T. Evstafyeva, U. Sperhake, T. Helfer, R. Croft, M. Radia, B.-X. Ge, and E. A. Lim, Unequal-mass boson-star binaries: initial data and merger dynamics, Class. Quant. Grav. 40, 085009 (2023), arXiv:2212.08023 [gr-qc] .
  * Baumgarte and Shapiro [2010] T. W. Baumgarte and S. L. Shapiro,  _Numerical relativity : solving Einstein’s equations on the computer_ (Cambridge University Press, Cambridge, 2010).
  * Berlin  _et al._ [2022] A. Berlin, D. Blas, R. Tito D’Agnolo, S. A. R. Ellis, R. Harnik, Y. Kahn, and J. Schütte-Engel, Detecting high-frequency gravitational waves with microwave cavities, Phys. Rev. D 105, 116011 (2022), arXiv:2112.11465 [hep-ph] .
  * Kahn  _et al._ [2023] Y. Kahn, J. Schütte-Engel, and T. Trickle, Searching for High Frequency Gravitational Waves with Phonons, (2023), arXiv:2311.17147 [hep-ph] .
  * Escudero  _et al._ [2023] M. Escudero, C. K. Pooni, M. Fairbairn, D. Blas, X. Du, and D. J. E. Marsh, Axion Star Explosions: A New Source for Axion Indirect Detection, (2023), arXiv:2302.10206 [hep-ph] .
  * Du  _et al._ [2023] X. Du, D. J. E. Marsh, M. Escudero, A. Benson, D. Blas, C. K. Pooni, and M. Fairbairn, Soliton Merger Rates and Enhanced Axion Dark Matter Decay, (2023), arXiv:2301.09769 [astro-ph.CO] .
  * Vaquero  _et al._ [2019] A. Vaquero, J. Redondo, and J. Stadler, Early seeds of axion miniclusters, JCAP 04, 012, arXiv:1809.09241 [astro-ph.CO] .
  * Kasuya and Kawasaki [2009] S. Kasuya and M. Kawasaki, Axion isocurvature fluctuations with extremely blue spectrum, Phys. Rev. D 80, 023516 (2009), arXiv:0904.3800 [astro-ph.CO] .
  * Erickcek and Sigurdson [2011] A. L. Erickcek and K. Sigurdson, Reheating Effects in the Matter Power Spectrum and Implications for Substructure, Phys. Rev. D 84, 083503 (2011), arXiv:1106.0536 [astro-ph.CO] .
  * Bolliet  _et al._ [2021] B. Bolliet, J. Chluba, and R. Battye, Spectral distortion constraints on photon injection from low-mass decaying particles, Mon. Not. Roy. Astron. Soc. 507, 3148 (2021), arXiv:2012.07292 [astro-ph.CO] .
  * Brubaker  _et al._ [2017] B. M. Brubaker, L. Zhong, S. K. Lamoreaux, K. W. Lehnert, and K. A. van Bibber, HAYSTAC axion search analysis procedure, Phys. Rev. D 96, 123008 (2017), arXiv:1706.08388 [astro-ph.IM] .
  * Dror  _et al._ [2021] J. A. Dror, H. Murayama, and N. L. Rodd, Cosmic axion background, Phys. Rev. D 103, 115004 (2021), [Erratum: Phys.Rev.D 106, 119902 (2022)], arXiv:2101.09287 [hep-ph] .
  * Brouwer  _et al._ [2022] L. Brouwer _et al._ (DMRadio), Proposal for a definitive search for GUT-scale QCD axions, Phys. Rev. D 106, 112003 (2022), arXiv:2203.11246 [hep-ex] .
  * Croft  _et al._ [2023] R. Croft, T. Helfer, B.-X. Ge, M. Radia, T. Evstafyeva, E. A. Lim, U. Sperhake, and K. Clough, The gravitational afterglow of boson stars, Class. Quant. Grav. 40, 065001 (2023), arXiv:2207.05690 [gr-qc] .

##  Appendix A Shutting down of resonance

Axion stars can experience an electromagnetic instability if their
characteristic frequency lies within the electromagnetic resonance band, as we
first demonstrated in [25], and further in this work for the case of a merged
star. In Fig. 9, we show how the parametric resonance process responsible for
the decay of the star into electromagnetic radiation does not dissipate the
whole star, but shuts down due to the dilution of the star.

Figure 9: Radial profiles of the scalar energy density $\rho_{\phi}$ of a
single star decay (top) and the merger process (bottom). The original set up
is shown in black (single star on the top and two stars on the bottom plot),
and the remnant after the electromagnetic decay has happened, in green. The
dashed green line is a magnification of the remnant energy density, to
demonstrate its energy distribution, which indeed has a larger radius and
broader shape than the corresponding original star, as the star has diluted in
the parametric resonance process. The magnification factors are 100 and 10,
for the single and merged stars, respectively. The merged star created trough
the collision of the two sigle stars in the bottom plot, which subsequently
decays through parametric resonance, is shown in orange. We note that these
radial profiles have been taken through the center of the simulation box, with
the other two coordinates at $512m^{-1}$, and we have chosen the $x$-axis for
the merger process, as the stars are originally aligned along $x$.

As the star enters parametric resonance, its shedding of energy causes the
star to dilute, increasing its size. Since the characteristic frequency of the
star is linked to its size, the change in frequency means that the star slides
off the resonance frequency band and hence the process stops, leaving behind a
less compact, dilute remnant. This is seen in the upper plot of Fig. 9 for the
single star, where the original star represented by the black line has decayed
into a significantly less energetic, less compact star, shown in green. The
magnification of the energy density of the remnant star, shown with a dashed
green line, demonstrates its broader energy distribution. This indicates a
larger radius, and hence different characteristic frequency to the original
star.

Similarly, in the bottom plot the original two stars are shown in black, the
merged star created by their collision in orange, and the remnant after the
electromagnetic decay in green, with the magnification of the remnant energy
with a dashed green line. The corresponding times, shown in the legend, can be
mapped to the process represented as a whole in Fig. 3. As the merged star is
created by the collision of the single stars, its energy distribution has some
additional features compared to the single star, which are discussed in
Section III. These features are also present in the remnant, which exhibits a
broader energy distribution and hence a larger radius than the original stars,
as expected.

##  Appendix B Convergence testing

###  B.1 General convergence

We demonstrate the convergence of our simulation by monitoring the ratio of
the difference in the scalar field $\phi$ value between the high and medium,
and the medium and low resolution runs, summing over a stencil of 7 points,
centered at the middle of our simulation box. The box had a size
$L=1024m^{-1}$ and a maximum of 7 refinement levels. The number of coarse grid
points used for each resolution was $N^{3}=256^{3}$, $N^{3}=384^{3}$, and
$N^{3}=512^{3}$ (from lower to higher resolution). The region of convergence,
defined by the ratio being below 1, is demonstrated by the grey area in Fig.
10, alongside the ratio values (shown with the black solid line), as a
function of time.

We note that there is a peak at around $t\sim 400m^{-1}$ which does not fall
in the region of convergence. This is due to the regridding forced on the
simulation by the gravitational wave extraction tool. Our numerical relativity
code GRChombo includes an adaptive mesh refinement system that automatically
adds refinement levels in specific regions of the simulation box, based on set
conditions and threshold values. However, the extraction of gravitational
waves at a set radius and refinement level overrides the original, smooth
refinement process by forcing a specific refinement layer on the simulation at
each time step. As the boundaries of each layer introduce numerical errors to
the simulation, we postulate that it is this regridding process that is
responsible for the peak in our convergence plot. This was shown by running
the same simulation without GW extraction, in which case the peak at around
$t\sim 400m^{-1}$ disappeared. Moreover, the differences in the $\phi$ value
as shown in Fig. 10 get smaller as the simulation progresses, which is
consistent with convergence.

Figure 10: Convergence testing using a stencil of 7 points centered in the
middle of the simulation box. We plot the ratio of the difference in the
scalar field $\phi$ value between the high and medium, and the medium and low
resolution runs, summing over the stencil. Hence, all data points within the
grey region indicate convergence. The simulation box size is $L=1024m^{-1}$,
and low, medium, and high resolution runs have $N^{3}=256^{3}$,
$N^{3}=384^{3}$, and $N^{3}=512^{3}$ number of coarse grid points,
respectively.

###  B.2 Gravitational wave convergence

In addition to testing the convergence of the bulk simulation, as explained
above, we monitored the convergence of the gravitational wave waveforms, which
is demonstrated in Fig. 11. The figure shows the imaginary part of the
$l=2,m=0$ mode of the Weyl scalar, but we note that the real part followed a
similar pattern. Convergence can be clearly seen by eye, as the difference
between the high and medium resolution waveforms is smaller than that between
the medium and the low resolutions.

Figure 11: Convergence testing of the gravitational wave waveform for the
decay of a single star. The number of coarse grid points used for each
resolution is shown in the legend. Convergence can be seen by eye as the
higher resolutions (in green and red dashed) are closer together than the
lower ones (red dashed and black). The waveforms were extracted at refinement
level 2 at $r=80m^{-1}$.

Here we demonstrate the convergence for the GW signal from the decay of a
single, stable star of mass $M_{s}=0.60m_{\mathrm{Pl}}^{2}/m$. The merger GW
signal (shown in Fig. 6) had better convergence. While Fig. 11 shows the
signal from the decay of a single (non-excited) star, we also confirmed the
convergence for the full combined merger and decay process. Due to the high
frequency of the gravitational wave waveform of the decay, it was challenging
to keep track of it, and indeed we found linear drift in the waveforms between
different resolution runs for the decay after merger, corresponding to 3-4
wavelengths in time. However, the frequencies of the different resolution
waveforms matched. Overall, despite the simulation and gravitational waveforms
converging, we found some variation in the waveform amplitudes.

###  B.3 Regridding

In order to extract an overall clean waveform, we had to tune our regridding
for the merger and the decay signal separately, using our adaptive mesh
refinement code. The goal of regridding is to add higher resolution meshes in
the simulation box areas where the interesting physics, and hence the need for
higher resolution, is happening. For our simulations, this means that the star
is resolved with a higher resolution grid than the empty space around it.
(However, as mentioned above, the extraction of gravitational waves away from
the star, closer to the boundaries of the box, means we also force a higher
resolution at the extraction radius). In practice, this tuning of the
regridding means that we had to change the threshold values for the variables
with respect to which the refinement is happening (in our case, the scalar
field $\phi$ and the conformal factor $\chi$).

To get the full signal as shown in Fig. 6, we had to decrease the thresholds
between the merger and the decay signal. The merger signal exhibited some
numerical noise that was reduced with lesser regridding, while the high
oscillation frequency of the decay signal required more regridding and higher
extraction levels. The decrease in the regridding thresholds (hence, more
frequent regridding), moving from the merger to the decay, was $18\%$ and
$23\%$, for $\phi$ and $\chi$, respectively. This regridding change had no
effect on the overall evolution of the system, and demonstrates the difficulty
of gravitational wave extraction from our oscillating scalar star.

Finally, we note that we found some linear drift in time in the energy fluxes
and gravitational wave signals from the simulation due to this change of
regridding in the middle of the run (between the merger and the decay of the
merged star). Hence, the gravitational signal in Fig. 8 has been corrected to
match the electromagnetic one, as we confirmed by running the simulation from
the middle that those two indeed line up (our full GW signal displayed in Fig.
6 has the regridding change in the middle, as explained above).

##  Appendix C Gravitational luminosity

The gravitational wave power (luminosity) was calculated according to the
well-known equation (see e.g. [62] for similar use):

$\frac{dE}{dt}=\lim_{r\rightarrow\infty}\frac{1}{16\pi}\sum_{l,m}\left|\int_{t_{0}}^{t}r\Psi_{4}^{lm}dt^{\prime}\right|^{2},$
(5)

where we used the modes up to $l=6$ for the decay signal of the star and $l=2$
for the merger signal of the star. To deal with the numerical error causing a
cumulative drift in the total gravitational energy of the system, the integral
of $r\Psi_{4}^{lm}$ was processed with a Butterworth high-pass filter (5th
order), as the numerical error introduces low-frequency modes.


# Gaussian Processes for Music Audio Modelling and Content Analysis

Pablo A. Alvarado, Dan Stowell
Queen Mary University of London

###### Abstract

Real music signals are highly variable, yet they have strong statistical
structure. Prior information about the underlying physical mechanisms by which
sounds are generated and rules by which complex sound structure is constructed
(notes, chords, a complete musical score), can be naturally unified using
Bayesian modelling techniques. Typically algorithms for Automatic Music
Transcription independently carry out individual tasks such as multiple-F0
detection and beat tracking. The challenge remains to perform joint estimation
of all parameters. We present a Bayesian approach for modelling music audio,
and content analysis. The proposed methodology based on Gaussian processes
seeks joint estimation of multiple music concepts by incorporating into the
kernel prior information about non-stationary behaviour, dynamics, and rich
spectral content present in the modelled music signal. We illustrate the
benefits of this approach via two tasks: pitch estimation, and inferring
missing segments in a polyphonic audio recording.

##  1 Introduction

In music information research, the aim of audio content analysis is to
estimate musical concepts which are present but hidden in the audio data [17].
With this purpose, different signal processing techniques are applied to music
signals for extracting useful information and descriptors related to the
musical concepts. Here, musical concepts refers to parameters related to
written music, such as pitch, melody, chords, onset, beat, tempo and rhythm.
Then, perhaps the most general application is one which involves the
prediction of several musical dimensions, that of recovering the score of a
music track given only the audio signal [10]. This is known as automatic music
transcription (AMT) [5].

AMT refers to extraction of a human readable and interpretable description
from a recording of a music performance. We refer to polyphonic AMT in cases
where more than a single musical pitch plays at a given time instant. The
general task of interest is to infer automatically a musical notation, such as
the traditional western music notation, listing the pitch values of notes,
corresponding timestamps and other expressive information in a given audio
signal of a performance [7]. Transcribing polyphonic music is a nontrivial
task, especially in its more unconstrained form when the task is performed on
an arbitrary acoustical input, and music transcription remains a very
challenging problem [4].

Real music signals are highly variable, but nevertheless they have strong
statistical structure. Prior information about the underlying structures, such
as knowledge of the physical mechanisms by which sounds are generated, and
knowledge about the rules by which complex sound structure is compiled (notes,
chords, a complete musical score), can be naturally unified using Bayesian
hierarchical modelling techniques. This allows the formulation of highly
structured probabilistic models [7]. On the other hand, typically, algorithms
for AMT are developed independently to carry out individual tasks such as
multiple-F0 detection, beat tracking and instrument recognition. The challenge
remains to combine these algorithms, to perform joint estimation of all
parameters [5].

We present the design, implementation, and results of experiments of an
alternative Bayesian approach for audio content analysis on monophonic, and
polyphonic music signals with the possibility of being used for AMT. We use
Gaussian process (GP) models for jointly uncovering music concepts from audio,
by introducing a direct connection between the music concepts and the model
hyper-parameters. The proposed methodology allows to incorporate in the model
prior information about physical or mechanistic behaviour, nonstationarity,
time dynamics (local periodicity, and non constant amplitude envelope),
spectral harmonic content, and musical structure, latent in the modelled music
signal. Specifically in the context of music informatics, we present kernels
that embody a probabilistic model of music notes as time-limited harmonic
signals with onsets and offsets. A comparison with related work is provided in
section 4.3. We illustrate the benefits of this approach via two tasks: pitch
estimation, and inferring missing segments in a polyphonic audio recording. As
part of working towards a high-resolution AMT system, out method correctly
estimates polyphonic pitch while performing these sample-level tasks.

##  2 Gaussian process regression for music signals

Gaussian process-based machine learning is a powerful Bayesian paradigm for
nonparametric nonlinear regression and classification [14]. Gaussian Processes
(GPs) can be defined as distributions over functions such that any finite
number of function evaluations $\textbf{f}=[f(t_{1}),\cdots,f(t_{N})]$, have a
jointly normal distribution [12]. A GP is completely specified by its mean
function $\mu(t)=\mathbb{E}[f(t)]$ (in this work it is assumed to be
$\mu(t)=0$), and its kernel or covariance function

$\displaystyle
k(t,t^{\prime})=\mathbb{E}\left[(f(t)-\mu(t))(f(t^{\prime})-\mu(t^{\prime}))\right],$
(1)

where $k(t,t^{\prime})$ has hyper-parameters $\boldsymbol{\theta}$. We write
the Gaussian process as

$\displaystyle f(t)\sim\mathcal{GP}(\mu(t),k(t,t^{\prime})).$ (2)

The regression problem concerns the prediction of a continuous quantity [12],
here a function $f(t)$, given a data set
$\mathcal{D}=\left\\{(t_{i},y_{i})\right\\}_{i=1}^{N}$, where $y_{i}$ are
assumed as noisy measurements of $f(t)$ at typically regularly-spaced time
instants $t_{i}$ (though GP regression framework allows for irregular sampling
or missing data), i.e. $y_{i}=f(t_{i})+\epsilon_{i}$, where
$\epsilon_{i}\sim\mathcal{N}(0,\sigma_{\text{noise}}^{2})$. In GP regression
for mono channel audio signals, instead of estimating parameters
$\boldsymbol{\eta}$ of fixed-form functions
$f(t,\boldsymbol{\eta}):\mathbb{R}\mapsto\mathbb{R}$ where the time input
variable $t\in\mathbb{R}$, we model the whole function $f(t)$ as a GP. That
is, instead of putting a prior over the function parameters
$\boldsymbol{\eta}$, we introduce a prior over the function $f(t)$ itself
[16]. Learning in GP regression corresponds to computing the posterior
distribution over the function $f(t)$ conditioned on the observed data
$\textbf{y}=\left[y_{1},\cdots,y_{N}\right]^{\top}$ [15, 13].

The underlying idea in GP regression is that the correlation function
introduces dependences between function $f(t)$ values at different inputs.
Thus, the function values at the observed points give information also of the
unobserved points [14]. The structure of the kernel (1) captures high-level
properties of the unknown function $f(t)$, which in turn determines how the
model generalizes or extrapolates to new test time instants [9]. This is quite
useful because we can introduce prior knowledge about what we believe the
proprieties of music signals are, by choosing a proper kernel that reflects
those characteristics. In section 2.2 we study in more detail the design of
kernels.

###  2.1 Model definition

Under a non-parametric Bayesian regression approach using Gaussian processes
we are interested in calculating the posterior distribution over a stochastic
function evaluated at test points $\textbf{t}_{*}$, that is, the joint
distribution of the vector f observed only via noisy measurements y, then

$\displaystyle p(\textbf{f}|\textbf{y})=\frac{p(\textbf{y}|\textbf{f})\times
p(\textbf{f}|\boldsymbol{\theta})}{p(\textbf{y})},$ (3)

where $p(\textbf{y}|\textbf{f})$ corresponds to the likelihood,
$p(\textbf{f}|\boldsymbol{\theta})$ to the prior, $\boldsymbol{\theta}$ are
the model hyper-parameters (prior parameters), $p(\textbf{y})$ is the evidence
or marginal-likelihood, and $p(\textbf{f}|\textbf{y})$ is the posterior or
conditional predictive distribution. We describe each of these four
expressions in the next sections.

####  2.1.1 Likelihood

Assuming that conditioned on the $f(t_{i})$ the signal observations $y_{i}$
are i.i.d. (independent and identically distributed), then the joint
probability distribution of all the observations y follows a Gaussian
distribution corresponding to

$\displaystyle
p(\textbf{y}|\textbf{f},\sigma^{2}_{\text{noise}})=\mathcal{N}(\textbf{y}|\textbf{f},\sigma^{2}_{\text{noise}}\textbf{I}_{N}),$
(4)

where $f_{i}=f(t_{i})$ and $\textbf{I}_{N}$ is an identity matrix of size $N$.

####  2.1.2 Prior

Using the definition of GPs introduced at the beginning of this section, and
knowing that we have a finite set of corrupted observations y, then the finite
set of GP function evaluation values $\mathbf{f}$ follows a normal marginal
distribution $p(\textbf{f}|\boldsymbol{\theta})$ conditioned on the hyper-
parameters $\boldsymbol{\theta}$, whose mean is zero and whose covariance is
defined by a Gram matrix $\textbf{K}_{f}$, this is

$\displaystyle
p(\textbf{f}|\boldsymbol{\theta})=\mathcal{N}(\textbf{f}|\boldsymbol{0},\textbf{K}_{f}),$
(5)

where the covariance matrix is calculated using (1), i.e.
$[\textbf{K}_{f}]_{i,j}=k(t_{i},t_{j})$ [6].

####  2.1.3 Marginal-Likelihood

The marginal-likelihood (or evidence) $p(\textbf{y})$ mentioned before in (3)
is the integral of the likelihood times the prior [12]

$\displaystyle p(\textbf{y})=\int
p(\textbf{y}|\textbf{f})p(\textbf{f}|\boldsymbol{\theta})\textrm{d}\textbf{f}.$
(6)

Since the likelihood $p(\textbf{y}|\textbf{f})$ and the prior
$p(\textbf{f}|\boldsymbol{\theta})$ are multivariate Gaussian distributions,
we can calculate directly the integral in (6). Using the properties of the
normal distribution [6] for marginal and conditional normal distributions, we
obtain

$\displaystyle
p(\textbf{y})=\mathcal{N}(\textbf{y}|\boldsymbol{0},\textbf{K}_{y}),$ (7)

where the values in the matrix
$\textbf{K}_{y}=\textbf{K}_{f}+\sigma^{2}_{\text{noise}}\textbf{I}$ depend on
the hyper-parameters $\boldsymbol{\theta}$ (we have included
$\sigma^{2}_{\text{noise}}$ in the hyper-parameters vector). The reason it is
called the marginal likelihood, rather than just likelihood, is because we
have marginalized out the latent Gaussian vector f [11].

####  2.1.4 Posterior

The computation of the posterior distribution of the Gaussian process
conditioned on the set of measurements y and estimation of the parameters
$\boldsymbol{\theta}$ of the covariance function of the process correspond to
learning in this non-parametric model [14]. Using the properties of Gaussian
distribution [6, 12], the posterior has the form

$\displaystyle
p(\textbf{f}|\textbf{y})=\mathcal{N}(\textbf{y}|\boldsymbol{\mu}_{\text{pos}},\textbf{K}_{\text{pos}}),$
(8)

where the posterior mean is
$\boldsymbol{\mu}_{\text{pos}}=\textbf{K}_{f}\textbf{K}_{y}^{-1}\textbf{f}$,
and the posterior covariance matrix is
$\textbf{K}_{\text{pos}}=\textbf{K}_{f}-\textbf{K}_{f}^{\top}\textbf{K}_{y}\textbf{K}_{f}$.

###  2.2 Kernel design

The covariance function (1) used for computing the prior distribution (5)
allows us to introduce in the model all the knowledge and beliefs we have
about the properties of the data. We are trying to model music signals, and
some of the broad properties of audio signals are non-stationarity, rich
spectral content, dynamics (locally periodic, non constant amplitude
envelope), mechanistic behaviour, and music structure. Therefore we seek
covariance functions that can describe or reflect these properties.

One powerful technique for constructing new kernels is to build them out of
simpler kernels as building blocks [19, 6]. Two useful properties we can use
to build valid kernels are as follows:
$k(t,t^{\prime})=\phi(t)k_{1}(t,t^{\prime})\phi(t^{\prime}),$ and
$k(t,t^{\prime})=k_{1}(t,t^{\prime})+k_{2}(t,t^{\prime}),$ where $\phi(\cdot)$
is any function. Other properties can be found in [6]. We use these properties
for building non-stationary covariance functions.

####  2.2.1 Kernels for describing non-stationarity

To construct non-stationary kernels we combine basic stationary covariance
functions. We use change-windows in order to be able to model notes or sound
events which are not continuously active but have a beginning and an ending in
the music signal. As in [9] we define a change-window by multiplying two
sigmoid functions, that is

$\displaystyle\phi(t)=\frac{1}{1+e^{-\varsigma(t-\alpha)}}\times\frac{1}{1+e^{-\varsigma(\beta-t)}},$
(9)

where $\varsigma$ determine how fast the change-window rises to its maximum
value or falls to zero, whereas $\alpha,\beta$ defines the onset and the
offset of the change-window respectively. The parameters of the change-windows
are directly related with the location, onset and offset of the sound events
we want to model. In the present work we will use manually-specified
onset/offset locations.

To illustrate the construction of a non-stationary kernel, an example model
with two change-windows and its corresponding covariance functions is
developed here. We assume a GP $f(t)\sim\mathcal{GP}(0,k_{f}(t,t^{\prime}))$
described by a linear combination of other two GPs
$f_{i}(t)\sim\mathcal{GP}(0,k_{i}(t,t^{\prime}))$ for $i=1,2$, each one
weighted by its corresponding change-window $\phi_{i}(t)$, i.e.
$f(t)=\phi_{1}(t)f_{1}(t)+\phi_{2}(t)f_{2}(t).$ Using the definition of
covariance function (1) we can calculate the kernel $k_{f}(t,t^{\prime})$ for
the complete process $f(t)$ as follows:
$k_{f}(t,t^{\prime})=\phi_{1}(t)k_{1}(t,t^{\prime})\phi_{1}(t^{\prime})+\phi_{2}(t)k_{2}(t,t^{\prime})\phi_{2}(t^{\prime}),$
where $k_{1}(t,t^{\prime})$ and $k_{2}(t,t^{\prime})$ are the covariance
functions for the GPs $f_{1}(\cdot)$ and $f_{2}(\cdot)$ respectively. This
kernel configuration was used for generating the two samples shown in Figure
1(b). In both cases $k_{1}(t,t^{\prime})$ and $k_{2}(t,t^{\prime})$ have the
harmonic kernel form (16) that we describe shortly (section 2.2.4). Figure
1(a) contains the shape of the two change-windows used ($\phi_{1}(t)$ blue,
$\phi_{2}(t)$ red). From Figure 1 we see that with the proposed methodology we
are able to describe non-stationary functions.

(a) Two change-windows.

(b) Two function realizations.

Figure 1: Example of GP model with nonstationary kernel.

####  2.2.2 Kernel general form for $M$ change-windows

The previous example could be generalized for modelling the complete process
$f(t)$ as a linear combination of $M$ random process, representing each one a
note or sound event. In this way

$\displaystyle f(t)=\sum_{m=1}^{M}\phi_{m}(t)f_{m}(t),$ (10)

where each Gaussian process $[f_{1},f_{2},\cdots,f_{M}]$ is independent with
respect to each other. It is important to highlight that $M$ is directly
related with the number of notes or sound events in the signal. On the other
hand, $\phi_{m}(t)$ are respectively the change-windows or weight functions
that allow a specific GP $f_{m}(t)$ to appear or vanish in certain parts of
the input space (time). In this way, the general expression for the covariance
function $k_{f}(t,t^{\prime})$ is given by

$\displaystyle
k_{f}(t,t^{\prime})=\sum_{m=1}^{M}\phi_{m}(t)k_{m}(t,t^{\prime})\phi_{m}(t^{\prime}).$
(11)

We see that the overall process has a kernel consisting of a linear
combination of the corresponding covariance functions of every subprocess.

####  2.2.3 Rich spectral content kernel property

In the previous section we described the general form of the proposed kernel.
Here we address the structure of the covariance functions
$k_{m}(t,t^{\prime})$ in (11). We assume that each GP $f_{m}(t)$ in (10) is
stationary. A random process is stationary (wide sense stationary WSS) if its
mean is constant, and its kernel is a covariance function of
$\tau=t-t^{\prime}$, then we can write $k(t,t^{\prime})=k(\tau)$ [18, 12]. It
can be shown that the spectral density or power spectrum $S(s)$ of a WSS
process corresponds to the Fourier transform (FT) of the covariance function,
that is

$\displaystyle S(s)=\int_{-\infty}^{\infty}k(\tau)e^{-js\tau}\textrm{d}\tau,$
(12)

thus

$\displaystyle
k(\tau)=\frac{1}{2\pi}\int_{-\infty}^{\infty}S(s)e^{js\tau}\textrm{d}s.$ (13)

This is known as the Wiener-Khintchine theorem [12, 18]. Taking that into
account, we can do frequency-domain analysis for several covariance functions
and decide which kernel is more appropriate for modelling the spectral content
of music signals. Again we illustrate these concepts by an example, we compare
the FT of the exponentiated quadratic covariance function
$k_{\text{EQ}}(\tau)$ with the FT of the exponentiated cosine kernel
$k_{\text{EC}}(\tau)$, that is

$\displaystyle k_{\text{EQ}}(\tau)$
$\displaystyle=\sigma^{2}\exp\left(-\frac{\tau^{2}}{2l^{2}}\right),$ (14)
$\displaystyle k_{\text{EC}}(\tau)$
$\displaystyle=\sigma^{2}\exp\left[z\cos(\omega\tau)\right],$ (15)

where to keep values up to one we set $\sigma^{2}=\exp(-z)$ in (15), as well
as $z=2$, and $\omega=2\pi 6$. For (14) $l=0.01$. covariance function (14) is
probably the most widely-used kernel within the kernel machines field, because
the GP with a exponentiated-quadratic covariance function is very smooth [12].
Figures 2(a) and 2(d) represent the form of these kernels. The spectral
density of the GP with kernel (14) contains only low frequency components and
does not have any harmonic structure (Figure 2(b)). Figure 2(c) shows two
different realizations sampled from a GP with $k_{\text{EQ}}(\tau)$ kernel,
these functions evolve smoothly without any periodic or harmonic properties.
On the other hand, the spectral density (Figure 2(e)) of the exponentiated
cosine kernel (15) presents a DC part, as well as components at the natural
frequency $6$Hz and harmonics (integer numbers of $6$Hz). Figure 2(f) shows
two functions sampled from a GP with covariance function (15). These
realizations present constant amplitude-envelope and periodic properties with
a fundamental frequency together with several harmonics.

(a) $k_{\text{EQ}}(\tau)$

(b) FT of $k_{\text{EQ}}(\tau)$

(c) Samples using $k_{\text{EQ}}$

(d) $k_{\text{EC}}(\tau)$

(e) FT of $k_{\text{EC}}(\tau)$

(f) Samples using $k_{\text{EC}}$

(g) $k_{\text{ECQ}}(\tau)$

(h) FT of $k_{\text{ECQ}}(\tau)$

(i) Samples using $k_{\text{ECQ}}$

Figure 2:  Frequency analysis of kernels (14), (15), (16) depicted in (a),
(d), (g) respectively. (b), (e), (h) are their corresponding FT. (c), (f), (i)
show sampled functions.

####  2.2.4 Time dynamics kernel property

In order to allow periodic kernels to describe functions where the amplitude
envelope changes in time, we introduce a modification of expression (15). To
do so, we multiply the periodic kernel with (14). The resulting covariance
function, called exponentiated-cosine-quadratic corresponds to

$\displaystyle
k_{\text{ECQ}}(\tau)=\sigma^{2}\exp\left[z\cos(\omega\tau)-\frac{\tau^{2}}{2l^{2}}\right],$
(16)

where we assume $\sigma^{2}=\exp(-z)$, $\omega=2\pi 6$, $z=5$ and $l=4$.
Figure 2(h) depicts the FT of (16). We see that its spectral density keeps
similar to the one obtained for (15) (see Figure 2(e)). But the realizations
sampled from a GP with this covariance function (Figure 2(i)) show a smooth
variation in the amplitude envelope, and also maintain the properties
described by the previous kernel (15), i.e. a periodic structure with natural
frequency and harmonics. This covariance function (16) seems to be more
appropriate for modelling music signals in comparison with the two kernels
presented previously ((14)-(15)).

##  3 Empirical Evaluation

Experiments were done over real audio. We evaluated different kernel
configurations on a pitch estimation task, and on a missing data imputation
task. All experiments assume we previously know the number of change-windows
and its locations. In the pitch estimation task all the parameters of the
covariance function are known, except those related with the fundamental
frequency of each sound event, i.e. the value of $\omega_{m}$ in (15) and (16)
when using these kernels in the general model (10). Thus, we focus on
optimizing only these model hyperparameters from the data. In the missing data
imputation task the score of the modelled piece of music audio is used for
tuning manually the model hyperparameters.

###  3.1 Data

In this study we used two short audio excerpts, in order to explore the
method, so that we can efficiently fit models and search in the hyperparameter
space. The excerpt used for pitch estimation experiments corresponds to $0.7$
seconds of the song Black Chicken 37 by Buena Vista Social Club. This segment
of audio contains three notes of a bass melody (Figure 3(a)). In the missing
data imputation task we used polyphonic audio corresponding to $1.14$ seconds
of Chopin’s Nocturne Op. 15 No. 1, where more than one note occur at the same
time. The segments of signal in red in Figure 3(b) represent gaps of missing
data. We reduced the sample frequency of both audio excerpts from $44.1$KHz to
$8$KHz.

(a) Signal used for pitch estimation.

(b) Signal used for filling missing-data gaps.

Figure 3: (a) analysed audio (blue line), change-windows (dashed lines). (b)
observed data (blue line), missing-data gaps (red line), change-windows
(dashed lines).

###  3.2 Hyper-parameter estimation

To infer the hyper-parameters we consider an empirical Bayes approach, which
allow us to use continuous optimization methods. We maximize the marginal
likelihood. This moves us up one level of the Bayesian hierarchy, and reduces
the chances of overfitting [11]. Given an expression for the log marginal
likelihood and its partial derivatives, we can estimate the kernel parameters
using any standard gradient-based optimizer [11]. A gradient descent method
was used for optimization.

##  4 Results and Discussion

###  4.1 Pitch estimation

For the pitch estimation task we tested two different models with kernels
(15), and (16) respectively. We performed hyperparameters learning using all
the observed signal. This is because in this experiment rather than evaluating
the prediction of the trained models, we were interested in the accuracy of
pitch estimation. Covariance function (14) does not have any parameter we can
link to the fundamental frequency of each sound event, that is why we omitted
it here.

####  4.1.1 Results using $k_{\text{EC}}(\tau)$

We performed regression on the signal shown in Figure 3(a) using the kernel
(15). Figure 4(a) shows the posterior mean of the predictive distribution
after training (blue continuous line). The black circle points correspond to
observed data. We see the trained model is able to estimate the pitch for each
sound event with a RMS error of $0.6282$ semitones (Table 1). On the other
hand, the amplitude-envelope evolution of the signal is beyond the scope of
the structure that this kernel can model. This is because this covariance
function can only describe constant amplitude-envelope, periodic signals, with
a fundamental frequency and several harmonics (see Figure 2(f)).

(a) Observations (dots), and posterior mean (continuous line) using (15).

(b) Observations (dots), and posterior mean (continuous line) using (16).

Figure 4: Posterior mean for the pitch estimation experiments. (a) using
$k_{\text{EQ}}(\tau)$, and (b) using $k_{\text{EQC}}(\tau)$.

####  4.1.2 Results using $k_{\text{ECQ}}(\tau)$

To face the issue of modelling time dynamics we modified the previous
covariance function (15), by multiplying it with an exponentiated quadratic
kernel (14). This allows to “smooth” the strictly periodic behaviour of (15).
The resulting kernel corresponds to (16). From Figure 4(b) we see that
although the posterior mean of the predictive distribution does not exactly
fit the data, the model is able to learn the pitch of each of the three sound
events with a smaller RMS error (Table 1), as well as the time dynamics or
variations in the amplitude envelope of the signal.

Table 1: Pitch estimation RMS error (semitones). kernel | RMS
---|---
$k_{\text{EQ}}(\tau)$ | $--$
$k_{\text{EC}}(\tau)$ | $0.6282$
$k_{\text{EQC}}(\tau)$ | $0.1075$

###  4.2 Filling gaps of missing data in audio

We compared three different models predicting missing-data gaps. We studied
kernels (14), (15), and (16). In Figure 3(b) first gap (red segment) contains
the transient (onset and attack [3]) of a sound event, whereas the second gap
is located in a more stable segment of the data (smooth decay). Figures
5(a)-5(b) depict the prediction using (14). These figures correspond to zoom
in small sections of the signal where the gaps occur (Figure 3(b)). We see
that the model using this kernel overfits the data, i.e. the posterior mean
(blue line) fits all the observed data (black dots) with high confidence (grey
shaded area), but the confidence decreases and the prediction is quite poor in
the input space zones where the data is not available (red dots). Also, we see
that the model using (14) does not expect any periodic behaviour in the gaps.
The RMS error for both gaps is presented in Table 2.

(a) Prediction on transient gap using $k_{\text{EQ}(\tau)}$.

(b) Prediction on smooth decay gap using $k_{\text{EQ}(\tau)}$.

(c) Prediction on transient gap using $k_{\text{EC}(\tau)}$.

(d) Prediction on smooth decay gap using $k_{\text{EC}(\tau)}$.

(e) Prediction on transient gap using $k_{\text{EQC}(\tau)}$.

(f) Prediction on smooth decay gap using $k_{\text{EQC}(\tau)}$.

Figure 5: Zoom in a portion of missing-data gaps. In each figure the
continuous blue line represent the posterior mean, grey shaded areas
correspond to the posterior variance, red dots are missing data, whereas black
dots are observed data.

Figures 5(c)-5(d) show the prediction using covariance function (15). In the
transient gap (Figure 5(c)) the posterior mean (blue line) does not follows
the data, this is because transients are short intervals during which the
signal evolves in a nonstationary, nontrivial and unpredictable way [3].
opposite to this, the model using kernel (15) can only describe the behaviour
of constant amplitude-envelope periodic stochastic functions. In the second
gap (Figure 5(d)) the posterior mean describes properly the periodic behaviour
of the data, but it does not follow the amplitude-envelope of the
observations. This is because this covariance function is able to describe
periodic functions that have several harmonic components. The drawback of this
kernel is that it assumes constant the amplitude of the periodic stochastic
functions that describes. These different performance on the prediction is
reflected on the RMS error obtained for each gap (Table 2).

Results using (16) are presented in Figures 5(e)-5(f). We see that in Figure
5(f) the posterior mean describes properly the periodic behaviour and
amplitude envelope smooth evolution of the modelled signal. We observe that
prediction on the decay gap using (16) is closer to the actual data (red dots)
than the results obtained with (15) as well as (14). This is reflected in the
smallest RMS error in table 2. This is because (16) allows to describe
periodic functions that have several harmonic components and time-varying
amplitude envelope. On the other hand, the prediction performance reduces for
the transient gap (Figure 5(e)). In order to model the onset, attack and decay
of a sound event, covariance function (16) could be modified for modelling
nonstationary amplitude envelope evolution.

Table 2: Filling gaps prediction RMS error. kernel | RMS transient gap | RMS decay gap
---|---|---
$k_{\text{EQ}}(\tau)$ | $0.2265$ | $0.3172$
$k_{\text{EC}}(\tau)$ | $0.2143$ | $0.0964$
$k_{\text{EQC}}(\tau)$ | $0.0912$ | $0.0355$

###  4.3 Related work

In [20] GPs are used for time-frequency analysis as probabilistic inference.
Natural signals are assumed to be formed by the superposition of distinct
time-frequency components, with the analytic goal being to infer these
components by applying Bayes’ rule [20]. GPs have also been used for
underdetermined audio source separation. In [8] the mixture signal is modelled
as a linear combination of independent convolved versions of latent GPs or
sources. The model splits the mixture signal in frames also considered
independent, by using weight-functions. Thus each source is modelled as a
series of concatenated locally stationary frames, each one with its
corresponding covariance function. With this assumption the resulting signal
is supposed to be non-stationary [8]. On the other hand, despite the approach
we present also assumes the latent GPs $f_{m}$ in (10) as non-correlated, the
observed signal is not framed into independent segments. Instead of using
weight-functions that act over the observed data, we introduce change-windows
$\phi_{m}$ influencing each latent GP ending up with latent processes
representing specific sound events that happen at certain segments of time.
Therefore the proposed model keeps the correlation between the observations
throughout all the signal. That is what allows to make prediction in gaps of
missing data (section 4.2). GPs have been used also for estimating spectral
envelope and fundamental frequency of singing voice [21], and for time-domain
audio source separation [22].

##  5 Conclusions

In this article we discussed a Gaussian processes regression framework for
modelling music audio. We compared different models in pitch estimation as
well as in prediction of missing data. We showed which kernels were more
appropriate for describing properties of music signals, specifically:
nonstationarity, dynamics, and spectral harmonic content. The advantage of
this approach is that by designing a proper kernel we can introduce prior
knowledge and beliefs about the properties of music signals, and use all that
prior information to improve prediction. The presented work could be extended
using efficient representations of GPs in order to model larger audio signals.
Other kernels could be studied, as the spectral mixture for modelling harmonic
content [2], and Latent Force models [1] for describing mechanistic
characteristics of the signal.

## References

  * [1] M.A. Alvarez, D. Luengo, and N.D. Lawrence.  Linear latent force models using gaussian processes.  Pattern Analysis and Machine Intelligence, IEEE Transactions on, 35(11):2693–2705, Nov 2013.
  * [2] Ryan Prescott Adams Andrew Gordon Wilson.  Gaussian process kernels for pattern discovery and extrapolation.  International Conference on Machine Learning (ICML), 2013., 2013\.
  * [3] J. P. Bello, L. Daudet, S. Abdallah, C. Duxbury, M. Davies, and M. B. Sandler.  A tutorial on onset detection in music signals.  IEEE Transactions on Speech and Audio Processing, 13(5):1035–1047, Sept 2005.
  * [4] E. Benetos et al.  Automatic music transcription: Breaking the glass ceiling.  In Proceedings of the 13th International Society for Music Information Retrieval Conference, pages 379–384, 2012.
  * [5] E. Benetos et al.  Automatic music transcription: challenges and future directions.  Journal of Intelligent Information Systems, 41(3):407–434, 2013\.
  * [6] C. Bishop.  Pattern Recognition and Machine Learning.  Springer-Verlag, New York, 2006.
  * [7] A. Cemgil et al.  Bayesian statistical methods for audio and music processing.  The Oxford Handbook of Applied Bayesian Analysis, pages 1–45, 2010\.
  * [8] A. Liutkus, R. Badeau, and G. Richard.  Gaussian processes for underdetermined source separation.  IEEE Transactions on Signal Processing, 59(7):3155–3167, July 2011\.
  * [9] J. Lloyd et al.  Automatic construction and Natural-Language description of nonparametric regression models.  In Proceedings of the 28th AAAI Conference on Artificial Intelligence, pages 1242–1250, 2014.
  * [10] M. Muller et al.  Signal processing for music analysis.  IEEE Journal of Selected Topics in Signal Processing, 5(6):1088–1110, 2011.
  * [11] K. Murphy.  Machine Learning: A Probabilistic Perspective.  MIT Press, 2012.
  * [12] C. Rasmussen and C. Williams.  Gaussian Processes for Machine Learning.  MIT Press, 2005.
  * [13] M. Ebden S. Reece N. Gibson S. Aigrain. S. Roberts, M. Osborne.  Gaussian processes for timeseries modelling.  2012\.
  * [14] S. Sarkka, A. Solin, and J. Hartikainen.  Spatiotemporal learning via infinite-dimensional bayesian filtering and smoothing: A look at gaussian process regression through kalman filtering.  IEEE Signal Processing Magazine, 30(4):51–61, July 2013.
  * [15] Simo Särkkä.  Artificial Neural Networks and Machine Learning – ICANN 2011: 21st International Conference on Artificial Neural Networks, Espoo, Finland, June 14-17, 2011, Proceedings, Part II, chapter Linear Operators and Stochastic Partial Differential Equations in Gaussian Process Regression, pages 151–158.  Springer Berlin Heidelberg, Berlin, Heidelberg, 2011.
  * [16] Simo Särkkä and Arno Solin.  Image Analysis: 18th Scandinavian Conference, SCIA 2013, Espoo, Finland, June 17-20, 2013. Proceedings, chapter Continuous-Space Gaussian Process Regression and Generalized Wiener Filtering with Application to Learning Curves, pages 172–181.  Springer Berlin Heidelberg, Berlin, Heidelberg, 2013.
  * [17] X. Serra et al.  Roadmap for Music Information Research.  Creative Commons BY-NC-ND 3.0 license, 2013.
  * [18] K. Shanmugan and A. Breipohl.  Random Signals: Detection, Estimation and Data Analysis.  Wiley, 1988.
  * [19] J. Shawe and N. Cristianini.  Kernel Methods for Pattern Analysis.  Cambridge University Press, New York, 2004.
  * [20] R. E. Turner and M. Sahani.  Time-frequency analysis as probabilistic inference.  IEEE Transactions on Signal Processing, 62(23):6171–6183, Dec 2014\.
  * [21] K. Yoshii and M. Goto.  Infinite kernel linear prediction for joint estimation of spectral envelope and fundamental frequency.  In 2013 IEEE International Conference on Acoustics, Speech and Signal Processing, 2013.
  * [22] K. Yoshii, R. Tomioka, D. Mochihashi, and M. Goto.  Beyond nmf: time-domain audio source separation without phase reconstruction.  In Proceedings of the 14th International Society for Music Information Retrieval Conference, November 4-8 2013.