# Small-scale structure in vector dark matter

Mustafa A. Amin  Mudit Jain  Rohith Karur  Philip Mocz

###### Abstract

We investigate the differences in the small-scale structure of vector dark
matter (VDM) and scalar dark matter (SDM) using 3+1 dimensional simulations of
single/multicomponent Schrödinger-Poisson system. We find that the amount of
wave interference, core-to-halo mass ratio (and its scatter), spin of the
core, as well as the shape of the central regions of dark matter halos can
distinguish VDM and SDM. Starting with a collection of idealized halos (self-
gravitating solitons) as an initial condition, we show that the system
dynamically evolves to an approximately spherically symmetric configuration
that has a core surrounded by a halo of interference patterns in the mass
density. In the vector case, the central soliton in less dense and has a
smoother transition to an $r^{-3}$ tail compared to the scalar case. As
compared to SDM, wave interference in VDM is $\sim 1/\sqrt{3}$ times smaller,
resulting in fewer low and high density regions, and more diffuse granules in
the halo. The ratio of VDM core mass to the total halo mass is lower than that
in SDM, with a steeper dependence on the total energy of the system and a
slightly larger scatter. Finally, we also initiate a study of the evolution of
intrinsic spin angular momentum in the VDM case. We see a positive correlation
between the total intrinsic spin in the simulation and the spin of the final
central core, with significant scatter. We see large intrinsic spin in the
core being possible even with vanishing amounts total angular momentum in the
initial conditions (at least instantaneously). Our results point towards the
possibility of distinguishing VDM from SDM using astrophysical and terrestrial
observations.

##  1 Introduction

The identity of particles/fields that make up dark matter remains uncertain.
Observational evidence [1] suggests that dark matter is very weakly
interacting with the Standard model, is non-relativistic in the contemporary
universe, and has clumped efficiently under the influence of gravity (for a
historical overview, see [2]). However, the mass and spin of the fundamental
quanta of dark matter are not known. The mass of the fundamental particles
making up dark matter can range from $\sim 10^{-21}\rm eV$ [3] to $\sim
m_{\mathrm{pl}}$. The bounds are softened further if dark matter is multi-
component, or composite [4]. We have no robust constraints on the spin of the
particles/fields that make up dark matter either.

When dark matter is light, $m\ll 10\,\textrm{eV}$, and bosonic, the typical
occupation number of the field and the overlap of the particle de Broglie
wavelengths is high enough in astrophysical and cosmological settings that a
classical wave description becomes appropriate (as opposed to discrete point-
like particle as in CDM)[5]. This “wave dark matter” results in novel
phenomenon resulting from coherence and interference of the effective
classical field, common to all wave dynamics. If the mass of the underlying
boson is sufficiently small, these wave effects can manifest themselves on
macroscopic and even astrophysical scales. Significant effort has been devoted
in recent years to study such wave phenomenon, including suppression of
structure on small scales [6, 7], solitons [8, 9, 10], interference patterns
[8], turbulence [11], vortices [12], superradiance [13] etc. in the context of
spin-0 dark matter (scalar dark matter or SDM) in the nonrelativistic limit.
See [14, 5, 15, 16] for reviews and references therein. This is natural since
the well-motivated QCD axion, as well as ultralight axions motivated by the
String Axiverse [17, 18] construction, generically predict the existence of
such light, weakly interacting scalars.

There is no evidence yet, however, that dark matter is a light scalar field.
It is natural to ask about the intrinsic spin of the underlying wave dark
matter. Gravitational production of light (but not ultralight) vector dark
matter was explored in detail in [19] (also see [20, 21, 22]). Vector dark
matter production was explored in a misaligned axion/Higgs scenario in [23,
24, 25, 26], with the misaligned scenarios allowing for ultralight masses.
Numerical simulations of vector dark matter production in the misaligned
scenario were carried out in [23]. However, exploring the astrophysical
implications of vector wave dark matter in the contemporary universe on
nonlinear scales has been restricted to analytic calculations [27, 28, 29].

With this in mind, we carry out 3+1 dimensional simulations of spin-$1$ wave
dark matter (which we refer to as vector dark matter or VDM) in the
nonrelativistic limit. In this limit, the dynamics of vector dark matter can
be simulated by a multicomponent Schrödinger-Poisson system [27, 28]. This
essentially allows us to use the existing numerical framework used for scalar
dark matter, but now with three (equal mass) components. We restrict our
attention to the mergers of idealized halos, in the form of solitons, to
explore how structure formation might proceed on nonlinear scales in VDM. We
characterize the core and halo resulting from these mergers, and compare the
results with the corresponding SDM case. For identical initial conditions in
terms of their initial mass density, the final configurations for SDM and VDM
differ from each other. Even if we ignore the knowledge of initial conditions,
and simply compare the shape of the final configuration, we are able to
distinguish the VDM from the SDM case. We also track the spin density of VDM
(which is trivially zero in SDM) for the first time, which could provide
another novel handle on VDM when coupled to (for example) the electromagnetic
field, or when approaching the relativistic limit [28]. Thus, SDM and VDM can
be potentially distinguished in an astrophysical setting giving us hope of
observationally probing a fundamental degree of freedom (spin) of ultralight
dark matter. The key to understanding the differences between SDM and VDM is
that the wave-interference effects are smaller in VDM compared to SDM.111We
note that in the non-relativistic limit and without additional interactions,
we cannot distinguish between three scalars with identical masses and a vector
field with the same boson mass (although effects related to the production
mechanism might hold clues to the underlying nature of the field). The spin
discussed here can then be thought of as an isospin [28]. Separately, for
phenomenology of multiple axion species with different masses in the late
universe, see for example, [17, 30, 31, 32, 33].

The rest of the paper is organized as follows. In Sec. 2, we introduce our
model for VDM along with its nonrelativistic limit. We also provide an
understanding of interference in VDM waves, as well as solitons in VDM. We
explore binary soliton mergers in Sec. 3, and calculate the fraction of total
mass that remains bound in the final soliton. In Sec. 4, we consider the
merger of $N=\mathcal{O}(10)$ solitons. We compare the results of the merger
in VDM and SDM, including core mass, density profiles, size of interference
granules, as well as spin angular momentum density. In Sec. 5, we briefly
discuss observational implications including dynamical heating of stars, cores
of dwarf galaxies, and DM substructure. We summarize our main results, as well
a future outlook in 6. Details of the numerical simulation, as well as some
details of our analytic calculations are deferred to the Appendix.

##  2 Preliminaries

###  2.1 Model and equations of motion

A (dark) massive spin-$1$ field $W_{\mu}$ minimally coupled to gravity and
without non-gravitational self-interactions is described by the following
action:

$\displaystyle S$
$\displaystyle=\int\mathrm{d}^{4}x\,\sqrt{-g}\Bigl{[}-\frac{1}{4}g^{\mu\alpha}g^{\nu\beta}\,\mathcal{G}_{\mu\nu}\mathcal{G}_{\alpha\beta}+\frac{1}{2}\frac{m^{2}c^{2}}{\hbar^{2}}\,g^{\mu\nu}W_{\mu}W_{\nu}+\frac{c^{3}}{16\pi
G}R+...\Bigr{]},$ (2.1)

where $\mathcal{G}_{\mu\nu}=\partial_{\mu}W_{\nu}-\partial_{\nu}W_{\mu}$. The
‘$...$’ in (2.1) represents the Standard Model Lagrangian and other possible
dark sector(s). Here, $m$ is the mass of the vector boson. We can represent
the spatial part of the (real-valued) vector field ${\bm{W}}$ in terms of a
complex vector ${\bm{\mathsf{\Psi}}}$ as

$\displaystyle{\bm{W}}(t,\bm{x})\equiv\frac{\hbar}{\sqrt{2mc}}\Re\left[{\bm{\mathsf{\Psi}}}(t,\bm{x})e^{-imc^{2}t/\hbar}\right],$
(2.2)

where ${\bm{\mathsf{\Psi}}}$ has dimensions of $[\textrm{length}]^{-3/2}$.
Similarly,
$W_{0}(t,\bm{x})\equiv\hbar/\sqrt{2mc}\,\Re\left[\psi_{0}(t,\bm{x})e^{-imc^{2}t/\hbar}\right]$.
We are interested in the non-relativistic behaviour of the vector field where
the spatial variation in the field is slow compared to the Compton scale
$\lambda_{m}=\hbar/mc$ and we are in the Newtonian gravity regime. We focus on
sufficiently subhorizon dynamics, and hence ignore Hubble expansion. In this
case, the dynamics are described by the non-relativistic action for the
complex vector field ${\bm{\mathsf{\Psi}}}$ and the Newtonian gravitational
potential $\Phi$:

$\displaystyle\mathcal{S}_{nr}$
$\displaystyle=\int\mathrm{d}t\mathrm{d}^{3}x\,\Biggl{[}\frac{i\hbar}{2}{\bm{\mathsf{\Psi}}}^{\dagger}\dot{{\bm{\mathsf{\Psi}}}}+\mathrm{c.c.}-\frac{\hbar^{2}}{2m}\nabla{\bm{\mathsf{\Psi}}}^{\dagger}\cdot\nabla{\bm{\mathsf{\Psi}}}+\frac{1}{8\pi
G}\,\Phi\nabla^{2}\Phi-m\,\Phi\,{\bm{\mathsf{\Psi}}}^{\dagger}{\bm{\mathsf{\Psi}}}\Biggr{]},$
(2.3)

and the corresponding multi-component Schrödinger-Poisson (SP) system of
equations of motion:222To include the effects of Hubble expansion, simply
replace $\nabla\rightarrow\nabla/a$ and
$\partial_{t}\rightarrow\partial_{t}+3H/2$ where $a$ is the scalefactor and
$H=\dot{a}/a$.

$\displaystyle i\hbar\frac{\partial}{\partial t}{\bm{\mathsf{\Psi}}}$
$\displaystyle=-\frac{\hbar^{2}}{2m}\nabla^{2}{\bm{\mathsf{\Psi}}}+m\,\Phi\,{\bm{\mathsf{\Psi}}}\,,\qquad\qquad\nabla^{2}\Phi=4\pi
Gm\,{\bm{\mathsf{\Psi}}}^{\dagger}{\bm{\mathsf{\Psi}}}.$ (2.4)

This is our master equation that we work with throughout this work. We re-
iterate that ${\bm{\mathsf{\Psi}}}$ a complex 3-tuple with components
$[{\bm{\mathsf{\Psi}}}]_{i}=\psi_{i}$ with $i=1,2,3$ and
${\bm{\mathsf{\Psi}}}^{\dagger}{\bm{\mathsf{\Psi}}}=\sum_{i=1}^{3}|\psi_{i}|^{2}$.
For scalar dark matter, we have a single component field (which leads to the
“usual” Schrödinger-Poisson system). The nonrelativistic equations above were
derived earlier in [27, 28]. For a generalization to the spin-$s$ case, also
see [28]. In the future, it might also be interesting to systematically
explore the relativistic corrections to this multicomponent system in the
nonlinear regime [34].

####  2.1.1 Conserved Quantities

Note that in our convention, the number density, mass density, and spin
density are

$\displaystyle\mathcal{N}(t,\bm{x})={\bm{\mathsf{\Psi}}}^{\dagger}{\bm{\mathsf{\Psi}}},\qquad\rho(t,\bm{x})=m{\bm{\mathsf{\Psi}}}^{\dagger}{\bm{\mathsf{\Psi}}},\quad{\textrm{and}}\quad\bm{s}=i\hbar{\bm{\mathsf{\Psi}}}\times{\bm{\mathsf{\Psi}}}^{\dagger}.$
(2.5)

The conserved quantities associated with our non-relativistic VDM are:

$\displaystyle N$
$\displaystyle=\int\mathrm{d}^{3}x{\bm{\mathsf{\Psi}}}^{\dagger}{\bm{\mathsf{\Psi}}},\quad\textrm{and}\quad
M=mN,\qquad\qquad\qquad\qquad\,\,\,\textrm{(particle number and rest mass)}$
(2.6) $\displaystyle E$
$\displaystyle=\int\mathrm{d}^{3}x\Bigl{[}\frac{\hbar^{2}}{2m}\,\nabla{\bm{\mathsf{\Psi}}}^{\dagger}\cdot\nabla{\bm{\mathsf{\Psi}}}-\frac{Gm^{2}}{2}{\bm{\mathsf{\Psi}}}^{\dagger}{\bm{\mathsf{\Psi}}}\int\frac{\mathrm{d}^{3}y}{4\pi|{\bm{x}}-{\bm{y}}|}{\bm{\mathsf{\Psi}}}^{\dagger}({\bm{y}}){\bm{\mathsf{\Psi}}}(\bm{y})\Bigr{]}\,,\qquad\qquad\quad\textrm{(energy)}$
(2.7) $\displaystyle{\bm{S}}$
$\displaystyle=\hbar\int\mathrm{d}^{3}x\,i{\bm{\mathsf{\Psi}}}\times{\bm{\mathsf{\Psi}}}^{\dagger}\,,\quad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\quad\,\textrm{(spin
angular momentum)}$ (2.8) $\displaystyle\bm{L}$
$\displaystyle=\hbar\int\mathrm{d}^{3}x\,\Re\left(i\,{\bm{\mathsf{\Psi}}}^{\dagger}\nabla{\bm{\mathsf{\Psi}}}\times\bm{x}\right).\,\,\qquad\qquad\qquad\qquad\qquad\qquad\textrm{(orbital
angular momentum)}$ (2.9)

Note that the spin and orbital angular momentum are separately conserved in
the non-relativistic system. Importantly, by definition, spin angular momentum
is identically zero for SDM (but not VDM). For details of the non-relativistic
action and conserved quantities for a general spin-$s$ bosonic field
(including VDM) see [28].

####  2.1.2 Fluid equations

We can also transform our multicomponent SP system eq. (2.4) into a set of
three, coupled fluid equations (following the Madelung transform commonly used
in SDM [35]). With the following field re-definition,
$\psi_{j}=\sqrt{\rho_{j}/m}\,e^{iS_{j}}$, and the velocity
$\bm{u}_{i}\equiv\hbar\nabla S_{i}/m$, we have

$\displaystyle\frac{\partial\rho_{j}}{\partial
t}+\nabla\cdot(\rho_{j}\bm{u}_{j})=0\,,\qquad\frac{\partial\bm{u}_{j}}{\partial
t}+(\bm{u}_{j}\cdot{\nabla})\bm{u}_{j}=\frac{1}{m}\nabla(Q_{j}-m\Phi),\quad\textrm{where}\quad
j=1,2,3$ (2.10)

where $Q_{j}=(\hbar^{2}/2m)\nabla^{2}\sqrt{\rho_{j}}/\sqrt{\rho_{j}}$. In
terms of the Madelung variables, the spin density
$s_{i}=i(\hbar/m^{2})\epsilon_{ijk}\sqrt{\rho_{j}\rho_{k}}e^{i(S_{j}-S_{k})}$.
The vorticity for each of the three fluids
$\bm{\omega}_{j}=\nabla\times\bm{u}_{j}=0$ if $\rho_{j}\neq 0$. Note that zero
vorticity does not imply zero spin density. If $\bm{\omega}_{i}\neq 0$ for
some fixed $i$ (with ${\bm{\omega}}_{j\neq i}=0$), then $\bm{s}=s_{i}\hat{i}$.

We numerically solve eq. (2.4), but the conservation/fluid equations can be
useful in gaining physical intuition for the behaviour of the system
(including for example, vortices [12] in three fluids.).

####  2.1.3 Scales and scalings

The SP system (Eq. (2.4)) has certain scaling symmetries, which significantly
increases the generality of the results. Specifically, if we have a solution
${\bm{\mathsf{\Psi}}}(t,\bm{x})$ for a system with total mass $M$ and vector
boson mass $m$, then
${\bm{\mathsf{\Psi}}}_{\beta\gamma}(t,\bm{x})=\beta^{5/2}\gamma^{2}{\bm{\mathsf{\Psi}}}(\gamma^{2}\beta^{3}t,\gamma\beta^{2}\bm{x})$
is a solution for a system with mass $\gamma M$ and $\beta m$. Moreover,

$\\{m,M\\}\rightarrow\\{\beta m,\gamma
M\\}\Longrightarrow\\{t,r,\rho,{\bm{S}}\\}\rightarrow\\{t/(\gamma^{2}\beta^{3}),r/(\gamma\beta^{2}),\gamma^{4}\beta^{6}\rho,(\gamma/\beta){\bm{S}}\\}.$
(2.11)

For ease of comparison with astrophysical scales, let us define

$\displaystyle m_{20}\equiv\frac{mc^{2}}{10^{-20}\,{\rm
eV}},\quad\mathcal{M}_{5}\equiv\frac{M}{2.3\times 10^{5}M_{\odot}}\,.$ (2.12)

The Compton length and time scales are then given by

$\displaystyle\lambda_{m}=\hbar/(mc)=6.4\times
10^{-7}\,\textrm{kpc}/m_{20},\quad\tau_{m}=\hbar/(mc^{2})=2.1\times
10^{-3}\,\textrm{yr}/m_{20}.$ (2.13)

###  2.2 Wave Interference

Consider the density resulting from the superposition of two unit amplitude
plane waves in a spin-$s$ field ($s=0$ for SDM and $s=1$ for VDM):
${\bm{\mathsf{\Psi}}}_{a}(\bm{x})=V^{-1/2}\,{\bm{\mathsf{\epsilon}}}^{(s)}_{a}e^{i\bm{k}_{a}\cdot\bm{x}}$,
where $V$ is the volume, $a=1,2$, and ${\bm{\mathsf{\epsilon}}}^{(s)}_{a}$ is
a unit complex vector:

$|{\bm{\mathsf{\Psi}}}_{a}(\bm{x})+{\bm{\mathsf{\Psi}}}_{b}(\bm{x})|^{2}=2\,V^{-1}\left(1+\Re\left[{\bm{\mathsf{\epsilon}}}^{(s)\dagger}_{a}\cdot{\bm{\mathsf{\epsilon}}}^{(s)}_{a}e^{-i(\bm{k}_{a}-\bm{k}_{b})\cdot\bm{x}}\right]\right)=2\,V^{-1}\left(1+{\rm
int}_{(s)})\right).$ (2.14)

Here, the factor of $2$ is the number of waves and ${\rm int}_{(s)}$ is the
interference term. Without loss of generality, we set $\bm{x}=0$. The
interference term is simply the cosine of the angle between the two waves
(unit vectors in $2$ dimensions for SDM and $6$ for VDM), with their heads
lying on a unit $4s+1$ sphere. Assuming a uniform distribution on the sphere,
the cosine of the angle between these waves ${\rm
int}_{(s)}\equiv\eta=\cos\theta$ has the following distribution
$p^{(s)}(\eta)=\pi^{-1/2}\left\\{\Gamma(2s+1)/\Gamma(2s+1/2)\right\\}(1-x^{2})^{(4s-1)/2}$.
While the mean is zero, the standard deviation

$\sqrt{\langle{\rm
int}_{(s)}^{2}\rangle}=\frac{1}{\sqrt{2(2s+1)}},\qquad\textrm{with}\qquad\frac{\sqrt{\langle{\rm
int}_{(1)}^{2}\rangle}}{\sqrt{\langle{\rm
int}_{(0)}^{2}\rangle}}=\frac{1}{\sqrt{3}}.$ (2.15)

That is, interference decreases for higher spin fields. This is a reflection
of the intuitive fact that in a larger component field, orthogonal components
do not interfere. This simple fact has important implications for differences
between VDM and SDM.

###  2.3 Solitons

Ground state solitons in VDM are characterized by the “chemical potential”
$\mu$ and a unit complex 3-vector ${\bm{\mathsf{\epsilon}}}$ (with
${\bm{\mathsf{\epsilon}}}^{\dagger}{\bm{\mathsf{\epsilon}}}=1$):

${{\bm{\mathsf{\Psi}}}}_{\rm sol}({t},{\bm{x}})={\psi}_{\rm
sol}({\mu},{r})e^{i{\mu c^{2}}{t}/\hbar}{\bm{\mathsf{\epsilon}}},$ (2.16)

where ${\psi}_{\rm sol}$ is a real valued spherically symmetric function that
satisfies

$-{\mu c^{2}}{\psi}_{\rm sol}=-\hbar^{2}\frac{{\nabla}^{2}}{2m}{\psi}_{\rm
sol}+m{\Phi}{\psi}_{\rm sol},\quad{\nabla}^{2}{\Phi}=4\pi Gm{\psi}_{\rm
sol}^{2}.$ (2.17)

Note that the profile for a VDM soliton satisfies the same time-independent
equation as a SDM soliton. The mass and spin of this soliton are given by

${M}_{\rm sol}\approx
60.7\frac{m_{\mathrm{pl}}^{2}}{m}\sqrt{\frac{\mu}{m}},\quad{{\bm{S}}}_{\rm
sol}\approx
i({\bm{\mathsf{\epsilon}}}\times{\bm{\mathsf{\epsilon}}}^{\dagger})60.7\frac{m_{\mathrm{pl}}^{2}}{m^{2}}\sqrt{\frac{\mu}{m}}\hbar.$
(2.18)

The special cases of maximally polarized solitons configurations are given by
${\bm{\mathsf{\epsilon}}}_{(0)}=\hat{z}$ and
${\bm{\mathsf{\epsilon}}}_{(1)}=(\hat{x}+i\hat{y})/\sqrt{2},$ along with their
spatial rotations [28]. Configurations with
${\bm{\mathsf{\epsilon}}}={\bm{\mathsf{\epsilon}}}_{(0)}$ are linearly
polarized, with zero total spin angular momentum [27]. Whereas, configurations
with ${\bm{\mathsf{\epsilon}}}={\bm{\mathsf{\epsilon}}}_{(1)}$ have a maximal
spin angular momentum $|{{\bm{S}}}_{\rm sol}|=\hbar{M}_{\rm sol}/m$. For all
other solitons, we expect the spin angular momentum to lie between these
maximal values. That is $0\leq|{{\bm{S}}}_{\rm sol}|\leq\hbar{M}_{\rm sol}/m$
[28]. Note that these solitons (with and without spin) have lower energies
than those with spherically symmetric vector fields vanishing at the origin
[36, 37], at least in the nonrelativistic limit.

In [8], the scalar soliton profile was parameterized by a characteristic width
$r_{c}$, so that the density and mass can be characterized as

${\rho}_{\rm sol}(r)\approx 1.9\times 10^{7}m_{22}^{-2}\frac{\left({\rm
kpc}/{r}_{c}\right)^{4}}{(1+0.091(r/{r}_{c})^{2})^{8}}\frac{M_{\odot}}{\textrm{kpc}^{3}},\qquad
M_{\rm sol}\approx 2.2\times
10^{8}\left(\frac{\textrm{kpc}}{r_{c}}\right)m_{22}^{-2}M_{\odot}.$ (2.19)

Using eq. (2.19) and eq. (2.17), we have $r_{c}=6.8\times
10^{-5}m_{22}^{-1}\sqrt{m/\mu}\,\textrm{kpc}$. For the solitons in our
simulations we typically have ${\mu/m}\sim 10^{-12}$.

We re-iterate that the soliton profile is characterized by the same function
in VDM and SDM. The analysis is identical with
${\bm{\mathsf{\epsilon}}}\rightarrow e^{i\varphi}$ (a phase) for SDM. Also
note that the vector solitons discussed above are gravitationally bound; it is
also possible to have vector solitons (vector oscillons [38]) which are bound
by attractive self-interactions. In such self-interaction supported solitons,
linearly polarized solitons have a lower energy compared to the circularly
polarized ones (at fixed particle number). As in the case of scalars (see for
example, [39, 40]), it will be possible to construct non-relativistic solitons
with a combination of self-interactions and gravity [41].

##  3 Two soliton mergers

In this section we explore the merger of 2 solitons, as a warm up to the
$N$-soliton case. For simplicity, we restrict ourselves to head-on collisions
only.

We begin with two identical VDM solitons with a characteristic radius
$r_{c}\approx 1\,{\rm kpc}\times(\mathcal{M}_{4}m_{20}^{2})^{-1}$ and
separation of $\approx 10{r}_{c}$, with masses $M_{\rm sol,1}=M_{\rm sol,2}$.
We define $\mathcal{M}_{4}=(M/2.2\times 10^{4}M_{\odot})$. We give each of
them a small $v/c\approx 3.3\times 10^{-7}\,\mathcal{M}_{4}m_{20}$ velocity
towards each other (the typical velocity expected in our $N$ soliton
simulations at this distance). Each VDM soliton has its own complex unit
vector ${\bm{\mathsf{\epsilon}}}_{1,2}$ (see eq. (2.16)). After the collision,
if a new soliton forms, it has a mass $M_{\rm sol,f}=f_{\rm v}(M_{\rm
sol,1}+M_{\rm sol,2})$, where $1-f_{\rm v}$ is mass fraction that does not end
up in the final soliton. Here, we are partly repeating similar analyses
carried out for solitons in SDM [42, 43]. We note that in our analysis, the
merged configuration is often highly excited, and it is not always clear that
it is a settled soliton within the duration of the simulation. Nevertheless,
we will continue to refer to the merged object as a soliton for simplicity in
this section.333Since the merged object is not spherically symmetric, and can
have large excitations, we take the non-spherical density profile at a given
instant and rotate and average it (6 times) to create an approximately
spherically symmetric profile. This procedure is repeated for many time
slices, and time-averaged to get a new radial profile. This angular and time
averaged profile is then fit to a soliton profile to define the width $r_{c}$.
The merged soliton mass is calculated from this $r_{c}$. The ratio of this
final soliton mass to the total initial mass is what we define as $f$.

Along with the mass loss fraction, we also investigate the evolution of
intrinsic spin for the VDM case. As with mass, while the total spin angular
momentum is conserved, the final soliton need not carry all of the initial
spin angular momentum. This naturally generates the possibility of creating a
merged object with intrinsic spin angular momentum, even if the two initial
solitons had none.

We carried out about $\sim 10$ runs for the case with
${\bm{\mathsf{\epsilon}}}_{1}={\bm{\mathsf{\epsilon}}}_{2}e^{i\varphi}$, and
$10$ runs for
${\bm{\mathsf{\epsilon}}}_{1}\neq{\bm{\mathsf{\epsilon}}}_{2}e^{i\varphi}$.
Our simulations have periodic boundary conditions, so instead of calculating
the mass that escapes, we focus on the merged object. The merged configuration
typically shows large oscillations, which we average over and fit with a
single parameter profile (2.19), to obtain a value for $r_{c}$ and $M_{\rm
sol,f}$. We also kept track of the intrinsic spin angular momentum, which is
calculated using eq. (2.6) in the region of interest. Our key findings are as
follows:

  1. 1.

If ${\bm{\mathsf{\epsilon}}}_{1}={\bm{\mathsf{\epsilon}}}_{2}e^{i\varphi}$,
the collision is identical to the collision of two solitons in SDM with a
phase difference $\varphi$ between them (confirming the expectation in [28]).
In this case the solitons merge to form another soliton $0\leq\varphi<\pi$
with $f_{\rm v}=f_{\rm s}\approx 0.61\pm 0.01$.

  2. 2.

If ${\bm{\mathsf{\epsilon}}}_{1}\neq{\bm{\mathsf{\epsilon}}}_{2}e^{i\varphi}$,
then in general, such a collision cannot be mimicked by binary soliton
collisions in SDM [28]. Restricting ourselves to the cases where
${\bm{\mathsf{\epsilon}}}^{\dagger}_{1}\cdot{\bm{\mathsf{\epsilon}}}_{2}>0$
(so that the solitons can merge relatively quickly due to constructive
interference), we find that $f_{\rm v}\approx 0.56\pm 0.03$.

  3. 3.

To investigate the evolution of spin angular momentum, we consider the case
where ${\bm{\mathsf{\epsilon}}}_{1}=\hat{x}$ and
${\bm{\mathsf{\epsilon}}}_{2}=i\hat{y}$. In this case, both initial solitons
have zero spin angular momentum (they are linearly polarized). At the end of
the merger, we find that the final core has a spin angular momentum
${\bm{S}}_{\rm c,f}\propto\hbar\hat{z}$. The field which is not in the core
carries an equal amount of (opposite signed) spin angular momentum. We note
that the core is not quite a soliton for this particular case; density of
different components are not spherically symmetric around a common
origin.444If instead of ${\bm{\mathsf{\epsilon}}}_{2}=i\hat{y}$ we had chosen
${\bm{\mathsf{\epsilon}}}_{2}=\hat{y}$, the dynamics of the fields would
remain identical. However, the resulting spin density would not. More
generally, if we have a solution ${\bm{\mathsf{\Psi}}}$ of our system, then
$U{\bm{\mathsf{\Psi}}}$ is also a solution if $U^{\dagger}U=1$. Moreover the
number density $\mathcal{N}$ remains the same. However, the spin density
$\bm{s}$ does not remain the same in the general case. How does the spin
density change? For $U=\mathcal{O}$, with $\mathcal{O}^{T}\mathcal{O}=1$, we
have $\bm{s}\rightarrow\mathcal{O}\bm{s}$ as expected. However for
$U_{ij}\propto\delta_{ij}$, with $U_{jj}=ie^{\alpha_{j}}$, we have
$i{\bm{\mathsf{\Psi}}}\times{\bm{\mathsf{\Psi}}}^{\dagger}\rightarrow
2\Im[\epsilon_{ijk}e^{i\alpha_{ij}}\psi_{i}\psi_{j}^{*}]\neq U\bm{s}$ where
$\alpha_{ij}\equiv\alpha_{i}-\alpha_{j}$.

This is a preliminary investigation of 2 soliton mergers in VDM, and the
dependence on impact parameters, angular momentum, initial spin, energy etc.
needs to be investigated further. We note that that our $f_{\rm s}\approx
0.61$ is less than $\approx 0.7$ quoted in the literature [42], which could be
due to different ways in which the mass loss fraction is calculated as well as
the initial relative velocities used.

Figure 1: Starting with a collection of idealized halos (solitons), we
eventually evolve to an an approximately spherically symmetric configuration
with a central core surrounded by a halo. Top row is vector dark matter (VDM),
bottom row is scalar dark matter (SDM). For identical initial conditions in
density, the VDM final central core is less dense and its halo shows less
interference, as compared to SDM. The core to halo transition is also smoother
in VDM compared to SDM. In the above images, the color represents the
projected mass density in the simulation volume. Lighter colors correspond to
higher mass densities.

##  4 Many soliton mergers

We begin with $N\sim\mathcal{O}[10]$ solitons whose positions are chosen
randomly within our simulation volume. As we let the system evolve,
gravitational interactions bring the solitons closer. Field interference and
nonlinear evolution leads to a complex transient phase, after which, the
density settles into an approximately spherically symmetric density
configuration. The typical time-scale of this transient phase is less than the
dynamical time scale $t_{\rm dyn}=1/\sqrt{G\bar{\rho}}$ of our systen. When
comparing SDM and VDM, the initial density at each point in the simulation
volume is always identical. Snapshots of the time evolution of SDM and VDM are
shown in Fig. 1.

We consider the case where all the solitons have the same initial radii, as
well as the case where we draw the radii from a Gaussian distribution. We also
consider the case where we change the total mass $M_{\rm tot}$ while fixing
the number of solitons, as well as the case where we fix the number of
solitons and change the total mass. For SDM, the initial phase for each
soliton is drawn from a uniform distribution between $0$ and $2\pi$. For VDM
solitons, a complex unit vector ${\bm{\mathsf{\epsilon}}}$ is a $6$
dimensional unit vector with its head uniformly distribution on a unit
5-sphere.

Figure 2: The ratio of the final core mass to the total mass is plotted
against the invariant quantity $\Xi$ for $\sim 80$ SDM and $\sim 80$ VDM
configurations. The numerical power law fits, $M_{\rm core}/M_{\rm
tot}=a\,\Xi^{\alpha}$, are also shown in the above panel. The fitted exponent
$\alpha$ (with errors $\delta\alpha\\}\lesssim 0.01$ obtained from the
numerics are in agreement with expectations of simple analytic arguments
related to repeated binary mergers. The bottom plot shows the deviation of
$M_{\rm core}/M_{\rm tot}$ from the best fit lines. The average scatter,
standard deviation of $\Delta M_{\rm core}/M_{\rm tot}$ divided by the fit, in
VDM is $\sim 15\%$ larger than SDM.

(a)

(b)

Figure 3: Left panel (3a): Angle-averaged, late-time central core+halo
profiles for $\sim 80+80$ simulations spanning a range of initial conditions
including different total mass, initial number of solitons, locations of
solitons, phases and spins of solitons (i.e. $\Xi$ spans an order of
magnitude). The radial coordinate and density are normalized by $r_{c}$ and
$\rho(r=0)$ to highlight the differences in profile shape of VDM and SDM
coalesced cores independent of the initial conditions. Solid lines indicate
average over different simulations, the shaded regions indicate the spread in
all profiles. A marker at $r/r_{c}\approx 3.5$ shows a general transition
between core/halo regions in both SDM and VDM scenarios. Right panel (3b):
Final radial density from 10 simulations (time averaged over roughly 1 period
of radial oscillations of the core), where the initial mass is narrowly
distributed around $M_{\rm tot}=2.3\times
10^{5}\;M_{\odot}\times\mathcal{M}_{5}$, the size of the simulation volume is
$L=100\,{\rm{kpc}}\times(\mathcal{M}_{5}m_{20}^{2})^{-1}$ and the number of
initial solitons was fixed at 21.

###  4.1 Core-halo mass

In our simulations, we find a tight correlation between the mass of the final
core and the total mass of the system: $M_{\rm core}/M_{\rm
tot}\propto\Xi^{\alpha}$, where $\Xi$ is a measure of the total energy of the
system, and $\alpha$ is different for VDM and SDM. This correlation is shown
in Fig. 2. Below, we present an explanation for the observed
relationship.555Note that from simulations, $M_{\rm core}$ is obtained by
first fitting the central density region with (2.19) to obtain an $r_{c}$, and
then including mass withing a sphere of radius $2r_{c}$. For a soliton,
$M_{\rm core}=M(r<2r_{c})=(3/4)M_{\rm sol}$. Also note that $r_{c}$ for the
core is not fixed in time, and towards the end of the simulations shows
oscillatory time variations of roughly $20\%$.

Beginning with $N$ solitons of mass $M_{\rm sol}^{\rm i}$ each, and
distributed randomly throughout the box, the total energy is (scaled to yield
a dimensionless scale-invariant measure $\Xi$)

$\displaystyle\Xi\equiv\frac{|E_{\rm tot}|}{M_{\rm tot}^{3}(Gm/\hbar)^{2}}$
$\displaystyle\approx\frac{1}{M_{\rm
tot}^{3}(Gm/\hbar)^{2}}\left[N\frac{G(M_{\rm sol}^{\rm i})^{2}}{2R_{\rm
sol}^{\rm i}}+(1.88)N(N-1)\frac{G(M_{\rm sol}^{\rm i})^{2}}{L}\right],$ (4.1)
$\displaystyle\approx\frac{1}{20N^{2}}.$ (4.2)

In the first line, $L$ is the box size and $R_{\rm sol}^{\rm i}\ll L$ is the
initial solitons’ radius. In the last equality, we have assumed that the first
term in eq. (4.1) dominates over the second.666Note that $R_{\rm sol}^{\rm
i}\equiv 9.95\hbar^{2}/(GM_{\rm sol}^{\rm i}m^{2})$ contains $99\%$ of the
soliton’s mass, and we also include gradient contributions to the individual
soliton’s energy. The factor of $1.88$ in the second term arises from an
ensemble average over initial locations of the solitons. We have ignored
corrections due to overlap of solitons. For our system, the second term is
only marginally smaller than the first.

To relate the final soliton mass to the initial one, we follow a simplified
version of the arguments used in [42]. See our Appendix A.2 for details.
Suppose that for $N$ initial solitons, $N_{\rm major}\equiv\gamma N$ are major
binary mergers, and in each major merger a fraction $f$ of the progenitors
mass is contained in the resulting merged soliton. In this case, the final
soliton mass after $\gamma N$ major mergers is given by $M_{\rm sol}^{\rm
f}\approx(\gamma N)^{\log_{2}(2f)}{M_{\rm sol}^{\rm i}}=(\gamma
N)^{\log_{2}(f)}M_{\rm tot}$. Solving for $N$ and using this result in eq.
(4.1) we arrive at

$\displaystyle\frac{M_{\rm core}}{M_{\rm tot}}\approx
a\,\Xi^{-\log_{2}(\sqrt{f})}\qquad\textrm{where}\qquad
a=\frac{3}{4}\frac{\gamma}{f}\left(\frac{\gamma}{\sqrt{5}}\right)^{\log_{2}(f)}\,.$
(4.3)

Here, we used $M_{\rm core}=M_{\rm sol}(r<2r_{c})=(3/4)M_{\rm sol}^{\rm f}$.
The fraction $f$ depends on whether we have VDM or SDM. From the previous
section on two-soliton collisions we found that for VDM, $f_{\rm v}\approx
0.56\pm 0.03$, and for SDM, $f_{\rm s}\approx 0.61\pm 0.01$. These $f$s
translate to the following exponents of $\Xi$ for VDM and SDM:

$\displaystyle-\log_{2}(\sqrt{f_{\rm v}})\approx 0.42\pm
0.04\qquad-\log_{2}(\sqrt{f_{\rm s}})\approx 0.36\pm 0.02,$ (4.4)

which match the numerical obtained exponents well (see Fig. 2). We have not
calculated the fraction of major mergers $\gamma$, so we do not try to
estimate the coefficient $a$. Nevertheless, if we assume that $\gamma$ is the
same for SDM and VDM (with $1/4\lesssim\gamma\lesssim 1$), then the analytic
estimate captures the result that $M_{\rm core}/M_{\rm tot}$ for VDM should be
lower than SDM for a given $\Xi$ (in the range shown in the figure), and that
the exponent of $\Xi$ for VDM should be higher than that for SDM.

We expect that the tight correlation between core and halo mass in VDM will be
disrupted for cosmological simulations (as opposed to our idealized halo
mergers), as was the case for SDM [44, 45].

###  4.2 Late-time density profiles

We will focus on the (angular averaged) radial density profile of the late
time configuration in our simulation volume. We find that the central region
is well described by a soliton-like core, which eventually transitions into a
power law tail at larger radii. Starting with identical initial density
distributions in VDM and SDM, we typically find that the central solitonic
core is less dense for VDM.

Fitting for the cores with a soliton profile (see (2.19)), we extract a core
radius $r_{c}$ for each run. For a soliton, this $r_{c}$ uniquely fixes the
central density $\rho_{c}\propto r_{c}^{-4}$. For ease of comparison between
simulations with different initial conditions (which include different (total)
initial masses, number of solitons, distribution of radii, locations of
solitons etc.) we scale the density and radius of the final configuration by
$\rho_{c}$ and $r_{c}$. The resulting collection of density profiles are shown
in Fig. 3a.

The key observations common to VDM and SDM based on these normalized density
profiles are as follows:

  * •

For both VDM and SDM, a soliton-like core is clearly visible for
$r/r_{c}\lesssim 1$. At $r/r_{c}\gtrsim 1$, the profile starts dropping
rapidly. A transition from core to an $r^{-3}$ tail, occurs between $1\lesssim
r/r_{c}\lesssim 10$. For $r/r_{c}\gtrsim 10$, we see
$\rho/\rho_{c}\propto(r/r_{c})^{-3}$.777The periodic box makes it difficult to
trust the detailed power law when the radii become comparable to the size of
the box. The $r^{-3}$ is not robust at radii comparable to the size of the
simulation volume. The transition region is qualitatively delineated by
$r/r_{c}=3.5$ in Fig. 3a.

The key distinguishing feature between VDM and SDM is that

  * •

the transition from the soliton-like profile to $r^{-3}$ profile occurs a lot
more smoothly in VDM compared to SDM. This shape information is relatively
independent of our initial conditions.

The power law regime joins the soliton profile for $r/r_{c}\sim 1$ in case of
VDM. For SDM, the soliton-like profile persists for $r/r_{c}\gtrsim 1$, after
which there is a transitory power law (shallower than $r^{-3}$), before
joining the $r^{-3}$ tail at $r/r_{c}\sim 10$.

(a)

(b)

Figure 4: The normalized histogram of density (probability density function)
from a SDM and a VDM simulation, each with $M_{\rm tot}=2.2\times
10^{5}\;M_{\odot}\times\mathcal{M}_{5}$. At low densities compared to the mean
($\rho/\bar{\rho}\ll 1$), we see important qualitative differences between SDM
and VDM (4a) (left panel), with a dearth of ultra-low density regions in VDM
due to reduced interference. In the right panel, (4b), we see a lack of ultra-
high densities in VDM, again due to reduced interference. The general PDF
shapes of SDM and VDM densities, especially at low densities, are found to be
robust for a wide range of initial conditions.

###  4.3 Density PDF

As discussed earlier, interference effects are reduced in VDM compared to SDM.
As a result, we expect fewer deviations from the average density at low and
very high densities for VDM compared to SDM. This is seen in the PDF of
density extracted from our simulations (see Fig. 4a and Fig. 4b). The
horizontal axis is normalized by the average density in the box.

To gain further (still qualitative) understanding of the low-density shape
difference (near $\rho/\bar{\rho}\rightarrow 0$) between VDM and SDM, we
consider the density resulting from the superposition of a large number of
unit amplitude waves with random phases. In that case the density pdf for a
spin-$s$ field is a Gamma distribution
$\mathcal{P}^{(s)}(\rho/\bar{\rho})\sim{(2s+1)^{(2s+1)}}/{\Gamma[2s+1]}\left({\rho}/{\bar{\rho}}\right)^{2s}e^{-(2s+1)\rho/\bar{\rho}}$.
In particular, for $s=0$ (SDM) and $s=1$(VDM) and for $\rho/\bar{\rho}\ll 1$,
we have

$\mathcal{P}^{(0)}(\rho/\bar{\rho})\sim
1-\frac{\rho}{\bar{\rho}}+\ldots,\qquad\mathcal{P}^{(1)}(\rho/\bar{\rho})\sim\left(\frac{\rho}{\bar{\rho}}\right)^{\\!2}+\ldots$
(4.5)

which explains the qualitative behavior of the pdf at low densities in VDM and
SDM seen Fig. 4a. We caution that assuming unit amplitudes for all waves is
not justified, and we do not expect the numerical pdf to agree with our
analysis above quantitatively. Nevertheless, this qualitative understanding
and shape will already be useful for observational implications (see Section
5).

At high densities (Fig. 4b), we see that VDM has a shorter tail compared to
SDM which is again a direct consequence of reduced interference in VDM. The
extreme high density regions are found in the central core for both VDM and
SDM.

(a)

(b)

Figure 5: In Fig. 5a, we show the two point correlation function of fractional
density fluctuations in the halo (c.f. (4.6)) for a single simulation with
$M_{\rm tot}\approx 2.2\times 10^{5}\,M_{\odot}\times\mathcal{M}_{5}$ and
$N=21$ at late times. The central core is removed before calculating the
correlation function. The lower amplitude of the VDM curve reflects the
smaller variance of the density fluctuations, whereas the similarity in width
of the correlation functions (inset) shows that the physical size of the
interference granules is similar in SDM and VDM. In Fig. 5b, we show the
potential energy in the system as a function of time. The potential energy
(and kinetic energy) in SDM shows more temporal fluctuations compared to VDM.
The kinetic energy can be obtained from $K/|E_{\mathrm{tot}}|=1-U/|E_{\rm
tot}|$.

###  4.4 Two point correlation functions

To understand the typical amplitude and characteristic scale of density
fluctuations in the halo surrounding the central core, we construct the two
point correlation function

$\displaystyle\xi(\bm{l})=V^{-1}\int_{V}d^{3}r\,\delta(\bm{r})\delta(\bm{r}+\bm{l})\,,$
(4.6)

where we have removed the central high density core to reduce the influence of
the central core, so that $V$ is the volume enclosed between $r\sim 8r_{c}$
and $r\sim 20r_{c}$ from the center of the core. Here,
$\delta(\bm{r})=(\delta\rho/\rho)(\bm{r})\equiv[\rho(\bm{r})-\bar{\rho}(r)]/\bar{\rho}(r)$
where $\bar{\rho}(r)$ is the angular averaged density. Fig. 5a shows the
angular averaged $\xi(\bm{l})$.888This figure has changed from our previous
version where we used the total density rather than the fractional deviation
of the density in the correlation function. The VDM correlation function at
$l=0$ is a measure of the variance of $\delta(\bm{r})$. We find $[\xi(0)]_{\rm
vdm}/[\xi(0)]_{\rm sdm}\sim 1/3$ in agreement with our analysis in Section
2.2, indicating a reduced interference in VDM compared to SDM. When the
correlation functions in SDM and VDM are scaled by the $[\xi(0)]_{\rm vdm}$
and $[\xi(0)]_{\rm sdm}$ respectively, we find that they overlap, and have the
same width. We take this to be an indication that the spatial size of the
interference granules in SDM and VDM are the same in the outer region of the
halo, $r\gtrsim 10r_{c}$, where the mass enclosed as a function of radius
becomes identical in SDM and VDM.999Using the radial profile of mass density,
we find that the circular velocity, and hence approximate de Broglie scale in
the outer regions of the halo is $\lambda_{\rm
dB}\sim(\mathcal{M}_{5}^{2}m_{20}^{2})^{-1}1{\rm\,kpc}$.

Another measure of the density structure in VDM and SDM is the power spectrum,
which is provided in A.4.

Figure 6: The left panel shows the magnitude of the spin density at the end of
the simulation, whereas the zoomed inset shows the vector spin density in the
central core. The time-averaged spin vector per boson in the core (within
$2r_{c}$) is shown in the middle, along with its typical precession around the
mean over a de-Broglie time scale. We take the smallness of the variation to
be a sign that we have a soliton-like core in this cases, however, there are
can be larger variations in magnitude and direction of spin in other cases.
The top panel of the rightmost plot shows a correlation between the initial
spin per boson in our simulation (which is conserved) and the final spin per
boson in the core. The red points are ensemble mean of the magnitude of the
time-averaged vector spin in the core, where the ensemble consists of similar
initial spin/boson simulations. The error bars show a 90% confidence interval
within this ensemble. The bottom panel shows the ensemble mean and standard
deviation of the precession of the core spin. We caution that there might be a
core, but not necessarily a soliton present at the centre in some of the
cases. Note that a significant spin density in the core can be generated even
at small initial values of the total spin.

###  4.5 Spin distribution

We are interested in the intrinsic spin of the final configuration. To better
understand how this intrinsic spin can be generated, we consider two classes
of initial conditions. In one case, we start with initial conditions where
each soliton has zero intrinsic spin, with ${\bm{\mathsf{\epsilon}}}_{i}$ for
each soliton a real random unit vector. In another class, each soliton has a
random ${\bm{\mathsf{\epsilon}}}_{i}$ complex vector (with unit norm), which
leads to some residual non-zero intrinsic spin in the simulation volume (even
though the spin directions for different solitons are not aligned).

In Fig. 6, we show the final result of one such simulation (from the latter
set), as well as the correlation between the initial spin in the box and the
final spin in the core. At late times, we find that the final central core is
fractionally polarized (with a mixture of circular and linear polarization)
with a non-zero intrinsic spin. We notice a strong correlation between the
initial spin in the box and the final spin of the core. For low initial spins,
we find that the final spin in the core has a large variance as a function of
time even at late times. This might be an indication that there is no settled
soliton in the core for very low initial spins. For higher spins, there is a
distinct direction of spin in the core, with relatively small variations with
time (an example of which can be seen in Fig. 6). Wave interference effects in
SDM create oscillations in the core radius, and induce a random walk of the
core [46, 47]. A similar analysis in VDM could shed light on the time
dependence of spin in the core.

Since intrinsic spin for the entire simulation is conserved (separately from
the orbital angular momentum) in the non-relativistic limit, the halo carries
the rest of the spin (with an opposite sign).

##  5 Observational Implications

We discuss three application areas of VDM to astrophysical observations of
interest, for possible future study.101010For terrestrial experiments
searching for light VDM, see [48] and references therein. Also see [49] for
the relevance of stochastic gradients in axion fields on detection strategies.
Related analysis can be repeated for the VDM case taking into account both
interference patterns and solitons.

###  5.1 Dark matter substructure and dynamical heating

Density fluctuations resulting from wave interference in ultralight dark
matter can dynamically heat the old stellar population in the Milky Way,
thickening the scale height of its disk-like distribution [50]. The heating
rate $\mathcal{H}\propto(\delta\rho/\rho)^{2}\lambda_{\rm
dB}^{3}\propto(\delta\rho/\rho)^{2}(\hbar/m)^{3}$ where
$(\delta\rho/\rho)^{2}$ is the variance of fractional density fluctuations
around the smooth density profile, and $\lambda_{\rm dB}\sim\hbar/mv$ is the
characteristic size of the interference granules.111111The heating rate:
$\mathcal{H}\sim v^{-1}{(GM_{\rm gran})^{2}}n_{\rm gran}\ln\Lambda$, where
$M_{\rm gran}\sim\delta\rho\lambda_{\rm dB}^{3}$ is the mass of a typical
interference granule, $n_{\rm gran}\sim\lambda_{\rm dB}^{-3}$ is the number
density of granules and $\Lambda$ is the usual Coulomb logarithm. We are
essentially treating granules as scatterers which change the velocity of test
particles (stars). Using these approximations, we arrive at $\mathcal{H}\sim
v^{-1}{(G\rho)^{2}}\left({\delta\rho}/{\rho}\right)^{2}\lambda_{\rm
gran}^{3}\ln\Lambda$. The smooth density profile $\rho$ and characteristic $v$
are assumed to be constrained observationally (though additional care might be
needed here). For interference granules, $\delta\rho/\rho\sim 1$. See [50] and
[51] for more details. Using the scale height (and hence heating rate) as an
observable for our galaxy, the authors in [50] argued for an SDM boson mass of
$m\gtrsim 6\times 10^{-23}~{}{\rm eV}$.

VDM would similarly lead to a dynamical heating of the old stellar disk.
However, as we discussed in Section 2.2, interference effects are reduced in
VDM compared to SDM; $[(\delta\rho/\rho)^{2}]_{\rm
vdm}/[(\delta\rho/\rho)^{2}]_{\rm sdm}\sim 1/3$. Furthermore, we expect the de
Broglie scale in SDM and VDM to be the same for the same boson mass. These are
confirmed by Fig. 5a. In case the boson masses are different, the size of the
interference granules changes as $[\lambda_{\rm dB}]_{\rm vdm}/[\lambda_{\rm
dB}]_{\rm sdm}\sim m_{\rm{sdm}}/m_{\rm vdm}$, while we still expect the same
reduced interference in VDM compared to SDM. As a result, the constraint on
$m_{\rm{vdm}}$ is relaxed by $3^{1/3}$ compared to the SDM case, ie. $m_{\rm
vdm}>4.2\times 10^{-23}\,{\rm eV}$.121212We thank Neal Dalal for pointing out
an error in the previous version regarding the constraint on the vector boson
mass. A more careful analysis which takes into account the spectral
information of the interference patterns and granules in VDM compared to SDM
would be interesting to pursue. This analysis can be further generalized to
spin-$s$ bosonic fields where the
$[(\delta\rho/\rho)^{2}]_{s}/[(\delta\rho/\rho)^{2}]_{s=0}\sim 1/(2s+1)$.

A recent study has looked at constraining the SDM particle mass by considering
the effect of dynamical heating on the stellar kinematics in ultra-faint dwarf
galaxies [51], and found a tight constraint of $m\gtrsim 3\times
10^{-19}~{}{\rm eV}$. We expect the constraint to be similarly relaxed by
$\sim 3^{1/3}$ for VDM. Furthermore, as discussed in [51], the property of the
core might have a significant impact on this constraint (making it stronger).
It would be interesting to consider the difference in the central core
properties between VDM and SDM for constraining their respective masses.

###  5.2 Dark matter cores in dwarf galaxies

While standard CDM predicts signature $r^{-1}$ cuspy centers of dark matter
halos, at least without the presence of strong baryonic feedback,
observational evidence strongly points to large low-density cores in dwarf
galaxies (but see [52]). Ultralight DM naturally points to a cored DM halo,
and also predicts core sizes that are inversely proportional to mass. However,
[53, 54, 55] point out that a simple soliton profile given by its 1-parameter
family solutions and a fixed boson mass has difficulty matching the diversity
and scatter in observations of cores of halos with virial masses below
$10^{11}~{}M_{\odot}$, despite the fact that the soliton profile can be
matched well to dark-matter dominated dwarf spheroidal systems Sculptor and
Fornax [14]. An outstanding issue is that when examining a broad sample of low
mass galaxies, the inferred dark matter cores tend to have central densities
$\rho_{\rm core}$ scaling inversely with core radius, $r_{\rm core}^{-1}$ [56,
57], while a soliton profile has a scaling $\rho_{\rm core}\propto r_{\rm
core}^{-4}$ in both SDM and VDM. An investigation of whether VDM may naturally
introduce more scatter (and change in observed slope) into its core shape and
the core-halo mass relation, as our idealized simulations suggest they may, is
needed. In more realistic simulations, there could be additional scatter due
to tidal stripping as suggested by [45] for the SDM case. It is conceivable
that ultralight VDM might ameliorate the existing tensions with observational
data further.

###  5.3 DM Substructure and Lensing

Multiply imaged quasars at high redshift probe dark matter substructure in
halos; in CDM theory populations of $10^{9}~{}M_{\odot}$ subhalos can alter
the flux ratios of the images without significantly changing their positions
[58]. In VDM and SDM, while these theories predict a reduced population of
subhalos, they also feature wave interference that causes additional lensing.
The projected density and its power spectra can be used to estimate the size
of the effects of VDM and SDM (e.g. see plane wave perturbation analysis in
[59]). The reduced interference effects in VDM compared to SDM can change the
lensing predictions, potentially probing the nature of the underlying field
and changing constraints on the boson mass. Again, a more careful analysis is
warranted for lensing due to substructure in VDM, including interference
patterns and solitons.

##  6 Summary

Our goal in this paper was to understand whether the intrinsic spin of light
dark matter has an impact on its small scale structure assuming only
gravitational interactions. To this end, we have provided the first 3+1
dimensional simulation probing the small scale structure of vector dark matter
(VDM) and compared the results with those of scalar dark matter (SDM). See
Fig. 1. The key differences arise from the reduced interference of the fields
in VDM compared to SDM.

Similarities: Starting with an identical initial conditions in terms of their
mass density (in the form of $N$ solitons), we found that in both cases, we
form an approximately spherically symmetric density configuration with a
central core and a surrounding halo. The central core (after angular
averaging) can be fit well by a single parameter (core radius $r_{c}$) soliton
profile, and the surrounding halo eventually transitions to a power-law
profile beyond $r\gtrsim 10r_{c}$. The ratio of the mass of the central core
to the total mass can be expressed as a power law (with scatter, and different
for VDM and SDM) of the initial energy of the system, $M_{\rm core}/M_{\rm
tot}\propto\Xi^{\alpha}$).

Differences
Radial density profile: For identical initial conditions, the central core is
less dense and broader in VDM compared to SDM. If we normalize for the central
density (or core radius), there still remains a distinct difference in shape
of density distribution at the transition between the core and the halo: VDM
has a smoother transition. This shape information is independent of the
initial conditions. See Fig. 3.

Core-halo mass relation: For VDM, we were able to analytically estimate
$M_{\rm core}/M_{\rm tot}\propto\,\Xi^{\alpha}$ with
$\alpha=-\log_{2}(\sqrt{f})$. The difference in the exponent betweem SDM and
VDM can be attributed to different mass loss fractions, $f$, in individual
binary soliton mergers (obtained numerically). Furthermore, there is $\sim
15\%$ more scatter of $M_{\rm core}/M_{\rm tot}$ in VDM simulations compared
to SDM ones around the power-law fit. See Fig. 2.

Interference: The reduced interference in VDM compared to SDM leads to fewer
extreme (underdense and overdense) regions in VDM compared to SDM. Once the
region containing the core is removed, the amplitude of the two-point
correlation function of the fractional density contrast is smaller in VDM
compared to SDM, while widths of the correlation functions are the same.
Heuristically, the interference granules in VDM are of similar size, but less
overdense than those in SDM. See Fig. 4 and 5.

Spin: We initiated the study spin angular momentum during structure formation
in VDM. See Fig. 6. While total spin is conserved, spin density undergoes a
rich evolution. We found a positive correlation between the initial spin/boson
in the simulation, and the final spin of the central core. Even for a small
initial spin/boson, the final spin per boson in the core can be significant.
The detailed dynamics of spin needs to be explored further.

Near Future Directions: Much more remains to be done in the context of light,
higher spin DM. The growing number of analyses of observational implications
of ultralight SDM can be adapted for ultralight VDM, including for example
effects on dynamical friction [60, 61, 62], vortices [12], Ly$\alpha$ forest
[3, 63], CMB and galaxy surveys [64, 65]. Again, we expect the difference in
interference effects to hold the key to delineating effects of VDM and SDM.
For a more careful analysis of observational constraints, and for additional
handles on structure formation in VDM, the production mechanism for VDM is
relevant (e.g. [19, 23, 24, 66]). The production mechanisms can influence the
initial shape of the power spectrum. A linear/quasi-linear analysis could then
be used to provide initial conditions for a late time, large scale simulation
of VDM (comparable to [8, 67, 44] done for ultralight SDM). Finally, non-
gravitational (electromagnetic) effects of the potentially large intrinsic
spin of solitons and VDM cores will be worth exploring as a novel probe of the
spin of DM in an astrophysical setting [28].

Note: Our manuscript was submitted on the day [68] appeared on the arXiv. In
that paper, the authors investigate the formation of dark photon stars (vector
solitons) in the early universe and their implication for VDM substructure. We
will explore the connection between their results and ours in the future.

##  7 Acknowledgements

We would like to thank Neal Dalal for an insightful discussion on constraints
on the boson mass, as well Ray Hagimoto, Siyang Ling, Andrew Long and HongYi
Zhang for helpful conversations. MA and MJ are partly supported by a NASA
grant 80NSSC20K0518. The numerical simulations were carried out on the NOTS
cluster supported by the Center for Research Computing at Rice University.

## References

  * [1] Planck Collaboration, N. Aghanim et al., “Planck 2018 results. VI. Cosmological parameters,” Astron. Astrophys. 641 (2020) A6, arXiv:1807.06209 [astro-ph.CO].
  * [2] G. Bertone and D. Hooper, “History of dark matter,” Rev. Mod. Phys. 90 no. 4, (2018) 045002, arXiv:1605.04909 [astro-ph.CO].
  * [3] V. Iršič, M. Viel, M. G. Haehnelt, J. S. Bolton, and G. D. Becker, “First constraints on fuzzy dark matter from Lyman-$\alpha$ forest data and hydrodynamical simulations,” Phys. Rev. Lett. 119 no. 3, (2017) 031302, arXiv:1703.04683 [astro-ph.CO].
  * [4] D. M. Jacobs, G. D. Starkman, and B. W. Lynn, “Macro Dark Matter,” Mon. Not. Roy. Astron. Soc. 450 no. 4, (2015) 3418–3430, arXiv:1410.2236 [astro-ph.CO].
  * [5] L. Hui, “Wave Dark Matter,” Ann. Rev. Astron. Astrophys. 59 (2021) 247–289, arXiv:2101.11735 [astro-ph.CO].
  * [6] T. Matos and L. A. Urena-Lopez, “A Further analysis of a cosmological model of quintessence and scalar dark matter,” Phys. Rev. D 63 (2001) 063506, arXiv:astro-ph/0006024.
  * [7] W. Hu, R. Barkana, and A. Gruzinov, “Cold and fuzzy dark matter,” Phys. Rev. Lett. 85 (2000) 1158–1161, arXiv:astro-ph/0003365.
  * [8] H.-Y. Schive, T. Chiueh, and T. Broadhurst, “Cosmic Structure as the Quantum Interference of a Coherent Dark Wave,” Nature Phys. 10 (2014) 496–499, arXiv:1406.6586 [astro-ph.GA].
  * [9] M. A. Amin and P. Mocz, “Formation, gravitational clustering, and interactions of nonrelativistic solitons in an expanding universe,” Phys. Rev. D 100 no. 6, (2019) 063507, arXiv:1902.07261 [astro-ph.CO].
  * [10] A. Arvanitaki, S. Dimopoulos, M. Galanis, L. Lehner, J. O. Thompson, and K. Van Tilburg, “Large-misalignment mechanism for the formation of compact axion structures: Signatures from the QCD axion to fuzzy dark matter,” Phys. Rev. D 101 no. 8, (2020) 083014, arXiv:1909.11665 [astro-ph.CO].
  * [11] P. Mocz, M. Vogelsberger, V. H. Robles, J. Zavala, M. Boylan-Kolchin, A. Fialkov, and L. Hernquist, “Galaxy formation with BECDM – I. Turbulence and relaxation of idealized haloes,” Mon. Not. Roy. Astron. Soc. 471 no. 4, (2017) 4559–4570, arXiv:1705.05845 [astro-ph.CO].
  * [12] L. Hui, A. Joyce, M. J. Landry, and X. Li, “Vortices and waves in light dark matter,” JCAP 01 (2021) 011, arXiv:2004.01188 [astro-ph.CO].
  * [13] R. Brito, V. Cardoso, and P. Pani, “Superradiance: New Frontiers in Black Hole Physics,” Lect. Notes Phys. 906 (2015) pp.1–237, arXiv:1501.06570 [gr-qc].
  * [14] D. J. E. Marsh and A.-R. Pop, “Axion dark matter, solitons and the cusp–core problem,” Mon. Not. Roy. Astron. Soc. 451 no. 3, (2015) 2479–2492, arXiv:1502.03456 [astro-ph.CO].
  * [15] J. C. Niemeyer, “Small-scale structure of fuzzy and axion-like dark matter,” arXiv:1912.07064 [astro-ph.CO].
  * [16] E. G. M. Ferreira, “Ultra-light dark matter,” Astron. Astrophys. Rev. 29 no. 1, (2021) 7, arXiv:2005.03254 [astro-ph.CO].
  * [17] A. Arvanitaki, S. Dimopoulos, S. Dubovsky, N. Kaloper, and J. March-Russell, “String Axiverse,” Phys. Rev. D 81 (2010) 123530, arXiv:0905.4720 [hep-th].
  * [18] M. Cicoli, V. Guidetti, N. Righi, and A. Westphal, “Fuzzy Dark Matter Candidates from String Theory,” arXiv:2110.02964 [hep-th].
  * [19] P. W. Graham, J. Mardon, and S. Rajendran, “Vector Dark Matter from Inflationary Fluctuations,” Phys. Rev. D 93 no. 10, (2016) 103520, arXiv:1504.02102 [hep-ph].
  * [20] Y. Ema, K. Nakayama, and Y. Tang, “Production of purely gravitational dark matter: the case of fermion and vector boson,” JHEP 07 (2019) 060, arXiv:1903.10973 [hep-ph].
  * [21] E. W. Kolb and A. J. Long, “Completely dark photons from gravitational particle production during the inflationary era,” JHEP 03 (2021) 283, arXiv:2009.03828 [astro-ph.CO].
  * [22] A. Ahmed, B. Grzadkowski, and A. Socha, “Gravitational production of vector dark matter,” JHEP 08 (2020) 059, arXiv:2005.01766 [hep-ph].
  * [23] P. Agrawal, N. Kitajima, M. Reece, T. Sekiguchi, and F. Takahashi, “Relic Abundance of Dark Photon Dark Matter,” Phys. Lett. B 801 (2020) 135136, arXiv:1810.07188 [hep-ph].
  * [24] R. T. Co, A. Pierce, Z. Zhang, and Y. Zhao, “Dark Photon Dark Matter Produced by Axion Oscillations,” Phys. Rev. D 99 no. 7, (2019) 075002, arXiv:1810.07196 [hep-ph].
  * [25] J. A. Dror, K. Harigaya, and V. Narayan, “Parametric Resonance Production of Ultralight Vector Dark Matter,” Phys. Rev. D 99 no. 3, (2019) 035036, arXiv:1810.07195 [hep-ph].
  * [26] M. Bastero-Gil, J. Santiago, L. Ubaldi, and R. Vega-Morales, “Vector dark matter production at the end of inflation,” JCAP 04 (2019) 015, arXiv:1810.07208 [hep-ph].
  * [27] P. Adshead and K. D. Lozanov, “Self-gravitating Vector Dark Matter,” Phys. Rev. D 103 no. 10, (2021) 103501, arXiv:2101.07265 [gr-qc].
  * [28] M. Jain and M. A. Amin, “Polarized solitons in higher-spin wave dark matter,” Phys. Rev. D 105 no. 5, (2022) 056019, arXiv:2109.04892 [hep-th].
  * [29] N. Blinov, M. J. Dolan, P. Draper, and J. Shelton, “Dark Matter Microhalos From Simplified Models,” Phys. Rev. D 103 no. 10, (2021) 103514, arXiv:2102.05070 [astro-ph.CO].
  * [30] J. Eby, M. Leembruggen, L. Street, P. Suranyi, and L. C. R. Wijewardhana, “Galactic Condensates composed of Multiple Axion Species,” JCAP 10 (2020) 020, arXiv:2002.03022 [hep-ph].
  * [31] L. Street, N. Y. Gnedin, and L. C. R. Wijewardhana, “Testing multiflavored ULDM models with SPARC,” arXiv:2204.01871 [astro-ph.CO].
  * [32] H.-K. Guo, K. Sinha, C. Sun, J. Swaim, and D. Vagie, “Two-scalar Bose-Einstein condensates: from stars to galaxies,” JCAP 10 (2021) 028, arXiv:2010.15977 [astro-ph.CO].
  * [33] D. Cyncynates, T. Giurgica-Tiron, O. Simon, and J. O. Thompson, “Resonant nonlinear pairs in the axiverse and their late-time direct and astrophysical signatures,” Phys. Rev. D 105 no. 5, (2022) 055005, arXiv:2109.09755 [hep-ph].
  * [34] B. Salehian, H.-Y. Zhang, M. A. Amin, D. I. Kaiser, and M. H. Namjoo, “Beyond Schrödinger-Poisson: nonrelativistic effective field theory for scalar dark matter,” JHEP 09 (2021) 050, arXiv:2104.10128 [astro-ph.CO].
  * [35] E. Madelung, “Quantentheorie in hydrodynamischer Form,” Zeitschrift fur Physik 40 no. 3-4, (Mar., 1927) 322–326.
  * [36] R. Brito, V. Cardoso, C. A. R. Herdeiro, and E. Radu, “Proca stars: Gravitating Bose–Einstein condensates of massive spin 1 particles,” Phys. Lett. B 752 (2016) 291–295, arXiv:1508.05395 [gr-qc].
  * [37] F. Di Giovanni, N. Sanchis-Gual, C. A. R. Herdeiro, and J. A. Font, “Dynamical formation of Proca stars and quasistationary solitonic objects,” Phys. Rev. D 98 no. 6, (2018) 064044, arXiv:1803.04802 [gr-qc].
  * [38] H.-Y. Zhang, M. Jain, and M. A. Amin, “Polarized Vector Oscillons,” arXiv:2111.08700 [astro-ph.CO].
  * [39] P.-H. Chavanis, “Mass-radius relation of newtonian self-gravitating bose-einstein condensates with short-range interactions. i. analytical results,” Physical Review D 84 no. 4, (Aug, 2011) . http://dx.doi.org/10.1103/PhysRevD.84.043531.
  * [40] P. H. Chavanis and L. Delfini, “Mass-radius relation of Newtonian self-gravitating Bose-Einstein condensates with short-range interactions: II. Numerical results,” Phys. Rev. D 84 (2011) 043532, arXiv:1103.2054 [astro-ph.CO].
  * [41] M. Jain, “Yang-Mills stars in Higgsed non-abelian dark matter,” to appear.
  * [42] X. Du, C. Behrens, and J. C. Niemeyer, “Substructure of fuzzy dark matter haloes,” Mon. Not. Roy. Astron. Soc. 465 no. 1, (2017) 941–951, arXiv:1608.02575 [astro-ph.CO].
  * [43] M. P. Hertzberg, Y. Li, and E. D. Schiappacasse, “Merger of Dark Matter Axion Clumps and Resonant Photon Emission,” JCAP 07 (2020) 067, arXiv:2005.02405 [hep-ph].
  * [44] S. May and V. Springel, “Structure formation in large-volume cosmological simulations of fuzzy dark matter: impact of the non-linear dynamics,” Mon. Not. Roy. Astron. Soc. 506 no. 2, (2021) 2603–2618, arXiv:2101.01828 [astro-ph.CO].
  * [45] H. Y. J. Chan, E. G. M. Ferreira, S. May, K. Hayashi, and M. Chiba, “The diversity of core–halo structure in the fuzzy dark matter model,” Mon. Not. Roy. Astron. Soc. 511 no. 1, (2022) 943–952, arXiv:2110.11882 [astro-ph.CO].
  * [46] X. Li, L. Hui, and T. D. Yavetz, “Oscillations and Random Walk of the Soliton Core in a Fuzzy Dark Matter Halo,” Phys. Rev. D 103 no. 2, (2021) 023508, arXiv:2011.11416 [astro-ph.CO].
  * [47] J. L. Zagorac, I. Sands, N. Padmanabhan, and R. Easther, “Schrödinger-Poisson Solitons: Perturbation Theory,” arXiv:2109.01920 [astro-ph.CO].
  * [48] D. Antypas et al., “New Horizons: Scalar and Vector Ultralight Dark Matter,” arXiv:2203.14915 [hep-ex].
  * [49] M. Lisanti, M. Moschella, and W. Terrano, “Stochastic properties of ultralight scalar field gradients,” Phys. Rev. D 104 no. 5, (2021) 055037, arXiv:2107.10260 [astro-ph.CO].
  * [50] B. V. Church, P. Mocz, and J. P. Ostriker, “Heating of Milky Way disc stars by dark matter fluctuations in cold dark matter and fuzzy dark matter paradigms,” MNRAS 485 no. 2, (May, 2019) 2861–2876, arXiv:1809.04744 [astro-ph.GA].
  * [51] N. Dalal and A. Kravtsov, “Not so fuzzy: excluding FDM with sizes and stellar kinematics of ultra-faint dwarf galaxies,” arXiv e-prints (Mar., 2022) arXiv:2203.05750, arXiv:2203.05750 [astro-ph.CO].
  * [52] K. Hayashi, M. Chiba, and T. Ishiyama, “Diversity of Dark Matter Density Profiles in the Galactic Dwarf Spheroidal Satellites,” ApJ 904 no. 1, (Nov., 2020) 45, arXiv:2007.13780 [astro-ph.GA].
  * [53] H. Deng, M. P. Hertzberg, M. H. Namjoo, and A. Masoumi, “Can Light Dark Matter Solve the Core-Cusp Problem?,” Phys. Rev. D 98 no. 2, (2018) 023513, arXiv:1804.05921 [astro-ph.CO].
  * [54] A. Burkert, “Fuzzy Dark Matter and Dark Matter Halo Cores,” ApJ 904 no. 2, (Dec., 2020) 161, arXiv:2006.11111 [astro-ph.GA].
  * [55] M. Safarzadeh and D. N. Spergel, “Ultra-light Dark Matter Is Incompatible with the Milky Way’s Dwarf Satellites,” ApJ 893 no. 1, (Apr., 2020) 21, arXiv:1906.11848 [astro-ph.CO].
  * [56] F. Donato, G. Gentile, P. Salucci, C. F. Martins, M. I. Wilkinson, G. Gilmore, E. K. Grebel, A. Koch, and R. Wyse, “A constant dark matter halo surface density in galaxies,” Mon. Not. Roy. Astron. Soc. 397 (2009) 1169–1176, arXiv:0904.4054 [astro-ph.CO].
  * [57] E. V. Karukes and P. Salucci, “The universal rotation curve of dwarf disc galaxies,” Mon. Not. Roy. Astron. Soc. 465 no. 4, (2017) 4703–4722, arXiv:1609.06903 [astro-ph.GA].
  * [58] R. B. Metcalf and P. Madau, “Compound Gravitational Lensing as a Probe of Dark Matter Substructure within Galaxy Halos,” ApJ 563 no. 1, (Dec., 2001) 9–20, arXiv:astro-ph/0108224 [astro-ph].
  * [59] S. Mao and P. Schneider, “Evidence for substructure in lens galaxies?,” MNRAS 295 no. 3, (Apr., 1998) 587–594, arXiv:astro-ph/9707187 [astro-ph].
  * [60] L. Hui, J. P. Ostriker, S. Tremaine, and E. Witten, “Ultralight scalars as cosmological dark matter,” Phys. Rev. D 95 no. 4, (2017) 043541, arXiv:1610.08297 [astro-ph.CO].
  * [61] Y. Wang and R. Easther, “Dynamical friction from ultralight dark matter,” Phys. Rev. D 105 no. 6, (2022) 063523, arXiv:2110.03428 [gr-qc].
  * [62] D. Traykova, K. Clough, T. Helfer, E. Berti, P. G. Ferreira, and L. Hui, “Dynamical friction from scalar dark matter in the relativistic regime,” Phys. Rev. D 104 no. 10, (2021) 103014, arXiv:2106.08280 [gr-qc].
  * [63] M. Nori, R. Murgia, V. Iršič, M. Baldi, and M. Viel, “Lyman $\alpha$ forest and non-linear structure characterization in Fuzzy Dark Matter cosmologies,” Mon. Not. Roy. Astron. Soc. 482 no. 3, (2019) 3227–3243, arXiv:1809.09619 [astro-ph.CO].
  * [64] R. Hlozek, D. Grin, D. J. E. Marsh, and P. G. Ferreira, “A search for ultralight axions using precision cosmological data,” Phys. Rev. D 91 no. 10, (2015) 103512, arXiv:1410.2896 [astro-ph.CO].
  * [65] A. Laguë, J. R. Bond, R. Hložek, K. K. Rogers, D. J. E. Marsh, and D. Grin, “Constraining ultralight axions with galaxy surveys,” JCAP 01 no. 01, (2022) 049, arXiv:2104.07802 [astro-ph.CO].
  * [66] A. J. Long and L.-T. Wang, “Dark Photon Dark Matter from a Network of Cosmic Strings,” Phys. Rev. D 99 no. 6, (2019) 063529, arXiv:1901.03312 [hep-ph].
  * [67] P. Mocz et al., “Galaxy formation with BECDM – II. Cosmic filaments and first galaxies,” Mon. Not. Roy. Astron. Soc. 494 no. 2, (2020) 2027–2044, arXiv:1911.05746 [astro-ph.CO].
  * [68] M. Gorghetto, E. Hardy, J. March-Russell, N. Song, and S. M. West, “Dark Photon Stars: Formation and Role as Dark Matter Substructure,” arXiv:2203.10100 [hep-ph].
  * [69] B. Schwabe, J. C. Niemeyer, and J. F. Engels, “Simulations of solitonic core mergers in ultralight axion dark matter cosmologies,” Phys. Rev. D 94 no. 4, (2016) 043513, arXiv:1606.05151 [astro-ph.CO].

##  Appendix A Appendix

###  A.1 Simulation setup

We simulate $\psi_{i}$ ($i=1,2,3$) governed by Eq. (2.4) using a split-step
Fourier method [11] (as was done for the scalar case, but now we have 3
components instead of 1). Our fields live on a discretized three dimensional
grid with resolution $\Delta x$. This method has a temporal error of
$\mathcal{O}(\Delta t^{2})$, where $\Delta t$ is chosen by the Courant-
Friedrichs-Lewy (CFL) condition [69]: $\Delta
t\leq\max\Big{[}\frac{m}{6\hbar}\Delta
x^{2},\frac{\hbar}{m|\Phi|_{\max}}\Big{]}$, every iteration.

Quantities like mass and spin that are calculated with norm-preserving
operations are conserved up to machine precision. We also monitor fractional
energy loss $\Delta E/E_{\mathrm{tot}}$ as a function of time as a measure of
our simulations fidelity. We have performed our simulations with $320^{3}$
grid points for $N$ soliton simulations, where we have observed $\Delta
E/E_{\mathrm{tot}}\lesssim\mathcal{O}(10^{-5})$ for VDM and $\Delta
E/E_{\mathrm{tot}}\lesssim\mathcal{O}(10^{-4})$ for SDM. We have used
$256^{3}$ grid-points for two soliton simulations, where we have observed
$\Delta E/E_{\mathrm{tot}}\lesssim\mathcal{O}(10^{-4})$ for both VDM and SDM.

###  A.2 Estimating the final core mass

This calculation essentially follows the arguments in [42].131313Note that
they refer to core as only the mass within $r_{c}$. This is $1/4$ of the mass
of the soliton. Consider the case where we have $N$ solitons at time $t_{0}$,
each of mass $M_{\rm sol}^{\rm i}=M_{\rm sol}^{(0)}$. At time $t_{1}$, they
merge pair wise, to yield $N/2$ solitons, each with mass $M_{\rm
sol}^{(1)}=2fM_{\rm sol}^{(0)}$. Here, $1-f$ is the fraction of the progenitor
mass that does not contribute to the merged soliton. At time $t_{2}$, the
$N/2$ solitons merge again to yield $N/4$ solitons, each with mass $M_{\rm
sol}^{(2)}=2fM_{\rm sol}^{(1)}=(2f)^{2}M_{\rm sol}^{(0)}$. This process
continues until we have a single final soliton, this happen at time $t_{n}$
where $n$ is determined from $N/2^{n}=1$. That is, the merger process is
complete after $n=\log_{2}(N)$ merger generations. With this, the final mass
of the soliton $M_{\rm sol}^{\rm f}$ is

$M_{\rm sol}^{\rm f}/M_{\rm sol}^{\rm i}=(2f)^{\log_{2}(N)}.$ (A.1)

Of course this is all very crude, and mergers will likely not proceed with
just equal mass solitons merging. Mergers with highly unequal mass are
expected to not contribute increasing the mass of the merged soliton. This
will reduce the effective $N$ to $\gamma N$.

How can the vector and scalar cases be different? First, the $(1-f)$ fraction
during merger might be different between VDM and SDM. We expect $f_{\rm
v}<f_{\rm s}$. It is unclear to us whether the reduced interference might also
lead to fewer mergers in VDM.

###  A.3 Continuity equations for mass and spin density

For our purposes, we were be particularly interested in following the mass
density, and the spin density in the system. The equations of motion
immediately lead to the following local conservation laws:

$\displaystyle\partial_{t}\rho+\nabla\cdot\bm{\mathsf{j}}_{\rho}=0\,,\qquad\bm{\mathsf{j}}_{\rho}=i\frac{\hbar}{2}\left[{\bm{\mathsf{\Psi}}}\nabla{\bm{\mathsf{\Psi}}}^{\dagger}-{\bm{\mathsf{\Psi}}}^{\dagger}\nabla{\bm{\mathsf{\Psi}}}\right],$
(A.2) $\displaystyle\partial_{t}\bm{s}+\nabla\cdot\bm{\mathsf{j}}_{\rm
s}=0\,,\qquad\bm{\mathsf{j}}_{\rm
s}=\frac{\hbar^{2}}{2m}\left[\nabla{\bm{\mathsf{\Psi}}}\times{\bm{\mathsf{\Psi}}}^{\dagger}-{\bm{\mathsf{\Psi}}}\times\nabla{\bm{\mathsf{\Psi}}}^{\dagger}\right].$
(A.3)

where $\bm{\mathsf{j}}_{\rho}$ and $\bm{\mathsf{j}}_{\rm s}$ are the mass
current vector, and spin current tensor respectively. Similar equations can be
written for other conserved quantities also.

Figure 7: $k^{3}$ scaled power spectrum of the density field for an ensemble
of 10 simulations concentrated around $M=2.3\times
10^{5}M_{\odot}\times\mathcal{M}_{5}$ with $N=21$ solitons each. The solid
lines represent the mean, whereas the shaded bands include variation from
different simulations.

###  A.4 Power Spectrum

We construct the power spectrum of the density field:

$\displaystyle P_{\rho}(k)=\langle|\rho_{\bm{k}}|^{2}\rangle\,,$ (A.4)

with $\rho_{\bm{k}}\equiv
V^{-1/2}\int_{V}d^{3}x\rho(\bm{x})e^{i\bm{k}\cdot\bm{x}}$ and
$\langle\ldots\rangle$ representing an average over all $|\bm{k}|=k$. Power
spectra for an ensemble of 10 runs are shown in Fig. 7 (with $M=2.3\times
10^{5}M_{\odot}\times\mathcal{M}_{5}$ and $N=21$ in each run).

From our simulations, we find that for VDM there is a characteristic peak at
$k_{\rm peak}\approx 2.5\,\textrm{kpc}^{-1}\times(\mathcal{M}_{5}m_{20}^{2})$,
whereas for SDM, the peak in the power spectrum is at $k_{\rm peak}\approx
3.2\,\textrm{kpc}^{-1}\times(\mathcal{M}_{5}m_{20}^{2})$. Note that for
identical initial conditions, we find that typically the power in VDM is
smaller, and peaks at smaller wave numbers compared to SDM. This is tightly
correlated with the fact that the central soliton is less dense in VDM
compared to SDM.