capbtabboxtable[][]

# Tidal Deformability of Fermion-Boson Stars:
Neutron Stars Admixed with Ultra-Light Dark Matter

Robin Fynn Diedrichs<EMAIL_ADDRESS>Institute for Theoretical
Physics, Goethe University, 60438 Frankfurt am Main, Germany  Niklas Becker
<EMAIL_ADDRESS>Institute for Theoretical Physics, Goethe
University, 60438 Frankfurt am Main, Germany  Cédric Jockel<EMAIL_ADDRESS>frankfurt.de Institute for Theoretical Physics, Goethe University, 60438
Frankfurt am Main, Germany  Jan-Erik Christian<EMAIL_ADDRESS>Institute for Theoretical Physics, Goethe University, 60438 Frankfurt am Main,
Germany  Laura Sagunski<EMAIL_ADDRESS>Institute for
Theoretical Physics, Goethe University, 60438 Frankfurt am Main, Germany
Jürgen Schaffner-Bielich<EMAIL_ADDRESS>Institute for
Theoretical Physics, Goethe University, 60438 Frankfurt am Main, Germany

###### Abstract

In this work we investigate the tidal deformability of a neutron star admixed
with dark matter, modeled as a massive, self-interacting, complex scalar
field. We derive the equations to compute the tidal deformability of the full
Einstein-Hilbert-Klein-Gordon system self-consistently, and probe the
influence of the scalar field mass and self-interaction strength on the total
mass and tidal properties of the combined system. We find that dark matter
core-like configurations lead to more compact objects with smaller tidal
deformability, and dark matter cloud-like configurations lead to larger tidal
deformability. Electromagnetic observations of certain cloud-like
configurations would appear to violate the Buchdahl limit. The self-
interaction strength is found to have a significant effect on both mass and
tidal deformability. We discuss observational constraints and the connection
to anomalous detections. We also investigate how this model compares to those
with an effective bosonic equation of state and find the interaction strength
where they converge sufficiently.

##  I Introduction

Neutron stars are highly compact remnants of massive stars. Due to the high
densities inside of neutron stars, they allow us to probe nuclear matter at
high densities, a region that is not readily accessible with analytic
techniques.

The equation of state (EoS) describes the interplay between density and
pressure, which is needed to close the Tolman-Oppenheimer-Volkoff (TOV)
equations [1, 2] that describe the density profile of a spherically symmetric
star and the curvature of space-time that is produced self-consistently. A
significant constraint on the EoS is the mass value of the most massive known
compact star. If an EoS is not able to generate a star of this mass, it cannot
describe reality. There are multiple pulsars with masses at or above $2$ M⊙
[3, 4, 5, 6, 7]. Recently even a $2.35\pm 0.17$ M⊙ neutron star was reported
by Romani et al. [8]. There is also some speculation that the lighter
companion of the GW190814 gravitational wave event [9] was the most massive
neutron star ever observed, with a mass of about 2.6 M⊙. However, there is
some evidence that the object should be considered the lightest observed black
hole instead [10, 11, 12, 13, 14, 15]. Such high masses require stiff EoSs,
where the energy density strongly rises with increasing pressure. This
constraint is supported by the NICER measurements of the pulsars J0030+0451
[16, 17, 18] and J0740+6620 [19, 20, 21], which report quite large radii. The
contrary is true for the neutron star merger event GW170817 detected by
LIGO/Virgo [22, 23, 24], which favors more compact configurations generated by
soft EoSs.

It is additionally possible that neutron stars accumulate dark matter (DM) in
a sufficient abundance to modify their observables, such as the mass, radius,
and tidal deformability. These quantities have been measured in recent
observations made by, e.g., NICER [16, 17, 18, 19, 20, 21] and the
LIGO/Virgo/Kagra(LVK) collaborations, which thus allow to constraint the
properties of DM. DM is an integral part of the $\Lambda$CDM model, which is
the concordant model of cosmology [25]. Despite decades of searches, its
nature and properties are still largely unknown [26, 27]. A possible contender
is DM being made up of an additional scalar field in the universe [28, 29].
The connection between neutron stars – where the highest densities of matter
are expected – and DM has also been explored in numerous publications and is
an active area of research [30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41,
42, 43, 44, 45, 46, 47, 48, 49, 50]. DM as a scalar field could be around
neutron stars as a cloud or inside neutron stars as a core. Neutron stars with
DM cores could form 1) from a DM ‘seed’ through accretion of baryonic matter
[51], 2) through mergers of neutron stars and boson stars, 3) through
accretion and subsequent accumulation of DM inside the neutron star [30, 31,
32, 34, 36, 46, 47] or 4) through the decay of standard model particles inside
the neutron star into DM [43, 44, 45, 48, 49]. The presence of DM clouds and
cores in and around neutron stars will affect the observable properties of the
neutron stars, thus making them indirect laboratories for DM properties. It
was previously shown that even large dark matter fractions of up to 20% are
not excluded from current observations [52]. Present and future gravitational
wave detectors have the potential to detect the possible presence of DM in
merging neutron stars and to constrain the properties of DM, such as its mass
and its self-interaction strength [51, 53, 54, 55, 56, 57, 58, 59, 59, 60, 61,
62, 63, 64, 65, 66, 30, 41, 67, 68].

In this work, we model DM as a minimally coupled complex scalar field that
only interacts with the standard model (SM) via gravity and study its impact
on the neutron star observables. To this end, we construct equilibrium
solutions and their first-order perturbations and solve the coupled Einstein-
Hilbert-Klein-Gordon (EHKG) system of equations. Such systems, termed fermion-
boson stars (FBS), were first introduced by Henriques et al. [69] and
subsequently analytically studied in terms of stability under radial
perturbations [70]. In [71] they were connected to current constraints on the
mass and radii of NSs and their dynamical properties were explored in [72, 73,
74, 75]. In all of these cases, these systems were investigated using a
perfect fluid for the nuclear matter and a classical scalar field for the
bosonic DM, which is an approach that we will also follow in this work. The
described system is closely related to boson stars [76, 77, 78], as it can be
seen as a boson star that coexists with a neutron star at the same location in
space.

The tidal deformability of such systems was first investigated in [53], where
the authors considered scalar bosonic DM with masses in the MeV to GeV range,
which is gauged by a U(1) vector boson and focused on the parameter space that
results in the formation of a dark halo. They further constructed an EoS for
the bosonic sector by using mean field theory. In order to obtain solutions
for their system, they thus extended the TOV equations to account for two
fluids at the same time. This model was subsequently further investigated
first in [52] in terms of detectability prospects and in [79], where the
resulting tidal deformability was presented for a wider range of parameters
that also include scenarios in which the DM form a core. Similarly, in [57,
80], scalar DM that self-interacts via a quartic coupling was considered.
Here, the authors used an effective EoS that was first derived in [78] and
then also used the two-fluid approach. The utilized EoS is however only valid
if the self-interactions are sufficiently strong.

Our method is also applicable to scalar fields with weak to no self-
interactions. First, we review the equilibrium solutions and show the effects
of the ultralight scalar field on the mass-radius relations. Then, we derive
the relevant equations for the tidal deformability, and show the results of
our numerical investigation. We find two classes of solutions, DM cores inside
the neutron star and DM clouds enveloping the neutron star matter. Core
solutions have higher compactness and lower tidal deformabilities, while cloud
solutions can have large tidal deformabilities and lower the compactness of
the overall object. For large DM fractions, a neutron star inside these clouds
could appear to violate the Buchdahl limit, if the bosonic component is not
observed, as for example in NICER observations. The large tidal
deformabilities of the cloud solutions would be observable in LVK observations
even for small DM fractions, while DM core solutions only have small effects
on the tidal deformability, which would be difficult to discern from EoS
effects.

This paper is structured as follows: In section II we present the construction
of equilibrium solutions and further extend these equations in section III to
also include the first-order perturbations. In section IV we present the
resulting tidal deformability and compare it to observational constraints. In
section V we compare the EHKG solutions to the two-fluid model. Finally, in
section VI we summarize our findings. Throughout this work, we use units in
which $G={\rm M}_{\odot}=c=1$. See also appendix A for information on the unit
conversion.

##  II Equilibrium Solutions

Figure 1: Left panel: Density plot that displays the total gravitational mass
as a function of the central rest-mass density ($\rho_{c}$) and central value
of the scalar field ($\phi_{c}$). Additionally, it displays the stability
curve as the solid black line calculated using Eq. (18), i.e. all
configurations that lie within the bottom left parameter region that is
bordered by the black line are stable against radial perturbations. Right
panel: Mass-radius diagram displaying the fermionic radius (the radius of the
fermionic component) vs the total gravitational mass for configurations that
are within the stable region displayed in the left panel. Each point
corresponds to a single configuration and is color-coded according to the rest
mass fraction of the dark matter component. The solid black line shows the
mass-radius curve for pure fermionic matter. For both plots a massive scalar
field with no self-interactions and the mass set to $m=1.3\times 10^{-10}$ eV
was considered in addition to the DD2 EoS.

In this section, we review the construction of equilibrium solutions of FBS,
which was first presented in [69]. We model dark matter as a massive and
complex scalar field that only interacts with the SM via gravity, such that
the action of the combined system is given by

$\displaystyle S=\int
d^{4}x\sqrt{-g}\left[\frac{R}{16\pi}-\nabla_{\alpha}\bar{\Phi}\nabla^{\alpha}\Phi-V(\bar{\Phi}\Phi)+\mathcal{L}_{m}\right],$
(1)

where $\mathcal{L}_{m}$ is the Lagrangian describing nuclear matter and
$V(\bar{\Phi}\Phi)$ is the scalar field’s potential. The scalar field is
invariant under a global U(1) symmetry that gives rise to a conserved Noether
current

$j_{\mu}=i\left(\bar{\Phi}\nabla_{\mu}\Phi-\Phi\nabla_{\mu}\bar{\Phi}\right),$
(2)

which allows to generally define the total number of bosons in the system as

$N_{\text{b}}\equiv\int d^{3}x\sqrt{-g}g^{0\mu}j_{\mu}.$ (3)

The energy-momentum tensor for the scalar part is given by

$\displaystyle\begin{split}T_{\mu\nu}^{(\Phi)}=-g_{\mu\nu}\left(\partial_{\alpha}\bar{\Phi}\partial^{\alpha}\Phi+V(\bar{\Phi}\Phi)\right)\\\
+\partial_{\mu}\bar{\Phi}\partial_{\nu}\Phi+\partial_{\mu}\Phi\partial_{\nu}\bar{\Phi}.\end{split}$
(4)

Varying the action with respect to the scalar field results in the Klein-
Gordon equation,

$\nabla_{\mu}\nabla^{\mu}\Phi=\Phi
V^{\prime}(\bar{\Phi}\Phi),\;\text{with}\;V^{\prime}(\bar{\Phi}\Phi):=\frac{dV}{d|\Phi|^{2}}.$
(5)

This equation directly implies that the energy-momentum tensor of the scalar
field is separately conserved from the perfect fluid energy-momentum tensor.

The energy-momentum tensor for nuclear matter is assumed to be of the perfect
fluid form:

$T_{\mu\nu}^{\text{(NS)}}=[\rho(1+\epsilon)+P]u_{\mu}u_{\nu}+Pg_{\mu\nu},$ (6)

where $\rho$ is the rest-mass energy density and $\epsilon$ is the internal
energy density, such that $\rho(1+\epsilon)$ describes the total energy
density $e$. Requiring that the Noether current is conserved, i.e.
$\nabla_{\mu}(\rho u^{\mu})=0$, allows to define the total number of baryons
in the system generally as

$N_{\text{f}}=\int d^{3}x\sqrt{-g}g^{0\mu}\rho u_{\mu}.$ (7)

We consider the system (for now) to be in spherical symmetric equilibrium,
such that the metric can be written as

$g_{\mu\nu}=\text{diag}\left(-e^{v(r)},e^{u(r)},r^{2},r^{2}\sin^{2}\theta\right).$
(8)

We further consider a static perfect fluid, such that
$u_{\mu}=(e^{v/2},0,0,0)$ and write the scalar field as

$\Phi(t,r)=\phi_{0}(r)e^{-i\omega t},$ (9)

Using the spherical symmetric ansatz together with the Klein-Gordon equation
results in an equation describing the radial dependence of the bosonic field

$\phi_{0}^{\prime\prime}=e^{u}\left(V^{\prime}(\phi_{0}^{2})-\omega^{2}e^{-v}\right)\phi_{0}+\left(\frac{u^{\prime}-v^{\prime}}{2}-\frac{2}{r}\right)\phi_{0}^{\prime}.$
(10)

Additionally, the Einstein equations simplify to the following two equations
regarding the metric functions $u(r)$ and $v(r)$:

$\displaystyle\begin{split}u^{\prime}&=8\pi
re^{u}\Big{[}\omega^{2}\phi_{0}^{2}e^{-v}+V(\phi\bar{\phi})\\\ &\hskip
71.13188pt+e^{-u}\phi_{0}^{\prime
2}+\rho(1+\epsilon)\Big{]}-\frac{e^{u}-1}{r},\end{split}$ (11)
$\displaystyle\begin{split}v^{\prime}&=8\pi
re^{u}\Big{[}\omega^{2}\phi_{0}^{2}e^{-v}-V(\phi\bar{\phi})\\\ &\hskip
71.13188pt+e^{-u}\phi_{0}^{\prime 2}+P\Big{]}+\frac{e^{u}-1}{r}.\end{split}$
(12)

Also, the conservation of the energy-momentum tensor of nuclear matter
$\nabla_{\mu}T^{\mu\nu\,\text{(NS)}}=0$ provides a differential equation for
$P$:

$P^{\prime}=-[\rho(1+\epsilon)+P]\frac{v^{\prime}}{2}$ (13)

This system of equations is closed by providing an EoS $P(\rho,\epsilon)$ (or
$P(e)$) for the nuclear matter part.

Further, for the considered system, the expressions for the total number of
fermions (nuclear matter) and bosons (dark matter) simplify to

$\displaystyle N_{\text{b}}$
$\displaystyle=8\pi\int_{0}^{\infty}drr^{2}e^{(u-v)/2}\omega\phi_{0}^{2},$
(14) $\displaystyle N_{\text{f}}$
$\displaystyle=4\pi\int_{0}^{R_{\text{f}}}drr^{2}e^{u/2}\rho,$ (15)

where $R_{\text{f}}$ denotes the fermionic radius, which is determined by the
radial position at which the fermionic pressure $P$ vanishes. The total
gravitational mass of the system is given by

$M_{\text{tot}}=\lim_{r\rightarrow\infty}\frac{r}{2}\left(1-e^{-u(r)}\right).$
(16)

In order to integrate these equations, it is still necessary to provide
suitable initial conditions. We do this by enforcing asymptotic flatness and
regularity at the origin, i.e.

$\displaystyle\begin{split}\lim_{r\rightarrow\infty}v(r)&=0,\quad\quad
v(0)=v_{c},\\\ \lim_{r\rightarrow\infty}u(r)&=0,\quad\quad u(0)=0,\\\
\lim_{r\rightarrow\infty}\phi_{0}(r)&=0,\quad\quad\hskip-4.83694pt\phi_{0}(0)=\phi_{c},\\\
\phi_{0}^{\prime}(0)&=0,\quad\quad\rho(0)=\rho_{c}.\end{split}$ (17)

Asymptotic flatness generally requires fine-tuning $v_{c}$ to some non-zero
value. However, as was also discussed in [69], it is possible to absorb a
constant shift in $v(r)$ (e.g. $v\rightarrow v^{\prime}=v-v_{c}$) by rescaling
the above set of equations by $\omega\rightarrow\omega^{\prime}=\omega
e^{-v_{c}/2}$. This rescaling leaves the set of equations invariant and has
the advantage that we automatically have $v_{c}=0$. After integrating to
obtain a solution, we can retrieve the physical values of $\omega$ and $v$ by
doing the inverse transformation using the asymptotic value of $v(r)$.

For given $\rho_{c},\phi_{c}$ it is necessary to find the value of $\omega$,
such that the boundary conditions at infinity (Eq. (17)) are fulfilled, i.e.
the eigenvalues. There are infinitely many eigenvalues, which are
characterized by how many nodes (i.e. radial positions with $\phi(r)=0$) are
present in the scalar field profile. We find the lowest eigenvalue, such that
there are no nodes. The bosonic ODE system is such that it will always diverge
at finite radii, due to finite numerical precision. We employ this in order to
efficiently find solutions. We use the fact that the scalar field profile
either diverges towards $+\infty$ or $-\infty$ and changes its direction of
divergence when $\omega$ passes an eigenvalue. This provides us with a binary
criterion and thus allows us to implement a bisection for $\omega$ which
converges exponentially fast.

Once a sufficiently accurate $\omega$ is found, we modify the integration,
such that $\phi_{0}$ is set to zero at a finite radius $r_{B}^{*}$. This
radius $r_{B}^{*}$ is found by the condition
$\phi_{0}(r_{B}^{*})/\phi_{c}<10^{-4}$. This is necessary because otherwise,
the numerical integration diverges at finite radii. Since we have the
additional neutron matter component, in some part of the parameter space, the
integration would diverge before $P$ has converged to $0$. For example, in
compact dark matter core configurations, the integration could diverge while
still inside the neutron star component. Therefore, we artificially set
$\phi_{0}=0$ for $r>r_{B}^{*}$, which allows us to circumvent the divergence
and accurately resolve the rest of the neutron star component. The condition
was chosen such that the remaining contribution of the scalar field to the
other quantities (i.e. the metric components) would be minimized. We have
checked for lower thresholds and the extracted results are the same.

We integrate the system until a radius is reached at which both the scalar
field $\phi_{0}$ and the fermionic component parameterized by $P$ have
converged to zero. There, we can extract the properties of interest, such as
the total mass $M_{\text{tot}}$ and number densities $N_{f},N_{b}$.

For some configurations, due to numerical precision limits, the scalar field
convergence condition cannot be fulfilled. This generally happens for small
initial field values $\phi_{c}\lesssim 10^{-4}$, where the bosonic cloud
extends far outside the neutron star. In these cases, we extract the mass
$M_{\text{tot}}$ at the point where its derivative has a global minimum. When
the scalar field diverges, also the metric components do, and with it the
calculated mass of the system. By taking the point where the derivative of the
mass has a global minimum, which roughly corresponds to where the scalar field
and its derivative is closest to zero, we get the best estimate of the mass of
the system before the divergence.

Once we have a Fermion-Boson-Star solution for given $\rho_{c},\phi_{c}$, the
stability of the physical system is of importance. To this end, we need to
calculate a whole family of solutions and use the stability criterion derived
in [70], which is a generalization of the stability criterion for neutron
stars. The idea is to find extrema in the total number of particles for fixed
mass, depending on the central values $\rho_{c},\phi_{c}$. At these lines must
be the transition between stable and unstable configurations

$\frac{dN_{\text{f}}}{d\sigma}=\frac{dN_{\text{b}}}{d\sigma}=0,$ (18)

where $d/d\sigma$ denotes the derivative in the direction of constant total
mass, i.e. up to a normalization factor

$\frac{dN_{\text{f}}}{d\sigma}\propto-\frac{\partial
M_{\text{tot}}}{\partial\rho_{c}}\frac{\partial
N_{\text{f}}}{\partial\phi_{c}}+\frac{\partial
M_{\text{tot}}}{\partial\phi_{c}}\frac{\partial
N_{\text{f}}}{\partial\rho_{c}}.$ (19)

Figure 1 shows what configurations are stable depending on the central value
of the rest mass density and the central value of the scalar field according
to the above condition for the case of a massive scalar field with no self-
interactions and the mass set to $m=1.34\times 10^{-10}$ eV. Additionally, the
resulting mass and radii for the stable configurations are also displayed.

##  III Tidal Deformability

In order to obtain the tidal deformability, we will follow the same procedure
that was used in [81] to obtain the tidal deformability of pure neutron stars
and subsequently also applied to pure boson stars in [82, 83]: We are
expanding the matter and gravitational field around a static, spherically
symmetric configuration and then insert this expansion into the linearized
Einstein equations to obtain a system of differential equations that allows
solving for the linear perturbations, from which we then extract the tidal
deformability.

Applying an external quadrupolar tidal field $\mathcal{E}_{ij}$ to a
spherically symmetric star results in it developing a quadrupolar moment
$Q_{ij}$ as a response. At linear order, this response is proportional to the
applied tidal field, such that
$Q_{ij}=-\lambda_{\textrm{tidal}}\mathcal{E}_{ij}$, where
$\lambda_{\textrm{tidal}}$ is the tidal deformability. The induced quadrupolar
moment modifies the $g_{tt}$ metric component, such that at leading order in
the asymptotic rest frame at large radii [84]

$g_{tt}=-1+\frac{2M_{\rm
tot}}{r}-\mathcal{E}_{ij}x^{i}x^{j}\left(1+\frac{3\lambda_{\textrm{tidal}}}{r^{5}}\right),$
(20)

where $x^{i}$ define a Cartesian coordinate system with
$r^{2}=\delta_{ij}x^{i}x^{j}$.

We now turn to explicitly deriving the equations governing the linear
perturbations from the linearized Einstein equations. We focus on static,
even-parity, and quadrupolar ($l=2$) metric perturbations, which we denote by
$h_{\mu\nu}$. Further, we choose to work in the Regge-Wheeler gauge, in which
$h_{\mu\nu}$ takes the form

$\displaystyle\begin{split}&h_{\mu\nu}=Y_{20}(\theta,\varphi)\times\\\
&\text{diag}\left(-e^{v(r)}H_{0}(r),e^{u(r)}H_{2}(r),r^{2}K(r),r^{2}K(r)\sin^{2}\theta\right),\end{split}$
(21)

where $H_{0}$, $H_{2}$ and $K$ describe the radial dependence of each
perturbed metric component and $Y_{20}$ is the $(l,m)=(2,0)$ spherical
harmonic. At the same time, we expand the scalar field. We denote the first-
order perturbation as $\delta\Phi$ and use the same ansatz as [83], such that

$\displaystyle\begin{split}\delta\Phi(t,r,\theta,\varphi)=\phi_{1}(r)\frac{e^{-i\omega
t}}{r}Y_{20}(\theta,\varphi),\end{split}$ (22)

where the same time dependence was chosen for the perturbations in order to
ensure that the energy-momentum tensor remains static. We can obtain a set of
differential equations that relate the perturbations to the background
solutions by expanding the Einstein equations to first order in $h_{\mu\nu}$
and $\phi_{1}$.

Inserting this expansion into the Klein-Gordon equation (Eq. (5)) and only
keeping terms linear in the perturbations results in

$\displaystyle\begin{split}\phi_{1}^{\prime\prime}&=\frac{u^{\prime}-v^{\prime}}{2}\phi_{1}^{\prime}+\left[-2\phi_{0}^{\prime}-r\phi_{0}^{\prime\prime}+\frac{v^{\prime}+u^{\prime}}{2}r\phi_{0}^{\prime}+\omega^{2}r\phi_{0}e^{u-v}\right]H_{0}\\\
&+\left[\frac{6e^{u}}{r^{2}}+\frac{v^{\prime}-u^{\prime}}{2r}+16\pi\phi_{0}^{\prime
2}+e^{u}\bigg{(}V^{\prime}(\phi_{0}^{2})+2\phi_{0}^{2}V^{\prime\prime}(\phi_{0}^{2})-\omega^{2}e^{-v}\right)\bigg{]}\phi_{1}.\end{split}$
(23)

Similarly, we expand the Einstein equations, i.e. we look at $\delta
G_{\mu\nu}=8\pi\delta T_{\mu\nu}$. The perturbed energy-momentum tensor of the
fermionic part is written as $\delta
T_{\nu}^{\mu\text{(NS)}}=\text{diag}(-\delta P/c_{s}^{2},\delta P,\delta
P,\delta P)$, where we used $\delta e=\delta P\,\partial e/\partial P=\delta
P/c_{s}^{2}$, with $c_{s}$ the sound speed. The perturbed energy-momentum
tensor of the scalar field is computed by expanding Eq. (4). Subtracting the
$\theta\theta$ from the $\phi\phi$ component of the perturbed Einstein
equations reveals $H_{2}(r)=-H_{0}(r)$. Adding the $\theta\theta$ component to
the $\phi\phi$ component allows to obtain an expression for $\delta P$, which
can be substituted into the $tt$ minus the $rr$ component to obtain a
differential equation for $H_{0}$:

$\displaystyle\begin{split}&H^{\prime\prime}_{0}+\left[\frac{v^{\prime}-u^{\prime}}{2}+\frac{2}{r}\right]H^{\prime}_{0}\\\
&+\left[-8\pi\frac{1+3\,c_{s}^{2}}{c_{s}^{2}}\phi_{0}^{\prime
2}+8\pi\omega^{2}e^{u-v}\frac{c_{s}^{2}-1}{c_{s}^{2}}\phi_{0}^{2}-\frac{u^{\prime}v^{\prime}+v^{\prime
2}}{2}+v^{\prime\prime}+\frac{3u^{\prime}+7v^{\prime}}{2r}+\frac{u^{\prime}+v^{\prime}}{2rc_{s}^{2}}-\frac{6}{r}e^{u}\right]H_{0}\\\
&=\left[-\frac{16\pi}{r}\frac{1+3\,c_{s}^{2}}{c_{s}^{2}}\phi_{0}^{\prime\prime}+\frac{8\pi}{r}\left(3u^{\prime}+v^{\prime}+\frac{u^{\prime}-v^{\prime}}{c_{s}^{2}}-\frac{4}{r}\frac{1+3\,c_{s}^{2}}{c_{s}^{2}}\right)\phi_{0}^{\prime}+\frac{16\pi}{r}e^{u}\left(V^{\prime}(\phi_{0}^{2})\frac{c_{s}^{2}+1}{c_{s}^{2}}+\omega^{2}e^{-v}\frac{c_{s}^{2}-1}{c_{s}^{2}}\right)\phi_{0}\right]\phi_{1}.\end{split}$
(24)

Here primes denote derivatives with respect to the coordinate radius $r$. The
above equation contains a term depending on $v^{\prime\prime}$, which is
explicitly given by

$\displaystyle\begin{split}v^{\prime\prime}={}&8\pi
e^{u}\left(rPu^{\prime}+rP^{\prime}+P\right)+16\pi
r\phi_{0}^{\prime}\phi_{0}^{\prime\prime}+8\pi\phi_{0}^{\prime 2}+16\pi
re^{u}\left(-V^{\prime}(\phi_{0}^{2})+\omega^{2}e^{-v}\right)\phi_{0}\phi_{0}^{\prime}\\\
&-8\pi
e^{u}V(\phi_{0}^{2})\left(ru^{\prime}+1\right)+8\pi\omega^{2}e^{u}\left[r(u^{\prime}-v^{\prime})e^{-v}+e^{-v}\right]\phi_{0}^{2}+e^{u}\frac{ru^{\prime}-1}{r^{2}}+\frac{1}{r^{2}}.\end{split}$
(25)

As mentioned in [83], for radii larger than the typical size of the combined
system, the differential equation for $H_{0}$ reduces to

$H_{0}^{\prime\prime}+\left(\frac{2}{r}+e^{u}\frac{2M}{r^{2}}\right)H_{0}^{\prime}-\left(\frac{6e^{u}}{r^{2}}+e^{2u}\frac{4M^{2}}{r^{4}}\right)H_{0}=0,$
(26)

which has a solution in terms of associate Legendre functions

$H_{0}\approx
c_{1}Q_{2}^{2}\left(\frac{r}{M}-1\right)+c_{2}P_{2}^{2}\left(\frac{r}{M}-1\right).$
(27)

Expanding this equation in $r/M$ and matching to Eq. (20) results in

$\displaystyle\lambda_{\text{tidal}}={}$
$\displaystyle\frac{16}{15}M^{5}(1-2\mathcal{C})^{2}[2+2\mathcal{C}(y-1)-y]$
(28)
$\displaystyle\times\\{3(1-\mathcal{C})^{2}[2-y+2\mathcal{C}(y-1)]\log(1-2\mathcal{C})$
$\displaystyle+2\mathcal{C}[6-3y+3\mathcal{C}(5y-8)]+4\mathcal{C}^{3}[13-11y$
$\displaystyle+\mathcal{C}(3y-2)+2\mathcal{C}^{2}(1+y)]\\}^{-1},$

where $y\equiv
r_{\text{ext}}H_{0}^{\prime}(r_{\text{ext}})/H_{0}(r_{\text{ext}})$,
$\mathcal{C}\equiv M_{\text{ext}}/r_{\text{ext}}$ and $r_{\text{ext}}$ denotes
radial position at which $\lambda_{\text{tidal}}$ is calculated. The
dimensionless tidal deformability is defined as
$\Lambda_{\text{tidal}}:=\lambda_{\text{tidal}}/M_{\text{tot}}^{5}$.

In order to determine the behavior of $\phi_{1}$ and $H_{0}$ at the origin
(and thus determine the initial conditions we have to impose), we expand all
quantities around the origin as

$\displaystyle\begin{split}\phi_{1}(r)=\sum_{i}\phi_{1,i}r^{i},\end{split}$
(29)

where $\phi_{1,i}$ are the expansion coefficients that do not have any
dependency on the radius. Similarly, we also expand $H_{0}$ with $H_{0,i}$ as
the coefficients. After plugging this expansion into Eqs. (23) and (24) and
solving the resulting polynomial equations order by order results in

$\displaystyle\begin{split}\phi_{1}(r)&=\phi_{1,3}\,r^{3}+\mathcal{O}\big{(}r^{5}\big{)},\\\
H_{0}(r)&=H_{0,2}\,r^{2}+\mathcal{O}\big{(}r^{4}\big{)}.\\\ \end{split}$ (30)

We additionally impose the boundary conditions

$\displaystyle\lim_{r\to\infty}\phi_{1}=0.$ (31)

Now, we can use the fact that Eq. (24) is invariant under a simultaneous
rescaling of $\phi_{1}$ and $H_{0}$. Due to this, we can rescale the equations
to automatically have $H_{0,2}=1$. Similarly to the procedure for $\omega$, we
use a bisection algorithm to then find the initial $\phi_{1,3}$ such that the
above conditions are fulfilled. $\phi_{1}$ converges to $0$ just as
$\phi_{0}$, so we also set $\phi_{1}(r)=0$ for $r>r_{B}^{*}$. This allows us
to circumvent the divergence of the perturbations, while having no effect on
the tidal deformability, since the equations for $\phi_{1}$, $H_{0}$ decouple
with $\phi_{0}\equiv 0$. Then, the tidal deformability is constant for any
$r>r_{B}^{*}$ and can easily be extracted.

In case the convergence condition cannot be fulfilled, we follow the procedure
in [83] and extract $y$ at $r_{\text{ext}}$ such that it is a local maximum.
Since there are two components in the neutron star at play, there can be
multiple local maxima, of which we choose the one at the largest radius.

The code is publicly available along with examples and the procedures to
obtain the results.111github.com/DMGW-Goethe/FBS-Solver

##  IV Results

Figure 2: The relation between total gravitational mass $M$ and the fermionic
radius $R_{f}$ for the FBSs with different DD2 forms mass fractions. The rows
correspond to three different bosonic masses $m=\\{1,10,0.1\\}\cdot
1.34\crossproduct 10^{-10}$ eV, while the columns correspond to three
different $\Lambda_{\text{int}}=\\{0,10,100\\}$. The EoS we employ for the
fermionic part is the DD2. Notice the different scale of the bottom plots.
Observing only the fermionic radius of these systems would appear to violate
the Buchdahl limit, even though the whole FBS does not.

Figure 3: The relation between total gravitational mass $M$ and the effective
gravitational radius $R_{g}$ for the FBSs with different DD2 forms mass
fractions. The effective gravitational radius is the radius at which $99\%$ of
the rest mass is contained. The rows correspond to three different bosonic
masses $m=\\{1,10,0.1\\}\cdot 1.34\crossproduct 10^{-10}$ eV, while the
columns correspond to three different $\Lambda_{\text{int}}=\\{0,10,100\\}$.
The EoS we employ for the fermionic part is the DD2. In the case of pure
neutron stars, the crust has comparatively low density, which makes this
effective gravitational radius significantly smaller than the fermionic one.
Notice the different scales of the bottom plots. For low masses, the bosonic
component forms a core and the total compactness of the object increases. For
higher masses, the bosonic component forms a cloud and can significantly
decrease the compactness of the object.

Figure 4:  The relation between dimensionless tidal deformability
$\Lambda_{\text{tidal}}=\lambda_{\text{tidal}}/M^{5}$ and total gravitational
mass $M$ for the FBSs with different DM mass fraction. The rows correspond to
three different bosonic masses $m=\\{1,10,0.1\\}\cdot 1.34\crossproduct
10^{-10}$ eV, while the columns correspond to three different interactions
strengths $\Lambda_{\text{int}}=\\{0,10,100\\}$. For the fermionic part, we
employ the DD2 EoS.

We now specialize to a potential that is quartic in the field:

$V(\bar{\Phi}\Phi)=m^{2}\bar{\Phi}\Phi+\frac{\lambda}{2}(\bar{\Phi}\Phi)^{2},$
(32)

where $m$ is the particle mass and $\lambda$ is the self-interaction
parameter. To allow for easy comparison with previous works, we use the
effective interaction parameter $\Lambda_{\text{int}}=\lambda/(8\pi m^{2})$.
This was originally introduced in [78] to quantify the self-interaction
strength, i.e. for $\Lambda_{\text{int}}\ll 1$ the total gravitational mass of
a pure boson stars scales as $M\propto 1/m$, while for
$\Lambda_{\text{int}}\gg 1$ we have $M\propto 1/m^{2}$. Also, in this regime
the stress-energy tensor becomes approximately isotropic, meaning that an EoS
might be used to model this case (see sec. V below). It is important to keep
in mind that $\Lambda_{\text{int}}$ was introduced in the context of pure
boson stars and thus the scaling relations of the total mass are not generally
valid for the mixed system, i.e. FBSs. Nonetheless, we still find it
convenient to use it as a general measure to compare different choices of the
mass and self-interaction strength.

We investigate nine different models with $m=\\{0.1,1,10\\}\cdot
1.34\crossproduct 10^{-10}$ eV and $\Lambda_{\text{int}}=\\{0,10,100\\}$. This
mass range is chosen such that the Compton wavelength of the bosonic field is
half the Schwarzschild radius of the sun, see the explanation in Appendix A.
The range of self-interaction is well within bullet cluster constraints for
dark matter, since [85, 86]

$\displaystyle\pi\Lambda_{\text{int}}^{2}m=\frac{\lambda^{2}}{64\pi
m^{3}}=\frac{\sigma}{m}<{}$ $\displaystyle 1\frac{\text{cm}^{2}}{\text{g}}$
(33) $\displaystyle\iff\Lambda_{\text{int}}<{}$ $\displaystyle
10^{50}\sqrt{\frac{1.34\cdot 10^{-10}\text{eV}}{m}},$

where $\sigma/m$ is the effective cross-section.

For the fermionic component, we employ the DD2 EoS (with electrons) from the
CompOSE database [87, 88]. The DD2 EoS is based on a relativistic mean-field
model with density dependent coupling constants which has been fitted to the
properties of nuclei and results from Brueckner-Hartree-Fock calculations for
dense nuclear matter. Thereby, the EoS describes also the EoS of pure neutron
matter from chiral effective field theory, see [89]. For the purpose of our
investigations the particular choice of the nuclear equation of state is not
of importance and does not change our conclusions.

###  IV.1 Mass-Radius Relations and Tidal Deformability

First, the mass-radius relations are plotted in Fig. 2 and Fig. 3.

We use a grid of $\rho_{c},\phi_{c}$ to populate the plots, selecting only the
stable configurations as explained in section II. Each point is colored by the
resulting DM mass fraction $N_{B}/(N_{B}+N_{F})$. Instead of a mass-radius
curve, this gives a mass-radius region for the FBSs with different fermionic
and bosonic content. Important to note is that in Fig. 2 we plot the fermionic
radius, the radius where the fermionic component vanishes. The bosonic radius
can be orders of magnitudes larger or smaller, depending on the mass and self-
interaction parameter. To better understand these objects, we also plot the
effective gravitational radius – the radius at which $99\%$ of the rest mass
is contained – in Fig. 3. Here, the compactness of the FBS can be inferred.
For pure neutron stars with the DD2, the crust has comparatively low density,
which makes this effective gravitational radius smaller than the fermionic
one. Which radius is more relevant for a given problem depends on the
observation, e.g. the fermionic radius would be crucial for electromagnetic
signatures, such as those observed by the NICER telescope. The effective
gravitational radius would be more relevant for the inspiral in binary mergers
and enters through the compactness and the tidal deformability.

Some general trends can be seen in the figures. Stars dominated by the
fermionic part are close to the pure DD2 solution, as expected. For stars
dominated by the bosonic component, the pure boson star solutions are
recovered. For $m=\\{1,10\\}\cdot 1.34\crossproduct 10^{-10}$ eV, the regions
in Fig. 2 extend to lower masses with similar apparent compactness. These
results are consistent with the lines shown in [75]. A look at Fig. 3 reveals
the behavior of these solutions. For $m=1.34\crossproduct 10^{-9}$ eV, the
bosonic component is predominantly inside the fermionic one as a DM core. For
$m=1.34\crossproduct 10^{-10}$ eV, the bosonic and fermionic distributions
have a similar extent, for low DM mass fraction the compactness is increased,
while for higher DM mass fraction the compactness decreases as the DM forms a
cloud. This is similar to the behavior seen in [90] for a different mass
range, where increasing the DM mass fraction leads to cloud formation. For
$m=1.34\crossproduct 10^{-11}$ eV, the bosonic component completely envelops
the fermionic one in a cloud and can significantly decrease the compactness of
the object (notice the different scales on the x-axis). The apparent
compactness of the fermionic part increases on the other hand. Here, only
observing the fermionic radius as in Fig. 2 would seem like a violation of GR,
as the apparent compactness exceeds the Buchdahl limit of $4/9$.

The relation between tidal deformability and total gravitational mass is
plotted in Fig. 4. Here, we show the dimensionless tidal deformability
$\Lambda_{\text{tidal}}=\lambda_{\text{tidal}}/M^{5}$. In blue-bordered lines,
the tidal deformability of the DD2 EoS is shown, while the tidal deformability
of a pure boson star is shown in yellow-bordered lines. The latter agrees with
the trend lines shown in [83].

For $m=1.34\crossproduct 10^{-9}$ eV, the DM is mostly confined to the inner
part of the neutron star as a core and therefore does not affect the tidal
deformability significantly. Only for stars completely dominated by DM, the
results are close to the pure boson star solutions. For larger interactions
$\Lambda_{\text{int}}\approx 100$, the tidal deformability is decreased.

Figure 5: Left panel: Resulting mass and radii of FBS for the two cases of
$m=1.34\times 10^{-10}$ eV and $m=1.34\times 10^{-11}$ eV shown together with
constraints from HESS J1731-347 [91], PSR J0030+0451 [17], PSR J0740+6620
[20], PSR J1311-3430 [92] and J0952-0607 [8]. In both cases, the self-
interaction was set to zero and the percentage number denotes the DM mass
fraction. For the NICER and HESS measurements, dashed lines display the
$1\sigma$, while straight lines show the $2\sigma$ regions. Right panel:
Dimensionless tidal deformability for the same set of parameters shown
together with the constraint coming from the GW170817 event [24].

For $m=1.34\crossproduct 10^{-11}$ eV on the other hand, where the bosonic
component forms a cloud, there is a significant effect on the tidal
deformability. The tidal deformability of boson stars is much higher than the
one of purely fermionic ones, so even small amounts of DM can significantly
increase the tidal deformability of the FBS. For constant $\rho_{c}$, the
tidal deformability increases orders of magnitude as $\phi_{c}$ increases.
Then, there is a turning point where the tidal deformability decreases while
increasing total gravitational mass and converges to the purely bosonic
solutions. Overall, this opens up a vast new parameter space, even for small
DM mass fractions. While the presence of these bosonic clouds in small
quantities would barely be observable in the mass-radius plane, it would
clearly affect the tidal deformability even in small quantities, as visible in
Fig. 5.

For $m=1.34\crossproduct 10^{-10}$ eV, the behavior is more dependent on the
interaction strength $\Lambda_{\text{int}}$. For weaker interactions, the
tidal deformability stays roughly in the same order of magnitude for constant
$\rho_{c}$, while slowly converging to the pure bosonic solution for
increasing $\phi_{c}$. For stronger interactions, the tidal deformability
actually increases as it converges to the bosonic solution, as the bosonic
component starts to form a cloud. This behavior is consistent with the
observations of [57], where an effective EoS was used for modeling the bosonic
component.

###  IV.2 Comparison to Observational Constraints

There are measurements of the (fermionic) radius of neutron stars by the NICER
telescope, tracking hot spots on their surface with X-ray observations. For
the millisecond pulsar PSR J0030+0451 they derive the constraints on the mass
$M=1.34^{+0.15}_{-0.16}$ M⊙ (68%) and radius $R=12.71^{+1.14}_{-1.19}$ (68%)
[17]. A second, heavier millisecond pulsar PSR J0740+6620 has been measured at
$M=2.07^{+0.07}_{-0.07}$ M⊙ (68%) with radius $R=12.39^{+1.30}_{-0.98}$ km
(68%) [20]. These measurements constitute only two single points on the mass-
radius curve (in the neutron star case) or region (in the FBS case), but it
can show which curves/regions would support the existence of such stars.

We plot the posterior distributions of these measurements in Fig. 5 which
should be compared to the regions in Fig. 2, where the fermionic radius is
plotted.

The FBS solutions with a core become more compact depending on the DM
fraction. For higher DM fraction, they are not able to produce the maximum
mass required by the PSR measurements. DM cloud solutions on the other hand
can easily reach higher maximum masses. This is in accordance with [90], who
modeled the FBS with an effective EoS and also included the changing photon
geodesics due to the DM cloud, and [52] who performed a Bayesian analysis with
the effective EoS.

Another measurement comes from the supernova remnant HESS J1731-347. Modeling
the X-ray spectrum with accurate distance information from GAIA, they report a
mass of $M=0.77^{+0.20}_{-0.17}$ M⊙ (68%) with radius $R=10.4^{+0.86}_{-0.78}$
km (68%) [91]. This is an unusually light neutron star, which standard star
evolution theory struggles to explain, see e.g. [93]. The authors of [91]
propose it to be a strange star, but looking at Fig. 2, this region is also
well populated by DM core solutions. Of course, one would have to repeat their
analysis with an actual bosonic component to get accurate constraints, which
we leave for future work.

Lastly, there is the observation of GW170817, a binary neutron star merger.
Reference [24] has derived constraints with minimal assumptions on the nature
of the compact objects. They use a mass-weighted linear combination of the
individual tidal deformabilities and cite an upper limit of $630$.
Alternatively, assuming neutron stars with the same EoS, ref. [23] has derived
constraints on the tidal deformability with the help of universal relations
[94, 95]. These constraints are not perfectly applicable to our case, as the
I-Love-Q relations are not necessarily applicable (although they might be [59]
– we leave this for future work) and our two FBS stars might have the same EoS
but different DM mass fractions. Nevertheless, we can make some initial
guesses. It can be seen that the measurements generally favor lower tidal
deformabilities. Extrapolating this to Fig. 4, this would mean that the DM
cloud scenarios with larger tidal deformability are disfavored. Favored on the
other hand are DM core situations, which can lower the tidal deformability. A
more thorough analysis might place quantitative constraints on these models,
which we leave for future work.

Previous studies using an effective EoS description for the bosonic component
reach similar conclusions and have placed initial constraints on different
mass ranges, such as [57, 79, 96].

Overall, the different measurements seem complementary, and combining them in
a proper analysis might significantly constrain the parameter space. DM cloud
solutions can have large tidal deformabilities, even for small DM fractions,
these would most likely be observable in LVK measurements. DM core effects on
the tidal deformability on the other hand can hide inside neutron stars and
only slightly change the properties even for larger DM fractions. But core
solutions can explain the HESS measurement, and assuming different DM
fractions for different stars, they are not ruled out by the maximum mass
measurement. Of course, these effects are somewhat degenerate with the neutron
star EoS. Breaking these degeneracies requires other methods, such as looking
at correlations in the galactic DM distribution with the neutron star (FBS)
mass distribution [79, 96].

##  V Comparison with An Effective EoS

Figure 6: Distribution of the relative error (e.g. $|\Lambda_{\rm
tidal,full}-\Lambda_{\rm tidal,eff}|/\Lambda_{\rm tidal,full}$) of the
dimensionless tidal deformability $\Lambda_{\rm tidal}$ (upper panel) and
total mass $M_{\mathrm{tot}}$ (lower panel) as a function of $\Lambda_{\rm
int}$. For example, the straight line shows the boundary, below which half of
the FBS configurations lie, meaning that half of them have a relative error of
less than the shown value for a given $\Lambda_{\rm int}$. Subscripts of full
and eff denote quantities obtained from the full system and from the effective
EoS, respectively. Only stable FBS were considered for the relative error at a
given $\Lambda_{\mathrm{int}}$ and computations were performed for
$m=6.7\times 10^{-11}$ eV. The agreement between the full system and the
effective EoS becomes generally better for large $\Lambda_{\mathrm{int}}$,
however, at some point, numerical inaccuracies in the full system dominate the
relative error, which starts to be problematic for $\Lambda_{\rm int}\gtrsim
400$.

Due to the significant numerical effort associated with solving the full
system of equations (eqs. 10 \- 13) self-consistently, earlier studies [78,
80] have used an effective EoS $P(e)$ for the scalar field, treating it like a
perfect fluid with pressure $P$ and total energy density $e$. The effective
EoS was originally derived in [78] for the cases where
$\Lambda_{\mathrm{int}}=\lambda/8\pi m^{2}>0$ is large (strong self-
interactions). It models exclusively the ground state of the scalar field and
assumes an isotropic energy-momentum tensor (which is only valid in the given
limit). The EoS has the advantage that the scalar field must not be solved for
directly, and the evolution equations simplify to the default TOV-equations.
The effective EoS is given by

$\displaystyle
P=\frac{4}{9}\rho_{0}\left[\left(1+\frac{3}{4}\frac{e}{\rho_{0}}\right)^{1/2}-1\right]^{2},$
(34)

where $\rho_{0}=m^{4}/2\lambda$. Note that our expressions for $\rho_{0}$ and
$\Lambda_{\mathrm{int}}$ deviate from [78, 80] by a factor of two due to the
different normalization of the scalar field $\Phi$ and the self-interaction
parameter $\lambda$ in the potential (32). The authors of [80] used the
effective EoS in a two-fluid system of perfect fluids, which interact only
gravitationally, to compute the tidal deformability of FBS. In the following,
we compare the results obtained from integrating the two-fluid model (see [80]
for details) and from solving the full system (10-13). In addition, the tidal
deformability is computed as one would for a single-fluid system (details in
[80]) for the two-fluid model, and as described in section III for the full
system.

For the initial conditions of the two-fluid model we choose the same
conditions as in [80]. For better comparability between the full system and
the effective EoS, we first want to find an expression relating the scalar
field $\phi$ to the energy density $e_{\textrm{eff}}$ of the effective fluid.
To derive this relation, we set the $T_{tt}$ component of Eq. (4) equal to the
$T_{tt}$-component of a perfect fluid (therefore
$T^{(\Phi)}_{tt}\stackrel{{\scriptstyle!}}{{=}}e_{\textrm{eff}}\cdot e^{v}$,
and use the approximations used in [78] (i.e. neglecting spatial derivatives).
We obtain an expression that depends only on the scalar field value $\phi$
(see Eq. (17)), the scalar field mass $m$ and the self-interaction parameter
$\lambda$

$e_{\textrm{eff}}(\phi)=2m^{2}\phi^{2}+\frac{3}{2}\lambda\phi^{4},$ (35)

where $e^{-v}\omega^{2}=m^{2}+\lambda\phi^{2}$ was substituted using the
Klein-Gordon equation (10). Equation (35) holds for all radii (under the
approximations stated above). To get the initial conditions for
$e_{\textrm{eff},c}$, one simply plugs in the corresponding central value of
the scalar field $\phi_{c}$.

Figure 7: Left panel: Tidal deformability $\Lambda_{\rm tidal}$ plotted
against the total gravitational mass $M_{\rm tot}$ for pure BS and various
self-interaction strengths $\Lambda_{\mathrm{int}}$. The boson mass is
$m=6.7\times 10^{-11}$ eV in all cases. The solid lines are the values
obtained using the full system eqs. (10-13) and the dashed lines are the
corresponding solutions using the effective bosonic EoS Eq. (34). Right panel:
Tidal deformability $\Lambda_{\rm tidal}$ with respect to the total
gravitational mass $M$ of different FBS for different self-interaction
strengths $\Lambda_{\mathrm{int}}$. The boson mass is $m=6.7\times 10^{-11}$
eV in all cases. All lines have a constant central value of the scalar field
$\phi_{c}=0.02$, but different central densities $\rho_{c}$. Only stars within
the stability region are shown. The solid lines are the values obtained using
the full system eqs. (10-13) and the dashed lines are the corresponding
solutions using the effective bosonic EoS Eq. (34).

Figure 6 shows the relative error $\epsilon_{\mathrm{rel}}$ for the quantities
$M_{\mathrm{tot}}$ and the tidal deformability $\Lambda_{\mathrm{tidal}}$,
computed using the full system and the effective two-fluid system, with
respect to $\Lambda_{\mathrm{int}}$. It can be seen that the errors (the
shaded regions) generally decrease for increasing $\Lambda_{\mathrm{int}}$.
This is consistent with the assumption that the effective EoS Eq. (34) becomes
exact only in the limit of strong self-interactions. For small
$\Lambda_{\mathrm{int}}$ the relative error reaches $100\%$ for the total mass
and diverges for the tidal deformability. This is to be expected since the
total mass converges to zero for pure boson stars when using the effective EoS
in the limit $\Lambda_{\mathrm{int}}\rightarrow 0$ (see fig. 2 in [78]), while
it reaches a constant value when computing the mass using the full system.
Likewise, due to the definition of the tidal deformability (see above), a
diverging error is to be expected. For $\Lambda_{\mathrm{int}}\approx 100$ the
maximal error of the total mass (tidal deformability) is on the order of
$88\,\%$ ($>10^{4}\,\%$), whereas the lower 95-th percentiles of errors are
noticeably smaller at around $<47\,\%$ ($<240\,\%$). This means that only
$5\%$ of the computed configurations have relative errors higher than $47\,\%$
($250\,\%$). The median error denoted by the solid blue line is around $1\,\%$
($2\,\%$). At $\Lambda_{\mathrm{int}}=300$ the maximal error reaches $85\,\%$
($>10^{4}\,\%$) and the median error reaches $0.4\,\%$ ($0.8\,\%$).
Asymptotically, the error is constrained by floating-point precision and the
inherent error of the effective EoS as compared to the full system.

To gain a better understanding how the effective EoS and the full system
compare, we compute the tidal deformability $\Lambda_{\mathrm{tidal}}$ using
both systems. The left panel of figure 7 shows the tidal deformability of pure
boson stars calculated for different self-interaction strengths
$\Lambda_{\mathrm{int}}=\\{10,100,200,400\\}$. The solid lines show the
solutions using the full system and the dashed lines are the values obtained
using the effective EoS. The effective EoS can qualitatively reproduce the
solution of the full system, even for small $\Lambda_{\mathrm{int}}$. With
increasing lambda, the agreement between full and effective system becomes
better. At around $\Lambda_{\mathrm{int}}=400$, the quantitative agreement
reaches a few % relative difference.

Next, we consider the case for mixed configurations with nonzero scalar field-
and central density. The right panel of Fig. 7 shows the tidal deformability
with respect to the FBS mass. Several curves of constant central scalar field
$\phi_{c}$ were calculated at different
$\Lambda_{\mathrm{int}}=\\{10,100,200,400\\}$. The choice of constant
$\phi_{c}$ is per se arbitrary but was made for the sake of simpler
comparability with future works. The solid lines show the solutions obtained
using the full system and the dashed lines were computed with the effective
EoS (all other values being equal). With increasing $\Lambda_{\mathrm{int}}$,
the solutions using the effective EoS agree with the full system with
increasing accuracy. Even though at lower $\Lambda_{\mathrm{int}}<200$ the
deviations are quite large, the qualitative trend is correctly recovered. At
$\Lambda_{\mathrm{int}}=400$, both systems produce reasonably similar results
(within a few % of relative difference). This supports the usage of the
effective EoS for large $\Lambda_{\mathrm{int}}\gtrsim 400$ also for the
computation of the tidal deformability $\Lambda$.

A few notes on the usefulness of the effective EoS Eq. (34) and the two-fluid
system: We were able to verify the general notion, that the effective EoS
becomes asymptotically more accurate, for most configurations. However, a
significant percentage of FBS configurations with high relative errors remain,
especially when considering the tidal deformability, where the relative error
surpasses 200% for roughly five percent of all configurations. This is due to
the different low mass limits and the definition of the dimensionless tidal
deformability. Nevertheless, we conclude that the usage of the effective EoS
is justified in the cases where $\Lambda_{\mathrm{int}}\gtrsim 400$, as the
errors are acceptable for most (massive) configurations. Of course, solving
the full system eqs. (10-13) will always yield the exact results in theory. In
practice, it can be numerically difficult to integrate the full system at high
$\Lambda_{\mathrm{int}}\gtrsim 400$ because (1) the frequency $\omega$ must be
tuned up to higher accuracy than what is possible using 64-bit floating-point
numbers and (2) increasingly small step-sizes are needed, to solve the
equations correctly. During our tests, we could determine that the more
relevant constraining factor is the high needed accuracy for $\omega$, rather
than the step-size. Smaller initial $\phi_{c}$ lead to larger bosonic radii
$\gg 10$ km, for which the numerical integration becomes problematic. This
concerns $5\,\%$ of the considered configurations. In contrast, the two-fluid
system together with the effective bosonic EoS is numerically robust and does
not require numerical root-finding for $\omega$, and can manage well with
larger numerical step-sizes. With equal step-sizes and initial conditions, the
two-fluid system takes around two orders of magnitude less computation time
than solving the full system. The speedup can be increased further when
considering that the two-fluid system also tolerates larger step-sizes while
staying numerically accurate.

##  VI Conclusions

In this work, we considered the impact of a complex scalar field on the mass
and tidal deformability of neutron stars. The scalar field was assumed to be
massive and self-interacting, but to only interact gravitationally with the
fermionic neutron star matter. We derived the equations describing the linear
perturbations of the combined FBS system induced by the presence of an
external gravitational tidal field and numerically solved them to obtain the
tidal deformability of the combined system. We found that the scalar field
masses $m$ and self-interaction strengths $\lambda$ which result in the core-
like configurations of the dark matter lead to objects with higher compactness
and reduced tidal deformability. This is the case for masses $m\gtrsim
1.34\times 10^{-10}$ eV. However, large self-interactions $\lambda$ allow for
higher FBS masses or can in some cases result in cloud-like configurations. In
some of these cases, observing only the fermionic radius would appear to
violate the Buchdahl limit.

When comparing the results to available observational data of pulsars, it
becomes clear that their uncertainties are currently too large (apart from the
pulsar mass measurements) to derive quantitative constraints on the dark
matter component. The degeneracy of the effects of DM in the FBS with the EoS
poses an additional challenge. As certain DM masses can increase the total
mass of the system while leaving the fermionic radius roughly constant, this
makes previously excluded EoS possible again, if they appear in a mixed
configuration of NS matter and DM. Likewise, the unusually light neutron star
HESS J1731-347 is difficult to reconcile with known high-mass pulsar
measurements, using a regular EoS.

The relatively weak constraint from GW170817 on the tidal deformability
($\Lambda_{\mathrm{tidal}}\leq 800$ at $M_{\mathrm{tot}}\approx
1.4\,\mathrm{M}_{\odot}$) is currently also not strong enough to significantly
narrow down the dark matter properties. With the upcoming joint run of LIGO,
Virgo and KAGRA, we expect more observational data, which will enable us to
derive quantitative constraints. We plan to investigate how to constrain dark
matter properties using these observations in the future.

In addition to solving for the scalar field explicitly, we also utilized an
effective EoS to describe its contribution to the stress-energy-tensor and
reduce the complexity of this model to a two-fluid system. This approach was
recently used by [80] to compute the tidal deformability. In this work, we
compared the result of using the effective EoS to solving the full system of
equations. We found that for $m=6.7\times 10^{-11}$ eV and interactions
strengths $\Lambda_{\mathrm{i}nt}>300-400$ with
$\Lambda_{\mathrm{int}}=\lambda/(8\pi m^{2})$, the usage of the effective EoS
is typically justified. We do not expect this conclusion to be dependent on
the value of the mass $m$, but rather only on $\Lambda_{\mathrm{int}}$. Still,
even for large values of $\Lambda_{\mathrm{int}}$, we find a significant
number of configurations with relative errors of $>\mathcal{O}(10^{2})$.

Finally, it would be interesting to study the exact impact the additional
scalar field has on binary merger dynamics. In [97] this was initially studied
for a non-self-interacting scalar field. In general, it will be necessary to
extend this study to also account for self-interactions, as this can
drastically modify the FBS properties and thus impact the observed
gravitational wave signal. We will study this in detail in the future.

###### Acknowledgements.

The authors acknowledge support by the Deutsche Forschungsgemeinschaft (DFG,
German Research Foundation) through the CRC-TR 211 ‘Strong-interaction matter
under extreme conditions’– project number 315477589 – TRR 211.

##  Appendix A Units

In this work, we considered units in which $c=G=M_{\odot}=1$. As a direct
consequence, distances are measured in units of $\approx 1.48$ km,
$\hbar\approx 1.2\times 10^{-76}\neq 1$ and $m_{\text{planck}}=\sqrt{\hbar
c/G}\approx 1.1\times 10^{-38}$.

We describe the Boson star using the Klein-Gordon equation, which in SI units
and flat spacetime reads as $(\square-(mc/\hbar)^{2})\phi=0$. The term
$mc/\hbar$ is the inverse of the reduced Compton-wavelength
$\lambda_{c}=\hbar/mc$, which sets the typical length scale for the system
even in the self-gravitating case. Setting it equal to the gravitational
radius $GM/c^{2}$, which in the case of mass-scales of $\sim M_{\odot}$ is
approximately $1.48$ km, leads to $m=\hbar/c\lambda_{c}$, which corresponds to
$1.34\times 10^{-10}$ eV, which then also automatically results in Boson stars
with masses $\sim 1$ M⊙. Previous works such as e.g. [71] therefore specify
the mass of the scalar particle in units of $1.34\times 10^{-10}$ eV.

##  Appendix B Alternative Conventions

There is an alternative convention for the metric used in some publications,
e.g. [71], where the spherically symmetric, stationary metric is described by

$ds^{2}=-\alpha(r)^{2}dt^{2}+a(r)^{2}dr^{2}+r^{2}(d\theta^{2}+\sin^{2}\theta
d\phi^{2})$ (36)

The ODEs for the equilibrium solution Eq. (10)-(13) are then given by

$\displaystyle\phi_{0}^{\prime\prime}={}$
$\displaystyle\left[-\frac{\omega^{2}a^{2}}{\alpha^{2}}+a^{2}V^{\prime}\right]\phi_{0}+\left[\frac{a^{\prime}}{a}-\frac{\alpha^{\prime}}{\alpha}-\frac{2}{r}\right]\phi_{0}^{\prime}$
(37) $\displaystyle a^{\prime}={}$
$\displaystyle\frac{a}{2}\left[\frac{1-a^{2}}{r}+8\pi
ra^{2}\left(\frac{\omega^{2}\phi^{2}}{\alpha^{2}}+V+\frac{{\phi_{0}^{\prime}}^{2}}{a^{2}}+\rho(1+\epsilon)\right)\right]$
(38)
$\displaystyle\alpha^{\prime}={}$
$\displaystyle\frac{\alpha}{2}\left[\frac{a^{2}-1}{r}+8\pi
ra^{2}\left(\frac{\omega^{2}\phi^{2}}{\alpha^{2}}-V+\frac{{\phi_{0}^{\prime}}^{2}}{a^{2}}+P\right)\right]$
(39) $\displaystyle P^{\prime}={}$
$\displaystyle-\left[\rho(1+\epsilon)+P\right]\frac{\alpha^{\prime}}{\alpha}$
(40)

where $V=V(\phi\bar{\phi})$. Plugging in the quartic potential Eq. (32), this
would give the same equations as in [71], except for a different normalization
of the field $\phi_{0}$, which differs by a factor of $\sqrt{2}$.

Making an ansatz for the pertubations as in Eq. (21)

$\displaystyle\begin{split}&h_{\mu\nu}=Y_{20}(\theta,\varphi)\times\\\
&\text{diag}\left(-\alpha(r)^{2}H_{0}(r),a(r)^{2}H_{2}(r),r^{2}K(r),r^{2}K(r)\sin^{2}\theta\right),\end{split}$
(41)

and performing the same steps leads to the perturbation equations for
$H_{0},\phi_{1}$ gives

$\displaystyle
H_{0}^{\prime\prime}-\left[\frac{a^{\prime}}{a}-\frac{\alpha^{\prime}}{\alpha}-\frac{2}{r}\right]H_{0}^{\prime}$
(42)
$\displaystyle-\left[8\pi\omega^{2}\phi_{0}^{2}\frac{a^{2}}{\alpha^{2}}\frac{1-c_{s}^{2}}{c_{s}^{2}}+8\pi{\phi_{0}^{\prime}}^{2}\frac{1+3c_{s}^{2}}{c_{s}^{2}}-2\frac{\alpha^{\prime\prime}}{\alpha}+2\frac{\alpha^{\prime}a^{\prime}}{\alpha
a}+4\frac{\alpha^{\prime
2}}{\alpha^{2}}-\frac{a^{\prime}}{ra}\frac{1+3c_{s}^{2}}{c_{s}^{2}}-\frac{\alpha^{\prime}}{r\alpha}\frac{1+7c_{s}^{2}}{c_{s}^{2}}+6\frac{a^{2}}{r^{2}}\right]H_{0}$
$\displaystyle=16\pi\left[\omega^{2}\phi_{0}\frac{a^{2}}{r\alpha^{2}}\frac{c_{s}^{2}-1}{c_{s}^{2}}+\phi_{0}V^{\prime}\frac{a^{2}}{r}\frac{1+c_{s}^{2}}{c_{s}^{2}}-\frac{\phi_{0}^{\prime\prime}}{r}\frac{1+3c_{s}^{2}}{c_{s}^{2}}+\phi_{0}^{\prime}\frac{a^{\prime}}{ra}\frac{1+3c_{s}^{2}}{c_{s}^{2}}+\phi_{0}^{\prime}\frac{\alpha^{\prime}}{r\alpha}\frac{c_{s}^{2}-1}{c_{s}^{2}}-2\frac{\phi_{0}^{\prime}}{r^{2}}\frac{1+3c_{s}^{2}}{c_{s}^{2}}\right]\phi_{1}$

$\displaystyle\phi_{1}^{\prime\prime}=\left[\frac{a^{\prime}}{a}-\frac{\alpha^{\prime}}{\alpha}\right]\phi_{1}^{\prime}$
$\displaystyle+\left[-\omega^{2}\frac{a^{2}}{\alpha^{2}}+32\pi{\phi_{0}^{\prime}}^{2}+2\phi_{0}^{2}a^{2}V^{\prime\prime}+a^{2}V^{\prime}-\frac{a^{\prime}}{ra}+\frac{\alpha^{\prime}}{r\alpha}+6\frac{a^{2}}{r^{2}}\right]\phi_{1}$
(43)
$\displaystyle+\left[\omega^{2}r\phi_{0}\frac{a^{2}}{\alpha^{2}}-r\phi_{0}^{\prime\prime}+\left(r\frac{a^{\prime}}{a}+r\frac{\alpha^{\prime}}{\alpha}-2\right)\phi_{0}^{\prime}\right]H_{0}$

which need to be complemented with the explicit description for
$\alpha^{\prime\prime}$

$\displaystyle\alpha^{\prime\prime}={}$ $\displaystyle
4\pi\omega^{2}\left[2r\phi_{0}^{2}aa^{\prime}+2r\phi_{0}a^{2}\phi_{0}^{\prime}+\phi_{0}^{2}a^{2}\right]\frac{1}{\alpha}+\left[4\pi
ra^{2}\left(-\frac{\omega^{2}\phi_{0}^{2}}{\alpha^{2}}+P-V+\frac{{\phi_{0}^{\prime}}^{2}}{a^{2}}\right)+\frac{a^{2}-1}{2r}\right]\alpha^{\prime}$
(44) $\displaystyle+\left[4\pi
r(2Paa^{\prime}-2Vaa^{\prime}-2\phi_{0}a^{2}\phi_{0}^{\prime}V^{\prime}+a^{2}P^{\prime}+2\phi_{0}^{\prime}\phi_{0}^{\prime\prime})+4\pi
a^{2}(P-V)+4\pi{\phi_{0}^{\prime}}^{2}+\frac{aa^{\prime}}{r}+\frac{1-a^{2}}{2r^{2}}\right]\alpha$

We have implemented both conventions into our code and checked that they give
the same results.

## References

  * Tolman [1939] R. C. Tolman, Phys. Rev. 55, 364 (1939).
  * Oppenheimer and Volkoff [1939] J. R. Oppenheimer and G. M. Volkoff, Phys. Rev. 55, 374 (1939).
  * 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, M. H. van Kerkwijk, M. Kramer, C. Bassa, V. S. Dhillon, T. Driebe, J. W. T. Hessels, V. M. Kaspi, V. I. Kondratiev, N. Langer, T. R. Marsh, M. A. McLaughlin, T. T. Pennucci, S. M. Ransom, I. H. Stairs, J. van Leeuwen, J. P. W. Verbiest, and D. G. Whelan, Science 340, 6131 (2013), arXiv:1304.6875 [astro-ph.HE] .
  * Fonseca  _et al._ [2016] E. Fonseca, T. T. Pennucci, J. A. Ellis, and other, Astrophys. J. 832, 167 (2016), arXiv:1603.00545 [astro-ph.HE] .
  * Cromartie  _et al._ [2019] H. T. Cromartie, S. M. R. E. Fonseca, and e. a. P. B. Demorest, Nat. Astron. 4, 72 (2019), arXiv:1904.06759 [astro-ph.HE] .
  * Nieder  _et al._ [2020] L. Nieder _et al._ , Astrophys. J. 902, L46 (2020), arXiv:2009.01513 [astro-ph.HE] .
  * Romani  _et al._ [2022] R. W. Romani, D. Kandel, A. V. Filippenko, T. G. Brink, and W. Zheng, Astrophys. J. Lett. 934, L18 (2022), arXiv:2207.05124 [astro-ph.HE] .
  * Abbott  _et al._ [2020] R. Abbott, T. D. Abbott, S. Abraham, F. Acernese, and e. a. K. Ackley, The Astrophysical Journal 896, L44 (2020).
  * Most  _et al._ [2020] E. R. Most, L. J. Papenfort, L. R. Weih, and L. Rezzolla, Mon. Not. Roy. Astron. Soc. 499, L82 (2020), arXiv:2006.14601 [astro-ph.HE] .
  * Fattoyev  _et al._ [2020] F. J. Fattoyev, C. J. Horowitz, J. Piekarewicz, and B. Reed, Phys. Rev. C 102, 065805 (2020), arXiv:2007.03799 [nucl-th] .
  * Dexheimer  _et al._ [2021] V. Dexheimer, R. O. Gomes, T. Klähn, S. Han, and M. Salinas, Phys. Rev. C 103, 025808 (2021), arXiv:2007.08493 [astro-ph.HE] .
  * Tews  _et al._ [2021] I. Tews, P. T. H. Pang, T. Dietrich, M. W. Coughlin, S. Antier, M. Bulla, J. Heinzel, and L. Issa, Astrophys. J. Lett. 908, L1 (2021), arXiv:2007.06057 [astro-ph.HE] .
  * Blaschke and Cierniak [2021] D. Blaschke and M. Cierniak, Astron. Nachr. 342, 227 (2021), arXiv:2012.15785 [astro-ph.HE] .
  * Nathanail  _et al._ [2021] A. Nathanail, E. R. Most, and L. Rezzolla, Astrophys. J. Lett. 908, L28 (2021), arXiv:2101.01735 [astro-ph.HE] .
  * Miller  _et al._ [2019] M. C. Miller, F. K. Lamb, A. J. Dittmann, _et al._ , Astrophys. J. Lett. 887, L24 (2019), arXiv:1912.05705 [astro-ph.HE] .
  * Riley  _et al._ [2019] T. E. Riley, A. L. Watts, S. Bogdanov, _et al._ , Astrophys. J. Lett. 887, L21 (2019), arXiv:1912.05702 [astro-ph.HE] .
  * Raaijmakers  _et al._ [2019] G. Raaijmakers, T. E. Riley, A. L. Watts, _et al._ , Astrophys. J. Lett. 887, L22 (2019), arXiv:1912.05703 [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] .
  * Raaijmakers  _et al._ [2021] G. Raaijmakers, S. K. Greif, K. Hebeler, T. Hinderer, S. Nissanke, A. Schwenk, T. E. Riley, A. L. Watts, J. M. Lattimer, and W. C. G. Ho, Astrophys. J. Lett. 918, L29 (2021), arXiv:2105.06981 [astro-ph.HE] .
  * Abbott  _et al._ [2017] B. P. Abbott, R. Abbott, T. D. Abbott, _et al._ (Virgo, LIGO Scientific), Phys. Rev. Lett. 119, 161101 (2017), arXiv:1710.05832 [gr-qc] .
  * Abbott  _et al._ [2018] B. P. Abbott _et al._ (LIGO Scientific, Virgo), Phys. Rev. Lett. 121, 161101 (2018), arXiv:1805.11581 [gr-qc] .
  * Abbott  _et al._ [2019] B. P. Abbott, R. Abbott, T. D. Abbott, _et al._ (LIGO Scientific, Virgo), Phys. Rev. X9, 011001 (2019), arXiv:1805.11579 [gr-qc] .
  * Aghanim  _et al._ [2020] N. Aghanim _et al._ (Planck), Astron. Astrophys. 641, A6 (2020), [Erratum: Astron.Astrophys. 652, C4 (2021)], arXiv:1807.06209 [astro-ph.CO] .
  * Bertone  _et al._ [2005] G. Bertone, D. Hooper, and J. Silk, Phys. Rept. 405, 279 (2005), arXiv:hep-ph/0404175 .
  * Workman  _et al._ [2022] R. L. Workman _et al._ (Particle Data Group), PTEP 2022, 083C01 (2022).
  * Khlopov  _et al._ [1985] M. Khlopov, B. A. Malomed, and I. B. Zeldovich, Mon. Not. Roy. Astron. Soc. 215, 575 (1985).
  * Ferreira [2021] E. G. M. Ferreira, Astron. Astrophys. Rev. 29, 7 (2021), arXiv:2005.03254 [astro-ph.CO] .
  * Goldman and Nussinov [1989] I. Goldman and S. Nussinov, Phys. Rev. D 40, 3221 (1989).
  * Kouvaris [2008] C. Kouvaris, Phys. Rev. D 77, 023006 (2008), arXiv:0708.2362 [astro-ph] .
  * Kouvaris and Tinyakov [2010] C. Kouvaris and P. Tinyakov, Phys. Rev. D 82, 063531 (2010), arXiv:1004.0586 [astro-ph.GA] .
  * Sandin and Ciarcelluti [2009] F. Sandin and P. Ciarcelluti, Astropart. Phys. 32, 278 (2009), arXiv:0809.2942 [astro-ph] .
  * Ciarcelluti and Sandin [2011] P. Ciarcelluti and F. Sandin, Phys. Lett. B 695, 19 (2011), arXiv:1005.0857 [astro-ph.HE] .
  * Leung  _et al._ [2011] S. C. Leung, M. C. Chu, and L. M. Lin, Phys. Rev. D 84, 107301 (2011), arXiv:1111.1787 [astro-ph.CO] .
  * Güver  _et al._ [2014] T. Güver, A. E. Erkoca, M. Hall Reno, and I. Sarcevic, JCAP 05, 013 (2014), arXiv:1201.2400 [hep-ph] .
  * Li  _et al._ [2012] A. Li, F. Huang, and R.-X. Xu, Astropart. Phys. 37, 70 (2012), arXiv:1208.3722 [astro-ph.SR] .
  * Xiang  _et al._ [2014] Q.-F. Xiang, W.-Z. Jiang, D.-R. Zhang, and R.-Y. Yang, Phys. Rev. C 89, 025803 (2014), arXiv:1305.7354 [astro-ph.SR] .
  * Tolos  _et al._ [2015] L. Tolos, J. Schaffner-Bielich, and Y. Dengler, Phys. Rev. D 92, 123002 (2015), [Erratum: Phys.Rev.D 103, 109901 (2021)], arXiv:1507.08197 [astro-ph.HE] .
  * Mukhopadhyay  _et al._ [2017] S. Mukhopadhyay, D. Atta, K. Imam, D. N. Basu, and C. Samanta, Eur. Phys. J. C 77, 440 (2017), [Erratum: Eur.Phys.J.C 77, 553 (2017)], arXiv:1612.07093 [nucl-th] .
  * Ellis  _et al._ [2018a] J. Ellis, G. Hütsi, K. Kannike, L. Marzola, M. Raidal, and V. Vaskonen, Phys. Rev. D 97, 123007 (2018a), arXiv:1804.01418 [astro-ph.CO] .
  * McKeen  _et al._ [2018] D. McKeen, A. E. Nelson, S. Reddy, and D. Zhou, Phys. Rev. Lett. 121, 061802 (2018), arXiv:1802.08244 [hep-ph] .
  * Baym  _et al._ [2018] G. Baym, D. H. Beck, P. Geltenbort, and J. Shelton, Phys. Rev. Lett. 121, 061801 (2018), arXiv:1802.08282 [hep-ph] .
  * Motta  _et al._ [2018a] T. F. Motta, P. A. M. Guichon, and A. W. Thomas, J. Phys. G 45, 05LT01 (2018a), arXiv:1802.08427 [nucl-th] .
  * Motta  _et al._ [2018b] T. F. Motta, P. A. M. Guichon, and A. W. Thomas, Int. J. Mod. Phys. A 33, 1844020 (2018b), arXiv:1806.00903 [nucl-th] .
  * Ivanytskyi  _et al._ [2020] O. Ivanytskyi, V. Sagun, and I. Lopes, Phys. Rev. D 102, 063028 (2020), arXiv:1910.09925 [astro-ph.HE] .
  * Bell  _et al._ [2021] N. F. Bell, G. Busoni, T. F. Motta, S. Robles, A. W. Thomas, and M. Virgato, Phys. Rev. Lett. 127, 111803 (2021), arXiv:2012.08918 [hep-ph] .
  * Husain  _et al._ [2022] W. Husain, T. F. Motta, and A. W. Thomas, JCAP 10, 028 (2022), arXiv:2203.02758 [hep-ph] .
  * Berryman  _et al._ [2022] J. M. Berryman, S. Gardner, and M. Zakeri, Symmetry 14, 518 (2022), arXiv:2201.02637 [hep-ph] .
  * Cassing  _et al._ [2022] M. Cassing, A. Brisebois, M. Azeem, and J. Schaffner-Bielich, (2022), arXiv:2210.13697 [gr-qc] .
  * Ellis  _et al._ [2018b] J. Ellis, A. Hektor, G. Hütsi, K. Kannike, L. Marzola, M. Raidal, and V. Vaskonen, Phys. Lett. B 781, 607 (2018b), arXiv:1710.05540 [astro-ph.CO] .
  * Rutherford  _et al._ [2022] N. Rutherford, G. Raaijmakers, C. Prescod-Weinstein, and A. Watts, (2022), arXiv:2208.03282 [astro-ph.HE] .
  * Nelson  _et al._ [2019] A. Nelson, S. Reddy, and D. Zhou, JCAP 07, 012 (2019), arXiv:1803.03266 [hep-ph] .
  * Horowitz and Reddy [2019] C. J. Horowitz and S. Reddy, Phys. Rev. Lett. 122, 071102 (2019), arXiv:1902.04597 [astro-ph.HE] .
  * Bauswein  _et al._ [2020] A. Bauswein, G. Guo, J.-H. Lien, Y.-H. Lin, and M.-R. Wu, (2020), arXiv:2012.11908 [astro-ph.HE] .
  * Dengler  _et al._ [2022] Y. Dengler, J. Schaffner-Bielich, and L. Tolos, Phys. Rev. D 105, 043013 (2022), arXiv:2111.06197 [astro-ph.HE] .
  * Karkevandi  _et al._ [2022] D. R. Karkevandi, S. Shakeri, V. Sagun, and O. Ivanytskyi, Phys. Rev. D 105, 023001 (2022), arXiv:2109.03801 [astro-ph.HE] .
  * Cardoso  _et al._ [2016] V. Cardoso, S. Hopper, C. F. B. Macedo, C. Palenzuela, and P. Pani, Phys. Rev. D 94, 084031 (2016), arXiv:1608.08637 [gr-qc] .
  * Maselli  _et al._ [2017a] A. Maselli, P. Pnigouras, N. G. Nielsen, C. Kouvaris, and K. D. Kokkotas, Phys. Rev. D 96, 023005 (2017a), arXiv:1704.07286 [astro-ph.HE] .
  * Maselli  _et al._ [2017b] A. Maselli, S. H. Völkel, and K. D. Kokkotas, Phys. Rev. D 96, 064045 (2017b), arXiv:1708.02217 [gr-qc] .
  * Mark  _et al._ [2017] Z. Mark, A. Zimmerman, S. M. Du, and Y. Chen, Phys. Rev. D 96, 084002 (2017), arXiv:1706.06155 [gr-qc] .
  * Gresham and Zurek [2019] M. I. Gresham and K. M. Zurek, Phys. Rev. D 99, 083008 (2019), arXiv:1809.08254 [astro-ph.CO] .
  * Toubiana  _et al._ [2021] A. Toubiana, S. Babak, E. Barausse, and L. Lehner, Phys. Rev. D 103, 064042 (2021), arXiv:2011.12122 [gr-qc] .
  * Wystub  _et al._ [2021] S. Wystub, Y. Dengler, J.-E. Christian, and J. Schaffner-Bielich, (2021), arXiv:2110.12972 [astro-ph.HE] .
  * Emma  _et al._ [2022] M. Emma, F. Schianchi, F. Pannarale, V. Sagun, and T. Dietrich, Particles 5, 273 (2022), arXiv:2206.10887 [gr-qc] .
  * Hippert  _et al._ [2022] M. Hippert, E. Dillingham, H. Tan, D. Curtin, J. Noronha-Hostler, and N. Yunes, (2022), arXiv:2211.08590 [astro-ph.HE] .
  * Das  _et al._ [2019] A. Das, T. Malik, and A. C. Nayak, Phys. Rev. D 99, 043016 (2019), arXiv:1807.10013 [hep-ph] .
  * Kain [2021] B. Kain, Phys. Rev. D 103, 043009 (2021), arXiv:2102.08257 [gr-qc] .
  * Henriques  _et al._ [1989] A. Henriques, A. R. Liddle, and R. Moorhouse, Physics Letters B 233, 99 (1989).
  * Henriques  _et al._ [1990] A. Henriques, A. R. Liddle, and R. Moorhouse, Physics Letters B 251, 511 (1990).
  * Di Giovanni  _et al._ [2022] F. Di Giovanni, N. Sanchis-Gual, P. Cerdá-Durán, and J. A. Font, Phys. Rev. D 105, 063005 (2022), arXiv:2110.11997 [gr-qc] .
  * Valdez-Alvarado  _et al._ [2013] S. Valdez-Alvarado, C. Palenzuela, D. Alic, and L. A. Ureña-López, Physical Review D 87 (2013), 10.1103/physrevd.87.084040.
  * Valdez-Alvarado  _et al._ [2020] S. Valdez-Alvarado, R. Becerril, and L. A. Ureña-López, Physical Review D 102 (2020), 10.1103/physrevd.102.064038.
  * Giovanni  _et al._ [2020] F. D. Giovanni, S. Fakhry, N. Sanchis-Gual, J. C. Degollado, and J. A. Font, Physical Review D 102 (2020), 10.1103/physrevd.102.084063.
  * Giovanni  _et al._ [2021] F. D. Giovanni, S. Fakhry, N. Sanchis-Gual, J. C. Degollado, and J. A. Font, Classical and Quantum Gravity 38, 194001 (2021).
  * Kaup [1968] D. J. Kaup, Phys. Rev. 172, 1331 (1968).
  * Ruffini and Bonazzola [1969] R. Ruffini and S. Bonazzola, Phys. Rev. 187, 1767 (1969).
  * Colpi  _et al._ [1986] M. Colpi, S. L. Shapiro, and I. Wasserman, Phys. Rev. Lett. 57, 2485 (1986).
  * Giangrandi  _et al._ [2022] E. Giangrandi, V. Sagun, O. Ivanytskyi, C. Providência, and T. Dietrich, (2022), arXiv:2209.10905 [astro-ph.HE] .
  * Leung  _et al._ [2022] K.-L. Leung, M.-c. Chu, and L.-M. Lin, Phys. Rev. D 105, 123010 (2022), arXiv:2207.02433 [astro-ph.HE] .
  * Hinderer [2008] T. Hinderer, Astrophys. J. 677, 1216 (2008), arXiv:0711.2420 [astro-ph] .
  * Mendes and Yang [2017] R. F. P. Mendes and H. Yang, Class. Quant. Grav. 34, 185001 (2017), arXiv:1606.03035 [astro-ph.CO] .
  * Sennett  _et al._ [2017] N. Sennett, T. Hinderer, J. Steinhoff, A. Buonanno, and S. Ossokine, Phys. Rev. D 96, 024002 (2017), arXiv:1704.08651 [gr-qc] .
  * Thorne [1998] K. S. Thorne, Phys. Rev. D 58, 124031 (1998), arXiv:gr-qc/9706057 .
  * Eby  _et al._ [2016] J. Eby, C. Kouvaris, N. G. Nielsen, and L. C. R. Wijewardhana, JHEP 02, 028 (2016), arXiv:1511.04474 [hep-ph] .
  * Sagunski  _et al._ [2021] L. Sagunski, S. Gad-Nasr, B. Colquhoun, A. Robertson, and S. Tulin, JCAP 01, 024 (2021), arXiv:2006.12515 [astro-ph.CO] .
  * Typel  _et al._ [2010] S. Typel, G. Ropke, T. Klahn, D. Blaschke, and H. H. Wolter, Phys. Rev. C 81, 015803 (2010), arXiv:0908.2344 [nucl-th] .
  * Typel  _et al._ [2015] S. Typel, M. Oertel, and T. Klähn, Phys. Part. Nucl. 46, 633 (2015), arXiv:1307.5715 [astro-ph.SR] .
  * Krüger  _et al._ [2013] T. Krüger, I. Tews, K. Hebeler, and A. Schwenk, Phys. Rev. C 88, 025802 (2013), arXiv:1304.2212 [nucl-th] .
  * Shakeri and Karkevandi [2022] S. Shakeri and D. R. Karkevandi, (2022), arXiv:2210.17308 [astro-ph.HE] .
  * Doroshenko  _et al._ [2022] V. Doroshenko, V. Suleimanov, G. Pühlhofer, and A. Santangelo, Nature Astronomy 6, 1444–1451 (2022).
  * Kandel and Romani [2023] D. Kandel and R. W. Romani, Astrophys. J. 942, 6 (2023), arXiv:2211.16990 [astro-ph.HE] .
  * Stockinger  _et al._ [2020] G. Stockinger _et al._ , Mon. Not. Roy. Astron. Soc. 496, 2039 (2020), arXiv:2005.02420 [astro-ph.HE] .
  * Yagi and Yunes [2013a] K. Yagi and N. Yunes, Science 341, 365 (2013a), arXiv:1302.4499 [gr-qc] .
  * Yagi and Yunes [2013b] K. Yagi and N. Yunes, Phys. Rev. D 88, 023009 (2013b), arXiv:1303.1528 [gr-qc] .
  * Sagun  _et al._ [2022] V. Sagun, E. Giangrandi, O. Ivanytskyi, C. Providência, and T. Dietrich, EPJ Web Conf. 274, 07009 (2022), arXiv:2211.10510 [astro-ph.HE] .
  * Bezares  _et al._ [2019] M. Bezares, D. Viganò, and C. Palenzuela, Phys. Rev. D 100, 044049 (2019), arXiv:1905.08551 [gr-qc] .


# Basic syntax from speech: Spontaneous concatenation in unsupervised deep
neural networks

Gašper Beguš
UC Berkeley
<EMAIL_ADDRESS>
Thomas Lu
UC Berkeley
<EMAIL_ADDRESS>
Zili Wang
UC Berkeley
<EMAIL_ADDRESS>
###### Abstract

Computational models of syntax are predominantly text-based. Here we propose
that basic syntax can be modeled directly from raw speech in a fully
unsupervised way. We focus on one of the most ubiquitous and basic properties
of syntax—concatenation. We introduce spontaneous concatenation: a phenomenon
where convolutional neural networks (CNNs) trained on acoustic recordings of
individual words start generating outputs with two or even three words
concatenated without ever accessing data with multiple words in the input.
Additionally, networks trained on two words learn to embed words into novel
unobserved word combinations. To our knowledge, this is a previously
unreported property of CNNs trained on raw speech in the Generative
Adversarial Network setting and has implications both for our understanding of
how these architectures learn as well as for modeling syntax and its evolution
from raw acoustic inputs.

##  1 Introduction

Concatenation (or compounding/conjoining elements) is one the most basic
operations in human syntax. Many animal communication systems use simple
symbols (call/sign$\sim$meaning pairs) that are not concatenated (termed
“elementary signals” by Nowak and Komarova 2001). In human syntax, on the
other hand, individual elements such as words can combine into “compound
signals” (Nowak and Komarova, 2001) with compositional meaning. The evolution
of concatenation (Progovac, 2015) as well as the existence of related
operations that are presumably uniquely human and domain-specific (such as the
proposed Merge; Chomsky 2014) have been the focus of debates in linguistics
and cognitive science. Models of human syntax are predominantly text-based,
but hearing human learners acquire syntax from acoustic inputs. Modeling
syntax from raw speech also has engineering applications: speech processing
increasingly by-passes text (Lakhotia et al., 2021). Understanding syntactic
capabilities and limitations of spoken language models can inform
architectural choices.

Here, we model how compound signals or concatenated words can arise
spontaneously in deep neural networks trained on raw speech in a fully
unsupervised manner. The sounds of human speech are a measurable, physical
property, but they also encode abstract linguistic information such as
syntactic, phonological, morphological, and semantic properties. For these
reasons, we model basic syntactic dependencies from raw acoustic inputs with
CNNs. We train CNN in the Generative Adversarial Network (GAN) setting. CNNs
and GANs are uniquely appropriate for modeling linguistic dependencies from
raw speech without supervision. These models have been demonstrated to learn
disentangled near-categorical representations of linguistically meaningful
units at the phonetic, phonological, morphological, and lexical semantic
levels from raw acoustic inputs (Beguš, 2020, 2021a, 2021b).

To test whether CNNs can spontaneously concatenate words, we conduct two sets
of experiments. In the two one-word experiments, we train the networks on
single-word inputs. Because the networks are trained in the GAN setting, the
Generator never accesses the training data directly, but generates innovative
outputs. It has been shown that GANs innovate in highly interpretable ways
that produce novel words or sound sequences Beguš (2021a). Here we test
whether innovations can produce spontaneously concatenated words. In the
second experiment, we train the networks on one-word and two-word inputs (the
two-word experiment) and withhold a subset of two-word combinations. We then
test whether words can be embedded into novel unobserved combinations in the
output. Such a design also mimics one-word and two-word stages in language
acquisition (Berk and Lillo-Martin, 2012).

##  2 Methods

###  2.1 The model

We train the ciwGAN and modified fiwGAN models (Beguš, 2021a). ciwGAN/fiwGAN
models are information-theoretic extensions of GANs (based on InfoGAN, WaveGAN
and DCGAN; Chen et al. 2016; Donahue et al. 2019; Radford et al. 2015)
designed to learn from audio inputs.

The ciwGAN/fiwGAN architectures involve three networks (Figure 1): the
Generator that takes latent codes $c$ (either one-hot vectors or binary codes)
and random latent space variable $z$ ($z\sim\mathcal{U}(-1,1)$) and through
six upconvolutional layers generates 2.048s audio (32,768 samples). The audio
is then fed to the Discriminator, which evaluates realness of the output via
the Wasserstein loss Arjovsky et al. (2017). The unique aspect of the
ciwGAN/fiwGAN architecture is a separate Q-network, which is trained to
estimate the Generator’s hidden code $c$.

During training the Generator learns to generate data such that it increases
the Discriminator’s error rate and decreases the Q-network’s error rate. In
other words, the Generator needs to learn to encode unique information into
its acoustic inputs, such that the Q-network is able to decode unique
information from its generated sounds. The training between the Generator and
the Q-network mimics the production-perception loop in speech communication:
after the training, the Generator learns to generate individual words given a
latent code $c$ and the Q-network learns to classify unobserved words with the
same corresponding codes (Beguš and Zhou, 2022). Since learning is completely
unsupervised, the Generator could in principle encode any information about
speech into its latent space, but the requirement to be maximally informative
causes it to encode linguistically meaningful properties (both lexical and
sublexical information; Beguš and Zhou 2022). Such a setting not only
replicates the production-perception loop, but is also one of the few
architectures featuring traces of communicative intent (between the Generator
and the Q-network). Unlike in generative models trained on next sequence
prediction or data replication where no communicative intent exists, the
training objective between the Generator and the Q-network is to increase
mutual information between the latent space and the data such that the
Q-network can retrieve the information (latent code) encoded into the speech
signal by the Generator.

Figure 1: The architecture of ciwGAN used in the two-second one-word
experiment.

CiwGAN and fiwGAN have been shown to be highly innovative in linguistically
interpretable ways. For example, the Generator produces new words or new sound
sequences that it never accesses during training. Crucially, the Generator
never directly accesses the data: it learns by generating data from noise such
that the Discriminator fails to distinguish real and generated data. In this
respect, it mimics learning by imitation in human language (rather than
replication as is the case with variational autoencoders).

###  2.2 Data

The training dataset consists of sliced lexical items from the TIMIT database
of spoken English Garofolo et al. (1993), such that each item is a single
spoken word.

In the first, one-second one-word experiment, we use 5 lexical items:, oily,
rag, suit, year and water. In this experiment based on a pre-trained model,
the 5-layer Generator outputs only 1.024s of audio and data is never left-
padded (only right padded) which controls for the effect of padding on
concatenation. We replicate the results with another one-word experiment
trained on box, greasy, suit, under, and water. Here, each item is randomly
padded with silence to a length of 2s to produce 100 distinct data points for
each class, for a total of 500 data points used in training (the two-second
one-word experiment).

In the third experiment (two-second two-word), we use 3 lexical items: greasy,
suit, and water. 100 data points each 2s in length are generated in an
analogous process to the first experiment, but for each combination of two
items (i.e. greasy, suit, and water alone, greasy followed by water, water
followed by greasy, and so on). However, we withhold the combination
suit/greasy, such that they do not appear together in the training set in any
order, to produce a final training set of 700 data points.

For the one-word experiments, we use the ciwGAN model with one-hot encoding
and five levels, such that each of the five unique words can be represented
with a unique one-hot vector. In the two-word experiment, we use a modified
fiwGAN (binary code). The binary code is limited to three bits, but each code
can have up to two values of 1 (e.g. [1,0,0] and [1,1,0]). We also train an
additional two-word two-second model with the same data but with 6-level one-
hot $c$ in the ciwGAN architecture.

##  3 Results

To test whether the models can spontaneously concatenate, we train the
networks for 8,011 (pretrained one-second one-word), 8,956 (two-second one-
word), 9,166 (two-second two-word fiwGAN), and 18,247 steps (two-second two-
word ciwGAN) and analyze generated data. We use the technique proposed in
Beguš (2020) to analyze the relationship between linguistically meaningful
units and latent space. According to this technique, setting individual latent
space variables to values outside of the training range reveals the underlying
linguistic value of each variable.

###  3.1 One-word model

In the one-second one-word model, the Generator learns to associate each
unique one-hot code with a unique lexical item in a fully unsupervised and
unlabeled manner. The Generator’s input during training is a one-hot vector
with values 0 or 1. For example, the network learns to represent suit with
$[1,0,0,0,0]$. To test this observation, we set the one-hot vector to values
outside the training range (e.g. $[5,0,0,0,0]$), which overrides lower-level
interactions in the latent space. This causes the Generator to output suit at
near categorical levels (9 times out of 10), revealing the underlying value of
the code. This further reveals that $[0,1,0,0,0]$ encodes year (8 times out of
10) and $[0,0,1,0,0]$ encodes water (10 times out of 10), since setting latent
codes to values greater than 1 results in the model almost categorically
outputting the associated word.

In addition to lexical learning, we observe a robust pattern: the networks
trained on one-word inputs generate two-word outputs when the one-hot values
are set to negative values outside of the training range. For example, when
the latent code is set to $[0,-2,-2,-2,0]$, the Generator consistently outputs
a two-word output suit year (8 times out of 10). For $[-3,-2,-3,-2,2]$, the
network consistently outputs rag year (8 times out of 10; Figure 3). These
concatenations occur despite the fact that the training data is always left-
aligned, the Generator never accesses the data directly and the Discriminator
only sees single words.

To show that this outcome is not an idiosyncratic property of one model and
that it is indeed the negative values that encode concatenated outputs, we
also analyze a separately trained two-second one-word model. The inputs to the
Discriminator in this case are also single words only, but they are longer
(2s) and randomly padded with silence on the left and right. While high
positive values occasionally yield two-word outputs in this model, negative
values are consistently associated with two-word outputs. For example,
$[-50,-50,0,-50,0]$ (with extreme values) consistently encodes box greasy (9
times out of 10), and $[-50,-50,-50,0,0]$ consistently encodes greasy under
(10 times out of 10). Positive values of the same codes produce completely
unintelligible outputs (noise).

In addition to several two-word concatenated outputs, the network even
occasionally generates a three-word concatenated output box under water for
the latent code $c$ with all negative values $[-3,-1,-1,-1,-1]$ (2 times out
of 10). Figure 2 illustrates the three-word sequence.

Figure 2: The three-word concatentated output box under water. Independently,
the second word (under) is somewhat difficult to analyze, but given only five
training words, it is clearly the closest output to under.

###  3.2 Two-word model

In the two-word experiment, the models get one-word and two-word inputs (thus
mimicking the two-word stage). The models are only trained on three words and
their combinations, except for the suit/greasy combination withheld.

In the ciwGAN two-word model, the Generator consistently outputs the
unobserved greasy suit for $[15,0,0,0,0,0]$ (17 times out of 20), which
suggests the network learned this unobserved combination as one of the
possible sequences and encoded it with a one-hot value. For the code
$[-1,4,4]$ (modified fiwGAN), the Generator occasionally outputs the three-
word output suit greasy water (1 time out of 20; Fig. 4), which contains the
unseen suit greasy pair. It appears that the negative values of the latent
code $c$ again encode unobserved novel combinations. We also observe repeated
three-word outputs such as water water suit as a consistent output of
$[0,50,-50]$ (20 times out of 20).

###  3.3 Repetition

In addition to two-word concatenation and embedding of words into novel
combinations, we also observe outputs with repeated words in all our trained
models. The training data never includes the same word repeated, yet the
models frequently include repeated words. For example, the two-second one-word
model consistently outputs greasy greasy for $[0,0,-40,0,0]$ (7 times out of
10; Fig. 4). This is significant because repetition or reduplication is one of
the most common processes in human language and language acquisition (Berent
et al., 2016; Dolatian and Heinz, 2020). Additionally, full or total
reduplication (where the entire word is repeated) is among the most
computationally complex morphophonological processes (Dolatian and Heinz,
2020) because it represents unbound copying at the segmental level. It has
been shown elsewhere that deep convolutional neural networks can learn partial
reduplication (where only a part of the word is repeated) from raw speech and
extend it to novel unobserved data (Beguš, 2021b). Our results suggest that
total reduplication (or unbound copying; Dolatian and Heinz 2020) can arise
spontaneously in these models.

###  3.4 Why negative values?

The Generator makes use of the unutilized space in the latent code to encode
unobserved but linguistically interpretable outputs. During training, the
Generator is trained on only two values of the latent code: 0 and 1. In the
one-second one-word model, individual codes represent unique individual words,
which suggest lexical learning emerges in the positive values in these models.
The network never accesses two-word inputs and it never gets negative values
in the latent code $c$ during the training. It appears that the network is
biased to concatenate and that it uses the unobserved latent code space to
encode unobserved concatenated outputs.

##  4 Conclusion

Our results suggest that the Generator network in the ciwGAN architecture not
only learns to encode information that corresponds to lexical items in its
audio outputs, but also spontaneoulsy concatenates those lexical items into
novel unobserved two-word or three-word sequences.

The ability of unsupervised deep neural networks trained on raw speech to
concatenate words into novel unobserved combinations has far-reaching
consequences. It means that we can model basic syntactic properties directly
from raw acoustic inputs of spoken language, which opens up potential to model
several other syntactic properties directly from speech with deep
convolutional neural networks as well as from other architectures.

From the perspective of evolution of syntax, the results suggest that a deep
neural network architecture with no language-specific properties can
spontaneously begin generating concatenated signals from simple signals. The
step from one-word stage to two-word stage is necessary both in evolution of
human language as well as during language acquisition. Our second experiment
mimics the two-word stage. We argue that unsupervised deep learning models not
only concatenate single words into multi-word outputs, but are also able to
embed words into novel unobserved combinations once the model is trained on
multiple-word inputs.

Further research into the relationship between basic syntactic properties that
spontaneously emerge in these fully unsupervised models trained on raw speech
and the structure of the latent space has the potential to yield insights for
the study of syntactic theory, language acquisition, and language evolution.
By evaluating these models on syntactic properties of spoken language, we
should also get a better understanding of computational limits of unsupervised
CNNs.

## Limitations

This paper models concatenation of acoustic lexical items. Syntax is
substantially more complex than concatenation (Berwick and Chomsky, 2019).
Exploration of other syntactic properties as well as of compositionality in
these models is left for future work. We also train the network on a
relatively small number of lexical items (5) and a small number of tokens
(100). The small number of lexical items is representative of the earliest
stages of language acquisition when the number of lexical items is highly
limited (Bates et al., 1994).

## Ethics Statement

Two models are trained for the purpose of this paper, and one model is
pretrained. The three models were trained for 16hrs on a single GPU (NVIDIA
1080ti). We use the TIMIT Garofolo et al. (1993) database for training. The
number of parameters is given in the Appendix A. We take the standard
hyperparameters (from Donahue et al. 2019 and Beguš 2021a). Because the
outputs are salient and rarely ambiguous, all transcriptions are performed by
the authors. Generated audio files and models’ checkpoints are available at
the anonymous link:
https://osf.io/przuq/?view_only=9d19a26f0bb84a3ea4db8e6844b37985.

## Acknowledgements

## References

  * Arjovsky et al. (2017) Martin Arjovsky, Soumith Chintala, and Léon Bottou. 2017.  Wasserstein generative adversarial networks.  In _Proceedings of the 34th International Conference on Machine Learning_ , volume 70 of _Proceedings of Machine Learning Research_ , pages 214–223, International Convention Centre, Sydney, Australia. PMLR.
  * Bates et al. (1994) Elizabeth Bates, Virginia Marchman, Donna Thal, Larry Fenson, Philip Dale, J. Steven Reznick, Judy Reilly, and Jeff Hartung. 1994.  Developmental and stylistic variation in the composition of early vocabulary.  _Journal of Child Language_ , 21(1):85–123.
  * Beguš (2020) Gašper Beguš. 2020.  Generative adversarial phonology: Modeling unsupervised phonetic and phonological learning with neural networks.  _Frontiers in Artificial Intelligence_ , 3:44.
  * Beguš (2021a) Gašper Beguš. 2021a.  CiwGAN and fiwGAN: Encoding information in acoustic data to model lexical learning with generative adversarial networks.  _Neural Networks_ , 139:305–325.
  * Beguš (2021b) Gašper Beguš. 2021b.  Identity-based patterns in deep convolutional networks: Generative adversarial phonology and reduplication.  _Transactions of the Association for Computational Linguistics_ , 9:1180–1196.
  * Beguš and Zhou (2022) Gašper Beguš and Alan Zhou. 2022.  Modeling speech recognition and synthesis simultaneously: Encoding and decoding lexical and sublexical semantic information into speech with no direct access to speech data.  In _Proc. Interspeech 2022_ , pages 5298–5302.
  * Berent et al. (2016) Iris Berent, Outi Bat-El, Diane Brentari, Amanda Dupuis, and Vered Vaknin-Nusbaum. 2016.  The double identity of linguistic doubling.  _Proceedings of the National Academy of Sciences_ , 113(48):13702–13707.
  * Berk and Lillo-Martin (2012) Stephanie Berk and Diane Lillo-Martin. 2012.  The two-word stage: Motivated by linguistic or cognitive constraints? _Cognitive Psychology_ , 65(1):118–140.
  * Berwick and Chomsky (2019) Robert C. Berwick and Noam Chomsky. 2019.  All or nothing: No half-merge and the evolution of syntax.  _PLOS Biology_ , 17(11):1–5.
  * Boersma and Weenink (2015) Paul Boersma and David Weenink. 2015.  Praat: doing phonetics by computer [computer program]. version 5.4.06.  Retrieved 21 February 2015 from http://www.praat.org/.
  * Chen et al. (2016) Xi Chen, Yan Duan, Rein Houthooft, John Schulman, Ilya Sutskever, and Pieter Abbeel. 2016.  Infogan: Interpretable representation learning by information maximizing generative adversarial nets.  In D. D. Lee, M. Sugiyama, U. V. Luxburg, I. Guyon, and R. Garnett, editors, _Advances in Neural Information Processing Systems 29_ , pages 2172–2180. Curran Associates, Inc.
  * Chomsky (2014) Noam Chomsky. 2014.  _The Minimalist Program_.  The MIT Press.
  * Dolatian and Heinz (2020) Hossep Dolatian and Jeffrey Heinz. 2020.  Computing and classifying reduplication with 2-way finite-state transducers.  _Journal of Language Modelling_ , 8(1):179–250.
  * Donahue et al. (2019) Chris Donahue, Julian J. McAuley, and Miller S. Puckette. 2019.  Adversarial audio synthesis.  In _7th International Conference on Learning Representations, ICLR 2019, New Orleans, LA, USA, May 6-9, 2019_ , pages 1–16. OpenReview.net.
  * Garofolo et al. (1993) J. S. Garofolo, Lori Lamel, W M Fisher, Jonathan Fiscus, D S. Pallett, N L. Dahlgren, and V Zue. 1993.  TIMIT acoustic-phonetic continuous speech corpus.  _Linguistic Data Consortium_.
  * Lakhotia et al. (2021) Kushal Lakhotia, Eugene Kharitonov, Wei-Ning Hsu, Yossi Adi, Adam Polyak, Benjamin Bolte, Tu-Anh Nguyen, Jade Copet, Alexei Baevski, Abdelrahman Mohamed, and Emmanuel Dupoux. 2021.  On generative spoken language modeling from raw audio.  _Transactions of the Association for Computational Linguistics_ , 9:1336–1354.
  * Nowak and Komarova (2001) Martin A Nowak and Natalia L Komarova. 2001.  Towards an evolutionary theory of language.  _Trends in Cognitive Sciences_ , 5(7):288–295.
  * Progovac (2015) Ljiljana Progovac. 2015.  _Evolutionary syntax_ , first edition edition.  Oxford studies in the evolution of language ; 20. Oxford University Press, Oxford, United Kingdom ; New York, NY.
  * Radford et al. (2015) Alec Radford, Luke Metz, and Soumith Chintala. 2015.  Unsupervised representation learning with deep convolutional generative adversarial networks.  _arXiv preprint arXiv:1511.06434_.

##  Appendix A Appendix

Generator | Discriminator
---|---
Layer | Dimension | Layer | Dimension
$z$ | 100 $\times$ 1 | input | 32768 $\times$ 1
fc + reshape | 16 $\times$ 2048 | conv0 | 8192 $\times$ 64
upconv0 | 64 $\times$ 1024 | conv1 | 2048 $\times$ 128
upconv1 | 256 $\times$ 512 | conv2 | 512 $\times$ 256
upconv2 | 1024 $\times$ 256 | conv3 | 128 $\times$ 512
upconv3 | 4096 $\times$ 128 | conv4 | 32 $\times$ 1024
upconv4 | 16384 $\times$ 64 | conv5 | 16 $\times$ 2048
upconv5 | 32768 $\times$ 1 | flatten + logit | 1 $\times$ 1

Q-Network
---
Layer | Dimension
input | 32768 $\times$ 1
conv0 | 8192 $\times$ 64
qconv1 | 2048 $\times$ 128
qconv2 | 512 $\times$ 256
qconv3 | 128 $\times$ 512
qconv4 | 32 $\times$ 1024
qconv5 | 16 $\times$ 2048
flatten + logit | 1 $\times$ 1

Table 1: The structure of the Generator, the Discriminator, and the Q-network
in the two-second experiments (based on Donahue et al. 2019; Beguš 2021a).

Figure 3: The suit year (top) output and the rag year (bottom) from the one-
second one-word model. All spectrograms are created in Praat Boersma and
Weenink (2015).

Figure 4: The suit greasy water output from the two-second two-word model
(top). The greasy greasy output from the two-second one-word model (bottom).