# Probing Quarkyonic Matter in Neutron Stars with the Bayesian Nuclear-Physics
Multi-Messenger Astrophysics Framework

Peter T. H. Pang<EMAIL_ADDRESS>Nikhef, Science Park 105, 1098 XG Amsterdam,
The Netherlands  Institute for Gravitational and Subatomic Physics (GRASP),
Utrecht University, Princetonplein 1, 3584 CC Utrecht, The Netherlands  Lars
Sivertsen<EMAIL_ADDRESS>Department of Physics and Astronomy, Iowa State
University, Ames, IA 50010  Rahul Somasundaram  Theoretical Division, Los
Alamos National Laboratory, Los Alamos, NM 87545, USA  Department of Physics,
Syracuse University, Syracuse, NY 13244, USA  Tim Dietrich  Institute for
Physics and Astronomy, University of Potsdam, D-14476 Potsdam, Germany  Max
Planck Institute for Gravitational Physics (Albert Einstein Institute), Am
Mühlenberg 1, Potsdam 14476, Germany  Srimoyee Sen  Department of Physics and
Astronomy, Iowa State University, Ames, IA 50010  Ingo Tews  Theoretical
Division, Los Alamos National Laboratory, Los Alamos, NM 87545, USA  Michael
Coughlin  School of Physics and Astronomy, University of Minnesota,
Minneapolis, Minnesota 55455, USA  Chris Van Den Broeck  Nikhef, Science Park
105, 1098 XG Amsterdam, The Netherlands  Institute for Gravitational and
Subatomic Physics (GRASP), Utrecht University, Princetonplein 1, 3584 CC
Utrecht, The Netherlands

###### Abstract

The interior of neutron stars contains matter at the highest densities
realized in our Universe. Interestingly, theoretical studies of dense matter,
in combination with the existence of two solar mass neutron stars, indicate
that the speed of sound $c_{s}$ has to increase to values well above the
conformal limit ($c_{s}^{2}\sim 1/3$) before decreasing again at higher
densities. The decrease could be explained by either a strong first-order
phase transition or a cross-over transition from hadronic to quark matter. The
latter scenario leads to a pronounced peak in the speed of sound reaching
values above the conformal limit, naturally explaining the inferred behavior.
In this work, we use the Nuclear-Physics Multi-Messenger Astrophysics
framework NMMA to compare predictions of the quarkyonic matter model with
astrophysical observations of neutron stars, with the goal of constraining
model parameters. Assuming quarkyonic matter to be realized within neutron
stars, we find that there can be a significant amount of quarks inside the
core of neutron stars with masses in the two solar mass range, amounting to up
to $\sim 0.13M_{\odot}$, contributing $\sim 5.9\%$ of the total mass.
Furthermore, for the quarkyonic matter model investigated here, the radius of
a $1.4M_{\odot}$ neutron star would be
$13.44^{+1.69}_{-1.54}(13.54^{+1.02}_{-1.04})$ km, at $95\%$ credibility,
without (with) the inclusion of AT2017gfo.

††preprint: LA-UR-23-29118††preprint: NP3M-P2300017

##  I Introduction

Measurements of neutron stars’ masses, radii, and deformabilities serve as
probes of the nuclear equation of state (EOS) of neutron-rich matter at high
density and low temperature Ozel and Freire (2016). This part of the QCD phase
diagram can not be accessed in lattice QCD simulations due to the sign
problem, which plagues the importance sampling of gauge configurations. Thus,
neutron star measurements remain one of the few tools to explore the behavior
of dense cold matter. For instance, the feasibility of constraining the EOS
via measurement of neutron stars’ tidal deformabilities with gravitational
waves has been demonstrated Flanagan and Hinderer (2008); Del Pozzo  _et al._
(2013). Constraints on the radius $R$ of neutron stars obtained from
gravitational-wave data suggest that $R<13.5$ km Annala  _et al._ (2018);
Abbott  _et al._ (2018); De  _et al._ (2018); Tews  _et al._ (2018a), which in
turn implies that the EOS of matter is soft at nuclear densities, of the order
of the nuclear saturation density $n_{\rm sat}\simeq 0.16$ fm-3. On the other
hand, observations of heavy pulsars indicate that the EOS needs to be
sufficiently stiff at higher densities in order to support neutron stars with
masses of more than twice the solar mass Watts  _et al._ (2016); Demorest  _et
al._ (2010); Antoniadis  _et al._ (2013); Drischler  _et al._ (2021);
Raaijmakers  _et al._ (2020).

This soft-stiff behavior of the EOS also reveals an interesting feature of the
speed of sound ($c_{s}$) in dense matter as a function of density. The speed
of sound can be calculated with uncertainty estimates at nuclear densities
using chiral effective theory (EFT) Hebeler  _et al._ (2013); Tews  _et al._
(2018b); Drischler  _et al._ (2020); Keller  _et al._ (2023) and at very high
density using perturbative QCD Gorda  _et al._ (2018, 2021, 2023a). Thus, we
know that the speed of sound $c_{s}$ is small, $c_{s}^{2}\ll 1$, at low
densities, and approaches the conformal bound of $c_{s}=1/3$ from below at
very high densities. At intermediate densities, where neither chiral EFT nor
perturbative QCD are applicable, the speed of sound could, in principle, be a
monotonic or nonmonotonic function of baryon density. The soft-stiff EOS
obtained from neutron-star measurements reveals that the speed of sound, in
fact, is likely a non-monotonic function of baryon density, violating the
conformal bound at a few times the nuclear saturation density Bedaque and
Steiner (2015); Tews  _et al._ (2018b). More specifically, $c_{s}^{2}$
increases monotonically from nuclear saturation density to a few times the
nuclear saturation density, overshooting the conformal bound of $1/3$ and
reaching a peak at intermediate densities. With even larger densities,
$c_{s}^{2}$ then decreases below the conformal bound, reaching a minimum and
eventually increasing to asymptotically approach the conformal bound from
below. Given that we do not have the tools to obtain the EOS of dense matter
from direct calculations in QCD, it is important to assess various models
which can replicate this non-monotonic behavior of the speed of sound and
simultaneously describe mass and radius measurements of neutron stars.

One such model is “quarkyonic matter” McLerran and Reddy (2019); Jeong  _et
al._ (2020); Sen and Sivertsen (2021) in which a special arrangement of
nucleons and quarks in a combined Fermi sphere forces a peak in the sound
speed. In the quarkyonic model, the nucleon and quark degrees of freedom (DOF)
are described by a single Fermi distribution function, as illustrated in Fig.
1. At high densities, the Fermi momentum for the baryons is large and the DOF
deep within the Fermi sphere are Pauli blocked. Creating a particle-hole
excitation from deep within the Fermi sea requires large energy and momentum,
and so these DOF can be regarded as weakly interacting. Since QCD is
asymptotically free, we expect the existence of quark matter at high densities
where quarks behave as nearly free particles, thus motivating treating the DOF
deep within the Fermi Sphere as quarks. The DOF near the Fermi surface can,
however, be excited with low energy and momentum transfers. Thus, confining
forces should remain important. This motivates treating the DOF near the Fermi
surface as nucleons arising from quark correlations.

At low densities, the radius of the inner quark Fermi sphere in the quarkyonic
model is zero. As the density increases, the inner quark Fermi sphere starts
forming at some threshold density. The peak in the speed of sound arises right
at those densities where the inner quark Fermi sphere first appears.
Interestingly, the introduction of quarks in other models is typically
accompanied by a first-order phase transition. A first-order transition leads
to a sudden decrease in the pressure gradient, or equivalently, a sudden drop
in the speed of sound and a softening of the EOS. In the quarkyonic picture,
however, the transition to quark matter happens instead through a crossover
with a milder decrease in the pressure gradient, allowing for higher pressures
at high densities. Hence, in this model, the EOS remains stiff. The quarkyonic
model, therefore, differs from other models of quark matter in this regard.

Although quarkyonic models have been successful in producing a peak in the
speed of sound, there have not been any studies to infer the quarkyonic-matter
model parameters from astrophysical data. In particular, given an
unconstrained quarkyonic model, we would like to investigate what we can learn
about the model parameters from astrophysical data and what the quark content
is in neutron stars. In this paper, we answer these questions using Bayesian
analysis and adopting a quarkyonic excluded-volume model proposed in Ref. Sen
and Sivertsen (2021).

The astrophysical observations considered here are the gravitational-wave
detections of the binary neutron star mergers GW170817 Abbott  _et al._ (2017)
and GW190425 Abbott  _et al._ (2020) by Advanced LIGO Aasi  _et al._ (2015)
and Advanced Virgo Acernese  _et al._ (2015), the GW170817-associated
electromagnetic counterparts AT2017gfo GBM (2017), the NICER observations on
PSR J0030+0451Miller  _et al._ (2019); Riley  _et al._ (2019) and PSR
J0740+6620 Miller  _et al._ (2021); Riley  _et al._ (2021), and the radio
observations on PSR J0348+4032 Antoniadis  _et al._ (2013) and PSR J1614-2230
Arzoumanian  _et al._ (2018).

This paper is organized as follows: We start by describing the employed
quarkyonic-matter model Sen and Sivertsen (2021) in Section II.1. We then
describe our Bayesian inference approach in Section II.2, and give a short
description of our implementation in Section II.3. Finally, we present our
results followed by our conclusions in Sections III and IV.

##  II Methods

###  II.1 Quarkyonic-Matter EOS

A non-dynamic quarkyonic matter model was first proposed in Ref. McLerran and
Reddy (2019), followed by a dynamic model for quarkyonic matter introduced in
Ref. Jeong  _et al._ (2020). In this dynamic model, both hadrons and quarks
appear as quasi-particles and are described using a single Fermi distribution
which has an outer shell of nucleons and an inner sphere of quarks. The baryon
density of nucleons $n_{\rm B}^{\rm N}$ for a particular configuration of
quarkyonic matter is computed by evaluating the momentum space volume of the
outer Fermi shell, whereas the quark baryon density $n_{\rm B}^{\rm Q}$ is
computed by evaluating the volume of the inner sphere. The total baryon
density is given by the sum of the two: $n_{\rm B}=n_{\rm B}^{\rm N}+n_{\rm
B}^{\rm Q}$. Thus, for a fixed baryon density, one can vary the shell width
made of the nucleons (or the radius of the inner quark Fermi sphere) while
maintaining a constant $n_{\rm B}$, to probe different configurations of
quarkyonic matter for that particular total baryon density. The dynamical
model in Ref. Jeong  _et al._ (2020) also introduces an energy density
functional for the systems, defined as the sum of individual energy density
functionals of the nucleons and quarks. By minimizing the energy density
functional for the quarkyonic Fermi distribution with respect to the shell
radius of the quarkyonic phase space, one can find the equilibrium quarkyonic
matter configuration for a fixed baryon density. One can then construct the
EOS of dense matter by computing the energy density for the equilibrium
configuration as a function of density. Note that the dynamical model of Ref.
Jeong  _et al._ (2020) considers only symmetric nuclear matter. This model was
extended to account for the astrophysically more relevant system of pure
neutron matter in Ref. Sen and Sivertsen (2021). We use this extended model in
this work. In the next subsection, we describe the construction of the
relevant energy density functionals, starting with the neutrons, and then
moving on to the quarks.

Figure 1: Illustration of the Fermi sphere for quarkyonic matter at low baryon
number densities before the quark onset (left panel) and at high densities
after the quark onset (right panel). The quark Fermi momentum in dynamic
quarkyonic-matter models well below $n_{0}$ is negligible and not shown in
this figure.

####  II.1.1 Neutron part of the EOS

In the excluded volume model for quarkyonic neutron matter, we consider
neutrons to be hard spheres of volume $v_{h}$ exhibiting two spin DOFs. This
allows us to define an effective density for neutrons. If we insert $N$
neutrons in a volume $V$, then the unoccupied volume $V_{\textrm{ex}}$ is
given by

$V_{\textrm{ex}}=V-Nv_{h}\equiv V-\frac{N}{n_{0}}\,.$ (1)

In the above, we have defined the hard core density of neutrons as
$n_{0}\equiv v_{h}^{-1}$. If we now add another neutron to this $N$-neutron
system, this extra neutron can only occupy a volume $v_{h}$ within the
available volume $V_{\textrm{ex}}$. This gives rise to an increased effective
density for the neutrons, which we will denote as the neutron excluded volume
density $n^{\textrm{N}}_{\textrm{ex}}$:

$n^{\textrm{N}}_{\textrm{ex}}\equiv\frac{N}{V_{\textrm{ex}}}=\frac{Vn_{\textrm{B}}^{\textrm{N}}}{V-\frac{N}{n_{0}}}=\frac{n_{\textrm{B}}^{\textrm{N}}}{1-\frac{N}{Vn_{0}}}=\frac{n_{\textrm{B}}^{\textrm{N}}}{1-\frac{n_{\textrm{B}}^{\textrm{N}}}{n_{0}}}\,.$
(2)

This implies that the neutron density $n_{\text{B}}^{\text{N}}$ is constrained
to be lower than the hard-core density $n_{0}$ since otherwise, the volume
occupied by neutrons would be larger than $V$. Because of this, the energy
density associated with nucleons will increase faster than that of a free
Fermi gas as a function of baryon number density. If we then consider quarks
to be free point particles, as we will describe in the next subsection, the
sum of the energy density contributions from quarks and neutrons will
eventually be lower if some of the baryon density is stored in quark DOF
compared to a pure neutron phase. This is the main idea behind this excluded
volume model: as we increase the number of neutrons in a constrained volume,
it eventually becomes energetically favorable for some of those neutrons to
split into their constituent quarks, leading to a smooth crossover between a
hadronic phase and a quark matter phase.

Keeping in mind that the neutron excluded volume density is the effective
density for the neutrons, we can write $n_{\textrm{ex}}^{\textrm{N}}$ in terms
of neutron Fermi momentum $\Delta$,

$n_{\textrm{ex}}^{\textrm{N}}=2\int_{0}^{\Delta}\frac{\textrm{d}^{3}k}{(2\pi)^{3}},$
(3)

where the factor of two accounts for the neutron’s two spin DOF. From this, it
follows that the non-interacting pure neutron energy density functional for
hard-core neutrons is given by

$\epsilon_{\textrm{N}}=2\Big{(}1-\frac{n_{\text{B}}^{\textrm{N}}}{n_{0}}\Big{)}\int_{0}^{\Delta}\frac{\textrm{d}^{3}k}{(2\pi)^{3}}\sqrt{M_{\textrm{N}}^{2}+k^{2}},$
(4)

where $M_{\textrm{N}}$ is the neutron mass.

To properly describe neutron stars, we will have to include neutron
interactions in our energy density functional. For this, we use a potential
ansatz inspired by Ref. McLerran and Reddy (2019), and include a nuclear
potential of the form

$V(n_{\textrm{B}}^{\textrm{N}})=\tilde{a}n_{\textrm{B}}^{\textrm{N}}\Big{(}\frac{n_{\textrm{B}}^{\textrm{N}}}{\rho_{0}}\Big{)}+\tilde{b}n_{\textrm{B}}^{\textrm{N}}\Big{(}\frac{n_{\textrm{B}}^{\textrm{N}}}{\rho_{0}}\Big{)}^{2}\,,$
(5)

where $\rho_{0}$ is the saturation density, and $\tilde{a}$ and $\tilde{b}$
are free parameters. The parameters $\tilde{a}<0$ and $\tilde{b}>0$ are to
mimic the short-range repulsion and long-range attraction between neutrons,
respectively, at densities $n_{\textrm{B}}\lesssim 2\rho_{0}$.

However, this is not sufficient to reproduce the EOS of low-density neutron
matter. The excluded volume energy density for the neutrons in Eq. (4) is
desirable at high densities to force quark onset, but the excluded volume
model has the undesirable feature of significantly modifying the low-density
neutron EOS as well. In order to correct for this effect, Ref. Sen and
Sivertsen (2021) proposed the following approach: We introduce a modified
potential $\tilde{V}(n_{\rm B}^{\rm N})$, and the energy density functional

$\epsilon_{\textrm{N}}=2\Big{[}1-\Big{(}\frac{n_{\textrm{B}}^{\textrm{N}}}{n_{0}}\Big{)}\Big{]}\int_{0}^{\Delta}\frac{\textrm{d}^{3}k}{(2\pi)^{3}}\sqrt{M_{\textrm{N}}^{2}+k^{2}}+\tilde{V}(n_{\textrm{B}}^{\textrm{N}}),$
(6)

to minimize the effect of the excluded volume at low density. $\tilde{V}$ is
designed such that in a low-density expansion in powers of ${n_{\rm B}^{\rm
N}}/{n_{0}}$, the excluded volume energy density of Eq. 6 mimics the energy
density of low density neutron matter given by

$\epsilon=2\int_{0}^{k_{F}}\frac{d^{3}k}{(2\pi)^{3}}\sqrt{k^{2}+M_{N}^{2}}+V(n_{N}^{N})$
(7)

up to corrections that go as some high power in the expansion parameter. Note
that, $n_{\rm B}^{N}=2\int_{0}^{k_{F}}\frac{d^{3}k}{(2\pi)^{3}}$. The model
proposed in Ref. Sen and Sivertsen (2021) demands that the two energy
densities match upto corrections that go as the $(n_{\rm B}^{N}/n_{0})^{3}$.
The corresponding expression for $\tilde{V}$ is given by

$\displaystyle\tilde{V}(n_{\textrm{B}}^{\textrm{N}})=V(n_{\textrm{B}}^{\textrm{N}})$
$\displaystyle-\frac{2}{3}\frac{(3\pi^{2}n_{\textrm{B}}^{\textrm{N}})^{5/3}}{10\pi^{2}M_{\textrm{N}}}\Big{(}\frac{n_{\textrm{B}}^{\textrm{N}}}{n_{0}}\Big{)}$
$\displaystyle-\frac{5}{9}\frac{(3\pi^{2}n_{\textrm{B}}^{\textrm{N}})^{5/3}}{10\pi^{2}M_{\textrm{N}}}\Big{(}\frac{n_{\textrm{B}}^{\textrm{N}}}{n_{0}}\Big{)}^{2}\,.$
(8)

####  II.1.2 Free quark part of the EOS

We will now discuss the energy density functional for the quarks. We treat the
quarks as almost free point particles of mass $M_{\rm q}=M_{\rm
N}/N_{c}=M_{\rm N}/3$, where $N_{c}$ is the number of quark colors and fixed
at $N_{c}=3$. The only effect of interactions of quarks is taken into account
through the constituent quark mass, $M_{\rm q}$. Here, we include two species
of quarks, up and down, with two spin DOFs. We can then express the
contribution to the baryon density from quarks in terms of the quark Fermi
momentum $k_{\textrm{F}}^{\textrm{q}}$, so that

$n_{\textrm{B}}^{\textrm{q}}=2\int_{0}^{k_{\textrm{F}}^{\textrm{q}}}\frac{\textrm{d}^{3}k}{(2\pi)^{3}},$
(9)

where q = {u,d} denotes the quark flavor, and the factor of two again comes
from the spin DOF. The quark energy density functional is then given by

$\epsilon_{\textrm{q}}=2\int_{0}^{k_{\textrm{F}}^{\textrm{q}}}\frac{\textrm{d}^{3}k}{(2\pi)^{3}}\sqrt{M_{\textrm{q}}^{2}+k^{2}}.$
(10)

Using this formulation of the quark part to the EOS, it was previously found
that the quark onset in this model happens too abruptly, leading to a
violation of causality and negative speeds of sound in some cases Jeong  _et
al._ (2020). Ref. Jeong  _et al._ (2020) proposed to cure this issue by
scaling the quark density of states by a factor

$g(k)=\frac{\sqrt{\Lambda^{2}+k^{2}}}{k}\,,$ (11)

where $\Lambda$ is a free parameter. This modification is motivated by having
the correct high energy behavior, $g(k)\to 1$ as $k\to\infty$, which is
consistent with a free Fermi gas. The regulator $\Lambda$ can be viewed as the
energy scale at which non-perturbative effects become unimportant. We should,
therefore, expect $\Lambda$ to be of the same order of magnitude as the QCD
confinement scale $\Lambda_{\textrm{QCD}}$.

However, Ref. Sen and Sivertsen (2021) pointed out that this modification to
the quark density of states, while curing the unphysical behavior of the speed
of sound, introduces a non-negligible density of quarks at all baryon
densities, which itself is unphysical. Ref. Sen and Sivertsen (2021) proposed
a remedy to this effect by giving a density dependence to $\Lambda$,
$\Lambda=\Lambda(n_{\textrm{B}})$, with

$\displaystyle\Lambda(n_{\textrm{B}})=\begin{cases}\Lambda_{0}\big{[}-20\big{(}\frac{n_{\textrm{B}}}{n_{0}}\big{)}^{7}+70\big{(}\frac{n_{\textrm{B}}}{n_{0}}\big{)}^{6}-84\big{(}\frac{n_{\textrm{B}}}{n_{0}}\big{)}^{5}+35\big{(}\frac{n_{\textrm{B}}}{n_{0}}\big{)}^{4}\big{]}\quad&\textrm{for}\quad
n_{\textrm{B}}<n_{0}\,,\\\ \Lambda_{0}&\textrm{for}\quad n_{\textrm{B}}\geq
n_{0}\,.\end{cases}$ (12)

This choice is motivated as follows: We know that the excluded volume
quarkyonic matter model does not produce quarks at low baryon density when the
regulator $\Lambda$ is set to zero. The need for a nonzero $\Lambda$ arises
only to remedy the unphysical behavior of the speed of sound at high
densities. Thus, if we introduce a density dependence in $\Lambda$ such that
it is zero at low baryon density and close to the QCD scale near the quark
onset density, we will eliminate quarks at low densities while also
eliminating the unphysical behavior of the speed of sound at high densities.
If we then want all thermodynamic properties to stay unaffected by $\Lambda$
at $n_{\textrm{B}}=0$, we need
$\Lambda(0)=\textrm{d}^{n}\Lambda(0)/\textrm{d}(n_{\textrm{B}})^{n}=0$ for
$n\leq 3$. Furthermore, we want $\Lambda(n_{\textrm{B}})$ to be in full effect
once a significant amount of quarks is present. It is then natural to impose
that $\Lambda(n_{\textrm{B}})=\Lambda_{0}$ for $n_{\textrm{B}}\geq n_{0}$,
where $\Lambda_{0}$ is a parameter of order $\Lambda_{\textrm{QCD}}$. We
further demand that
$\textrm{d}^{n}\Lambda(n_{\textrm{B}})/(\textrm{d}n_{\textrm{B}})^{n}$ is
continuous for $n\leq 3$ to ensure that chemical potential, pressure, and
speed of sound are all smooth functions of $n_{\textrm{B}}$. The polynomial
(12) was chosen as it is the lowest order polynomial that satisfies all these
criteria. However, as long as $\Lambda(n_{\textrm{B}})$ follows the properties
discussed above, its exact functional form is not found to be of significant
importance Sen and Sivertsen (2021). The quark energy density functional then
becomes

$\epsilon_{\textrm{q}}=2N_{\textrm{c}}\int_{0}^{k_{\text{F}}^{\textrm{q}}}\frac{\textrm{d}^{3}k}{(2\pi)^{3}}\sqrt{M_{\textrm{q}}^{2}+k^{2}}\frac{\sqrt{\Lambda^{2}(n_{\textrm{B}})+k^{2}}}{k}\,,$
(13)

where

$k_{\textrm{F}}^{\textrm{q}}=\sqrt{(3\pi^{2}n_{\textrm{B}}^{\textrm{q}}+\Lambda^{3}(n_{\textrm{B}}))^{2/3}-\Lambda^{2}(n_{\textrm{B}})}\,.$
(14)

We note that the introduction of a density dependent $\Lambda(n_{\text{B}})$
smears out the onset density for the quarks. Hence, they appear gradually as a
function of baryon density Sen and Sivertsen (2021). However, the quark
density is essentially zero compared to the neutron density up to baryon
densities close to $n_{0}$, and so the picture in Fig. 1 still holds.

To allow for two flavors of quarks (up and down), we describe the total quark
baryon density as

$n_{\textrm{B}}^{\textrm{Q}}=n^{u}_{\textrm{B}}+n^{d}_{\textrm{B}}\,,$ (15)

and relate the up ($n^{u}_{\textrm{B}}$) and down ($n^{d}_{\textrm{B}}$) quark
baryon densities by enforcing charge neutrality,

$n^{d}_{\textrm{B}}=2n^{u}_{\textrm{B}}\,.$ (16)

Then, the total quark energy density functional is given by

$\epsilon_{\textrm{Q}}=2N_{\textrm{c}}\sum_{q=u,d}\int_{0}^{k_{\text{F}}^{\textrm{q}}}\frac{\textrm{d}^{3}k}{(2\pi)^{3}}\sqrt{M_{\textrm{q}}^{2}+k^{2}}\frac{\sqrt{\Lambda^{2}(n_{\textrm{B}})+k^{2}}}{k}\,.$
(17)

####  II.1.3 Full quarkyonic EOS

In the quarkyonic model, the free quarks and the neutrons have a single Fermi
distribution, where the bottom of the neutron Fermi sea is set by the Fermi
momentum of the d-quarks, $k_{\textrm{F}}^{\textrm{d}}$. That is,
$k_{\textrm{F}}\equiv N_{\textrm{c}}k_{\textrm{F}}^{\textrm{d}}$ as
illustrated in Fig. 1. Combining Eqs. (6) and (17), the total energy density
functional is given by

$\displaystyle\epsilon=$ $\displaystyle
2\Big{[}1-\Big{(}\frac{n_{\textrm{B}}^{\textrm{N}}}{n_{0}}\Big{)}\Big{]}\int_{0}^{\Delta}\frac{\textrm{d}^{3}k}{(2\pi)^{3}}\sqrt{M_{\textrm{N}}^{2}+k^{2}}+\tilde{V}(n_{\textrm{B}}^{\textrm{N}})$
(18)
$\displaystyle+2N_{\textrm{c}}\sum_{q=u,d}\int_{0}^{k_{\text{F}}^{\textrm{q}}}\frac{\textrm{d}^{3}k}{(2\pi)^{3}}\sqrt{M_{\textrm{q}}^{2}+k^{2}}\frac{\sqrt{\Lambda^{2}(n_{\textrm{B}})+k^{2}}}{k}\,,$

where

$\Delta=(3\pi^{2}n_{\textrm{ex}}^{\textrm{N}}+k_{\textrm{F}}^{3})^{1/3}-k_{\textrm{F}}\,.$
(19)

The equilibrium configuration at a given baryon density is then found by
minimizing $\epsilon$ with respect to either $n_{\textrm{B}}^{\textrm{Q}}$ or
$k_{\textrm{F}}$. Once the equilibrium configuration is found, it is
straightforward to obtain the energy density of this equilibrium configuration
from the energy density functional. Other interesting thermodynamic properties
can then be obtained using the thermodynamic relations

$\displaystyle\mu(n_{\textrm{B}})$
$\displaystyle=\frac{\textrm{d}\epsilon(n_{\textrm{B}})}{\textrm{d}n_{\textrm{B}}}\,,\quad\quad\textrm{(chemical
potential)}$ (20) $\displaystyle P(n_{\textrm{B}})$
$\displaystyle=n_{\textrm{B}}\mu(n_{\textrm{B}})-\epsilon(n_{\textrm{B}})\,,\quad\textrm{(pressure)}$
(21) $\displaystyle c_{S}^{2}(n_{\textrm{B}})$
$\displaystyle=\frac{\textrm{d}P(n_{\textrm{B}})}{\textrm{d}\epsilon(n_{\textrm{B}})}\,,\quad\quad\quad\textrm{(speed
of sound)}$ (22)

for the equilibrium configuration.

We would like to point out some of the limitations of the quarkyonic matter
model employed here that have not yet been discussed. First, we have assumed
the baryonic part to be purely made up of neutrons, and have ignored beta
equilibrium because it is not trivial to treat the Fermi sphere for quarkyonic
matter in a dynamic model when it involves more than one species of nucleons.
We note, whoever, that there are non-dynamic quarkyonic matter models
accounting for beta-equilibrium Margueron  _et al._ (2021); Zhao and Lattimer
(2020). It can nevertheless be argued that the fraction of protons to neutrons
is less than $\lesssim 10\%$ McLerran and Reddy (2019) and should, therefore,
not have a significant effect on neutron star masses and radii. Second, we
have only considered two flavors of quarks. In principle, strange quarks could
be included as well and there are three-flavor excluded volume models of
quarkyonic matter in the literature Duarte  _et al._ (2020a, b). However, the
inclusion of strange quarks introduces hyperons, for which the nuclear
interactions are not well understood. This leads to the introduction of
additional free parameters, comlicating the analysis in this paper. Recently,
there have been developments in mean-field approaches to quarkyonic matter
Kumar  _et al._ (2023); Xia  _et al._ (2023) as well. However, we note that
there are no such models that dynamically set the neutron shell width like the
excluded volume model we present here. We, therefore, leave the analysis of
such models for future work. The quarkyonic model described in this section
was shown to be in agreement with neutron star mass-radius data Sen and
Sivertsen (2021) for $\Lambda_{0}$ close to the QCD scale and $n_{0}$ of the
order of the nuclear saturation density. In the next section, we are going to
be agnostic about these and the other model parameters and extensively
investigate what astrophysical data tells us about their values, as well as
the macroscopic significance of the quarkyonic matter model, using a Bayesian
approach.

###  II.2 Bayesian Inference

By using Bayes’ theorem, the posterior $p(\vec{\theta}|d,\mathcal{H})$ on a
set of parameters $\vec{\theta}$ under the hypothesis $\mathcal{H}$ and with
data $d$ is given by

$\displaystyle
p(\vec{\theta}|d,\mathcal{H})=\frac{p(d|\vec{\theta},\mathcal{H})p(\vec{\theta}|\mathcal{H})}{p(d|\mathcal{H})}\to\mathcal{P}(\vec{\theta})=\frac{\mathcal{L}(\vec{\theta})\pi(\vec{\theta})}{\mathcal{Z}}\,,$
(23)

where $\mathcal{P}(\vec{\theta})$, $\mathcal{L}(\vec{\theta})$,
$\pi(\vec{\theta})$, and $\mathcal{Z}$ are the posterior, likelihood, prior,
and evidence, respectively. The prior describes our knowledge of the
parameters before any observations. The likelihood quantifies how likely the
hypothesis can describe the data at a given point $\vec{\theta}$ in the
parameter space. Finally, the evidence, also known as the marginalized
likelihood, marginalizes the likelihood over the whole parameter space with
respect to the prior, i.e,

$\mathcal{Z}=\int_{\mathcal{V}}d\vec{\theta}\mathcal{L}(\vec{\theta})\pi(\vec{\theta})\,.$
(24)

Two hypotheses, $\mathcal{H}_{1}$ and $\mathcal{H}_{2}$, can be compared by
calculating the odds ratio between them, $\mathcal{O}^{1}_{2}$, given by

$\mathcal{O}^{1}_{2}\equiv\frac{p(d|\mathcal{H}_{1})}{p(d|\mathcal{H}_{2})}\frac{p(\mathcal{H}_{1})}{p(\mathcal{H}_{2})}\equiv\mathcal{B}^{1}_{2}\Pi^{1}_{2}\,,$
(25)

where $\mathcal{B}^{1}_{2}$ and $\Pi^{1}_{2}$ are the Bayes factor and prior
odds, respectively. If $\mathcal{O}^{1}_{2}>1$, $\mathcal{H}_{1}$ is more
plausible than $\mathcal{H}_{2}$, and vice versa. Throughout this work, the
prior odds is set to $1$, in which case the Bayes factor is the same as the
odd ratio.

To combine multiple independent observations, we express the likelihood as

$\mathcal{L}(\vec{\theta})=\prod_{i}\mathcal{L}_{i}(\vec{\theta})\,,$ (26)

where $\mathcal{L}_{i}(\vec{\theta})$ is the likelihood corresponding to the
$i$-th observation. Therefore, the combined posterior
$\mathcal{P}(\vec{\theta})$ is given by

$\displaystyle\mathcal{P}(\vec{\theta})$
$\displaystyle=\frac{\pi(\vec{\theta})}{\mathcal{Z}}\mathcal{L}(\vec{\theta})=\frac{\pi(\vec{\theta})}{\mathcal{Z}}\prod_{i}\mathcal{L}_{i}(\vec{\theta})\,.$
(27)

###  II.3 Implementation

The Nuclear-Physics Multi-Messenger Astrophysics framework NMMA Dietrich  _et
al._ (2020); Pang  _et al._ (2022) has been known for its capability and
flexibility of including various multi-messenger astrophysical observations,
nuclear-theory calculations, and heavy-ion collision experiments Huth  _et
al._ (2021). Here, however, we only consider the astrophysical observations to
constrain our theoretical model so that we obtain an astrophysics-only
motivated result.

To calculate the likelihood for different observational data given the
quarkyonic EOS model, a similar approach as in previous works is used Lackey
and Wade (2015); Essick  _et al._ (2020); Wysocki  _et al._ (2020); Pang  _et
al._ (2020). In particular, the likelihood of the EOS parameters
$\vec{E}\equiv\\{\tilde{a},\tilde{b},n_{0},\Lambda_{0}\\}$ given the $i$-th
observation is given by

$\displaystyle\mathcal{L}(\vec{E})=\int
d\vec{\theta}_{\small\textrm{macro}}\frac{\pi(\vec{\theta}_{\small\textrm{macro}}|\vec{E})}{\pi(\vec{\theta}_{\small\textrm{macro}})}\mathcal{P}(\vec{\theta}_{\small\textrm{macro}})\,,$
(28)

where $\vec{\theta}_{\small\textrm{macro}}$,
$\pi(\vec{\theta}_{\small\textrm{macro}}|\vec{E})$, and
$\pi(\vec{\theta}_{\small\textrm{macro}})$ are the macroscopic parameters of
interest, and the priors on $\vec{\theta}_{\small\textrm{macro}}$ with and
without an EOS imposed, respectively. The inference logic is summarized in
Fig. 2. The joint posterior is explored using the Nested Sampling algorithm
Multinest Feroz  _et al._ (2009) implemented in PyMultinest Buchner, J.  _et
al._ (2014); Buchner . The details of the likelihood evaluation for each
observation are described in the following.

Figure 2: The inference logic of our multi-messenger analysis. Each shaded
node corresponds to parameter(s) to be sampled. The unshaded nodes refer to
the latent parameters calculated based on the sample parameters and required
for the likelihood evaluation. An arrow pointing from $A$ to $B$ refers to $B$
depending on $A$. The NICER likelihood $\mathcal{L}_{\rm
NICER}\equiv\mathcal{L}_{\rm J0740+6620}\times\mathcal{L}_{\rm J0030+0451}$,
and the PSR likelihood $\mathcal{L}_{\rm PSR}\equiv\mathcal{L}_{\rm
J0348+4032}\times\mathcal{L}_{\rm J1614-2230}$. In total, eight neutron stars
are taken into account for the analysis.

####  II.3.1 PSR J0348+4032 and PSR J1614-2230

The radio observations on PSR J0348+4042 Antoniadis  _et al._ (2013), and PSR
J1614-2230 Arzoumanian  _et al._ (2018) have provided a lower bound on the
maximum mass of a neutron star. For these two observations, the mass of the
pulsar is the macroscopic parameter of interest. The likelihood is given by

$\mathcal{L}_{{\rm PSR-}j}(\vec{E})=\int^{m_{\rm max}}_{0}dm\mathcal{P}(m|{\rm
PSR-}j)\,,$ (29)

where $\mathcal{P}(m|{\rm PSR-}j)$ is the posterior distribution of the pulsar
$j$’s mass, and $m_{\rm max}$ is the maximum mass supported by the EOS with
parameter $\vec{E}$. We approximate the posterior distribution of the two
pulsars to be Gaussian with the reported values and $1$-$\sigma$ uncertainty,
similar to Ref. Dietrich  _et al._ (2020).

####  II.3.2 PSR J0030+0451 and PSR J0740+6620

The Neutron Star Interior Composition Explorer (NICER) mission aims to measure
both the masses and the radii of pulsars. It has provided estimates for the
mass and radius for the pulsars J0030+0451 Miller  _et al._ (2019); Riley  _et
al._ (2019) and J0740+6620 Miller  _et al._ (2021); Riley  _et al._
(2021)111An update of the NICER analysis on J0740+6620 is presented in Ref.
Salmi  _et al._ (2022). Due to the agreement between the previous results and
the updated one, here we are using the previously original results.. For
J0740+6620, data from the XMM-Newton telescope Struder  _et al._ (2001);
Turner  _et al._ (2001) has been additionally used for improving the total
flux measurement.

The corresponding likelihood for the NICER measurement of pulsar $j$ is given
by

$\displaystyle\mathcal{L}_{{\rm NICER-}j}(\vec{E},m)$ (30) $\displaystyle=\int
d\\!R\ \mathcal{P}_{{\rm NICER}-j}(m,R)\frac{\pi(m,R|\vec{E})}{\pi(m,R|I)}$
$\displaystyle\propto\int d\\!R\ \mathcal{P}_{{\rm
NICER}-j}(m,R)\delta(R-R(m,\vec{E}))$ $\displaystyle\propto\mathcal{P}_{{\rm
NICER}-j}(m,R=R(m,\vec{E}))\,,$

where $\mathcal{P}_{{\rm NICER}-j}(m,R)$ is the joint-posterior distribution
of mass and radius as measured by NICER and we use the fact that i) the radius
is a function of mass for a given EOS, and ii) the prior for the mass and
radius is taken to be uniform in Refs. Miller  _et al._ (2019); Riley  _et
al._ (2019); Miller  _et al._ (2021); Riley  _et al._ (2021).222The mass
measurement for J0740+6620 using radio observations Fonseca  _et al._ (2021)
is taken as the prior on the mass for the NICER analysis on J0740+6620 Miller
_et al._ (2021); Riley  _et al._ (2021). Hence, we consider it to be part of
the NICER likelihood for that star. Therefore, the prior on the mass of
J0740+6620 can be considered uniform. Similar to Ref. Huth  _et al._ (2021),
we use the results of Ref. Miller  _et al._ (2019) for PSR J0030+0451, while
for PSR J0740+6620 we average over the results presented in Refs. Miller  _et
al._ (2021) and Riley  _et al._ (2021).

####  II.3.3 GW170817 and GW190425

By analyzing the gravitational-wave signals GW170817 Abbott  _et al._ (2017)
and GW190425 Abbott  _et al._ (2020), the masses $m_{i}$ and the tidal
deformability $\Lambda_{i}$ of the two neutron stars in the binary can be
estimated333Note that it was also proposed that GW190425 originated from a
neutron-star-black-hole merger Foley  _et al._ (2020); Han  _et al._ (2020);
Kyutoku  _et al._ (2020).. The corresponding likelihood is given by

$\displaystyle\mathcal{L}_{{\rm GW}-j}(\vec{E},m_{i})$ $\displaystyle=\int
d\Lambda_{i}\frac{\pi(m_{i},\Lambda_{i}|\vec{E})}{\pi(m_{i},\Lambda_{i}|I)}\mathcal{P}_{{\rm
GW}-j}(m_{i},\Lambda_{i})$ (31) $\displaystyle=\int
d\Lambda_{i}\frac{\prod_{i}\delta(\Lambda_{i}-\Lambda(m_{i};\vec{E}))}{\pi(\Lambda_{i}|m_{i},I)}$
$\displaystyle\times\frac{\pi(m_{i}|\vec{E})}{\pi(m_{i}|I)}\mathcal{P}_{{\rm
GW}-j}(m_{i},\Lambda_{i})$ $\displaystyle\propto\left.\frac{\mathcal{P}_{{\rm
GW}-j}(m_{i},\Lambda_{i})}{\pi(\Lambda_{i}|,m_{i},I)}\right|_{\Lambda_{i}=\Lambda(m_{i};\vec{E})}\,,$

where $\mathcal{P}_{{\rm GW}-j}(m_{i},\Lambda_{i})$ is the joint-posterior
distribution on mass and tidal deformability of the binary neutron star $j$
measured by its gravitational-wave signal. We use that the tidal deformability
is a function of mass for a given EOS. For both events, we use the publicly
available posterior samples Collaboration and the Virgo Collaboration (2019,
2020).

####  II.3.4 AT2017gfo

For analyzing the observed kilonova AT2017gfo, we have done the Bayesian
inference with a Gaussian likelihood function given by

$\mathcal{L}_{\rm
EM}(\vec{\theta})\propto\exp\left(-\frac{1}{2}\sum_{ij}\left(\frac{m^{j}_{i}-m^{j,\rm{est}}_{i}(\vec{\theta})}{\sigma^{j}_{i}}\right)^{2}\right)\,,$
(32)

where $m^{j}_{i}$ and $\sigma^{j}_{i}$ are the observed apparent magnitude and
its corresponding statistical uncertainties, at observation time $t_{i}$,
respectively. Moreover, $m^{j,\rm{est}}_{i}$ are the theoretically predicted
apparent magnitudes for a given filter $j$ (a passband for a particular
wavelength interval).

For this analysis, the model presented in Ref. Bulla (2019) is used. For this
model, the macroscopic parameters of interest are the dynamic ejecta mass
$m^{\rm dyn}_{\rm ej}$ and the wind ejecta mass $m^{\rm wind}_{\rm ej}$. The
dynamic ejecta refer to the material ejected during the merger via torque and
shocks while the wind ejecta refer to material released from the disk formed
during the merger.

In order to connect the ejecta masses with the EOS and the masses of the two
neutron stars in the binary, fits to numerical-relativity simulations are
used. For the dynamical ejecta mass $m^{\rm dyn}_{\rm ej}$, the fit formula is
given by Krüger and Foucart (2020)

$\frac{m^{\rm dyn}_{\rm
ej,fit}}{10^{-3}M_{\odot}}=\left(\frac{a}{C_{1}}+b\left(\frac{m_{2}}{m_{1}}\right)^{n}+c\
C_{1}\right)+(1\leftrightarrow 2)\,,$ (33)

where $m_{i}$ and $C_{i}$ are the masses and the compactness of the two
components of the binary and the best-fit coefficients are given by
$a=-9.3335$, $b=114.17$, $c=-337.56$, and $n=1.5465$. As presented in Ref.
Krüger and Foucart (2020), this relation provides an accurate estimation of
the ejecta mass with an error that is well-approximated by a zero-mean
Gaussian with a standard deviation of $0.004M_{\odot}$. Therefore, the dynamic
ejecta mass is described by

$m^{\rm dyn}_{\rm ej}=m^{\rm dyn}_{\rm ej,fit}+\alpha\,,$ (34)

where $\alpha\sim\mathcal{N}(\mu=0,\sigma=0.004M_{\odot})$.

For the wind ejecta mass, we assume it to be proportional to the disk mass,

$m^{\rm wind}_{\rm ej}=\zeta\times m_{\rm disk,fit}\,,$ (35)

where $\zeta$ is an independent parameter in $[0,1)$. To estimate the disk
mass $m_{\rm disk}$, we follow the relation presented in Ref. Dietrich  _et
al._ (2020),

$\displaystyle\log_{10}\left(\frac{m_{\rm disk,fit}}{M_{\odot}}\right)=$ (36)
$\displaystyle\textrm{max}\left(-3,a\left(1+b\tanh\left(\frac{c-(m_{1}+m_{2})M_{\rm
thr}^{-1}}{d}\right)\right)\right)\,,$

with $a$ and $b$ given by

$\displaystyle a$ $\displaystyle=a_{o}+\delta a\cdot\xi\,,$ (37)
$\displaystyle b$ $\displaystyle=b_{o}+\delta b\cdot\xi\,.$

The parameter $\xi$ is given by

$\xi=\frac{1}{2}\tanh\left(\beta\left(q-q_{\rm trans}\right)\right)\,,$ (38)

where $q\equiv m_{2}/m_{1}\leq 1$ is the mass ratio and $\beta$ and $q_{\rm
trans}$ are free parameters. The best-fitting parameters are $a_{o}=-1.725$,
$\delta a=-2.337$, $b_{o}=-0.564$, $\delta b=-0.437$, $c=0.958$, $d=0.057$,
$\beta=5.879$, $q_{\rm trans}=0.886$. The threshold mass $M_{\rm thr}$ is
estimated to follow the prediction presented in Ref. Bauswein  _et al._
(2020). Because a fitting error in the disk mass is degenerate with the
proportionality parameters $\zeta$, no additional error parameters, similar to
$\alpha$ introduced in Eq. (34), are included. We point out that the relations
Eqs. (34) and (36) are generally calibrated to a large set of numerical-
relativity simulations covering a large range of the parameter space, and
while these relations are found to be ‘quasi-universal’, i.e., applicable
independent of the exact EOS, they have not been tested against quarkyonic
binary neutron stars mergers, which could result in a hidden systematic
uncertainty in our analysis.

Parameter | Prior
---|---
$\tilde{a}$ | Uniform$(-50{\rm MeV},50{\rm MeV})$
$\tilde{b}$ | Uniform$(-50{\rm MeV},50{\rm MeV})$
$n_{0}$ | LogUniform$(10^{-2}n_{\rm sat},20n_{\rm sat})$
$\Lambda_{0}$ | LogUniform$(10^{-5}{\rm GeV},10{\rm GeV})$
Table 1: Priors imposed for the EOS parameters in the analysis.
“Uniform$(a,b)$” is a uniform distribution in the range of $[a,b)$ and
“LogUniform$(a,b)$” is a log-uniform distribution in the range of $[a,b)$.
Therefore, if $x\sim{\rm LogUniform}(a,b)$, $\log x\sim{\rm Uniform}(\log
a,\log b)$.

The likelihood for the EOS parameters $\vec{E}$ and the component masses of
the two neutron stars, $m_{i}$ is given by

$\displaystyle\mathcal{L}_{\rm AT2017gfo}(\vec{E},m_{i},\alpha,\zeta)$ (39)
$\displaystyle=\int dm^{\rm dyn}_{\rm ej}dm^{\rm wind}_{\rm
ej}\frac{\pi(m^{\rm dyn}_{\rm ej},m^{\rm wind}_{\rm
ej}|\vec{E},m_{i},\alpha,\zeta)}{\pi(m^{\rm dyn}_{\rm ej},m^{\rm wind}_{\rm
ej}|I)}$ $\displaystyle\times\mathcal{P}_{\rm AT2017gfo}(m^{\rm dyn}_{\rm
ej},m^{\rm wind}_{\rm ej})$ $\displaystyle\propto\mathcal{P}_{\rm
AT2017gfo}(m^{\rm dyn}_{\rm ej},m^{\rm wind}_{\rm ej})|_{m^{\rm dyn}_{\rm
ej}=m^{\rm dyn}_{\rm ej}(\vec{E},m_{i},\alpha),m^{\rm wind}_{\rm
ej}(\vec{E},m_{i},\zeta)}\,.$

##  III Results

Figure 3: Posterior on the microscopic parameters of the quarkyonic EOS model.
The 2D contours are shown at $68\%$ and $95\%$ levels. The priors on each of
the parameters are also shown (grey). In addition to the distributions shown
in Tab. 1, the priors shown additionally impose that the resulting EOS can
support a neutron star.

In the following section, we discuss the results of our analysis of
astrophysical multi-messenger observations when applying the quarkyonic EOS
model. As mentioned in Sec. II.3.4, the ‘quasi-universal’ relations needed for
analyzing AT2017gfo have not been tested against the quarkyonic model.
Therefore, we report results without the inclusion of AT2017gfo as our main,
conservative results and results with the inclusion of AT2017gfo in
parentheses. All quoted values are medians with their $95\%$ credible
intervals unless mentioned otherwise.

###  III.1 Estimation on EOS parameters

The priors for the EOS parameters $\tilde{a}$, $\tilde{b}$, $n_{0}$, and
$\Lambda_{0}$ are defined Tab. 1. Additionally, we impose that

  * •

each parameter sample results in a valid EOS, i.e.,

    * –

the crust and core EOSs intercept, and

    * –

the EOS explores densities beyond the neutron-star crust,

  * •

and the EOS produces stable neutron stars, indicated by a stable
$M$-$R$-$\Lambda$ curve.

In Fig. 3, we show the priors and the posteriors with their median values and
with the uncertainty quoted at $95\%$ credibility.

Figure 4: Posterior on the neutron star radius as a function of mass. The
darker/solid (lighter/dashed) band shows the 68% (95%) credible interval of
the radius at a given mass. The prior band is shown in the 95% band.

The posterior on $\tilde{a}$ skews towards positive values, while the
posterior on $\tilde{b}$ is centered around $0$, with both of them having
sizeable uncertainty. This shows that interactions with higher-order density
dependences between neutrons are not pronounced within this model. This
finding can be attributed to the relatively low $n_{0}$, which indicates the
characteristic density at which quark matter starts to appear. This low
$n_{0}$ causes the quark part of the EOS to dominate over the neutrons’
interactions. However, we note that there are large uncertainties in the
neutron matter parameters $\tilde{a}$ and $\tilde{b}$, mainly due to the fact
that we have only used astrophysical observations to constrain them.

More interestingly, we find that the parameter $n_{0}$ is favored to be on the
order of a few times the nuclear saturation density, which is comparable to
the densities reached inside the heaviest neutron stars. This implies that,
assuming our quarkyonic matter model describes supranuclear-dense matter,
there will be a significant amount of quarks in the core of neutron stars with
masses in the two-solar-mass range. We will investigate this in more detail in
the next section.

Finally, we comment on the parameter $\Lambda_{0}$. As mentioned in section
II.1.2, the inclusion of $\Lambda_{0}$ was motivated by regulating the
infrared divergence that occurs at low densities due to the non-perturbative
effects of QCD. We, therefore, expect that
$\Lambda_{0}\sim\Lambda_{\text{QCD}}$. If we take $\Lambda_{\text{QCD}}$ to be
defined by the Landau pole of the QCD strong coupling, one finds
$\Lambda_{\text{QCD}}\sim 200-400$MeV, depending on the renormalization scheme
used Kataev and Molokoedov (2015). This is indeed consistent with what we
obtain for $\Lambda_{0}$, as seen in Fig. 3. We emphasize that although this
is an expected result, the value of $\Lambda_{0}$ is extracted completely by
comparison with astrophysical data. It is interesting and reassuring that
observations of neutron stars can constrain this microscopic parameter of our
model in a way that is consistent with expectations from QCD. This is because
the confinement scale is important for the high-density behavior of the EOS,
which, in turn, has a high impact on the mass of a neutron star. Conversely,
this explains, in part, why we are not able to constrain $\tilde{a}$ and
$\tilde{b}$ very well, as they mostly describe the low-density behavior of the
EOS. These findings are consistent with similar results showing that the
correlation between the low- and high-density EOS can be broken in neutron
stars Essick  _et al._ (2021a, b).

The posterior of the radius of neutron stars as a function of mass is shown in
Fig. 4. Within this quarkyonic-matter model, the radius of a $1.4M_{\odot}$
neutron star, $R_{1.4}$, is estimated to be
$13.44^{+1.69}_{-1.54}(13.54^{+1.02}_{-1.04})$ km, at 95% credibility, without
(with) the inclusion of AT2017gfo; the kilonova observation significantly
tightens the radius constraint. The extracted value for $R_{1.4}$ is higher
than in previous estimations (e.g. Tab. 1 in Ref. Pang  _et al._ (2022)),
reflecting the difference of model-independent and model-dependent data
analyses of neutron-star observations, as well as the absence of the
constraint from chiral effective field theory at low densities Tews  _et al._
(2018b). The latter also contributes to the wide posteriors on $\tilde{a}$ and
$\tilde{b}$.

###  III.2 Presence of quark matter

Figure 5: Posterior on the maximum quark density observed within a neutron
star. The posterior and prior are re-weighted such that the prior on the
maximum quark density is log-uniform, i.e., uniform in $\log n^{\rm Q}_{\rm
B}$. Figure 6: Posterior on the maximum quark-to-baryon ratio observed within
a neutron star. The posterior and prior are re-weighted such that the prior on
the ratio is log-uniform, i.e., uniform in $\log n^{\rm Q}_{\rm B}/n_{\rm B}$.
Figure 7: Posterior on the ($\log_{10}$) maximum quark baryon density
$\max\\{n^{\rm Q}_{\rm B}\\}$, the ($\log_{10}$) quark core mass for maximal
mass neutron star $M^{\rm Q}_{\rm TOV}$ and the associated radius $R_{\rm
TOV}$. A strong positive correlation is observed between the quark core mass
and the radius, and between maximum quark baryon density and quark core mass.

To study the presence of quarks in the core of neutron stars, we estimate the
maximum quark baryon density $\max\\{n^{\rm Q}_{\rm B}\\}$ and the maximum
quark-to-baryon ratio $\max\\{n^{\rm Q}_{\rm B}/n_{\rm B}\\}$. That is, we
consider the density of quarks and quark-to-baryon ratio inside the center of
the neutron star with largest mass. The estimated distributions of
$\max\\{n^{\rm Q}_{\rm B}\\}$ and $\max\\{n^{\rm Q}_{\rm B}/n_{\rm B}\\}$ are
shown in Fig. 5 and Fig. 6, respectively. In both cases, we re-weighted the
distributions to respect a log-uniform prior to enhance the impact of
astrophysical observations. Based on the estimated distributions, the
quarkyonic model is predicting $\max\\{n^{\rm Q}_{\rm B}\\}\gtrsim 1n_{\rm
sat}$ and the ratio $\max\\{n^{\rm Q}_{\rm B}/n_{\rm B}\\}\gtrsim 10\%$,
showing a significant presence of quarks within the most massive neutron
stars.

| Without AT2017gfo | With AT2017gfo
---|---|---
| $m_{q}[M_{\odot}]$ | $m_{q}/m[\%]$ | $m_{q}[M_{\odot}]$ | $m_{q}/m[\%]$
GW170817-$m_{1}$ | $0.01^{+0.02}_{-0.01}$ | $0.41^{+1.45}_{-0.37}$ | $0.00^{+0.02}_{-0.00}$ | $0.34^{+1.36}_{-0.31}$
GW170817-$m_{2}$ | $0.00^{+0.01}_{-0.00}$ | $0.23^{+0.91}_{-0.20}$ | $0.00^{+0.01}_{-0.00}$ | $0.25^{+1.07}_{-0.22}$
GW190425-$m_{1}$ | $0.02^{+0.05}_{-0.02}$ | $0.98^{+2.88}_{-0.94}$ | $0.02^{+0.06}_{-0.02}$ | $0.98^{+3.16}_{-0.94}$
GW190425-$m_{2}$ | $0.01^{+0.03}_{-0.01}$ | $0.52^{+1.79}_{-0.48}$ | $0.01^{+0.03}_{-0.01}$ | $0.51^{+1.96}_{-0.47}$
PSR J0030+0451 | $0.01^{+0.03}_{-0.01}$ | $0.40^{+1.65}_{-0.36}$ | $0.01^{+0.03}_{-0.01}$ | $0.40^{+1.74}_{-0.37}$
PSR J0740+6620 | $0.05^{+0.11}_{-0.05}$ | $2.57^{+5.41}_{-2.48}$ | $0.05^{+0.12}_{-0.05}$ | $2.57^{+5.82}_{-2.47}$
Maximal mass | $0.13^{+0.14}_{-0.10}$ | $5.87^{+6.72}_{-4.72}$ | $0.13^{+0.14}_{-0.11}$ | $5.86^{+7.08}_{-4.87}$
Table 2: Estimated quark core mass $m_{q}$ within neutron stars and the
associated mass fraction $m_{q}/m$ for the stars.

To give a concrete picture of the abundance of quarks in neutron stars, we
further estimate the total mass that the free quarks contribute to neutron
stars, see Tab. 2. As expected, the quark core mass is low for most of the
neutron stars. However, PSR J0740+6620, which is one of the most massive
neutron stars observed, would carry about $\sim 0.05M_{\odot}$ of free quarks,
or $\sim 2.6\%$ of its mass. For a neutron star at maximum mass, the quark
core can grow to $\sim 0.13M_{\odot}$, or $\sim 5.9\%$ of the stellar mass.

Moreover, a strong presence of free quarks shows a correlation with
macroscopic parameters. In Fig. 7, we show a corner plot connecting the
maximum quark baryon density $\max\\{n^{\rm Q}_{\rm B}\\}$, the quark mass in
the maximal-mass neutron star $M^{\rm Q}_{\rm TOV}$ and the associated radius
$R_{\rm TOV}$. We observe a strong positive correlation between the quark core
mass and the radius for the maximal mass neutron star within the quarkyonic-
matter model. Although it may seem to be counter-intuitive that stars
containing quarks have large radii, this positive correlation can be
attributed to the high contribution of the nucleonic matter to the radius of
the neutron star.

###  III.3 Speed of sound and normalized trace anomaly

Figure 8: Posterior on the speed of sound $c^{2}_{s}$ as a function of
density. The darker/solid (lighter/dashed) band show the $68\%$ ($95\%$)
credible interval of the speed of sound at a given density. The prior band is
shown at $95\%$ confidence level. Moreover, we show the range predicted from
the chiral effective field theory calculation of Ref. Tews  _et al._ (2018b)
up to $1.5n_{\rm sat}$. Figure 9: Posterior on the maximum speed of sound
$\max\\{c^{2}_{s}\\}$ and the associated number density $n_{\rm B}|_{{\rm
max}\ c_{s}}$ and the quark-to-baryon ratio $n^{\rm Q}_{\rm B}/n_{\rm
B}|_{{\rm max}\ c_{s}}$. A strong negative correlation is observed between the
maximum speed of sound and the ratio $n^{\rm Q}_{\rm B}/n_{\rm B}|_{{\rm max}\
c_{s}}$.

We further examine the speed of sound within neutron stars and show it as a
function of density in Fig. 8. It is obvious that the quarkyonic-matter model
violates the conformal limit of $c^{2}_{s}=1/3$ at about 3 times the nuclear
saturation density, but the speed of sound drops below that limit at higher
densities. For comparison, the range predicted by chiral effective field
theory is also shown in Fig. 8, and is on the lower end of the quarkyonic-
matter prediction. The maximum speed-of-sound, the densities at which the
maximum occurs, and the quark-to-baryon ratio at this density are shown in the
corner plot in Fig. 9. The estimate for the maximum speed of sound and its
density agree with the result reported in Ref. Legred  _et al._ (2021), in
which a nonparametric EOS model is used. The maximum speed of sound and the
quark-to-baryon ratio show a strong negative correlation, reflecting the
intuition that quark matter is generally softer than nucleonic matter.

Figure 10: Posterior on the normalized trace anomaly $\Delta_{\rm ta}$ as a
function of density. The darker/solid (lighter/dashed) band show the $68\%$
($95\%$) credible interval of the speed of sound at a given density. The prior
band is shown in the $95\%$ band. Moreover, the predicted range based on the
chiral effective field theory calculation of Ref. Tews  _et al._ (2018b) is
shown up to $1.5n_{\rm sat}$.

Recently, the authors of Ref. Fujimoto  _et al._ (2022) suggested that a peak
in the speed of sound, as in Fig. 8, is a sign of strongly-coupled conformal
matter. This was demonstrated by examining the behaviour of the normalized
trace anomaly, which is defined as

$\Delta_{\rm ta}\equiv\frac{1}{3}-\frac{P}{\epsilon}\,.$ (40)

From the definition it is evident that, as matter approaches conformality,
$\Delta_{\rm ta}\to 0$. Even though the speed of sound and the trace anomaly
are closely related, the behavior at intermediate densities is different and
provides complementary information. By expressing the sound speed in terms of
the trace anomaly and its derivative, Ref. Fujimoto  _et al._ (2022) showed
that the peak in the sound speed is, in fact, not a violation of the conformal
bound, but it is a steep approach to the conformal limit. In this work, we
check this interpretation by computing the trace anomaly as a function of the
density, see Fig. 10, and comparing it with the behavior of the sound speed.
In Fig. 10, we show that at $n\sim 6n_{\rm sat}$, $\Delta_{\rm ta}\to 0$, as
can be inferred from the behavior of the speed of sound in Fig. 8. However, at
intermediate densities, $\Delta_{\rm ta}$ decreases towards its conformal
value, exactly at those densities where we find a peak in the sound speed.
These results are similar to what was observed in Ref. Fujimoto  _et al._
(2022). Note, however, that the quarkyonic-matter model employed in this work
does not explicitly include any interactions between quarks and therefore it
cannot be interpreted as strongly coupled matter. This suggests that inferring
the presence of strongly coupled conformal matter from neutron-star
observations is non-trivial since the corresponding characteristic behaviour
of the trace anomaly can also be produced by weakly interacting matter.

Given the presence of quarks within the neutron stars within the quarkyonic-
matter model, we further investigate the consistency between the model and the
perturbative QCD (pQCD) calculations of the EOS at $\sim 40n_{\rm sat}$ Gorda
_et al._ (2021). Following the approach of Ref. Somasundaram  _et al._ (2023),
we connect to the integral pQCD constraints of Ref. Komoltsev and Kurkela
(2022) at the maximum-mass configuration because we are interested only in the
EOS of stable neutron stars.444See Ref. Gorda  _et al._ (2023b) for a similar
analysis but with a matching density of $10n_{\rm sat}$. We use the parameters
associated with the maximally massive neutron star, i.e., the associated
central pressure $p_{\rm TOV}$, the baryon number density $n_{\rm TOV}$ and
the chemical potential $\mu_{\rm TOV}$ to check if a given EOS is consistent
with the high-density pQCD constraint by verifying

  1. 1.

$\Delta p_{\rm min}<p_{\rm pQCD}-p_{\rm TOV}<\Delta p_{\rm max}$, and,

  2. 2.

$n_{\rm TOV}<n_{\rm pQCD}$,

where $p_{\rm pQCD}$ and $n_{\rm pQCD}$ are functions of $\mu_{\rm pQCD}$,
which is set to $2.6{\rm GeV}$. Furthermore, $\Delta p_{\rm min}$ and $\Delta
p_{\rm max}$ are calculated using Eqs. (12) and (13) of Ref. Gorda  _et al._
(2023b) respectively. If the proposed EOS satisfies both of the above
criteria, it has a consistency of $100\%$ and $0\%$ otherwise. Since the pQCD
calculation Gorda  _et al._ (2021) depends on the renormalization scale
parameter $X$, we calculate the consistency of a given EOS with pQCD by
marginalizing over $X\sim{\rm LogUniform}(1,4)$. Given the posterior samples,
without (with) the inclusion of AT2017gfo, the average consistency, i.e., the
average of the aforementioned calculation over all EOSs, is $99.88\%$
($99.89\%$) with a minimum consistency of $93.9\%$ ($94.7\%$). Therefore, we
conclude that the quarkyonic-matter model is consistent with pQCD
calculations.

##  IV Conclusions

In this work, we have studied a dynamic quarkyonic-matter model for
supranuclear matter inside neutron stars. This model is employed within the
nuclear multi-messenger astrophysics framework, which includes radio and X-ray
(NICER) observations of pulsars, gravitational-wave observations from binary
neutron star mergers, in particular GW170817, and electromagnetic observations
of the kilonova associated with GW170817. Assuming that the quarkyonic-matter
model describes neutron stars, we constrain the properties of quarkyonic
matter using astrophysical observations. The main results can be summarized as
follows:

  * •

We constrained the four model parameters without (with) AT2017gfo at 95%
credibility to be

$\displaystyle\tilde{a}$
$\displaystyle=2.62^{+18.83}_{-12.00}(3.52^{+20.27}_{-12.70})\text{MeV},$
$\displaystyle\tilde{b}$
$\displaystyle=-1.10^{+8.20}_{-7.53}(-1.30^{+8.77}_{-8.12})\text{MeV},$
$\displaystyle n_{0}$
$\displaystyle=3.86^{+2.08}_{-1.38}(3.85^{1.98}_{-1.36})n_{\text{sat}},$
$\displaystyle\Lambda_{0}$
$\displaystyle=306.51^{+120.87}_{-92.72}(305.76^{+108.07}_{-105.03})\text{MeV}.$
(41)

  * •

The model is predicting a substantial presence of quarks within neutron stars,
in particular:

    * –

The maximum quark baryon density exceeds $1n_{\rm sat}$,

    * –

The maximum quark-to-baryon ratio exceeds $10\%$ and

    * –

The mass contribution due to quarks can reach $\sim 0.13M_{\odot}$,
attributing to $\sim 6\%$ of the neutron star’s mass.

  * •

The quarkyonic-matter model predicts a peak in the speed-of-sound in neutron
stars, and the conformal limit is subsequently restored at $\sim 6n_{\rm
sat}$.

  * •

The quarkyonic-matter model predicts
$R_{1.4}=13.44^{+1.69}_{-1.54}(13.54^{+1.02}_{-1.04})$ km at 95% credibility
level, without (with) the inclusion of AT2017gfo.

Our work demonstrates that the quarkyonic-matter model can explain the current
observations on neutron stars and allows for a significant presence of quarks
within them. It is, therefore, a viable alternative to other models, including
free quarks. Although it is impossible to establish the existence of
quarkyonic matter with current neutron-star observations, future observations
using next-generation gravitational-wave observatories, such as the Einstein
telescope Punturo  _et al._ (2010); Hild  _et al._ (2011) or Cosmic Explorer
Reitze  _et al._ (2019), will likely provide crucial new data. Together with
advancements in theoretical nuclear physics, this will allow us to better
understand the nature of dense matter in the future.

###### Acknowledgements.

P.T.H.P. and C.V.D.B. are supported by the research programme of the
Netherlands Organisation for Scientific Research (NWO). L.S. and S.S. are
supported by the U.S. Department of Energy, Nuclear Physics Quantum Horizons
program through the Early Career Award DE-SC0021892. R.S. acknowledges support
from the Nuclear Physics from Multi-Messenger Mergers (NP3M) Focused Research
Hub which is funded by the National Science Foundation under Grant Number
21-16686, and by the Laboratory Directed Research and Development program of
Los Alamos National Laboratory under project number 20220541ECR. T.D.
acknowledges funding from the EU Horizon under ERC Starting Grant, no.
SMArt-101076369, support from the Deutsche Forschungsgemeinschaft, DFG,
project number DI 2553/7, and from the Daimler and Benz Foundation for the
project ‘NUMANJI’. The work of I.T. was supported by the U.S. Department of
Energy, Office of Science, Office of Nuclear Physics, under contract No. DE-
AC52-06NA25396, by the Laboratory Directed Research and Development program of
Los Alamos National Laboratory under project numbers 20220541ECR and
20230315ER, and by the U.S. Department of Energy, Office of Science, Office of
Advanced Scientific Computing Research, Scientific Discovery through Advanced
Computing (SciDAC) NUCLEI program. This material is based upon work supported
by NSF’s LIGO Laboratory which is a major facility fully funded by the
National Science Foundation.

## References

  * Ozel and Freire (2016) F. Ozel and P. Freire, Ann. Rev. Astron. Astrophys. 54, 401 (2016), arXiv:1603.02698 [astro-ph.HE] .
  * Flanagan and Hinderer (2008) E. E. Flanagan and T. Hinderer, Phys.Rev. D77, 021502 (2008), arXiv:0709.1915 [astro-ph] .
  * Del Pozzo  _et al._ (2013) W. Del Pozzo, T. G. F. Li, M. Agathos, C. Van Den Broeck, and S. Vitale, Phys. Rev. Lett. 111, 071101 (2013), arXiv:1307.8338 [gr-qc] .
  * Annala  _et al._ (2018) E. Annala, T. Gorda, A. Kurkela, and A. Vuorinen, Phys. Rev. Lett. 120, 172703 (2018), arXiv:1711.02644 [astro-ph.HE] .
  * Abbott  _et al._ (2018) B. P. Abbott _et al._ (Virgo, LIGO Scientific), Phys. Rev. Lett. 121, 161101 (2018), arXiv:1805.11581 [gr-qc] .
  * De  _et al._ (2018) S. De, D. Finstad, J. M. Lattimer, D. A. Brown, E. Berger, and C. M. Biwer, Phys. Rev. Lett. 121, 091102 (2018), arXiv:1804.08583 [astro-ph.HE] .
  * Tews  _et al._ (2018a) I. Tews, J. Margueron, and S. Reddy, Phys. Rev. C 98, 045804 (2018a), arXiv:1804.02783 [nucl-th] .
  * Watts  _et al._ (2016) A. L. Watts _et al._ , Rev. Mod. Phys. 88, 021001 (2016), arXiv:1602.01081 [astro-ph.HE] .
  * Demorest  _et al._ (2010) P. Demorest, T. Pennucci, S. Ransom, M. Roberts, and J. Hessels, Nature 467, 1081 (2010), arXiv:1010.5788 [astro-ph.HE] .
  * Antoniadis  _et al._ (2013) J. Antoniadis, P. C. Freire, N. Wex, T. M. Tauris, R. S. Lynch, _et al._ , Science 340, 6131 (2013), arXiv:1304.6875 [astro-ph.HE] .
  * Drischler  _et al._ (2021) C. Drischler, S. Han, J. M. Lattimer, M. Prakash, S. Reddy, and T. Zhao, Phys. Rev. C 103, 045808 (2021).
  * Raaijmakers  _et al._ (2020) G. Raaijmakers, S. K. Greif, T. E. Riley, T. Hinderer, K. Hebeler, A. Schwenk, A. L. Watts, S. Nissanke, S. Guillot, J. M. Lattimer, and R. M. Ludlam, The Astrophysical Journal 893, L21 (2020).
  * Hebeler  _et al._ (2013) K. Hebeler, J. M. Lattimer, C. J. Pethick, and A. Schwenk, Astrophys. J. 773, 11 (2013), arXiv:1303.4662 [astro-ph.SR] .
  * Tews  _et al._ (2018b) I. Tews, J. Carlson, S. Gandolfi, and S. Reddy, Astrophys. J. 860, 149 (2018b), arXiv:1801.01923 [nucl-th] .
  * Drischler  _et al._ (2020) C. Drischler, R. J. Furnstahl, J. A. Melendez, and D. R. Phillips, Phys. Rev. Lett. 125, 202702 (2020), arXiv:2004.07232 [nucl-th] .
  * Keller  _et al._ (2023) J. Keller, K. Hebeler, and A. Schwenk, Phys. Rev. Lett. 130, 072701 (2023), arXiv:2204.14016 [nucl-th] .
  * Gorda  _et al._ (2018) T. Gorda, A. Kurkela, P. Romatschke, M. Säppi, and A. Vuorinen, Phys. Rev. Lett. 121, 202701 (2018), arXiv:1807.04120 [hep-ph] .
  * Gorda  _et al._ (2021) T. Gorda, A. Kurkela, R. Paatelainen, S. Säppi, and A. Vuorinen, Phys. Rev. Lett. 127, 162003 (2021), arXiv:2103.05658 [hep-ph] .
  * Gorda  _et al._ (2023a) T. Gorda, R. Paatelainen, S. Säppi, and K. Seppänen, (2023a), arXiv:2307.08734 [hep-ph] .
  * Bedaque and Steiner (2015) P. Bedaque and A. W. Steiner, Phys. Rev. Lett. 114, 031103 (2015), arXiv:1408.5116 [nucl-th] .
  * McLerran and Reddy (2019) L. McLerran and S. Reddy, Phys. Rev. Lett. 122, 122701 (2019), arXiv:1811.12503 [nucl-th] .
  * Jeong  _et al._ (2020) K. S. Jeong, L. McLerran, and S. Sen, Phys. Rev. C 101, 035201 (2020), arXiv:1908.04799 [nucl-th] .
  * Sen and Sivertsen (2021) S. Sen and L. Sivertsen, Astrophys. J. 915, 109 (2021), arXiv:2011.04681 [astro-ph.HE] .
  * Abbott  _et al._ (2017) B. P. Abbott _et al._ (Virgo, LIGO Scientific), Phys. Rev. Lett. 119, 161101 (2017), arXiv:1710.05832 [gr-qc] .
  * Abbott  _et al._ (2020) B. Abbott _et al._ (LIGO Scientific, Virgo), Astrophys. J. Lett. 892, L3 (2020), arXiv:2001.01761 [astro-ph.HE] .
  * Aasi  _et al._ (2015) J. Aasi _et al._ (LIGO Scientific), Class. Quant. Grav. 32, 074001 (2015), arXiv:1411.4547 [gr-qc] .
  * Acernese  _et al._ (2015) F. Acernese _et al._ (VIRGO), Class. Quant. Grav. 32, 024001 (2015), arXiv:1408.3978 [gr-qc] .
  * GBM (2017) Astrophys. J. 848, L12 (2017), arXiv:1710.05833 [astro-ph.HE] .
  * Miller  _et al._ (2019) M. C. Miller _et al._ , Astrophys. J. Lett. 887, L24 (2019), arXiv:1912.05705 [astro-ph.HE] .
  * Riley  _et al._ (2019) T. E. Riley _et al._ , Astrophys. J. Lett. 887, L21 (2019), arXiv:1912.05702 [astro-ph.HE] .
  * Miller  _et al._ (2021) M. C. Miller _et al._ , Astrophys. J. Lett. 918, L28 (2021), arXiv:2105.06979 [astro-ph.HE] .
  * Riley  _et al._ (2021) T. E. Riley _et al._ , Astrophys. J. Lett. 918, L27 (2021), arXiv:2105.06980 [astro-ph.HE] .
  * Arzoumanian  _et al._ (2018) Z. Arzoumanian _et al._ (NANOGrav), Astrophys. J. Suppl. 235, 37 (2018), arXiv:1801.01837 [astro-ph.HE] .
  * Margueron  _et al._ (2021) J. Margueron, H. Hansen, P. Proust, and G. Chanfray, Phys. Rev. C 104, 055803 (2021).
  * Zhao and Lattimer (2020) T. Zhao and J. M. Lattimer, Phys. Rev. D 102, 023021 (2020).
  * Duarte  _et al._ (2020a) D. C. Duarte, S. Hernandez-Ortiz, and K. S. Jeong, Phys. Rev. C 102, 025203 (2020a), arXiv:2003.02362 [nucl-th] .
  * Duarte  _et al._ (2020b) D. C. Duarte, S. Hernandez-Ortiz, and K. S. Jeong, Physical Review C 102 (2020b), 10.1103/physrevc.102.065202.
  * Kumar  _et al._ (2023) A. Kumar, D. Dey, S. Haque, R. Mallick, and S. K. Patra, “Quarkyonic model for neutron star matter: A relativistic mean-field approach,” (2023), arXiv:2304.08223 [nucl-th] .
  * Xia  _et al._ (2023) C.-J. Xia, H.-M. Jin, and T.-T. Sun, “Quarkyonic matter and quarkyonic stars in an extended rmf model,” (2023), arXiv:2307.03032 [hep-ph] .
  * Dietrich  _et al._ (2020) T. Dietrich, M. W. Coughlin, P. T. H. Pang, M. Bulla, J. Heinzel, L. Issa, I. Tews, and S. Antier, Science 370, 1450 (2020), arXiv:2002.11355 [astro-ph.HE] .
  * Pang  _et al._ (2022) P. T. H. Pang _et al._ , (2022), arXiv:2205.08513 [astro-ph.HE] .
  * Huth  _et al._ (2021) S. Huth _et al._ , arXiv:2107.06229 (2021).
  * Lackey and Wade (2015) B. D. Lackey and L. Wade, Phys.Rev. D91, 043002 (2015), arXiv:1410.8866 [gr-qc] .
  * Essick  _et al._ (2020) R. Essick, I. Tews, P. Landry, S. Reddy, and D. E. Holz, Phys. Rev. C 102, 055803 (2020), arXiv:2004.07744 [astro-ph.HE] .
  * Wysocki  _et al._ (2020) D. Wysocki, R. O’Shaughnessy, L. Wade, and J. Lange, arXiv:2001.01747 (2020), arXiv:2001.01747 [gr-qc] .
  * Pang  _et al._ (2020) P. T. H. Pang, T. Dietrich, I. Tews, and C. Van Den Broeck, Phys. Rev. Res. 2, 033514 (2020), arXiv:2006.14936 [astro-ph.HE] .
  * Feroz  _et al._ (2009) F. Feroz, M. P. Hobson, and M. Bridges, Monthly Notices of the Royal Astronomical Society 398, 1601–1614 (2009).
  * Buchner, J.  _et al._ (2014) Buchner, J., Georgakakis, A., Nandra, K., Hsu, L., Rangel, C., Brightman, M., Merloni, A., Salvato, M., Donley, J., and Kocevski, D., A&A 564, A125 (2014).
  * (49) J. Buchner, “Pymultinest 2.9 documentation,” .
  * Salmi  _et al._ (2022) T. Salmi _et al._ , Astrophys. J. 941, 150 (2022), arXiv:2209.12840 [astro-ph.HE] .
  * Struder  _et al._ (2001) L. Struder _et al._ , Astron. Astrophys. 365, L18 (2001).
  * Turner  _et al._ (2001) M. J. L. Turner _et al._ , Astron. Astrophys. 365, L27 (2001), arXiv:astro-ph/0011498 .
  * Fonseca  _et al._ (2021) E. Fonseca _et al._ , Astrophys. J. Lett. 915, L12 (2021), arXiv:2104.00880 [astro-ph.HE] .
  * Foley  _et al._ (2020) R. J. Foley, D. A. Coulter, C. D. Kilpatrick, A. L. Piro, E. Ramirez-Ruiz, and J. Schwab, Mon. Not. Roy. Astron. Soc. 494, 190 (2020), arXiv:2002.00956 [astro-ph.HE] .
  * Han  _et al._ (2020) M.-Z. Han, S.-P. Tang, Y.-M. Hu, Y.-J. Li, J.-L. Jiang, Z.-P. Jin, Y.-Z. Fan, and D.-M. Wei, Astrophys. J. Lett. 891, L5 (2020), arXiv:2001.07882 [astro-ph.HE] .
  * Kyutoku  _et al._ (2020) K. Kyutoku, S. Fujibayashi, K. Hayashi, K. Kawaguchi, K. Kiuchi, M. Shibata, and M. Tanaka, Astrophys. J. 890, L4 (2020), arXiv:2001.04474 [astro-ph.HE] .
  * Collaboration and the Virgo Collaboration (2019) T. L. S. Collaboration and the Virgo Collaboration, “Parameter estimation sample release for gwtc-1,”  (2019).
  * Collaboration and the Virgo Collaboration (2020) T. L. S. Collaboration and the Virgo Collaboration, “Parameter estimation sample release for gw190425,”  (2020).
  * Bulla (2019) M. Bulla, Mon. Not. Roy. Astron. Soc. 489, 5037 (2019), arXiv:1906.04205 [astro-ph.HE] .
  * Krüger and Foucart (2020) C. J. Krüger and F. Foucart, Phys. Rev. D 101, 103002 (2020), arXiv:2002.07728 [astro-ph.HE] .
  * Bauswein  _et al._ (2020) A. Bauswein, S. Blacker, V. Vijayan, N. Stergioulas, K. Chatziioannou, J. A. Clark, N.-U. F. Bastian, D. B. Blaschke, M. Cierniak, and T. Fischer, Phys. Rev. Lett. 125, 141103 (2020), arXiv:2004.00846 [astro-ph.HE] .
  * Kataev and Molokoedov (2015) A. L. Kataev and V. S. Molokoedov, Phys. Rev. D 92, 054008 (2015).
  * Essick  _et al._ (2021a) R. Essick, I. Tews, P. Landry, and A. Schwenk, (2021a), arXiv:2102.10074 [nucl-th] .
  * Essick  _et al._ (2021b) R. Essick, P. Landry, A. Schwenk, and I. Tews, (2021b), arXiv:2107.05528 [nucl-th] .
  * Legred  _et al._ (2021) I. Legred, K. Chatziioannou, R. Essick, S. Han, and P. Landry, arXiv:2106.05313 (2021).
  * Fujimoto  _et al._ (2022) Y. Fujimoto, K. Fukushima, L. D. McLerran, and M. Praszalowicz, Phys. Rev. Lett. 129, 252702 (2022), arXiv:2207.06753 [nucl-th] .
  * Somasundaram  _et al._ (2023) R. Somasundaram, I. Tews, and J. Margueron, Phys. Rev. C 107, L052801 (2023), arXiv:2204.14039 [nucl-th] .
  * Komoltsev and Kurkela (2022) O. Komoltsev and A. Kurkela, Phys. Rev. Lett. 128, 202701 (2022), arXiv:2111.05350 [nucl-th] .
  * Gorda  _et al._ (2023b) T. Gorda, O. Komoltsev, and A. Kurkela, Astrophys. J. 950, 107 (2023b), arXiv:2204.11877 [nucl-th] .
  * Punturo  _et al._ (2010) M. Punturo _et al._ , Class. Quant. Grav. 27, 084007 (2010).
  * Hild  _et al._ (2011) S. Hild _et al._ , Class. Quant. Grav. 28, 094013 (2011), arXiv:1012.0908 [gr-qc] .
  * Reitze  _et al._ (2019) D. Reitze _et al._ , Bull. Am. Astron. Soc. 51, 035 (2019), arXiv:1907.04833 [astro-ph.IM] .


# nnSpeech: Speaker-Guided Conditional Variational Autoencoder for Zero-shot
Multi-speaker Text-to-Speech

###### Abstract

Multi-speaker text-to-speech (TTS) using a few adaption data is a challenge in
practical applications. To address that, we propose a zero-shot multi-speaker
TTS, named nnSpeech, that could synthesis a new speaker voice without fine-
tuning and using only one adaption utterance. Compared with using a speaker
representation module to extract the characteristics of new speakers, our
method bases on a speaker-guided conditional variational autoencoder and can
generate a variable Z, which contains both speaker characteristics and content
information. The latent variable Z distribution is approximated by another
variable conditioned on reference mel-spectrogram and phoneme. Experiments on
the English corpus, Mandarin corpus, and cross-dataset proves that our model
could generate natural and similar speech with only one adaption speech.

Index Terms—  zero-shot, multi-speaker text-to-speech, conditional variational
autoencoder

##  1 Introduction

Recently, with the great success of neural network based text-to-speech (TTS)
techniques, such as Tacotron [1], FastSpeech 1/2/s [2, 3], and Hifi-GAN [4]
like neural vocoder, we can efficiently synthesize more natural voice. The
single-speaker TTS model could be extended to multi-speaker scenarios using
multi-speaker corpora [2, 5]. However, these models are only applied to the
fixed set of speakers and the training requires a certain amount of human
recordings.

Multi-speaker TTS for unseen speaker has many applications, such as news
broadcast, personal assistant and metaverse. Several studies have been done in
this field and one popular solution is using few data to fine-tune the base
model that was trained by a large amount of data [6, 7]. It achieves great
voice quality when fine-tuning the whole model [8] or decoder [9]. There are
also many studies only fine-tuning the speaker embedding [10]. However, fine-
tuning based methods have two distinctive challenges: 1) Fine-tuning the model
means the memory storage and serving cost increasing because that we have to
train many new models for each speaker. 2). Fine-tuning cost a huge amount of
computing when we have a large number of users. Zero-shot multi-speaker TTS
has become a popular field in recent years [11, 12, 13], which can synthesize
voices directly by the trained base model. In some studies, pre-trained
speaker verification models are used to extract speaker information [14]. We
can also optimize the speaker encoder and TTS system jointly [15] to make the
speaker representation more optimal for the TTS system. But it is hard to
extract the speaker representation to guide the TTS system and this could
decrease voice quality or similarity. Zero-shot TTS needs training a generic
model that could perform well on unseen data. However, due to the amount of
data and the high dimension of the speaker vector, we can not estimate the
distribution of speaker representation correctly. A solution is to use more
adaption voices to get the averaged speaker representations [12]. But there
are some applications requiring TTS model generate speech using a few adaption
samples, such as personalized speech generation.

To solve the above problem, we introduce conditional variational autoencoder
(CVAE) to this task. CVAE [16, 17] is one of the most popular conditional
directed graphical models whose input observations modulate the prior on
Gaussian latent variables that generate the outputs. In this paper, we propose
a new method using both CVAE and the prior knowledge [17], speaker embeddings,
for the zero-shot multi-speaker TTS. Our model is named of nnSpeech (no new
speech) because our method can synthesize new voices directly instead of fine-
tuning a new model. In our TTS system, the phoneme is the condition. We assume
that the distribution of latent variable $Z$ with the conditions of phoneme
and speaker in multi-speaker TTS model could be approximated by the another
variable conditioned on a reference mel-spectrogram and phoneme. To better
training CVAE, we modified the standard CVAE and denotes it as speaker-guided
CVAE (sgCVAE).

##  2 METHOD

###  2.1 VAE and CVAE

Fig. 1: Graphical models of CVAE and Speaker-Guided CVAE (sgCVAE).

VAEs [18] are one of the auto-encoder neural network architectures. The
encoder network denotes the conditional distribution $q_{\phi}(Z|X)$ for the
latent space variable $Z$ given input data $X$ and the decoder network means
the distribution $p_{\theta}(X|Z)$ for the generated data $X$ given input
vector $Z$. VAEs can learn the parameters from training data to make the
$q_{\phi}(Z|X)$ consistent with the $p_{\theta}(Z|X)\propto
p_{\theta}(X|Z)p(Z)$. The log marginal distribution of data $X$ can be lower-
bounded by the evidence lower bound (ELBO) as followed:

$\mathcal{L}(\theta,\phi;X)=\mathbb{E}_{q_{\phi}(Z|X)}[log(p_{\theta}(X|Z))]-KL[q_{\phi}(Z|X)||p(Z)],$
(1)

where the $KL(q_{\phi}|p_{\theta})$ denotes the Kullback-Leibler (KL)
divergence that the minimum value is zero when the two distributions are
consistent. So we can minimize the KL divergence between $p_{\theta}(Z|X)$ and
$q_{\phi}(Z|X)$ by maximizing the ELBO, which can be further reduced to
formula Eq.1. One of the typical way to modeling $p_{\theta}(Z|X)$ and
$q_{\phi}(Z|X)$ is to assume normal distributions with a diagonal covariance
matrix. As for the prior distribution $p(Z)$, it can be designed as a specific
form according to the assumption that we need. The first term of Eq.1 is the
autoencoder reconstruction error, and the second term is the KL divergence
$q_{\phi}(Z|X)=\mathcal{N}(Z|\mu_{\phi}(x),\sigma_{\phi}(x))$ and
$p(Z)=\mathcal{N}(Z|0,I)$. Using a reparameterization
$z=\mu_{\phi}(x)+\sigma_{\phi}\odot\epsilon,\epsilon\sim\mathcal{N}(\epsilon|0,I)$,
we can sample $\epsilon$ from its distribution to generate $z$ instead of
sampling from $q_{\phi}(Z|X)$ directly, so that we can compute the gradient of
the lower bound with respect to $\theta$.

Conditional VAEs are the extended version of VAEs and are trained to maximize
the conditional log-likelihood of X given C. As shown in Fig. 1(a), the
variational lower bound to be maximized becomes

$\begin{split}\mathcal{L}(\theta,\phi;X,C)&=\mathbb{E}_{q_{\phi}(Z|X,C)}[\mathrm{log}(p_{\theta}(X|Z,C))]\\\
&-KL[q_{\phi}(Z|X,C)||p_{\theta}(Z|C)]\\\
&\leq\mathrm{log}p_{\theta}(X|C).\end{split}$ (2)

Be similar to VAEs, the prior distribution $p(Z|C)$ can be designed for the
requirements. When we assume the latent Z is independent with C, the
$p_{\theta}(Z|C)=p(Z)$. As for the neural network architecture, the
differences between VAE and CVAE are that the encoder and decoder networks can
take an auxiliary variable C as an additional input.

###  2.2 Speaker-guided CVAE for Zero-shot Multi-speaker TTS

Most multi-speaker TTS model consists of a phoneme encoder to encode the
phoneme information to a latent variable $Z$, a decoder to generate the mel-
spectrogram from $Z$ and a speaker representation added to $Z$ to control the
timbre. Speaker representation can be a jointly-optimized looked-up embedding
table or pretrained speaker verification system.

Based on the multi-speaker model, we model the multi-speaker TTS using CVAE.
The CVAE for multi-speaker TTS is composed of mel-spectrogram, phoneme, and
the latent variable. To better fuse the features, we define $X$, $C$ simply as
the output of mel encoder and phoneme encoder. The mel encoder utilizes the
speaker encoder of AdaIN-VC [19]. The phoneme encoder and mel decoder are
based on FastSpeech2 [3].

We define the conditional distribution
$p(X,Z|C)=p_{\theta}(X|C,Z)p_{\theta}(Z|C)$, and $p_{\theta}(X|C,Z)$ and
$p_{\theta}(Z|C)$ are modeled by the mel decoder network and the prior encoder
network. As shown in Fig.1(a), we use a $q_{\phi}(Z|C,X)$ to approximate the
true posterior $p_{\theta}(Z|C)$. But this assumption is counterintuitive,
because the latent variable $Z$ should consist the speaker information if we
want reconstruct the mel-spectrogram. So we propose the speaker-guided CVAE
(as shown in Fig.1(b) and Fig.2(b)) that assume the prior of latent $Z$ is
based on the speaker representation $S$, $p_{\theta}(Z|C)=p_{\theta}(Z|C,S)$,
and then use $q_{\phi}(Z|C,X)$ to approximate $p_{\theta}(Z|C,S)$. We assume
that the latent $Z$ contains the whole information of $C$, so
$p_{\theta}(X|C,Z,S)=p_{\theta}(X|Z,S)$. In speaker-guided CVAE, we modify the
conditional distribution to

$p(X,Z|C,S)=p_{\theta}(X|Z,S)p_{\theta}(Z|C,S).$

We hypothesise the latent variable $z$ follows the Gaussian distribution with
a diagonal covariance matrix, so
$q_{\phi}(Z|C,X)\sim\mathcal{N}(\mu_{1},\sigma_{1}^{2}\textbf{I})$ and
$p_{\theta}(Z|C,S)\sim\mathcal{N}(\mu_{2},\sigma_{2}^{2}\emph{{I}})$. Then the
$q_{\phi}(Z|C,X)$, named recognition network, and the $p_{\theta}(Z|C,S)$
named prior network are represented as:

$\begin{split}&\left[\begin{array}[]{cc}\mu_{1}\\\
\ln{log}\sigma^{2}_{1}\end{array}\right]=\textbf{MLP}_{q_{\phi}}(\left[\begin{array}[]{cc}X\\\
C\end{array}\right]),\\\ &\left[\begin{array}[]{cc}\mu_{2}\\\
\ln{log}\sigma^{2}_{2}\end{array}\right]=\textbf{MLP}_{p_{\theta}}(\left[\begin{array}[]{cc}C\\\
S\end{array}\right]).\end{split}$ (3)

The reparametrization trick is used to sample $Z_{1},Z_{2}$ from the
distribution $\mathcal{N}(\mu_{1},\sigma_{1}^{2}\emph{{I}})$ and
$\mathcal{N}(\mu_{2},\sigma_{2}^{2}\emph{{I}})$. Then, we can generate mel-
spectrogram $X$ based on $p_{\theta}(X|Z,S)$, which is represent by mel
decoder network. Additionally, We utilize a network to predict
$\widehat{S}=\textbf{MLP}_{S}(Z)$. In the inference step, the speaker
embedding $S$ will be replaced by predicted $\widehat{S}$. Finally, we can use
Hifi-GAN to generate voice [4].

(a) Training

(b) Inference

Fig. 2: Neural network architectures of nnSpeech.

###  2.3 Loss Functions

nnSpeech model is trained by maximizing the modified ELBO:

$\begin{split}\mathcal{L}(\theta,\phi;X,C,S)&=\mathbb{E}_{q_{\phi}(Z|X,C)}[log(p_{\theta}(X|Z,S))]\\\
&+\mathbb{E}_{q_{\phi}(Z|X,C)}[log(p_{\theta}(S|Z))]\\\
&-KL[q_{\phi}(Z|X,C)||p_{\theta}(Z|C,S)],\end{split}$ (4)

where the first two terms are the expectation of generated mel-spectrogram and
speaker embedding. Maximizing them means minimizing the mean square error of
the estimations and the coefficient of the third term, KL divergence, is
negative. Therefore, Maximizing the ELBO can be implemented by minimizing
three loss functions, mel loss, speaker loss, and KL loss. The mel loss and
speaker loss are defined as:

$\begin{split}&\mathcal{L}_{mel}=\mathbb{E}[\lvert\lvert
X-\widehat{X}\rvert\rvert_{2}],\\\ &\mathcal{L}_{spk}=\mathbb{E}[\lvert\lvert
S-\widehat{S}\rvert\rvert_{2}],\end{split}$ (5)

where the X means mel-spectrogram and S denotes speaker embedding vector. The
KL loss is calculated by

$\begin{split}\mathcal{L}_{kl}&=\frac{1}{2}\mathrm{log}\frac{|\Sigma_{2}|}{|\Sigma_{1}|}-n+tr(\Sigma_{2}^{-1}\Sigma_{1})\\\
&+(\mu_{1}-\mu_{2})^{T}\Sigma^{-1}_{2}(\mu_{1}-\mu_{2}),\end{split}$ (6)

where the $\mu_{1},\Sigma_{1}$ and $\mu_{2},\Sigma_{2}$ denotes the Gaussian
distribution $p_{\theta}$ and $q_{\phi}$. Then the total loss is:

$\mathcal{L}=\alpha*\mathcal{L}_{mel}+\beta*\mathcal{L}_{spk}+\gamma*\mathcal{L}_{kl}$
(7)

$\alpha$ and $\beta$ are set as 1, and $\gamma$ is the hyper-parameter of our
model that will be analyzed in the following experiments. Besides, our model
has duration loss, energy loss, and pitch loss like Fastspeech2.

##  3 EXPERIMENTS

###  3.1 Experimental Setup

We trained our model on a multi-speaker English corpus (903 speakers, 191.29
hours, a subset ’train-clean-360’ of the whole LibriTTS [20]), and a multi-
speaker Mandarin corpus (218 speakers, 85 hours, AISHELL-3 [21]). We randomly
chose several male and female speakers from both the LibriTTS and AISHELL3 for
evaluation. The only single speaker dataset LJSpeech [22] was used to do a
cross-dataset evaluation.

The speech data were resampled to 22050Hz and then extracted the mel-
spectrogram using 256 hop size and 1024 window size for all datasets. We also
calculated the duration, pitch, and energy for variance adaptor modular as the
same as FastSpeech2 [3].

In the training stage, we first trained the model without speaker embedding
predictor module 100,000 steps, and then trained the whole model without
speaker embedding module 500,000 steps on NVIDIA V100 GPU with the batch size
16. We used the Adam optimizer with $\beta_{1}=0.9$, $\beta_{2}=0.98$,
$\epsilon=10^{-9}$. In the testing stage, one speech was randomly selected as
the reference for each speaker.

Table 1: Subjective Comparison of the fine-tune based TTS methods Methods | Params/spk | Time/spk | MOS | VSS
---|---|---|---|---
LibriTTS | AISHELL3 | LJSpeech | LibriTTS | AISHELL3 | LJSpeech
GT | – | – | 4.73 $\pm$ 0.49 | 4.60 $\pm$ 0.60 | 4.75 $\pm$ 0.60 | – | – | –
GT+Vocoder | – | – | 4.48 $\pm$ 0.69 | 4.52 $\pm$ 0.68 | 4.71 $\pm$ 0.54 | 4.84 $\pm$ 0.68 | 4.79 $\pm$ 0.57 | 4.92 $\pm$ 0.43
Finetune spk emb | 256 | 1410 $\pm$ 396 | 3.45 $\pm$ 0.92 | 4.10 $\pm$ 0.87 | 3.33 $\pm$ 0.90 | 3.56 $\pm$ 1.03 | 3.91 $\pm$ 0.97 | 2.88 $\pm$ 1.00
Finetune decoder | 14 M | 1690 $\pm$ 526 |  3.48 $\pm$ 0.88 |  4.25 $\pm$ 0.83 |  3.39 $\pm$ 0.91 |  4.03 $\pm$ 0.98 |  4.68 $\pm$ 0.77 |  4.00 $\pm$ 0.87
X-vector [23] | – | – | 3.27 $\pm$ 0.85 | 3.95 $\pm$ 0.84 | 3.03 $\pm$ 0.94 | 3.30 $\pm$ 0.81 | 3.86 $\pm$ 0.82 | 2.21 $\pm$ 0.88
StyleSpeech [13] | – | – | 3.33 $\pm$ 0.91 | 4.07 $\pm$ 0.91 |  3.37 $\pm$ 0.90 | 3.36 $\pm$ 0.98 | 3.97 $\pm$ 0.76 |  3.24 $\pm$ 0.78
nn-Speech | – | – |  3.42 $\pm$ 0.92 |  4.13 $\pm$ 0.86 | 3.13 $\pm$ 0.98 |  3.46 $\pm$ 0.93 |  4.17 $\pm$ 0.73 | 2.65 $\pm$ 0.89

###  3.2 Evaluation

We evaluated our model both on subjective and objective tests. For subjective
test, we focused on the adaption voice quality, which consists of the voice
naturalness mean opinion score (MOS), and the voice similarity score (VSS), as
shown in Table 1. Fourteen listeners judged each sentence, and the averaged
MOS scores and VSS score were the final result. In addition, we also
calculated the parameters that needs to be stored and the consumed time of
fine-tune based methods. As for the adaption utterances, 20 sentences were
used for retraining 2,000 steps for fine-tune-based methods, and one same
reference utterance was selected randomly for the one-shot multi-speaker
methods. Besides, we calculated the mel-cepstral distortion (MCD) score as the
objective test.

We compared our method with fine-tuning the different parts of a multi-speaker
TTS model. It could be found that fine-tuning the decoder could achieve a much
better performance than fine-tuning the speaker embedding and zero-shot
methods. For zero-shot multi-speaker methods, we compared our model with two
Fastspeech2-based methods: 1). Using the x-vector [23] as the speaker
embedding in Fastspeech2; 2). StyleSpeech [13]. Our method outperformed the
x-vector in three datasets and got better performance than StyleSpeech both on
the MOS and VSS on AISHELL3 and LibriTTS. However, our model performed poorly
under cross-dataset LJspeech. Both fine-tuning decoder and using the style
adaptive layer normalization (StyleSpeech) perform better than fine-tuning
speaker embedding and using x-vector.

Table 2: Objective Comparison of the fine-tune based methods (MCD) Metric | LibriTTS | AISHELL3 | LJSpeech
---|---|---|---
Finetune spk emb | 7.55 $\pm$ 0.64 | 8.29 $\pm$ 0.65 | 7.65 $\pm$ 0.66
Finetune decoder |  7.18 $\pm$ 0.57 |  7.74 $\pm$ 0.57 |  7.05 $\pm$ 0.50
X-vector | 7.95 $\pm$ 0.65 | 8.89 $\pm$ 0.54 | 8.73 $\pm$ 0.68
StyleSpeech | 7.82 $\pm$ 0.68 | 8.64 $\pm$ 0.67 |  7.38 $\pm$ 0.51
nnSpeech |  7.44 $\pm$ 0.59 |  8.38 $\pm$ 0.53 | 8.27 $\pm$ 0.54

###  3.3 Method Analysis

In this section, we first studied the voice quality with different adaption
data. Next, we analyzed the hyper-parameter $\gamma$, the weight of KL
divergence loss. In addition, we conduct the ablation studies.

MCD under vary adaption data. We studied the MCD score with different adaption
data on the AISHELL3. For our model, we replaced the output of mel encoder $X$
and predicted $\hat{S}$ by their averaged vectors. It can be seen that the MCD
score decreases slightly with the increase of the adaption data, as shown in
Fig.3. Our method could get better performance when using four adaption
utterances, which proves that our method could extracts the speaker
information from a few adaption data.

Fig. 3: Performance with varying number of adaption voices.

Hyper-parameter analysis. As shown in Table 3, we analyzed the hyper-parameter
$\gamma$ on the dataset AISHELL3, and the model got the best performance when
the $\gamma$ was set as 0.0005 compared with the set as 0.005 and 0.05. That
is because the KL divergence is always large than MSE loss, such as
$\mathcal{L}_{mel}$ and $\mathcal{L}_{spk}$. So the model will be hard to
learn something from other losses if the weight of KL loss is large.

Table 3: The MCD with different $\gamma$ setting on AISHELL3 Metric |  $\gamma$=0.05 |  $\gamma$=0.005 |  $\gamma$=0.0005
---|---|---|---
MCD | 8.90 $\pm$ 0.49 | 8.50 $\pm$ 0.51 |  8.38 $\pm$ 0.53
Table 4: The ablation experiment on AISHELL3 Metric | Content add | w/o spk pred | nnSpeech
---|---|---|---
MCD | 8.83 $\pm$ 0.45 | 8.41 $\pm$0.63 |  8.38 $\pm$ 0.53

Analyses on the model architecture. We further compared the model architecture
with two other settings: 1). Add the content information $C$ to the latent
variable $Z$, so our model will become the standard CVAE; 2). remove the
speaker embedding predictor module and not add the predicted speaker embedding
$C$ to the latent $Z$. As shown in Table 4, we found that both settings get
the higher MCD than our model. For setting 1, we believe that the latent
variable $Z$ involves the content information and speaker information, so
adding with the content $C$ again is not necessary and maybe cover up the
speaker information. For setting 2, the result proves that the speaker
predictor module is helpful.

##  4 CONCLUSION

We proposed the speaker-guided CVAE for zero-shot multi-speaker TTS named
nnSpeech. Our method could improve the generalization of the zero-shot TTS
model. It performs excellent when only using one adaption voice both in the
English and Chinese corpus. We have done some model analysis to find the best
hyper-parameter and proves the our model is robust and each module is helpful.
However, our model performs poorly on the cross-dataset, and there is still a
gap between zero-shot methods and fine-tune methods. For future work, we will
focus on the cross-dataset problem. The result in section 3.2 shows using
adaptive layer normalization to modulate the decoder is a possible solution.

##  5 acknowledgement

This paper is supported by the Key Research and Development Program of
Guangdong Province No. 2021B0101400003 and the National Key Research and
Development Program of China under grant No. 2018YFB0204403.

## References

  * [1] Yuxuan Wang, RJ Skerry-Ryan, Daisy Stanton, Yonghui Wu, Ron J Weiss, Navdeep Jaitly, Zongheng Yang, Ying Xiao, Zhifeng Chen, Samy Bengio, et al.,  “Tacotron: Towards end-to-end speech synthesis,”  INTERSPEECH, pp. 4006–4010, 2017.
  * [2] Yi Ren, Yangjun Ruan, Xu Tan, Tao Qin, Sheng Zhao, Zhou Zhao, and Tie-Yan Liu,  “Fastspeech: Fast, robust and controllable text to speech,”  in Advances in Neural Information Processing Systems, 2019, pp. 3165–3174.
  * [3] Yi Ren, Chenxu Hu, Xu Tan, Tao Qin, Sheng Zhao, Zhou Zhao, and Tie-Yan Liu,  “Fastspeech 2: Fast and high-quality end-to-end text to speech,”  in ICLR, 2020.
  * [4] Jungil Kong, Jaehyeon Kim, and Jaekyoung Bae,  “Hifi-gan: Generative adversarial networks for efficient and high fidelity speech synthesis,”  Advances in neural information processing systems, 2020.
  * [5] Xulong Zhang, Jianzong Wang, Ning Cheng, Edward Xiao, and Jing Xiao,  “CycleGEAN:cycle generative enhanced adversarial network for voice conversion,”  in ASRU. 2021, pp. 1–6, IEEE.
  * [6] Jian Cong, Shan Yang, Lei Xie, Guoqiao Yu, and Guanglu Wan,  “Data efficient voice cloning from noisy samples with domain adversarial training,”  in INTERSPEECH, 2020, pp. 811–815.
  * [7] Mingjian Chen, Xu Tan, Bohan Li, Yanqing Liu, Tao Qin, Sheng Zhao, and Tie-Yan Liu,  “Adaspeech: Adaptive text to speech for custom voice,”  in ICLR, 2021.
  * [8] Zvi Kons, Slava Shechtman, Alex Sorin, Carmel Rabinovitz, and Ron Hoory,  “High quality, lightweight and adaptable tts using lpcnet,”  in INTERSPEECH, 2019, pp. 176–180.
  * [9] Henry B Moss, Vatsal Aggarwal, Nishant Prateek, Javier González, and Roberto Barra-Chicote,  “Boffin tts: Few-shot speaker adaptation by bayesian optimization,”  in ICASSP, 2020, pp. 7639–7643.
  * [10] Sercan O Arik, Jitong Chen, Kainan Peng, Wei Ping, and Yanqi Zhou,  “Neural voice cloning with a few samples,”  in Advances in Neural Information Processing Systems, 2018.
  * [11] Tao Wang, Jianhua Tao, Ruibo Fu, Jiangyan Yi, Zhengqi Wen, and Chunyu Qiang,  “Bi-level speaker supervision for one-shot speech synthesis,”  in INTERSPEECH, 2020, pp. 3989–3993.
  * [12] Erica Cooper, Cheng-I Lai, Yusuke Yasuda, Fuming Fang, Xin Wang, Nanxin Chen, and Junichi Yamagishi,  “Zero-shot multi-speaker text-to-speech with state-of-the-art neural speaker embeddings,”  in ICASSP, 2020, pp. 6184–6188.
  * [13] Dongchan Min, Dong Bok Lee, Eunho Yang, and Sung Ju Hwang,  “Meta-stylespeech: Multi-speaker adaptive text-to-speech generation,”  ICML, pp. 7748–7759, 2021.
  * [14] Ye Jia, Yu Zhang, Ron J Weiss, Quan Wang, Jonathan Shen, Fei Ren, Zhifeng Chen, Patrick Nguyen, Ruoming Pang, Ignacio Lopez Moreno, et al.,  “Transfer learning from speaker verification to multispeaker text-to-speech synthesis,”  in Advances in Neural Information Processing Systems, 2018.
  * [15] Mengnan Chen, Minchuan Chen, Shuang Liang, Jun Ma, Lei Chen, Shaojun Wang, and Jing Xiao,  “Cross-lingual, multi-speaker text-to-speech synthesis using neural speaker embedding,”  in INTERSPEECH, 2019, pp. 2105–2109.
  * [16] Kihyuk Sohn, Honglak Lee, and Xinchen Yan,  “Learning structured output representation using deep conditional generative models,”  in Advances in neural information processing systems, 2015.
  * [17] Tiancheng Zhao, Ran Zhao, and Maxine Eskénazi,  “Learning discourse-level diversity for neural dialog models using conditional variational autoencoders,”  in ACL, 2017, pp. 654–664.
  * [18] Diederik P Kingma and Max Welling,  “Stochastic gradient vb and the variational auto-encoder,”  in ICLR, 2017, p. 121.
  * [19] Ju-chieh Chou, Cheng-chieh Yeh, and Hung-yi Lee,  “One-shot voice conversion by separating speaker and content representations with instance normalization,”  in INTERSPEECH, 2019, pp. 664–668.
  * [20] Heiga Zen, Viet Dang, Rob Clark, Yu Zhang, Ron J. Weiss, Ye Jia, Zhifeng Chen, and Yonghui Wu,  “Libritts: A corpus derived from librispeech for text-to-speech,”  in INTERSPEECH, 2019, pp. 1526–1530.
  * [21] Yao Shi, Hui Bu, Xin Xu, Shaoji Zhang, and Ming Li,  “Aishell-3: A multi-speaker mandarin tts corpus and the baselines,”  arXiv preprint arXiv:2010.11567, 2020.
  * [22] Keith Ito and Linda Johnson,  “The lj speech dataset,” 2017.
  * [23] David Snyder, Daniel Garcia-Romero, Gregory Sell, Daniel Povey, and Sanjeev Khudanpur,  “X-vectors: Robust dnn embeddings for speaker recognition,”  in ICASSP, 2018, pp. 5329–5333.