# Polarized solitons in higher-spin wave dark matter

Mudit Jain<EMAIL_ADDRESS>Mustafa A. Amin<EMAIL_ADDRESS>Department of Physics and Astronomy, Rice University
Houston, TX, 77005, U.S.A.

###### Abstract

We first show that the effective non-relativistic theory of gravitationally
interacting, massive integer-spin fields (spin-$0$, $1$, and $2$ in
particular) is described by a $2s+1$ component Schrödinger-Poisson action,
where $s$ is the spin of the field. We then construct $s+1$ distinct,
gravitationally supported solitons in this non-relativistic theory from
identically polarized plane waves. Such solitons are extremally polarized,
with macroscopically large spin, but no orbital angular momentum. These $s+1$
solitons form a basis set, out of which partially polarized solitons can be
constructed. All such solitons are ground states, have a spherically symmetric
energy density but not field configurations. We discuss how solitons in
higher-spin fields can be distinguished from scalar solitons, and potential
gravitational and non-gravitational probes of them.

††preprint: APS/123-QED

##  I Introduction

Dark matter makes up approximately $84\%$ of the non-relativistic matter in
our cosmos [1], but its detailed nature remains uncertain. For example, the
mass of the fundamental particles making up dark matter can range from $\sim
10^{-21}\rm eV$ [2] to $\sim m_{\mathrm{pl}}$. The upper bound is extended
further if dark matter is multi-component, or composite [3]. We also have no
robust constraints on the spin of the particles/fields that make up dark
matter. What we do know from a plethora of observations is 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 [4]).

With the aim of better understanding dark matter, we wish to explore how the
fundamental degrees of freedom, such as mass and spin of the dark field can
manifest themselves on macroscopic scales. As an extreme case for the mass,
making the dark matter particle ultra-light is one way to boost its Compton/de
Broglie scale to astrophysical scales, thus allowing us to probe astrophysical
effects sensitive to this mass [5, 6, 7].

Figure 1: A visual summary of some of the main results of our paper.

In this paper, we are interested in whether the spin of the dark fields can
also have implications on macroscopic scales. Since we are interested in this
field being all or part of dark matter, we first derive the non-relativistic
limit of massive, higher-spin bosonic field that is capable of clustering
under gravity. We show that the non-relativistic behaviour of a spin-$s$
bosonic field is described by a $2s+1$ component Schrödinger-Poisson system.
Our results are also applicable to fields in the early universe, for example,
at the end of inflation [8, 9] or later [10].

Fields with non-linear interactions (gravitational or self-interactions) often
allow for coherent, metastable field configurations called solitons. With this
non-relativistic limit of our higher-spin dark fields at hand, we derive the
lowest energy soliton solutions in spin-$0$, spin-$1$ and spin-$2$ fields. In
general for a spin-$s$ field, we have $s+1$ distinct and “extremally”
polarized solitons with equal energy, where all the waves making up a given
soliton have the same multiplicity of spin (which we refer to as
polarization). We show that the total intrinsic spin of such solitons can be
macroscopically (even astrophysically) large, while still having zero orbital
angular momentum and a spherical energy density. We also construct
fractionally polarized solitons from these extremally polarized ones. All such
solitons have a universal mass density profile (up to scalings). We call all
such ground-state solitons “$p$-solitons”. For a visual summary of the main
results of our paper, see Fig. 1. We also discuss higher energy solitons, and
compare them to our $p$-solitons.

The non-relativistic limit of scalar (spin-$0$) dark matter has of course been
explored before. For recent works, including relativistic corrections, see
[11, 12, 13]. The non-relativistic limit of vector (spin-$1$) dark matter in
an expanding universe was discussed in [14], whereas the tensor (spin-$2$)
case in an asymptotically flat spacetime was discussed in [15]. In our work,
we recover these results and provide a framework that allows for a
generalization to higher spins.111The massive spin-$2$ case comes with
subtleties related to having two metrics, self-interactions, limiting the
coupling to the Standard model for these fields to be dark, as well as having
a consistent theory at sufficiently high energies. Some of these are discussed
in the Appendix.

Similarly, different types of solitons in the non-relativistic limit of
massive spin-$0$ [16, 17], spin-$1$ [14] and spin-$2$ field [15] appear in the
literature. A distinguishing feature of our analysis is that we construct
solitons with a focus on the polarization of their constituent plane waves.
This lead us to uncover a novel class of extremally polarized solitons with
spin $S_{\rm tot}/\hbar=\lambda M/m$ which can be macroscopically large for
$M\gg m$. Here, $m$ is the mass of the field, $M$ is the mass of the soliton
and $\lambda$ is the spin multiplicity. These coherent solitons (along with
fractionally polarized ones mentioned earlier) might open up new avenues for
observationally probing higher-spin fields.

We find that even within Newtonian gravity it might be possible to distinguish
interacting solitons with different polarizations. Going beyond Newtonian
gravity, which we do not pursue here, might remove degeneracies between
different polarizations of the higher-spin fields even further. We also
discuss possibilities of probing higher-spin dark matter via non-gravitational
interactions astrophysically, taking advantage of the polarization state of
the solitons. In this case, various terrestial experiments [18, 19, 20] may be
used for detectable signatures.

The paper is organized as follows. In section II we discuss our model for the
case of dark scalar, vector, and tensor massive fields, leaving additional
details in Appendix A. In section III we provide the effective non-
relativistic action (which is the Schrödinger-Poisson system) for these dark
integer spin fields, and discuss the various symmetries of the action. In
section IV we discuss the gravitationally bound solitons. In section V we
discuss their distinguishability, primarily within Newtonian gravity, and also
mention other non-gravitational couplings that can probe the spin nature of
the fields. In section VI we summarize and also highlight some future
directions worth investigating.

##  II Models

Our matter Lagrangian consists of the usual Standard Model (SM) sector, along
with some dark sector that includes additional massive spin-$0$, spin-$1$, or
even spin-$2$ fields. We take these fields to be real valued.

Explicitly, our general action has the form

$\displaystyle S=S_{\rm EH}+S_{\rm dark}+S_{\rm vis}\,,$ (1)

where $S_{\rm EH}$ is the gravity sector, $S_{\rm dark}$ is some dark sector
(incluing dark integer spin fields), and $S_{\rm vis}$ is the visible sector
(comprising of the SM). Our focus is only on the gravity + dark sector in this
paper. We consider perturbations of different fields around some background
metric $\bar{g}_{\mu\nu}$ which leads to the usual massless spin-$2$
fluctuations: $h_{\mu\nu}$ (the graviton), along with other perturbations in
different fields. We will focus on a given spin-$s$ field + gravity, instead
of considering massive spin-$0$, $1$ and $2$ together, although our formalism
can accomodate the latter scenario as well.

For most part, we are interested in sub-horizon physics where length scales
associated with configurations of these dark fields are much smaller than the
Hubble horizon. As a result, we ignore Hubble expansion, and take the
background metric to be222We use
$\bar{g}_{\mu\nu}=\mathrm{diag}(1,-a^{2}(t),-a^{2}(t),-a^{2}(t))$ for an
expanding universe when needed. Here, $a(t)$ is the scale factor normalized to
unity today. $\bar{g}_{\mu\nu}=\eta_{\mu\nu}=\mathrm{diag}(1,-1,-1,-1)$. We
also take $\hbar=c=1$.

In the next three subsections, we provide the general action up-to quadratic
order in the fields of interest, along with leading order gravitational
interactions. For the non-relativistic limit that we are interested in, the
leading order actions provided here are sufficient. The full nonlinear actions
are discussed in the Appendix.

###  II.1 Spin-$0$

The quadratic (free) action for the spin-$0$ field $\phi$, and metric
fluctuations $h_{\mu\nu}$, along with their leading interactions is given by

$\displaystyle S_{\rm EH}+S_{\rm dark}$
$\displaystyle=\int\mathrm{d}^{4}x\Bigl{[}m_{\mathrm{pl}}^{2}\mathcal{L}^{(2)}_{\rm
GR}(h)+\mathcal{L}^{(2)}_{m,0}(\phi)$
$\displaystyle\qquad\qquad\quad-\frac{1}{2}h_{\mu\nu}\mathcal{T}^{\mu\nu}(\phi)+...\Bigr{]},$
(2)

where $\mathcal{L}^{(2)}_{\rm GR}(h)$ is the linearized Einstein-Hilbert
action:

$\displaystyle\mathcal{L}^{(2)}_{\rm GR}(h)$
$\displaystyle=\frac{1}{4}\Bigl{\\{}\eta^{\alpha\beta}\eta^{\lambda\sigma}\eta^{\kappa\rho}\partial_{\beta}h_{\sigma\rho}\partial_{\alpha}h_{\lambda\kappa}-\eta^{\alpha\beta}\partial_{\beta}h\partial_{\alpha}h$
$\displaystyle+2\eta^{\lambda\beta}\eta^{\kappa\alpha}\partial_{\beta}h\partial_{\alpha}h_{\lambda\kappa}-2\eta^{\kappa\rho}\eta^{\sigma\beta}\eta^{\lambda\alpha}\partial_{\beta}h_{\sigma\rho}\partial_{\alpha}h_{\lambda\kappa}\Bigr{\\}},$
(3)

and $\mathcal{L}^{(2)}_{m,0}(\phi)$ and $\mathcal{T}_{\mu\nu}(\phi)$ are the
usual Klein-Gordon Lagrangian density and leading order energy-momentum
tensor:

$\displaystyle\mathcal{L}^{(2)}_{m,0}(\phi)$
$\displaystyle=\frac{1}{2}\eta^{\mu\nu}\partial_{\mu}\phi\,\partial_{\mu}\phi-\frac{1}{2}m^{2}\phi^{2}\,,$
$\displaystyle\mathcal{T}_{\mu\nu}(\phi)$
$\displaystyle=\partial_{\mu}\phi\,\partial_{\nu}\phi-\eta_{\mu\nu}\,\mathcal{L}^{(2)}_{m,0}(\phi).$
(4)

The ‘$...$’ in (II.1) represents higher-order interaction terms between $\phi$
and $h$ along with other dark sector fields. Self-interactions can easily be
included by adding $V_{\rm
nl}(\phi)=\lambda_{3}\phi^{3}+\lambda_{4}\phi^{4}\ldots$ to the quadratic
potential , but we ignore them in this paper.

While the massive spin-$0$ field has no constraints, the gravitational sector
does. We discuss these constraints and the action for the physical degrees of
freedom in the gravitational sector at the end of this section.

###  II.2 Spin-$1$

For a massive spin-$1$ field $W_{\mu}$, we have the following leading order
action for $h_{\mu\nu}$ and $W_{\mu}$:

$\displaystyle S_{\rm EH}+S_{\rm dark}$
$\displaystyle=\int\mathrm{d}^{4}x\Bigl{[}m_{\mathrm{pl}}^{2}\mathcal{L}^{(2)}_{GR}(h)+\mathcal{L}^{(2)}_{m,1}(W)$
$\displaystyle\qquad\qquad\quad-\frac{1}{2}h_{\mu\nu}T^{\mu\nu}(W)+...\Bigr{]},$
(5)

where $G_{\mu\nu}=\partial_{\mu}W_{\nu}-\partial_{\nu}W_{\mu}$;
$\mathcal{L}^{(2)}_{m,1}(W)$ and $T_{\mu\nu}(W)$ are the usual Proca
Lagrangian density and leading order energy momentum tensor:

$\displaystyle\mathcal{L}^{(2)}_{m,1}(W)$
$\displaystyle=-\frac{1}{4}\eta^{\mu\alpha}\eta^{\nu\beta}\,G_{\mu\nu}G_{\alpha\beta}+\frac{1}{2}m^{2}\,\eta^{\mu\nu}W_{\mu}W_{\nu}$
$\displaystyle T_{\mu\nu}(W)$
$\displaystyle=-\eta^{\alpha\beta}G_{\mu\alpha}G_{\nu\beta}+m^{2}W_{\mu}W_{\nu}-\eta_{\mu\nu}\,\mathcal{L}^{(2)}_{m,1}(W).$
(6)

Similar to the case of the scalar field, we can include self-interactions via
$V_{\textrm{nl}}=-\lambda_{4}(W_{\mu}W^{\mu})^{2}+\lambda_{6}(W_{\mu}W^{\mu})^{3}+\ldots$.
As we shall show elsewhere [21], such a setup with $\lambda_{4},\lambda_{6}>0$
can arise as an effective theory, and admits vector oscillon solutions
supported by self-interactions. However, for our present purposes, we assume
no self-interactions.

#### Action for the physical d.o.f

To aid the transition to the non-relativistic limit, we eliminate the
constraints in the massive spin-$1$ sector and provide the action for the
three physical degrees of freedom (namely the three spin-multiplicity states).
Varying the action (II.2) with respect to $W_{0}$ and substituting back,
yields333This is a straightforward, but long exercise which we omit here to
reduce clutter.

$\displaystyle\mathcal{L}^{(2)}_{m,1}-\frac{1}{2}h_{\mu\nu}T^{\mu\nu}(W)$
$\displaystyle=\frac{1}{2}\dot{W}_{i}\,\mathcal{P}_{ij}(m)\,\dot{W}_{j}-\frac{1}{2}W_{i}\left(-\nabla^{2}+m^{2}\right)\,\mathcal{P}_{ij}(m)\,W_{j}$
$\displaystyle\quad-\frac{1}{2}h_{\mu\nu}\mathcal{T}^{\mu\nu},$ (7)

where $\mathcal{T}^{\mu\nu}$ is the energy momentum tensor (with $W_{0}$
substituted in (II.2)), and $\mathcal{P}_{ij}(m)$ is the projection operator

$\displaystyle\mathcal{P}_{ij}(m)\equiv\delta_{ij}+\frac{\partial_{i}\partial_{j}}{-\nabla^{2}+m^{2}}\,.$
(8)

As mentioned earlier, the gravitational sector constraints are discussed at
the end of this section.

###  II.3 Spin-$2$

In the case of a massive spin-$2$ field $H_{\mu\nu}$, we take the quadratic
actions for $h_{\mu\nu}$ and $H_{\mu\nu}$ \+ leading interactions to be:

$\displaystyle S_{h,H}$
$\displaystyle=\int\mathrm{d}^{4}x\Bigl{\\{}m_{\mathrm{pl}}^{2}\mathcal{L}^{(2)}_{GR}(h)+\mathcal{L}^{(2)}_{m,2}(H)$
$\displaystyle\qquad\qquad\qquad-\frac{1}{2}h_{\mu\nu}T^{\mu\nu}(H)+...\Bigr{\\}},$
(9)

where $\mathcal{L}^{(2)}_{m,2}(H)$ is the Fierz-Pauli Lagrangian density and
$T^{\mu\nu}(H)$ is the associated energy momentum tensor:

$\displaystyle\mathcal{L}^{(2)}_{m,2}(H)$
$\displaystyle=\frac{1}{2}\eta^{\alpha\beta}\eta^{\lambda\sigma}\eta^{\kappa\rho}\partial_{\beta}H_{\sigma\rho}\partial_{\alpha}H_{\lambda\kappa}-\frac{1}{2}\eta^{\alpha\beta}\partial_{\beta}H\partial_{\alpha}H$
$\displaystyle-\eta^{\kappa\rho}\eta^{\sigma\beta}\eta^{\lambda\alpha}\partial_{\beta}H_{\sigma\rho}\partial_{\alpha}H_{\lambda\kappa}+\eta^{\lambda\beta}\eta^{\kappa\alpha}\partial_{\beta}H\partial_{\alpha}H_{\lambda\kappa}$
$\displaystyle+\frac{1}{2}m^{2}\left[\eta^{\mu\nu}\eta^{\alpha\beta}H_{\mu\nu}H_{\alpha\beta}-\eta^{\sigma\lambda}\eta^{\rho\kappa}H_{\sigma\rho}H_{\lambda\kappa}\right];$
$\displaystyle T_{\mu\nu}(H)$
$\displaystyle=\Bigl{\\{}\partial_{\mu}H_{\sigma\rho}\partial_{\nu}H^{\sigma\rho}+2\eta^{\kappa\rho}\partial^{\alpha}H_{\nu\rho}\partial_{\alpha}H_{\mu\kappa}$
$\displaystyle-2\partial^{\beta}H_{\beta\nu}\partial^{\alpha}H_{\alpha\mu}-2\partial_{\nu}H_{\mu\rho}\partial_{\alpha}H^{\alpha\rho}+\partial_{\nu}H\partial^{\alpha}H_{\mu\alpha}$
$\displaystyle+\partial^{\alpha}H\partial_{\nu}H_{\alpha\mu}-2\partial_{\mu}H_{\nu\rho}\partial_{\alpha}H^{\alpha\rho}+\partial_{\mu}H\partial^{\alpha}H_{\nu\alpha}$
$\displaystyle+\partial^{\alpha}H\partial_{\mu}H_{\alpha\nu}-\partial_{\mu}H\partial_{\nu}H\Bigr{\\}}$
$\displaystyle-2m^{2}\left(\eta^{\alpha\beta}H_{\mu\alpha}H_{\nu\beta}-H_{\mu\nu}H\right)$
$\displaystyle-2\partial_{\xi}\Bigl{\\{}H_{\kappa\mu}\partial_{\nu}H^{\kappa\xi}+H_{\kappa\nu}\partial_{\mu}H^{\kappa\xi}+\eta^{\kappa\rho}H_{\kappa\mu}\partial^{\xi}H_{\nu\rho}$
$\displaystyle+\eta^{\kappa\rho}H_{\kappa\nu}\partial^{\xi}H_{\mu\rho}-H^{\alpha\xi}\partial_{\nu}H_{\mu\alpha}-H^{\alpha\xi}\partial_{\mu}H_{\nu\alpha}$
$\displaystyle-2\delta^{\xi}_{\nu}H_{\kappa\mu}\partial_{\beta}H^{\beta\kappa}-2\delta^{\xi}_{\mu}H_{\kappa\nu}\partial_{\beta}H^{\beta\kappa}-H_{\mu\nu}\partial_{\beta}H^{\beta\xi}$
$\displaystyle+H_{\mu\nu}\partial^{\xi}H+\delta^{\xi}_{\nu}H_{\alpha\mu}\partial^{\alpha}H+\delta^{\xi}_{\mu}H_{\alpha\nu}\partial^{\alpha}H$
$\displaystyle-\eta_{\mu\nu}H^{\alpha\xi}\partial_{\alpha}H\Bigr{\\}}-\eta_{\mu\nu}\,\mathcal{L}^{(2)}_{m,2}(H)\,.$
(10)

At this leading order, we obtained the energy momentum tensor by promoting
$\eta_{\mu\nu}$ to $g_{\mu\nu}=\eta_{\mu\nu}+h_{\mu\nu}$ in
$\mathcal{L}^{(2)}_{m,2}(H)$, along with covariant derivatives with respect to
$g$, and then simply reading out the terms that couple to $h_{\mu\nu}$. That
is, simply minimally coupling it to gravity just like any other matter field.
We also arrive at the above action from the bi-gravity theory [22, 23, 24, 25,
26, 27, 28] (which is a ghost-free, nonlinear completion of our quadratic
action) in the Appendix.

The massive spin-$2$ field scenario includes additional subtleties compared to
the lower spin cases. Unlike the spin-$0$ and spin-$1$ case, self-interactions
are unavoidable in the massive spin-$2$ sector within bi-gravity. However,
they are conveniently suppressed in the non-relativistic limit (see discussion
around eq. (65) in the Appendix). Furthermore, for the massive spin-$2$ field
to be dark, its direct couplings to the visible sector that arise in the bi-
gravity theory need to be suppressed. This is also possible using a particular
choice of an effective metric [29, 30] that couples to the SM. We discuss this
further in Sec. [A.3.1] of the Appendix. Within the bi-gravity context, the
extension to FLRW background is also non-trivial, and might have instabilities
in the early universe [31] when using the effective metric mentioned above.
For our purposes, we remain agnostic regarding a full non-linear interacting
theory of a massive spin-$2$ dark field in an FLRW background, and simply work
with the above quadratic action.

#### Actions for physical d.o.f

For massive tensor field $H_{\mu\nu}$, the constrained variables are $H_{00}$,
$H_{0i}$ and $\operatorname{Tr}[H_{ij}]$, leaving 5 physical degrees of
freedom (the 5 spin-multiplicity states). After solving for these constraints
(obtained by varying the action (II.3) with respect to these variables), and
plugging them back into the action, we get the following Lagrangian density

$\displaystyle\mathcal{L}^{(2)}_{m,2}-\frac{1}{2}h_{\mu\nu}T^{\mu\nu}$
$\displaystyle=\frac{1}{2}\dot{H}_{ij}\,\mathcal{P}_{ik}(m)\,\mathcal{P}_{jl}(m)\,\dot{H}_{kl}$
$\displaystyle-\frac{1}{2}H_{ij}\,\mathcal{P}_{ik}(m)\,(-\nabla^{2}+m^{2})\,\mathcal{P}_{jl}(m)\,H_{kl}$
$\displaystyle-\frac{1}{2}h_{\mu\nu}\mathcal{T}^{\mu\nu}.$ (11)

Here, $\mathcal{P}_{ij}$ is the projection matrix (8), and
$\mathcal{T}_{\mu\nu}(H_{ij})$ is the energy momentum tensor (dependent on the
physical degrees of freedom only).

###  II.4 The gravitational sector

Recall that the action has the general form:

$\displaystyle S_{\rm EH}+S_{\rm dark}$
$\displaystyle=\int\mathrm{d}^{4}x\Bigl{[}m_{\mathrm{pl}}^{2}\mathcal{L}^{(2)}_{\rm
GR}(h)+\mathcal{L}^{(2)}_{m,s}(\bm{\mathcal{F}})$
$\displaystyle\qquad\qquad\quad-\frac{1}{2}h_{\mu\nu}\mathcal{T}^{\mu\nu}(\bm{\mathcal{F}})+...\Bigr{]},$
(12)

where $\bm{\mathcal{F}}$ now represents the physical degrees of freedom in the
spin-$s$ field, and $\mathcal{T}^{\mu\nu}$ is the corresponding energy-
momentum tensor. Varying this action with respect to $h_{0i}$,
$\mathrm{Tr}[h_{ij}]$, and $h_{00}\equiv 2\Phi$, we obtain

$\displaystyle\nabla^{2}\,\mathcal{P}_{ij}(0)\,h_{0i}$
$\displaystyle=\partial_{i}\dot{h}_{ij}-\partial_{j}{\rm
Tr}[\dot{h}_{ij}]-\frac{\mathcal{T}_{0j}}{m_{\mathrm{pl}}^{2}}\,,$
$\displaystyle\nabla^{2}{\rm Tr}[h_{ij}]$
$\displaystyle=\partial_{i}\partial_{j}h_{ij}+\frac{\mathcal{T}_{00}}{m_{\mathrm{pl}}^{2}}\,,$
$\displaystyle 2\,\nabla^{2}\Phi$
$\displaystyle=-\frac{\partial_{i}\partial_{j}}{\nabla^{2}}\,\ddot{h}_{ij}-\Box\frac{1}{\nabla^{2}}\frac{\mathcal{T}_{00}}{m_{\mathrm{pl}}^{2}}\,.$
(13)

When these constraints are plugged back in the above action, we get the
quadratic action for the physical degrees of freedom in the gravitational
sector + interactions. The quadratic action is the same as the first two terms
on the right hand side of (II.3), with $m=0$ and $H_{ij}\rightarrow
h_{ij}$444The projection operator $P_{ij}(0)$, along with the constraints,
only allows for the 2 physical degrees of freedom in $h_{ij}$ (which are the
transverse and traceless gravitational waves)., and the interaction term is

$\displaystyle\frac{1}{2}\,h_{\mu\nu}\,\mathcal{T}^{\mu\nu}$
$\displaystyle=\frac{1}{2}h_{ij}\Biggl{\\{}\mathcal{T}_{ij}+\partial_{i}\partial_{j}\frac{1}{\nabla^{2}}\frac{1}{\nabla^{2}}\ddot{\mathcal{T}}_{00}-\partial_{j}\frac{1}{\nabla^{2}}\dot{\mathcal{T}}_{0i}$
$\displaystyle\qquad\quad-\partial_{i}\frac{1}{\nabla^{2}}\dot{\mathcal{T}}_{0j}+2\,\partial_{i}\partial_{j}\frac{1}{\nabla^{2}}\partial_{k}\frac{1}{\nabla^{2}}\dot{\mathcal{T}}_{0k}\Biggr{\\}}.$
(14)

##  III Effective non-relativistic theory

With the Lagrangian densities at hand from the previous section, we now
decompose $\bm{\mathcal{F}}$ (which represents the physical d.o.f of different
integer spin-fields and carries one and two spatial indices for spin-$1$ and
spin-$2$ cases respectively) in the following fashion

$\displaystyle\bm{\mathcal{F}}({\bf
x},t)=\frac{1}{\sqrt{2m}}\left[e^{-imt}\tilde{{\bm{\mathsf{\Psi}}}}({\bf
x},t)+\mathrm{h.c.}\right].$ (15)

To obtain the non-relativistic limit, we work with the slowly varying piece in
$\tilde{{\bm{\mathsf{\Psi}}}}$ that we denote as ${\bm{\mathsf{\Psi}}}$.
Essentially, we discard all the terms that carry the oscillating factor
$e^{\pm i2mt}$ and two time derivative terms in the action (which would be
suppressed by factors of $k/m$, where $k$ is a characteristic wave-number).
The projection operator (8) simplifies to $P_{ij}|_{\rm
nr}=\delta_{ij}+\mathcal{O}(k^{2}/m^{2})$. Upon making these approximations,
we arrive at the following free Schrödinger action for massive scalar, vector
and tensor fields:

$\displaystyle\mathcal{S}^{\mathrm{eff}}_{\mathrm{free}}$
$\displaystyle=\int\mathrm{d}^{4}x\,\Biggl{\\{}\frac{i}{2}\operatorname{Tr}\left[{\bm{\mathsf{\Psi}}}^{\dagger}\dot{{\bm{\mathsf{\Psi}}}}\right]+\mathrm{h.c.}$
$\displaystyle\qquad\qquad\qquad\qquad-\frac{1}{2m}\operatorname{Tr}[\nabla{\bm{\mathsf{\Psi}}}^{\dagger}\cdot\nabla{\bm{\mathsf{\Psi}}}]\Biggr{\\}}.$
(16)

For the case of spin-$0$, ${\bm{\mathsf{\Psi}}}$ carries no spatial index,
while for the case of spin-$1$ and spin-$2$, it carries one and two spatial
indices respectively. We shall refer to their components as

$\displaystyle\psi=[{\bm{\mathsf{\Psi}}}]\qquad$
$\displaystyle\textrm{spin-$0$}\,,$
$\displaystyle\psi_{i}=[{\bm{\mathsf{\Psi}}}]_{i}\qquad$
$\displaystyle\textrm{spin-$1$}\,,$ (17)
$\displaystyle\psi_{ij}=[{\bm{\mathsf{\Psi}}}]_{ij}\qquad$
$\displaystyle\textrm{spin-$2$}\,.$

Note that for the spin-$1$ and spin-$2$ cases,
$\operatorname{Tr}[{\bm{\mathsf{\Psi}}}^{\dagger}{\bm{\mathsf{\Psi}}}]=\psi_{i}^{\dagger}\psi_{i}$
and
$\operatorname{Tr}[{\bm{\mathsf{\Psi}}}^{\dagger}{\bm{\mathsf{\Psi}}}]=\psi_{ij}^{\dagger}\psi_{ij}$
respectively, where summation over indices in implicit.555Throughout the text,
we implicitly assume summation over repeated indices, unless mentioned
otherwise. Also for the massive spin-$2$ case, we have
$\operatorname{Tr}[{\bm{\mathsf{\Psi}}}]=\mathcal{O}(k^{2}{\bm{\mathsf{\Psi}}}/m^{2})$
and therefore can be neglected in the non-relativistic limit. Furthermore, we
obtain the following general structure for the non-relativistic energy
momentum tensor

$\displaystyle\mathcal{T}_{00}$
$\displaystyle=m\,\operatorname{Tr}\left[{\bm{\mathsf{\Psi}}}^{\dagger}{\bm{\mathsf{\Psi}}}\right]+\mathcal{O}\left(\frac{k^{2}}{m}{\bm{\mathsf{\Psi}}}^{\dagger}{\bm{\mathsf{\Psi}}}\right),$
$\displaystyle\mathcal{T}_{0i}$
$\displaystyle=\mathcal{O}\left(k\,{\bm{\mathsf{\Psi}}}^{\dagger}{\bm{\mathsf{\Psi}}}\right),$
$\displaystyle\mathcal{T}_{ij}$
$\displaystyle=\mathcal{O}\left(\frac{k^{2}}{m}{\bm{\mathsf{\Psi}}}^{\dagger}{\bm{\mathsf{\Psi}}}\right).$
(18)

With this, we can see from (II.4) that there are no gravitons (gravitational
waves) produced in the non-relativistic limit since the source term for
$h_{ij}$ is
$\mathcal{O}(k^{2}{\bm{\mathsf{\Psi}}}^{\dagger}{\bm{\mathsf{\Psi}}}/m)$ for
all the three cases. From (II.4), it is also clear that the vector constraint
$h_{0i}$ is not sourced at the leading order in either of the three cases. The
only constraint that survives is the Newtonian potential $\Phi$, which is
determined by

$\displaystyle\nabla^{2}\Phi=\frac{m}{2m_{\mathrm{pl}}^{2}}\,\operatorname{Tr}[{\bm{\mathsf{\Psi}}}^{\dagger}{\bm{\mathsf{\Psi}}}]+\mathcal{O}\left(\frac{k^{2}}{m}{\bm{\mathsf{\Psi}}}^{\dagger}{\bm{\mathsf{\Psi}}}\right).$
(19)

Note that the trace ${\rm Tr}[h_{ij}]$ is equal to twice the Newtonian
potential, ${\rm Tr}[h_{ij}]=2\Phi$. Putting it all together, using (III),
(III), and (19), we get the following Schrödinger Poisson action

$\displaystyle\mathcal{S}^{\mathrm{eff}}_{nr}$
$\displaystyle=\int\mathrm{d}^{4}x\,\Bigl{[}\frac{i}{2}\operatorname{Tr}\left[{\bm{\mathsf{\Psi}}}^{\dagger}\dot{{\bm{\mathsf{\Psi}}}}\right]+\mathrm{c.c.}-\frac{1}{2m}\operatorname{Tr}[\nabla{\bm{\mathsf{\Psi}}}^{\dagger}\cdot\nabla{\bm{\mathsf{\Psi}}}]$
$\displaystyle\qquad\qquad\quad+m_{\mathrm{pl}}^{2}\,\Phi\nabla^{2}\Phi-m\,\Phi\,\operatorname{Tr}[{\bm{\mathsf{\Psi}}}^{\dagger}{\bm{\mathsf{\Psi}}}]\Bigr{]}.$
(20)

The corresponding equation of motion is the Schrödinger-Poisson system

$\displaystyle i\frac{\partial}{\partial t}{\bm{\mathsf{\Psi}}}$
$\displaystyle=-\frac{1}{2m}\nabla^{2}{\bm{\mathsf{\Psi}}}+m\,\Phi\,{\bm{\mathsf{\Psi}}}\,,$
$\displaystyle\nabla^{2}\Phi$
$\displaystyle=\frac{m}{2m_{\mathrm{pl}}^{2}}\,\operatorname{Tr}[{\bm{\mathsf{\Psi}}}^{\dagger}{\bm{\mathsf{\Psi}}}].$
(21)

Extension to an FLRW universe can be achieved via $\nabla\rightarrow\nabla/a$
and $\partial/\partial t\rightarrow\partial/\partial t+3H/2$ where
$H=\dot{a}/a$.666For the spin-2 case, one needs to ensure that the choice of
FLRW background is consistent within the bi-gravity + matter theory. While
some aspects of this are discussed in the appendix, we leave a detailed
exploration to future work.

###  III.1 Conservation Laws

We now highlight various symmetries of the non-relativistic effective theory
(III). These symmetries will be helpful in understanding the space of soliton
solutions.

First, the action in invariant under ${\bm{\mathsf{\Psi}}}({\bf
x})\rightarrow\mathcal{M}({\bm{\mathsf{R}}})\,{\bm{\mathsf{\Psi}}}(\tilde{\bf
x})$, where ${\bm{\mathsf{R}}}$ is a rotation matrix and $\tilde{\bf
x}={\bm{\mathsf{R}}}^{-1}{\bf x}$. For scalars,
$\mathcal{M}({\bm{\mathsf{R}}})\,{\bm{\mathsf{\Psi}}}(\tilde{\bf
x})=\psi(\tilde{\bf x})$; for vectors
$\mathcal{M}({\bm{\mathsf{R}}})\,{\bm{\mathsf{\Psi}}}(\tilde{\bf
x})=R_{ij}\psi_{j}(\tilde{\bf x})$; and for tensors
$\mathcal{M}({\bm{\mathsf{R}}})\,{\bm{\mathsf{\Psi}}}(\tilde{\bf
x})=R_{ik}\psi_{kl}(\tilde{\bf x})R_{lj}$. The conserved charge density
associated with this, namely the total angular momentum density is the
following

$\displaystyle J_{k}$
$\displaystyle=s\,\Re\left(i\,\varepsilon_{ijk}[{\bm{\mathsf{\Psi}}}{\bm{\mathsf{\Psi}}}^{\dagger}]_{ij}\right)+\Re\left(i\,\varepsilon_{ijk}\mathrm{Tr}[{\bm{\mathsf{\Psi}}}^{\dagger}\partial_{i}{\bm{\mathsf{\Psi}}}]x^{j}\right)$
(22)

where we identify

$\displaystyle S_{k}$
$\displaystyle=s\,\Re\left(i\,\varepsilon_{ijk}[{\bm{\mathsf{\Psi}}}{\bm{\mathsf{\Psi}}}^{\dagger}]_{ij}\right),$
(23) $\displaystyle L_{k}$
$\displaystyle=\Re\left(i\,\varepsilon_{ijk}\mathrm{Tr}[{\bm{\mathsf{\Psi}}}^{\dagger}\partial_{i}{\bm{\mathsf{\Psi}}}]x^{j}\right)$

as the intrinsic spin and orbital angular momentum density
respectively.777Note that the total spin, and total orbital angular momentum
are independently conserved. In the above expressions, $s=\\{0,1,2\\}$ for
spin-$0$ (scalar), spin-$1$ (vector), and spin-$2$ (tensor) cases
respectively, and $\varepsilon_{a}$ are the totally anti-symmetric matrices
(generators of rotations). Also, $[..]_{ij}$ means the matrix obtained by
taking a tensor product of the elements within. For the vector case
$[{\bm{\mathsf{\Psi}}}{\bm{\mathsf{\Psi}}}^{\dagger}]_{ij}=\psi_{i}\psi^{\dagger}_{j}$,
while for the tensor case
$[{\bm{\mathsf{\Psi}}}{\bm{\mathsf{\Psi}}}^{\dagger}]_{ij}=\psi_{ik}\psi^{\dagger}_{kj}$.
For the scalar case this is zero.

The effective non-relativistic action (III) has a global $U(1)$ symmetry
(${\bm{\mathsf{\Psi}}}\rightarrow{\bm{\mathsf{\Psi}}}e^{i\alpha}$) leading to
a conserved particle number

$\displaystyle N=\int d^{3}x{\rm
Tr}[{\bm{\mathsf{\Psi}}}^{\dagger}{\bm{\mathsf{\Psi}}}]\,.$ (24)

Along with rotational invariance, the usual time-translation
invariance888There will also be linear momentum associated with invariance
under spatial translations. This however will not be of any direct use for us.
of (III) yields a conserved energy:

$\displaystyle E$
$\displaystyle=\int\mathrm{d}^{3}x\Bigl{[}\frac{1}{2m}\,\operatorname{Tr}[\nabla{\bm{\mathsf{\Psi}}}^{\dagger}\cdot\nabla{\bm{\mathsf{\Psi}}}]$
$\displaystyle+\frac{m^{2}}{4m_{\mathrm{pl}}^{2}}\operatorname{Tr}[{\bm{\mathsf{\Psi}}}^{\dagger}{\bm{\mathsf{\Psi}}}]\int\frac{\mathrm{d}^{3}y}{4\pi|{\bf
x}-{\bf y}|}\operatorname{Tr}[{\bm{\mathsf{\Psi}}}^{\dagger}({\bf
y}){\bm{\mathsf{\Psi}}}(\bf y)]\Bigr{]}$ (25)

There are additional conserved charges, which become apparent when we
decompose our field ${\bm{\mathsf{\Psi}}}$ into a polarization bases, which we
turn to next.

###  III.2 Decomposition into polarization basis

A massive spin-$s$ field admits $2s+1$ spin multiplicity states in some
particular direction $\hat{n}$, labelled by $\lambda\in\\{-s,...,s\\}$. These
states are characterized by the set
$\\{{\bm{\mathsf{\epsilon}}}^{(\lambda)}_{s,\hat{n}}\\}$, such that upon
substituting a plane wave ${\bm{\mathsf{\Psi}}}^{(\lambda)}=V^{-1/2}e^{i{\bf
k}\cdot{\bf x}}{\bm{\mathsf{\epsilon}}}^{(\lambda)}_{s,\hat{n}}$ in (23), we
get

$\displaystyle\hat{n}\cdot{\bm{S}}\left({\bm{\mathsf{\Psi}}}^{(\lambda)}\right)=\hat{n}\cdot{\bm{S}}\left({\bm{\mathsf{\epsilon}}}^{(\lambda)}_{s,\hat{n}}\right)=\frac{\lambda}{V}\qquad\forall\lambda\in\\{-s,..s\\}$
(26)

where $V$ is a spatial volume. For explicit forms of
${\bm{\mathsf{\epsilon}}}^{(\lambda)}_{s,\hat{n}}$, see (44) and (IV.1.3)
ahead for spin-$1$ and spin-$2$ cases, where we work with $\hat{n}=\hat{z}$
without loss of generality. The set
$\\{{\bm{\mathsf{\epsilon}}}^{(\lambda)}_{s,\hat{n}}\\}$ is orthogonal and
complete in the sense:

$\displaystyle\operatorname{Tr}[{{\bm{\mathsf{\epsilon}}}}^{(\lambda)\;\dagger}_{s,\hat{n}}{{\bm{\mathsf{\epsilon}}}}^{(\lambda^{\prime})}_{s,\hat{n}}]=\delta_{\lambda,\lambda^{\prime}},\quad\sum_{\lambda}\left[{{\bm{\mathsf{\epsilon}}}}^{(\lambda)}_{s,\hat{n}}\;{{\bm{\mathsf{\epsilon}}}}^{(\lambda)\dagger}_{s,\hat{n}}\right]_{ij}\propto\delta_{ij}.$
(27)

Using these $\\{{\bm{\mathsf{\epsilon}}}_{s,\hat{n}}^{(\lambda)}\\}$, the
field ${\bm{\mathsf{\Psi}}}$ admits the following decomposition

$\displaystyle{\bm{\mathsf{\Psi}}}({\bf
x},t)=\sum_{\lambda}\psi^{(\lambda)}_{s}({\bf
x},t){{\bm{\mathsf{\epsilon}}}}^{(\lambda)}_{s,\hat{n}}\,,$ (28)

where $\psi^{(\lambda)}_{s}$ is the field with polarization $\lambda$ in the
$\hat{n}$ direction. In terms of these different polarized fields
$\psi^{(\lambda)}_{s}$, the action has the following form

$\displaystyle\mathcal{S}^{\mathrm{eff}}_{nr}$
$\displaystyle=\sum_{\lambda}\int\mathrm{d}^{4}x\,\Bigl{[}\frac{i}{2}\psi^{(\lambda)\;\dagger}_{s}\dot{\psi}^{(\lambda)}_{s}+\mathrm{c.c.}$
$\displaystyle\qquad\qquad-\frac{1}{2m}\nabla\psi^{(\lambda)\;\dagger}_{s}\cdot\nabla\psi^{(\lambda)}_{s}+m_{\rm
pl}^{2}\,\Phi\nabla^{2}\Phi$
$\displaystyle\qquad\qquad\qquad\quad-m\,\Phi\,\psi^{(\lambda)\;\dagger}_{s}\psi^{(\lambda)}_{s}\Bigr{]}.$
(29)

Correspondingly, the equation of motion is the Schrödinger Poisson system
where we have a set of Schrödinger field equations for each
$\psi^{(\lambda)}_{s}$, plus the Newtonian Gauss’ law

$\displaystyle i\frac{\partial}{\partial t}\psi^{(\lambda)}_{s}$
$\displaystyle=-\frac{1}{2m}\nabla^{2}\psi^{(\lambda)}_{s}+m\,\Phi\,\psi^{(\lambda)}_{s}$
$\displaystyle\nabla^{2}\Phi$
$\displaystyle=\frac{m}{2m_{\mathrm{pl}}^{2}}\sum_{\lambda}\psi^{(\lambda)\;\dagger}_{s}\psi^{(\lambda)}_{s}.$
(30)

The orbital and spin angular momentum densities (22), in terms of the
$\psi_{s}^{(\lambda)}$ are

$\displaystyle S_{k}$
$\displaystyle=s\,\sum_{\lambda,\lambda^{\prime}}\Re\left(i\,\psi_{s}^{(\lambda^{\prime})\,\dagger}\psi^{(\lambda)}_{s}\,\varepsilon_{ijk}[{{\bm{\mathsf{\epsilon}}}}^{(\lambda)}_{s,\hat{n}}{{\bm{\mathsf{\epsilon}}}}^{(\lambda^{\prime})\;\dagger}_{s,\hat{n}}]_{ij}\right)$
$\displaystyle L_{k}$
$\displaystyle=\sum_{\lambda}\Re\left(i\,\psi_{s}^{(\lambda)\,\dagger}\varepsilon_{ijk}\,\partial_{i}\psi_{s}^{(\lambda)}\,x^{j}\right).$
(31)

Action (III.2), written in terms of the $\psi^{(\lambda)}_{s}$, is helpful to
identify another set of symmetries. We have global U(1) invariance for each of
the $2s+1$ degrees of freedom, giving $2s+1$ conserved particle numbers

$\displaystyle
N^{(\lambda)}=\int\mathrm{d}^{3}x\,\psi^{(\lambda)\;\dagger}_{s}\psi^{(\lambda)}_{s},$
(32)

where $N=\sum_{\lambda}N^{(\lambda)}$. The fact that each of the polarized
fields has an associated conserved particle number, will become important to
physically understand superpositions of extremally polarized solitons in order
to form fractionally polarized ones.

Finally, the conserved energy (III.1) written in terms of
$\psi^{(\lambda)}_{s}$ is

$\displaystyle E$
$\displaystyle=\sum_{\lambda}\Bigg{\\{}\int\mathrm{d}^{3}x\Bigl{[}\frac{1}{2m}\nabla\psi^{(\lambda)\;\dagger}_{s}\cdot\nabla\psi^{(\lambda)}_{s}$
$\displaystyle+\frac{m^{2}}{4m_{\mathrm{pl}}^{2}}\psi^{(\lambda)\;\dagger}_{s}\psi^{(\lambda)}_{s}\sum_{\lambda^{\prime}}\int\frac{\mathrm{d}^{3}y}{4\pi|{\bf
x}-{\bf y}|}\,\psi^{(\lambda^{\prime})\;\dagger}_{s}({\bf
y})\psi^{(\lambda^{\prime})}_{s}({\bf y})\Bigr{]}\Bigg{\\}}.$ (33)

On a conceptual note, the non-relativistic action (III.2) cannot distinguish
between (i) a set of $2s+1$ (spatial-)scalars (ii) a spin-$s$ field with
$2s+1$ spin multiplicity degrees of freedom. While the set of
(spatial-)scalars will have no spin angular momentum, they will still have a
conserved “isospin” (III.2).999The phenomenology with multiple scalar fields,
but with different masses, was explored in [32]. The non-relativistic theory
cannot distinguish between the two cases using gravitational physics alone.
This equivalence can be broken when we include relativistic corrections.

##  IV Polarized Solitons

Making use of the polarization basis (28), we classify different (lowest
energy) solitons based on their spin multiplicities. We first discuss
‘extremally polarized’ solitons (composed of identically polarized plane
waves), where only one of the $2s+1$ polarized fields $\psi_{s}^{(\lambda)}$
is non-zero. Thereafter, we discuss fractionally polarized solitons obtained
by linear superpositions of extremally polarized ones.

###  IV.1 Extremally polarized solitons

With only one of the polarization fields non-zero, we can assume the following
ansatz

$\displaystyle\psi^{(\lambda)}_{s}({\bf x},t)=\psi({\bf x})\,e^{i\mu t}\,.$
(34)

That is, $\psi^{(\lambda)}\neq 0$ for a particular $\lambda$, and zero
otherwise. Here $\mu=const.$ can be thought of as the ‘chemical potential’. In
this case, we have the usual (scalar) Schrödinger-Poisson system

$\displaystyle-\mu\,\psi$
$\displaystyle=-\frac{1}{2m}\nabla^{2}\psi+m\,\Phi\,\psi$
$\displaystyle\nabla^{2}\Phi$
$\displaystyle=\frac{m}{2m_{\mathrm{pl}}^{2}}\,|\psi({\bf x})|^{2},$ (35)

along with the following energy (III.1):

$\displaystyle E$
$\displaystyle=\int\mathrm{d}^{3}x\Bigl{[}\frac{1}{2m}\nabla\psi^{\dagger}\cdot\nabla\psi$
$\displaystyle\qquad\qquad\qquad+\frac{m^{2}}{4m_{\mathrm{pl}}^{2}}\psi^{\dagger}\psi\int\frac{\mathrm{d}^{3}y}{4\pi|{\bf
x}-{\bf y}|}\,\psi^{\dagger}({\bf y})\psi({\bf y})\Bigr{]}.$ (36)

The lowest energy solution to eq. (IV.1) constitutes the well known,
spherically symmetric, scalar soliton configuration [16]. The unique profile
(up to scaling of the amplitudes by $\mu/m$) is shown in Fig. 2.

Figure 2: The “universal” non-relativistic field profile and gravitational
potential for the soliton solutions. Note that both profiles must be
multiplied by $\mu/m$ to get the correct solution for each $\mu/m$. Figure 3:
A visualization of the two distinct extremally polarized vector solitons. The
left soliton has vanishing spin density ($\lambda=0$), and ${\bm{W}}$ is
oscillating along the $z$-axis. The right soliton has a spin density
${\bm{S}}=\lambda|\psi|^{2}\hat{z}$ with $\lambda=1$. The big arrows inside
the soliton represent the direction of the field ${\bm{W}}$, while the little
arrows on the circles represent their motion in time. The total spin
$|{\bm{S}}_{\rm tot}|=\lambda M/m\approx
60.7\,\lambda\,(m_{\mathrm{pl}}/m)^{2}(\mu/m)^{1/2}$, where $M$ is the total
mass of the soliton. Figure 4: The three distinct extremally polarized tensor
solitons. For visualization, we plot the eigenvectors of the traceless,
$3\times 3$ matrix ${\bm{\mathsf{H}}}^{(\lambda)}$ representing the polarized
solitons in the massive spin-$2$ field. The eignevectors are scaled by their
corresponding eigenvalues. The leftmost soliton has a vanishing spin density,
with each eigenvector (along the co-ordinate axes) of
${\bm{\mathsf{H}}}^{(\lambda)}$ oscillating in phase. The middle soliton has a
spin density ${\bm{S}}=\lambda|\psi|^{2}\hat{z}$ with $\lambda=1$. The
rightmost soliton has a spin density ${\bm{S}}=\lambda|\psi|^{2}\hat{z}$ with
$\lambda=2$. The total spin $|{\bm{S}}_{\rm tot}|=\lambda M/m\approx
60.7\,\lambda\,(m_{\mathrm{pl}}/m)^{2}(\mu/m)^{1/2}$ which can be
macroscopically large for $\lambda\neq 0$.

Since $\psi$ is real and $\mu$ spatially independent, the orbital angular
momentum in (III.2) vanishes:

$\displaystyle\bm{L}=\Re\left(i\,\psi\nabla\psi\times\bm{r}\right)=0\,.$ (37)

For higher-spin fields, the spin angular momentum density can be non-
vanishing. For the ansatz (34), it is equal to

$\displaystyle{\bm{S}}=\lambda\psi^{2}\hat{n}\,,\qquad\lambda=\\{-s,...,0,...,s\\}.$
(38)

This reflects the extremally polarized nature of the solitons. For each
$\lambda$, this is a coherent collection of plane waves, all polarized along
the $\hat{n}$ direction. Explicitly, the total spin angular momentum

$\displaystyle{{\bm{S}}}_{\rm tot}=\int d^{3}x{\bm{S}}=\lambda\int
d^{3}x|\psi|^{2}\hat{n}=\lambda N\hat{n}.$ (39)

where $N$ is the total particle number in the soliton.

For extremally polarized solitons, the particle number, energy and spin-
angular momentum are

$\displaystyle N$ $\displaystyle=\frac{M}{m}\approx
60.7\frac{m_{\mathrm{pl}}^{2}}{m^{2}}\left(\frac{\mu}{m}\right)^{\\!1/2}\,,$
(40) $\displaystyle E$
$\displaystyle\approx-19.2\frac{m_{\mathrm{pl}}^{2}}{m}\left(\frac{\mu}{m}\right)^{\\!3/2}\,,$
(41) $\displaystyle{\bm{S}}_{\rm tot}$ $\displaystyle\approx\lambda\times
60.7\frac{m_{\mathrm{pl}}^{2}}{m^{2}}\left(\frac{\mu}{m}\right)^{\\!1/2}\hat{n}\,,$
(42)

where $M$ is the total mass of the soliton. The numerical co-efficients are
obtained from the universal profile shown in Fig. 2. Heuristically
$(\mu/m)^{1/2}\sim(1/mL)$ where $L$ is the characteristic size of the soliton,
and $M\sim(m_{\mathrm{pl}}/m)^{2}/L$. Since in the non-relativistic regime
$\mu/m\ll 1$, we expect the maximal values of the above quantities (40) to be
bounded from above by the case $\mu\sim m$. Significant deviations from the
above expression can be expected as one approaches this limit [13].

It is also worth noting that these solitons are perfectly virialized, $E_{\rm
kin}/E_{\rm pot}=-1/2$, where $E_{\rm kin}$ is the term containing gradients,
and $E_{\rm pot}$ is the other term (gravitational potential energy) in
(III.1).

####  IV.1.1 Spin-$0$

For the case of spin-$0$, we have the following real field solution for
$\phi(=\bm{\mathcal{F}})$ in (15)

$\displaystyle\phi({\bf x},t)=\frac{\sqrt{2}\psi({\bf x})}{\sqrt{m}}\cos\omega
t\,,$ (43)

where $\omega\equiv m-\mu$, and there is of-course no intrinsic spin angular
momentum.

####  IV.1.2 Spin-$1$

For the massive spin-$1$ case, we have three distinct states corresponding to
$\pm 1$ and $0$ polarizations, which, for $\hat{n}=\hat{z}$, are conveniently
represented by the following orthonormal set of vectors:

$\displaystyle{\bm{\mathsf{\epsilon}}}^{(\pm
1)}_{1,\hat{z}}=\frac{1}{\sqrt{2}}\begin{pmatrix}1\\\ \pm i\\\
0\end{pmatrix}\;;$
$\displaystyle\qquad{\bm{\mathsf{\epsilon}}}^{(0)}_{1,\hat{z}}=\begin{pmatrix}0\\\
0\\\ 1\end{pmatrix}.$ (44)

satisfying (27). For ${\bm{\mathsf{\Psi}}}^{(\lambda)}=\psi e^{i\mu
t}{\bm{\mathsf{\epsilon}}}_{1,\hat{z}}^{(\lambda)}$, we have
${\bm{S}}=\lambda|\psi|^{2}\hat{z}$ where $\lambda=0,\pm 1$. Extremally
polarized solitons in terms of the real-valued vector field
${{\bm{W}}}(=\bm{\mathcal{F}})$ in (15) are

$\displaystyle{{\bm{W}}}^{(\pm 1)}({\bf x},t)$ $\displaystyle=\frac{\psi({\bf
x})}{\sqrt{m}}\begin{pmatrix}\cos\omega t\\\ \pm\sin\omega t\\\
0\end{pmatrix}\,,$ $\displaystyle{\bm{W}}^{(0)}({\bf x},t)$
$\displaystyle=\frac{\sqrt{2}\,\psi({\bf x})}{\sqrt{m}}\begin{pmatrix}0\\\
0\\\ \cos\omega t\end{pmatrix}.$ (45)

In fig. 3 we show these extremally polarized solitons.

We note that in [14], the authors provide ground state solitons
${\bm{\mathsf{\Psi}}}=\psi(r)e^{i\mu t}\\{w_{x},w_{y},w_{z}\\}$ where $w_{i}$
are components of a complex unit vector. For the extremally polarized cases,
this corresponds to the choices (44) for their $w_{i}$. Our focus on the spin
aspect of fields dictated this choice, which naturally leads to extremally
polarized solitons. More general solitons with arbitrary $w_{i}$ are discussed
in IV.2.

####  IV.1.3 Spin-2

For the massive spin-$2$ case we have 5 polarization states corresponding to
spin multiplicities $\pm 2$, $\pm 1$, and $0$. Again, for $\hat{n}=\hat{z}$,
these are represented by the following orthonormal (and trace free) set of
tensors101010These can be obtained through tensor products of the spin-1
polarization vectors: $\epsilon^{(\pm 2)}_{2,\hat{z}}=\epsilon^{(\pm
1)}_{1,\hat{z}}\odot\epsilon^{(\pm 1)}_{1,\hat{z}}$; $\epsilon^{(\pm
1)}_{2,\hat{z}}=(\epsilon^{(\pm
1)}_{1,\hat{z}}\odot\epsilon^{(0)}_{1,\hat{z}}+\epsilon^{(0)}_{1,\hat{z}}\odot\epsilon^{(\pm
1)}_{1,\hat{z}})/\sqrt{2}$; and
$\epsilon^{(0)}_{2,\hat{z}}=(2\epsilon^{(0)}_{1,\hat{z}}\odot\epsilon^{(0)}_{1,\hat{z}}-\epsilon^{(+1)}_{1,\hat{z}}\odot\epsilon^{(-1)}_{1,\hat{z}}-\epsilon^{(-1)}_{1,\hat{z}}\odot\epsilon^{(+1)}_{1,\hat{z}})/\sqrt{6}$.

$\displaystyle{\bm{\mathsf{\epsilon}}}^{(\pm 2)}_{2,\hat{z}}$
$\displaystyle=\dfrac{1}{2}\begin{pmatrix}1&&\pm i&&0\\\ \pm i&&-1&&0\\\
0&&0&&0\end{pmatrix}$ $\displaystyle{\bm{\mathsf{\epsilon}}}^{(\pm
1)}_{2,\hat{z}}$ $\displaystyle=\dfrac{1}{2}\begin{pmatrix}0&&0&&1\\\
0&&0&&\pm i\\\ 1&&\pm i&&0\end{pmatrix}$
$\displaystyle{\bm{\mathsf{\epsilon}}}^{(0)}_{2,\hat{z}}$
$\displaystyle=\dfrac{1}{\sqrt{6}}\begin{pmatrix}-1&&0&&0\\\ 0&&-1&&0\\\
0&&0&&2\end{pmatrix}$ (46)

For ${\bm{\mathsf{\Psi}}}^{(\lambda)}=\psi e^{i\mu
t}{\bm{\mathsf{\epsilon}}}_{2,\hat{z}}^{(\lambda)}$, we have the spin density
${\bm{S}}=\lambda|\psi|^{2}\hat{z}$ where $\lambda=0,\pm 1,\pm 2$.

The five extremally polarized solitons in the real-valued trace-free tensor
field are:

$\displaystyle{{\bm{\mathsf{H}}}}^{(\pm 2)}({\bf x},t)$
$\displaystyle=\frac{\psi({\bf x})}{\sqrt{2m}}\begin{pmatrix}\cos\,\omega
t&\pm\sin\,\omega t&0\\\ \pm\sin\,\omega t&-\cos\,\omega t&0\\\
0&0&0\end{pmatrix}$ $\displaystyle{{\bm{\mathsf{H}}}}^{(\pm 1)}({\bf x},t)$
$\displaystyle=\frac{\psi({\bf x})}{\sqrt{2m}}\begin{pmatrix}0&0&\cos\,\omega
t\\\ 0&0&\pm\sin\,\omega t\\\ \cos\,\omega t&\pm\sin\,\omega t&0\end{pmatrix}$
$\displaystyle{{\bm{\mathsf{H}}}}^{(0)}({\bf x},t)$
$\displaystyle=\frac{\psi({\bf x})}{\sqrt{3m}}\begin{pmatrix}-1&&0&&0\\\
0&&-1&&0\\\ 0&&0&&2\end{pmatrix}\,\cos\,\omega t,$ (47)

Fig. 4 shows these extremally polarized solitons (with spin density along the
$z$ axis).

###  IV.2 Fractionally polarized solitons

Here we construct non-extremal polarized solitons, obtained through linear
superpositions of the extremal ones. Since we have U(1) invariance for each
polarization field $\psi_{s}^{(\lambda)}$, we can superpose them to form new
solutions. That is to say that we can have

$\displaystyle{\bm{\mathsf{\Psi}}}({\bf x},t)$ $\displaystyle=\psi({\bf
x})\sum_{\lambda}c_{\lambda}e^{i(\mu
t-\varphi_{\lambda})}{{\bm{\mathsf{\epsilon}}}}^{(\lambda)}_{s,\hat{n}}$ (48)

with

$\displaystyle\qquad\sum_{\lambda}c^{2}_{\lambda}=1.$ (49)

The reason to enforce (49) is to ensure that we have the same particle number
density as in the case of extremally polarized solitons (and hence also the
same energy), which in turn guarantees that $\psi({\bf x})$ obeys the same
Schrödinger Poisson system (IV.1).

The corresponding real valued counterpart is

$\displaystyle\bm{\mathcal{F}}({\bf
x},t)=\sum_{\lambda}c_{\lambda}\bm{\mathcal{F}}^{(\lambda)}({\bf
x},t+\varphi_{\lambda}/\omega)\,,$ (50)

where $\bm{\mathcal{F}}^{(\lambda)}$ are extremally polarized solitons (equal
to ${{\bm{W}}}^{(\lambda)}$ for spin-$1$ (IV.1.2);
${{\bm{\mathsf{H}}}}^{(\lambda)}$ for spin-$2$ (IV.1.3)).

Similar to the case of extremally polarized solitons, the orbital angular
momentum density, obtained by substituting (48) into (III.2) is zero. On the
other hand the spin density is

$\displaystyle S_{k}=s|\psi({\bf
x})|^{2}\sum_{\lambda\lambda^{\prime}}\Re\left[ic_{\lambda}c_{\lambda^{\prime}}e^{i(\varphi_{\lambda}-\varphi_{\lambda^{\prime}})}\epsilon_{ijk}[{\bm{\mathsf{\epsilon}}}_{s,\hat{n}}^{(\lambda^{\prime})}{\bm{\mathsf{\epsilon}}}_{s,\hat{n}}^{(\lambda)\dagger}]_{ij}\right].$
(51)

The total spin need not be equal to $\lambda N$ where
$\lambda\in\\{-s,..,s\\}$ for general $\\{c_{\lambda},\varphi_{\lambda}\\}$.
However, for a class of $\\{c_{\lambda},\varphi_{\lambda}\\}$ which simply
amount to rotations of extremally polarized solitons, we will again get
$\lambda N$ for the total spin.

Superposing basis solutions (extremally polarized) to form new solutions,
while keeping the total particle number fixed, is just taking different
fractions of these extremally polarized solitons and putting them on top of
each other (c.f. (50)). The reason this is allowed is because there is a U(1)
invariance within each polarization sector. From a phenomenological point of
view, there could be extra spin induced interactions that favor same
polarization states for a collection of particles in some region (e.g. Ising
model). For such situations, it may be that there is a higher chance of
extremally polarized solitons to form over fractionally polarized ones due to
Bose-Einstein statistics. We leave such questions for future work.

###  IV.3 $(s/2,s/2)$ representation of polarized solitons

There exists a simple understanding of the space of all polarized solitons in
terms of the $(m,n)$ representation of the rotation group SO(3). Guided by the
U(1) and SO(3) invariance of the general action (III.2), and the fact that
there are $2s+1$ d.o.f. with $\lambda\rightarrow-\lambda$ obtainable via
simple rotation, we can represent our full space of soliton solutions via
$(s/2,s/2)$ representation of the SO(3) group. The representation is
reducible, meaning it admits a total of $s+1$ SO(3) invariant sub-spaces, each
containing the different (absolute value of-) polarizations
$|\lambda|\in\\{0,1,..s\\}$. These are our distinct extremally polarized
solitons. These $s+1$ extremally polarized solitons form a basis, and can give
rise to fractionally polarized solitons via appropriate superpositions (due to
separate U(1) within each polarization sector). The full space of solutions is
the product space of these $s+1$ subspaces, hosting a general soliton. A
physical soliton spontaneously breaks the SO(3) invariance of the action
(III), (or equivalently (III.2)). Table 1 shows the representation of
gravitationally bound polarized solitons in bosonic fields considered in this
paper.

Table 1: $(m,n)$ representation of the gravitationally bound soliton states in non-relativistic integer spin field theories  Soliton states |  Massive field (spin) |
---|---|---
(0,0) | scalar (0) |
(1/2,1/2) | vector (1) |
(1,1) | tensor (2) |
Figure 5: The left panel shows collisions between vector solitons that can be
replicated by solitons in a single scalar field. The right panel shows
examples of collisions which cannot be replicated by solitons in a single
scalar field.

###  IV.4 Beyond ground-state solitons

Spherically symmetric, single node solitons for the spin-$1$ and spin-$2$
fields are the ‘hedgehog’-like configurations, with Cartesian components [14,
15]111111For relativistic hedgehogs in complex-valued Proca-fields, see [33].

$\displaystyle W_{j}({\bf x},t)$ $\displaystyle=f(r)\frac{x^{j}}{r}\cos\omega
t\,,$ $\displaystyle H_{ij}({\bf x},t)$
$\displaystyle=g(r)\left(3\frac{x^{i}x^{j}}{r^{2}}-\delta_{ij}\right)\cos\omega
t\,,$ (52)

where $f(0)=g(0)=0$. That is, there is a node in the profile at the origin.
Both hedgehogs have higher energies (at a fixed particle number) compared to
the ground state solitons discussed earlier, and have zero spin and orbital
angular momentum. Explicitly, after fixing the particle number to be the same
as the polarized solitons ($N^{s}_{\rm hh}=N$), we have $E^{s=1}_{\rm
hh}\approx 0.33E$ and $E^{s=2}_{\rm hh}\approx 0.17E$ where $E<0$ and $N$ are
given in (40). Note that $\Delta E=E^{s}_{\rm hh}-E>0$. A linear stability
analysis was provided in [15] to argue that the hedgehogs in spin-$2$ case are
unstable and might transition to $p$-solitons. As with scalar solitons,
excited configurations with additional nodes and orbital angular momentum
might be possible with higher-spin fields, albeit with shorter lifetimes.

So far we have only allowed for spherically symmetric energy densities. It is
possible to construct non-spherically symmetric configurations, such as domain
walls, strings/vortices etc. [14]; the possible space of extended field
configurations with higher-spin fields is likely to be quite rich. The full
classification is beyond the scope of the present paper, but it is worth
pursuing since it might provide new avenues to probe these higher-spin fields.

##  V Distinguishability & probes of polarized solitons

Having shown that we have quite a rich space of soliton solutions, we briefly
discuss some of the phenomenological implications. Alongside these
implications, we address some conceptual questions: Can higher-spin solitons
be distinguished from scalar solitons? Can solitons with different
polarizations be distinguished using only gravitational interactions?

###  V.1 Gravitational interactions

Let us consider collisions between solitons $A$ and $B$ in a spin-$s$ field.
We show below that only if the two solitons differ by just an overall phase,
can the collision be mimicked by two scalar solitons. Otherwise, in general,
the higher-spin nature of the fields will leave an imprint in the observables
related to the collision of the two solitons.

For simplicity, let us consider a collision between two extremally polarized
solitons, initially far away from each other, such that the field admits the
following ansatz

$\displaystyle{\bm{\mathsf{\Psi}}}_{s}\Bigr{|}_{t=0}=\psi_{A}\left({\bf
x}+{\bf
x}_{0}\right){\bm{\mathsf{\epsilon}}}_{s,\hat{n}}^{(\lambda)}+e^{i\theta}\psi_{B}\left({\bf
x}-{\bf
x}_{0}\right){\bm{\mathsf{\epsilon}}}_{s,\hat{n}^{\prime}}^{(\lambda^{\prime})}\,.$
(53)

The corresponding initial number and current densities are

$\displaystyle\mathcal{N}_{s}\Big{|}_{t=0}$
$\displaystyle=\psi_{A}^{2}+\psi_{B}^{2}+\psi_{A}\psi_{B}\left(e^{i\theta}\operatorname{Tr}[{\bm{\mathsf{\epsilon}}}^{(\lambda^{\prime})}_{s,\hat{n}^{\prime}}{\bm{\mathsf{\epsilon}}}^{(\lambda)\dagger}_{s,\hat{n}}]+c.c.\right)$
$\displaystyle{\bm{j}}_{s}\Big{|}_{t=0}$
$\displaystyle=\frac{i}{2m}\Bigl{\\{}\left(e^{-i\theta}\operatorname{Tr}[{\bm{\mathsf{\epsilon}}}^{(\lambda)}_{s,\hat{n}}{\bm{\mathsf{\epsilon}}}^{(\lambda^{\prime})\dagger}_{s,\hat{n}^{\prime}}]-c.c.\right)\psi_{B}\nabla\psi_{A}$
$\displaystyle\qquad+\left(e^{i\theta}\operatorname{Tr}[{\bm{\mathsf{\epsilon}}}^{(\lambda^{\prime})}_{s,\hat{n}^{\prime}}{\bm{\mathsf{\epsilon}}}^{(\lambda)\dagger}_{s,\hat{n}}]-c.c.\right)\psi_{A}\nabla\psi_{B}\Bigr{\\}}.$
(54)

For the scalar case ($s=0$,
${\bm{\mathsf{\epsilon}}}_{0,\hat{n}}^{(\lambda)}=1$), we have

$\displaystyle\mathcal{N}_{0}\Big{|}_{t=0}$
$\displaystyle=\psi_{A}^{2}+\psi_{B}^{2}+2\psi_{A}\psi_{B}\cos\theta$ (55)
$\displaystyle\bm{j}_{0}\Big{|}_{t=0}$
$\displaystyle=\frac{\sin\theta}{m}\left(\psi_{B}\nabla\psi_{A}-\psi_{A}\nabla\psi_{B}\right).$
(56)

Now to compare the collision scenario with the scalar case, we wish to start
with equal number densities in the scalar and higher-spin cases. That is,
$\mathcal{N}_{s}=\mathcal{N}_{0}$ at $t=0$. To achieve this condition we have
two possibilities (1)
$\operatorname{Tr}[{\bm{\mathsf{\epsilon}}}^{(\lambda)}_{s,\hat{n}}{\bm{\mathsf{\epsilon}}}^{(\lambda^{\prime})\dagger}_{s,\hat{n}^{\prime}}]\neq
1$ and real, with $\theta=\pi/2$; (2)
$\operatorname{Tr}[{\bm{\mathsf{\epsilon}}}^{(\lambda)}_{s,\hat{n}}{\bm{\mathsf{\epsilon}}}^{(\lambda^{\prime})\dagger}_{s,\hat{n}^{\prime}}]=1$,
that is $\lambda=\lambda^{\prime}$ and $\hat{n}=\hat{n}^{\prime}$. For
possibility (1), ${\bm{j}}_{s}=\vec{0}$ while ${\bm{j}}_{0}\neq\vec{0}$ at
$t=0$, so the collision will proceed differently in the scalar vs. the higher-
spin case. Specifically, the densities, and hence the gravitational potential,
will also evolve differently in the two cases. Possibility (2) essentially
reduces to the scalar case as expected. This result is summarized in Fig. 5
(as an example for spin-$1$ fields).

Thus, using collisions we can tell whether two solitons in a given spin-$s$
field had the same spin density or not. This would not be possible with
individual solitons in isolation using the gravitational potential alone. We
leave a detailed analysis of how exactly the potential evolves and the
associated phenomenological avenues such as motion of test particles in this
dynamical potential [34] or merger rates of solitons [35] for future
work.121212With R. Karur and P. Mocz, we are investigating soliton
interactions using numerical simulations of the multi-component SP system.

More generally, many existing studies carried out for scalar dark matter, can
be repeated for the case of higher-spin dark matter including the spectra of
associated cosmological perturbations, soliton formation mechanisms [36], halo
formation [35], dynamical friction [37], generation of gravitational waves
[38], time delays caused in pulsar timing arrays by time-dependent pressures
and energy densities [39], black-hole superradiance and ‘gravitational
atoms’[40], transient vortices [41] etc. Some of these have already been
extended to higher-spin fields, see for example, [42, 43, 44, 45].

###  V.2 Non-gravitational interactions

One can envision many scenarios where the higher-spin fields couple to the
Standard Model such that it opens up channels for polarized solitons to be
detectable. One possibility is to have some high energy scale(s)
$g_{\mathcal{F}\gamma}^{-1}$, such that in the low-energy effective theory we
have:

$\displaystyle\mathcal{L}_{int}$
$\displaystyle\sim\begin{cases}g^{2}_{W\gamma}W_{\mu}W^{\mu}\,F_{\alpha\beta}\tilde{F}^{\alpha\beta}&{\rm
spin-}1\\\
g^{2}_{H\gamma}(H_{\mu\nu}H^{\mu\nu}-H^{2})\,F_{\alpha\beta}\tilde{F}^{\alpha\beta}&{\rm
spin-}2\end{cases}\;.$ (57)

In the non-relativistic limit, such interactions reduce to

$\displaystyle\mathcal{L}_{int}\sim
g^{2}_{\mathcal{F}\gamma}\operatorname{Tr}[{\bm{\mathcal{F}}}{\bm{\mathcal{F}}}]\,F_{\alpha\beta}\tilde{F}^{\alpha\beta}\,,$
(58)

where $\bm{\mathcal{F}}={\bm{W}}$ and ${\bm{\mathsf{H}}}$ for spin-$1$ and
spin-$2$ cases respectively. Here, $F_{\mu\nu}$ is the electromagnetic field
strength tensor, and $\tilde{F}_{\mu\nu}$ is its dual. Such couplings may have
similar phenomenological implications like the axion photon case [46, 47], but
with an important difference due to the polarization state of the soliton.
Specifically, $\operatorname{Tr}[{\bm{\mathcal{F}}}{\bm{\mathcal{F}}}]$ is
time-independent for maximally spinning solitons ($\lambda\neq 0)$, but is
time-dependent (periodic) for the $\lambda=0$ case and also for fractionally
polarized solitons. Besides emission from a single soliton, a collision
between any two non-zero spin solitons may result in radiation that has a
specific polarization pattern, depending upon their polarization (for the
spin-0 soliton collisions, see for example [48, 49, 50]).

The above interactions are CP odd. One can have CP even interactions as well,
by replacing $F\tilde{F}$ by $F^{2}$ or other contractions between $F$ and our
higher-spin fields. However, one needs to be careful to avoid issues of ghosts
and superluminality.131313We would like to thank M. P. Hertzberg for bringing
these issues to our attention. Another possibility for the massive spin-$1$
case could be to have kinetic mixing with the usual electromagnetism $\sim
F_{\mu\nu}G^{\mu\nu}$.141414We thank Andrew Long for this suggestion.

##  VI Summary and future directions

Massive, scalar (spin-$0$), vector (spin-$1$) or tensor fields (spin-$2$) can
make up all/part of the dark matter, or play a role in the early universe. In
this paper, we explored the non-relativistic limit of such fields, and derived
polarized ground-state solitons in such fields. We summarize our main results
below, and also discuss possible future directions.

###  VI.1 Summary

Non-relativistic Limit: Starting with the quadratic (free) action for massive
spin-$0$,$1$,$2$ fields + leading gravitational interactions, we derived an
effective action and equations of motion in the non-relativisitic regime. We
arrived at the non-relativistic action by first using the constraint equations
and deriving a quadratic action for physical degrees of freedom in the massive
spin-$s$ field. Our non-relativistic system for a field with spin $s$, is a
multi-component Schrödinger-Poisson (SP) system with $2s+1$ independent d.o.f.

$\bm{p}$-Solitons: Using a polarization basis, we enumerate the full space of
lowest energy soliton solutions in the SP system for nonzero-spin fields,
which we call $p$-solitons. These include extremally polarized solitons which
are coherent collections of identically polarized plane waves. The extremally
polarized solitons can be used as a basis set to construct fractionally
polarized solitons via superposition. All $p$-solitons have a spherically
symmetric (universal) energy density profile, although the field configuration
is not spherically symmetric.

The orbital angular momentum is zero for $p$-solitons, but spin angular
momentum need not be. In the extremal case, the spin angular momentum is
simply the spin multiplicity $\lambda$ times the macroscopically large
particle number of the soliton. Explicitly

$\displaystyle{\bm{S}}_{\rm tot}=\lambda\frac{M}{m}\,\,\hat{n}\,,$ (59)

where $M\gg m$ is the total mass of the soliton and $m$ is the mass of the
dark spin-$s$ field, and the above expression is in units of $\hbar$. The
instrinsic spin of our solitons can also have implications for substructure in
dark matter today, as well as for the formation of primordial black holes,
baryogenesis in the early Universe, etc.

For comparison, the average value of a dimensionless (orbital) angular
momentum parameter, $\Lambda\equiv J\sqrt{|E|}/GM^{5/2}$, which is used to
quantify the angular momentum of dark matter halos in N-body simulations, is
$\approx 0.045$ [51]. For our $p$-solitons, assuming a uniform distribution of
the magnitude of the total spin, we have $\langle\Lambda\rangle=\Lambda_{\rm
max}/2=0.12s$. Here, $\Lambda_{\rm max}$ is the value for an extremally
polarized soliton and $s$ is the spin of the field. For another comparison,
the ratio of the spin angular momentum of our extremal soliton to that of an
equal mass black hole is $S_{\rm tot}/J_{\rm bh}=8\pi\lambda
m_{\mathrm{pl}}^{2}/(aMm)$ with $a<1$. This ratio can be larger than unity
when $M$ is sufficiently small.151515We used black-hole angular momentum
$J_{\rm bh}=aGM^{2}$.

We also compared the energy, particle number and spin of $p$-solitons with
other higher energy, but zero total angular momentum solitons with hedgehog
configurations of the fields.

Distinguishability and Implications: We argued that it is possible to
distinguish two interacting $p$-solitons from their scalar counterparts.
Limiting ourselves to gravitational interactions in the non-relativistic
limit, we discussed how collisions between two higher-spin $p$-solitons cannot
be mimicked by two scalar solitons unless the $p$-solitons are identically
polarized. Such collisions can in principle be probed by test particles,
including photons. We highlighted how non-gravitational couplings to photons
can offer further detectability avenues via specific polarization patterns of
the outgoing electromagnetic radiation generated by, for example, $p$-soliton
collisions.

###  VI.2 Future Directions

Formation mechanisms: Formation mechanisms of solitons in scalar fields have
been studied before both with and without self-interactions [52, 53, 36]; for
higher-spin fields this remains an open question. Coupling the higher-spin
fields to a coherently oscillating scalar can lead to efficient production of
such fields (for example, see [54, 55, 56]). Higher-spin fields can also be
populated via gravitational particle production (albeit with large masses [57,
58]). Solitons could form via self-interactions, or via gravitational
clustering in these or in ultra-light fields. We hope to carry out a numerical
investigation of the cosmic history of dark, nonzero-spin fields and formation
mechanism of their solitons, which would allow a statistical prediction of
their their masses and spins in the early or contemporary universe. For a
recent analytic approach for the vector case, see [59].

Extension to higher spin fields: In general, an integer spin-$s$ degree of
freedom in embedded in a rank $s$ field (meaning an object carrying $s$ space-
time indices). Without mass, the fields are gauge fields (for $s\geq 1$) and
depending upon the spin, admit a rich gauge invariance (see [60, 61] and
references therein for a classic review). The gauge structure makes sure that
there are only two physical $\pm s$ spin multiplicity states. Another
equivalent way to construct the Lagrangian density is to add all the possible
Lorentz invariant kinetic terms in a way that leads to no ghosts, with only
two degrees of freedom surviving (similar to the case of spin-$1$ and
spin-$2$). Then, add all possible Lorentz invariant mass terms such that no
ghosts, but only the rest $2s-1$ spin multiplicity states now appear. Even
though we don’t carry out the exercise explicitly for a full Lorentz invariant
relativistic theory, we expect that in the non-relativistic and weak field
gravity regime, the action is still (III.2). We have kept most of our
expressions in this paper as general as possible, and only in IV have we
considered the specific cases of spin-$1$ and spin-$2$ solitons. We leave a
detailed study of a fully Lorentz invariant spin-3 and higher fields for
future work.

Relativistic corrections: In this paper we restricted ourselves to the non-
relativistic regime, and found $s+1$ distinct, but degenerate (in energy)
soliton configurations for a spin-$s$ bosonic field. The fate of these
solutions when relativistic corrections are added, especially for spin-$1$ and
spin-$2$, is an important question to consider (for the spin-0 case, see
[13]). In the non-relativistic limit, we essentially have $2s+1$ spin
multiplicity fields with a conserved particle number within each, and thus all
of the different $p$-solitons were degenerate and stable. However at sub-
leading order, the energy momentum tensor (especially for the cases of
spin-$1$ and higher fields) contains non-trivial components due to the spin
structure which we do not focus on in this paper. This might lift the
degeneracy between the $s+1$ solitons, and affect their stability. In
particular, gravitational waves (c.f. (II.4)) can be sourced at higher orders,
which, depending upon the spin of the $p$-soliton, may differ in their
polarization patterns.

Self-interactions: It is also possible to include non-gravitational self-
interactions within our framework. Attractive self-interactions (even without
gravity), are known to support metastable solitons (oscillons) in the scalar
case (see, for example, [62, 63]). For the vector case, we pursue the
investigation of oscillons held together by self-interactions alone elsewhere
[21], and find similar $s+1$ stable solitons. An important difference we find
is that there is no conserved particle number within each spin-multiplicity
sector. As a result, the extremally polarized solutions are not degenerate in
energy, and cannot be superposed to form new additional ones, as we were able
to do in the gravitationally supported case.

In the scalar and vector case, it is possible to find regimes in which either
self-interaction, or gravity alone, or a combination of both can be used to
support solitons. For the case of massive spin-$2$, self-interactions seem
unavoidable and the only ghost free possibility is the full non-linear bi-
gravity theory. These self-interactions may become relevant beyond the non-
relativistic regime, and must be studied carefully to analyze the fate and
lifetimes of our polarized solitons within the full (bi-gravity) theory. We
leave a detailed analysis of such questions for future work.

###### Acknowledgements.

We thank Siyang Ling, Andrew Long and Rachel Rosen for many useful
conversations regarding bi-gravity, as well as Peter Adshead, Mark Hertzberg,
Kaloian Lozanov, Zong-Gang Mou and Hong-Yi Zhang for helpful discussions and
comments on the paper. We are supported by a NASA ATP theory grant NASA-ATP
Grant No. 80NSSC20K0518We.

## Appendix A

In this appendix, we provide the full non-linear actions for spin-$0$,
spin-$1$ and spin-$2$ fields coupled to Einstein gravity. In particular,
additional details are provided for some of the issues relevant for a dark,
massive spin-$2$ field.

###  A.1 Spin-$0$

In this case, the full nonlinear actions are

$\displaystyle S_{\rm EH}+S_{\rm
vis}=\int\mathrm{d}^{4}x\,\sqrt{-g}\left[-\frac{m_{\mathrm{pl}}^{2}}{2}R+\mathcal{L}(g,\psi_{\rm
SM})\right],$ $\displaystyle S_{\rm
dark}=\int\mathrm{d}^{4}x\sqrt{-g}\left[\frac{1}{2}g^{\mu\nu}\nabla_{\mu}\phi\,\nabla_{\mu}\phi-V(\phi)\right]+...$
(60)

where $R$ is the Ricci scalar, and
$V(\phi)=(1/2)m^{2}\phi^{2}+\lambda_{3}\phi^{3}+\ldots$. Expanding the metric
around a flat background

$\displaystyle g_{\mu\nu}=\eta_{\mu\nu}+h_{\mu\nu}\,,$ (61)

and including only leading gravitational interactions, we arrive at the action
(II.1) for $h_{\mu\nu}$ and $\phi$ provided in the main text. There, we also
ignored terms beyond quadratic order in $V(\phi)$. These can be included, or
ignored, without much difficulty. We can also consider an FLRW background with
$\eta_{\mu\nu}$ replaced by and FLRW metric. In this case we must either allow
for a background $\phi$ (spatially-independent) or another homogeneous source
for the Friedmann equations.

###  A.2 Spin-$1$

Our actions $S_{\rm EH}$ and $S_{\rm vis}$ are same as in (A.1), whereas the
action for the spin-$1$ field is

$\displaystyle S_{\rm dark}$
$\displaystyle=\int\mathrm{d}^{4}x\sqrt{-g}\Bigl{[}-\frac{1}{4}g^{\mu\alpha}g^{\nu\beta}\,G_{\mu\nu}G_{\alpha\beta}$
$\displaystyle\qquad\qquad\qquad\qquad-V(W_{\mu}W^{\mu})\Bigr{]}+...$ (62)

where $G_{\mu\nu}=\partial_{\mu}W_{\nu}-\partial_{\nu}W_{\mu}$ and
$V(W_{\mu}W^{\mu})=(1/2)m^{2}W_{\mu}W^{\mu}-\lambda_{4}(W_{\mu}W^{\mu})^{2}+\ldots$
could be self interactions arising in the effective low energy theory.
Expanding the metric as in (61), and including only the quadratic part of $V$,
yields the action (II.2) for $h_{\mu\nu}$ and $W_{\mu}$ used in the main text.
Again, we can also include an FLRW background.

###  A.3 Spin-$2$

For the case of a massive spin-$2$ field, the only ghost-free possibility
seems to be the full non-linear bi-metric theory [22, 23, 24, 25, 26, 27, 28]
in which there are two metrics, $\mathfrak{g}$ and $\mathfrak{f}$. Composite
fluctuations in these metrics give rise to massless and massive spin-$2$
degrees of freedom. Denoting $S_{\mathfrak{g},\mathfrak{f}}=S_{\rm EH}+S_{\rm
dark}$, the bimetric action is

$\displaystyle S_{\mathfrak{g},f}$
$\displaystyle=-\frac{m_{\mathrm{pl}}^{2}}{2(1+\alpha^{2})}\int\mathrm{d}^{4}x\Bigl{\\{}\sqrt{-\mathfrak{g}}\,R(\mathfrak{g})+\alpha^{2}\,\sqrt{-\mathfrak{f}}\,R(\mathfrak{f})$
$\displaystyle\qquad\qquad\qquad\qquad-\frac{2\alpha^{2}}{1+\alpha^{2}}m^{2}\,\sqrt{-\mathfrak{g}}\,\mathcal{V}(X)\Bigr{\\}},$
(63)

where $\alpha$ is a constant and
$X^{\mu}_{\;\;\sigma}X^{\sigma}_{\;\;\nu}=\mathfrak{g}^{\mu\rho}\mathfrak{f}_{\rho\nu}$.
The potential
$\mathcal{V}=\beta_{0}+\beta_{1}\,X^{\mu}_{\;\;\mu}+\beta_{2}\,X^{\mu}_{\;\;[\mu}\,X^{\nu}_{\;\;\nu]}+\beta_{3}\,X^{\mu}_{\;\;[\mu}\,X^{\nu}_{\;\;\nu}\,X^{\alpha}_{\;\;\alpha]}+\beta_{4}\,X^{\mu}_{\;\;[\mu}\,X^{\nu}_{\;\;\nu}\,X^{\alpha}_{\;\;\alpha}\,X^{\beta}_{\;\;\beta]}$
where $\beta_{i}$ are constants related by $\beta_{1}+2\beta_{2}+\beta_{3}=1$.
We expand both metrics around the flat-spacetime background as161616Note that
the metric $\mathfrak{g}$ contains both the massless and massive spin-$2$
excitations (not to be confused with the usual GR metric $g$, appearing in the
previous two cases).

$\displaystyle\mathfrak{g}_{\mu\nu}$
$\displaystyle=\eta_{\mu\nu}+\frac{1}{m_{\mathrm{pl}}}h_{\mu\nu}-\frac{\alpha\sqrt{2}}{m_{\mathrm{pl}}}H_{\mu\nu}\,,$
$\displaystyle\mathfrak{f}_{\mu\nu}$
$\displaystyle=\eta_{\mu\nu}+\frac{1}{m_{\mathrm{pl}}}h_{\mu\nu}+\frac{\sqrt{2}}{\alpha\,m_{\mathrm{pl}}}H_{\mu\nu}\,,$
(64)

and obtain quadratic (free) actions for $h$ and $H$ \+ leading interactions,
see eq. (II.3) used in the main text.

We note that in order for $\mathfrak{g}$ and $\mathfrak{f}$ to both have the
flat background ansatz (eq. (A.3)), we must require the background Einstein
equations for both to be consistent with this choice. The Ricci terms
($\sqrt{-\mathfrak{g}}\,R(\mathfrak{g})$ and
$\sqrt{-\mathfrak{f}}\,R(\mathfrak{f})$) are consistent. However the potential
$\mathcal{V}$ produces cosmological constant terms with
$\Lambda_{\mathfrak{g}}=\alpha^{2}m_{\mathrm{pl}}^{2}(\beta_{0}+3\beta_{1}+3\beta_{2}+\beta_{3})/(1+\alpha^{2})$
and
$\Lambda_{\mathfrak{f}}=m_{\mathrm{pl}}^{2}(\beta_{1}+3\beta_{2}+3\beta_{3}+\beta_{4})/(1+\alpha^{2})$
respectively in the two background equations. We set these to zero to be
consistent with our choice (A.3).171717Within bi-gravity (without any matter
couplings) with equal backgrounds, the massless $h_{\mu\nu}$ sector alone re-
organizes to the usual GR + cosmological constant
$\Lambda=\Lambda_{\mathfrak{g}}=\Lambda_{\mathfrak{f}}$ [64].

In general there is also a whole tower of self interactions in the massive
spin-$2$ sector. Parametrically, the vertices go like

$\displaystyle\sim\frac{m^{2}}{\alpha^{n}m_{\mathrm{pl}}^{n-2}}H^{n}\quad\mathrm{and}\quad\sim\frac{m^{2}\alpha^{n}}{m_{\mathrm{pl}}^{n-2}}H^{n}\quad;\quad
n\geq 3.\,$ (65)

In particular, the trilinear interaction vertex ($n=3$ above), is the only
coupling that parametrically goes like the leading gravitational coupling
($h_{\mu\nu}\,T^{\mu\nu}(H)/m_{\mathrm{pl}}\sim
h\,H^{2}m^{2}/m_{\mathrm{pl}}$). All (non-zero) higher order couplings are
suppressed by extra factors of $m_{\mathrm{pl}}$. Notice that this coupling
cannot lead to on-shell decays of $H$ particle. Moreover in the non-
relativistic limit, this vertex contains an oscillatory factor $e^{\pm imt}$
(see (15) ahead). As a result, we ignore this 3-point interaction for this
paper.

####  A.3.1 Coupling to matter

We also need to worry about the interactions of the massive spin-$2$ field
with the SM. Since we have two metrics $\mathfrak{f}$ and $\mathfrak{g}$, and
both contain the massive and massless spin-$2$ degrees of freedom (E.q.
(A.3)), other SM fields will couple to both $h$ and $H$ if the metric that
appears in the visible action $S_{\rm vis}$ is only $\mathfrak{g}$. This leads
to decay channels of type $H\rightarrow{\rm SM}+{\rm SM}$, which puts strong
cosmological constraints on $\alpha$ [64]. However, there exists at-least one
broader class of metrics in the literature

$\displaystyle
g^{\mathrm{eff}}_{\mu\nu}=a^{2}\,\mathfrak{g}_{\mu\nu}+2\,ab\,\mathfrak{g}_{\mu\lambda}X^{\lambda}_{\;\;\nu}+b^{2}\,\mathfrak{f}_{\mu\nu},$
(66)

that couple to the usual matter sector, and are ghost free [29, 30].
Explicitly, we have $S_{\rm
vis}=\int\mathrm{d}^{4}x\sqrt{-g_{\mathrm{eff}}}\,\mathcal{L}(g_{\mathrm{eff}},\psi_{\rm
SM})$. Expanding $g^{\rm eff}$ using (A.3) yields

$\displaystyle g^{\mathrm{eff}}_{\mu\nu}$
$\displaystyle=(a+b)^{2}\left[\eta_{\mu\nu}+\frac{1}{m_{\mathrm{pl}}}h_{\mu\nu}+\frac{\sqrt{2}(b-a\alpha^{2})}{\alpha(a+b)m_{\mathrm{pl}}}H_{\mu\nu}\right]$
$\displaystyle\quad+\mathcal{O}(H^{2}).$ (67)

Then with $a=1/(1+\alpha^{2})$ and $b={\alpha^{2}}/{(1+\alpha^{2})}$, we can
kill the leading order $H_{\mu\nu}$ term in the effective metric (and hence
$S_{\rm vis}$). This removes the $H\rightarrow{\rm SM}+{\rm SM}$ channel,
alleviating the cosmological bounds on $\alpha$ derived in [64]. More
importantly for our purposes, this allows us to have $\alpha$ values that are
neither too large or small and therefore safely neglect the trilinear terms in
(65) in the non-relativistic limit. Furthermore this choice removes the
leading interaction terms between $H_{\mu\nu}$ and other dark sector fields.
As a result, one may even consider scenarios in which all different spin dark
fields are present in the same theory.

####  A.3.2 Cosmological Solutions

Unlike the spin-$0$ and spin-$1$ cases, including FLRW background comes with
additional complications. For example, in an FLRW background there might exist
gradient instabilities in the vector perturbations during the radiation
dominated era [31]. As discussed in that paper, it might be possible to cure
these by extending the bi-gravity model.

Another thing to consider is the Higuchi bound on the mass of any spin-s field
in De-Sitter, $m^{2}\geq s(s-1)H^{2}$ (see [65]). For very high scale
inflation, this requires the mass of a $s\geq 2$ spin field to be even higher,
resulting in ultrasmall solitons. This however is not an issue for
sufficiently small scale inflation, with the least conservative bound coming
from Big Bang Nucleosynthesis, $H\gtrsim
T_{\mathrm{BBN}}/m_{\mathrm{pl}}^{2}\sim 10^{-16}$ eV.

Furthermore, the cut-off of the bi-gravity theory
$\Lambda_{3}\sim(m^{2}m_{\mathrm{pl}})^{1/3}$ (obtained for example, by
considering unitarity bounds from $2\rightarrow 2$ scattering of various
helicity modes of the field) can be quite small when $m$ is sufficiently
small. While this may not be an issue in the non-relativistic limit where
$k\ll m$, a consistent, calculable set-up at higher energies which become
accessible in the early universe, might require further work.

## References

  * Aghanim  _et al._ [2020] N. Aghanim _et al._ (Planck), Planck 2018 results. VI. Cosmological parameters, Astron. Astrophys. 641, A6 (2020), arXiv:1807.06209 [astro-ph.CO] .
  * Iršič  _et al._ [2017] 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, 031302 (2017), arXiv:1703.04683 [astro-ph.CO] .
  * Jacobs  _et al._ [2015] D. M. Jacobs, G. D. Starkman, and B. W. Lynn, Macro Dark Matter, Mon. Not. Roy. Astron. Soc. 450, 3418 (2015), arXiv:1410.2236 [astro-ph.CO] .
  * Bertone and Hooper [2018] G. Bertone and D. Hooper, History of dark matter, Rev. Mod. Phys. 90, 045002 (2018), arXiv:1605.04909 [astro-ph.CO] .
  * Niemeyer [2019] J. C. Niemeyer, Small-scale structure of fuzzy and axion-like dark matter (2019), arXiv:1912.07064 [astro-ph.CO] .
  * Grin  _et al._ [2019] D. Grin, M. A. Amin, V. Gluscevic, R. Hlǒzek, D. J. E. Marsh, V. Poulin, C. Prescod-Weinstein, and T. L. Smith, Gravitational probes of ultra-light axions (2019), arXiv:1904.09003 [astro-ph.CO] .
  * Ferreira [2020] E. G. Ferreira, Ultra-Light Dark Matter (2020), arXiv:2005.03254 [astro-ph.CO] .
  * Amin  _et al._ [2014] M. A. Amin, M. P. Hertzberg, D. I. Kaiser, and J. Karouby, Nonperturbative Dynamics Of Reheating After Inflation: A Review, Int. J. Mod. Phys. D24, 1530003 (2014), arXiv:1410.3808 [hep-ph] .
  * Lozanov [2019] K. D. Lozanov, Lectures on Reheating after Inflation (2019), arXiv:1907.04402 [astro-ph.CO] .
  * Allahverdi  _et al._ [2020] R. Allahverdi _et al._ , The First Three Seconds: a Review of Possible Expansion Histories of the Early Universe (2020), arXiv:2006.16182 [astro-ph.CO] .
  * Namjoo  _et al._ [2018] M. H. Namjoo, A. H. Guth, and D. I. Kaiser, Relativistic Corrections to Nonrelativistic Effective Field Theories, Phys. Rev. D98, 016011 (2018), arXiv:1712.00445 [hep-ph] .
  * Salehian  _et al._ [2020] B. Salehian, M. H. Namjoo, and D. I. Kaiser, Effective theories for a nonrelativistic field in an expanding universe: Induced self-interaction, pressure, sound speed, and viscosity, JHEP 07, 059, arXiv:2005.05388 [astro-ph.CO] .
  * Salehian  _et al._ [2021] 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 (2021), arXiv:2104.10128 [astro-ph.CO] .
  * Adshead and Lozanov [2021] P. Adshead and K. D. Lozanov, Self-gravitating Vector Dark Matter, Phys. Rev. D 103, 103501 (2021), arXiv:2101.07265 [gr-qc] .
  * Aoki  _et al._ [2018] K. Aoki, K.-i. Maeda, Y. Misonoh, and H. Okawa, Massive Graviton Geons, Phys. Rev. D 97, 044005 (2018), arXiv:1710.05606 [gr-qc] .
  * RUFFINI and BONAZZOLA [1969] R. RUFFINI and S. BONAZZOLA, Systems of self-gravitating particles in general relativity and the concept of an equation of state, Phys. Rev. 187, 1767 (1969).
  * Chavanis [2011] P.-H. Chavanis, Mass-radius relation of newtonian self-gravitating bose-einstein condensates with short-range interactions. i. analytical results, Physical Review D 84, 10.1103/physrevd.84.043531 (2011).
  * Caldwell  _et al._ [2017] A. Caldwell, G. Dvali, B. Majorovits, A. Millar, G. Raffelt, J. Redondo, O. Reimann, F. Simon, and F. Steffen (MADMAX Working Group), Dielectric Haloscopes: A New Way to Detect Axion Dark Matter, Phys. Rev. Lett. 118, 091801 (2017), arXiv:1611.05865 [physics.ins-det] .
  * Baryakhtar  _et al._ [2018] M. Baryakhtar, J. Huang, and R. Lasenby, Axion and hidden photon dark matter detection with multilayer optical haloscopes, Phys. Rev. D 98, 035006 (2018), arXiv:1803.11455 [hep-ph] .
  * Chiles  _et al._ [2021] J. Chiles _et al._ , First Constraints on Dark Photon Dark Matter with Superconducting Nanowire Detectors in an Optical Haloscope (2021), arXiv:2110.01582 [hep-ex] .
  * Zhang  _et al._ [2021] H.-Y. Zhang, M. Jain, and M. A. Amin, Polarized Vector Oscillons (2021), arXiv:2111.08700 [astro-ph.CO] .
  * de Rham  _et al._ [2011] C. de Rham, G. Gabadadze, and A. J. Tolley, Resummation of Massive Gravity, Phys. Rev. Lett. 106, 231101 (2011), arXiv:1011.1232 [hep-th] .
  * Hassan and Rosen [2012a] S. F. Hassan and R. A. Rosen, Resolving the Ghost Problem in non-Linear Massive Gravity, Phys. Rev. Lett. 108, 041101 (2012a), arXiv:1106.3344 [hep-th] .
  * Hassan  _et al._ [2012] S. F. Hassan, R. A. Rosen, and A. Schmidt-May, Ghost-free Massive Gravity with a General Reference Metric, JHEP 02, 026, arXiv:1109.3230 [hep-th] .
  * Hassan and Rosen [2012b] S. F. Hassan and R. A. Rosen, Bimetric Gravity from Ghost-free Massive Gravity, JHEP 02, 126, arXiv:1109.3515 [hep-th] .
  * Hinterbichler [2012] K. Hinterbichler, Theoretical Aspects of Massive Gravity, Rev. Mod. Phys. 84, 671 (2012), arXiv:1105.3735 [hep-th] .
  * de Rham [2014] C. de Rham, Massive Gravity, Living Rev. Rel. 17, 7 (2014), arXiv:1401.4173 [hep-th] .
  * Schmidt-May and von Strauss [2016] A. Schmidt-May and M. von Strauss, Recent developments in bimetric theory, J. Phys. A 49, 183001 (2016), arXiv:1512.00021 [hep-th] .
  * Noller and Melville [2015] J. Noller and S. Melville, The coupling to matter in Massive, Bi- and Multi-Gravity, JCAP 01, 003, arXiv:1408.5131 [hep-th] .
  * Bonifacio  _et al._ [2018] J. Bonifacio, K. Hinterbichler, A. Joyce, and R. A. Rosen, Massive and Massless Spin-2 Scattering and Asymptotic Superluminality, JHEP 06, 075, arXiv:1712.10020 [hep-th] .
  * Comelli  _et al._ [2015] D. Comelli, M. Crisostomi, K. Koyama, L. Pilo, and G. Tasinato, Cosmology of bigravity with doubly coupled matter, JCAP 04, 026, arXiv:1501.00864 [hep-th] .
  * Eby  _et al._ [2020] J. Eby, M. Leembruggen, L. Street, P. Suranyi, and L. C. R. Wijewardhana, Galactic Condensates composed of Multiple Axion Species, JCAP 10, 020, arXiv:2002.03022 [hep-ph] .
  * Brito  _et al._ [2016] 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, 291 (2016), arXiv:1508.05395 [gr-qc] .
  * Boskovic  _et al._ [2018] M. Boskovic, F. Duque, M. C. Ferreira, F. S. Miguel, and V. Cardoso, Motion in time-periodic backgrounds with applications to ultralight dark matter haloes at galactic centers, Phys. Rev. D 98, 024037 (2018), arXiv:1806.07331 [gr-qc] .
  * Mocz  _et al._ [2017] 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, 4559 (2017), arXiv:1705.05845 [astro-ph.CO] .
  * Amin and Mocz [2019] M. A. Amin and P. Mocz, Formation, gravitational clustering, and interactions of nonrelativistic solitons in an expanding universe, Phys. Rev. D 100, 063507 (2019), arXiv:1902.07261 [astro-ph.CO] .
  * Lancaster  _et al._ [2020] L. Lancaster, C. Giovanetti, P. Mocz, Y. Kahn, M. Lisanti, and D. N. Spergel, Dynamical Friction in a Fuzzy Dark Matter Universe, JCAP 01, 001, arXiv:1909.06381 [astro-ph.CO] .
  * Helfer  _et al._ [2019] T. Helfer, E. A. Lim, M. A. G. Garcia, and M. A. Amin, Gravitational Wave Emission from Collisions of Compact Scalar Solitons, Phys. Rev. D 99, 044046 (2019), arXiv:1802.06733 [gr-qc] .
  * Khmelnitsky and Rubakov [2014] A. Khmelnitsky and V. Rubakov, Pulsar timing signal from ultralight scalar dark matter, JCAP 02, 019, arXiv:1309.5888 [astro-ph.CO] .
  * Brito  _et al._ [2015] R. Brito, V. Cardoso, and P. Pani, Superradiance: New Frontiers in Black Hole Physics, Lect. Notes Phys. 906, pp.1 (2015), arXiv:1501.06570 [gr-qc] .
  * Hui  _et al._ [2021] L. Hui, A. Joyce, M. J. Landry, and X. Li, Vortices and waves in light dark matter, JCAP 01, 011, arXiv:2004.01188 [astro-ph.CO] .
  * Cembranos  _et al._ [2017] J. A. R. Cembranos, A. L. Maroto, and S. J. Núñez Jareño, Perturbations of ultralight vector field dark matter, JHEP 02, 064, arXiv:1611.03793 [astro-ph.CO] .
  * Baryakhtar  _et al._ [2017] M. Baryakhtar, R. Lasenby, and M. Teo, Black hole superradiance signatures of ultralight vectors, Phys. Rev. D 96, 035019 (2017).
  * Baumann  _et al._ [2019] D. Baumann, H. S. Chia, J. Stout, and L. ter Haar, The Spectra of Gravitational Atoms, JCAP 12, 006, arXiv:1908.10370 [gr-qc] .
  * Armaleo  _et al._ [2020] J. M. Armaleo, D. López Nacir, and F. R. Urban, Pulsar timing array constraints on spin-2 ULDM, JCAP 09, 031, arXiv:2005.03731 [astro-ph.CO] .
  * Hertzberg and Schiappacasse [2018] M. P. Hertzberg and E. D. Schiappacasse, Dark Matter Axion Clump Resonance of Photons, JCAP 11, 004, arXiv:1805.00430 [hep-ph] .
  * Amin  _et al._ [2021] M. A. Amin, A. J. Long, Z.-G. Mou, and P. Saffin, Dipole radiation and beyond from axion stars in electromagnetic fields, JHEP 06, 182, arXiv:2103.12082 [hep-ph] .
  * Hertzberg  _et al._ [2020] M. P. Hertzberg, Y. Li, and E. D. Schiappacasse, Merger of Dark Matter Axion Clumps and Resonant Photon Emission, JCAP 07, 067, arXiv:2005.02405 [hep-ph] .
  * Levkov  _et al._ [2020] D. Levkov, A. Panin, and I. Tkachev, Radio-emission of axion stars, Phys. Rev. D 102, 023501 (2020), arXiv:2004.05179 [astro-ph.CO] .
  * Amin and Mou [2021] M. A. Amin and Z.-G. Mou, Electromagnetic Bursts from Mergers of Oscillons in Axion-like Fields, JCAP 02, 024, arXiv:2009.11337 [astro-ph.CO] .
  * Vitvitska  _et al._ [2002] M. Vitvitska, A. A. Klypin, A. V. Kravtsov, R. H. Wechsler, J. R. Primack, and J. S. Bullock, The origin of angular momentum in dark matter halos, Astrophys. J. 581, 799 (2002), arXiv:astro-ph/0105349 .
  * Kolb and Tkachev [1994] E. W. Kolb and I. I. Tkachev, Nonlinear axion dynamics and formation of cosmological pseudosolitons, Phys. Rev. D49, 5040 (1994), arXiv:astro-ph/9311037 [astro-ph] .
  * Amin  _et al._ [2012] M. A. Amin, R. Easther, H. Finkel, R. Flauger, and M. P. Hertzberg, Oscillons After Inflation, Phys. Rev. Lett. 108, 241302 (2012), arXiv:1106.3335 [astro-ph.CO] .
  * Agrawal  _et al._ [2020] P. Agrawal, N. Kitajima, M. Reece, T. Sekiguchi, and F. Takahashi, Relic Abundance of Dark Photon Dark Matter, Phys. Lett. B 801, 135136 (2020), arXiv:1810.07188 [hep-ph] .
  * Co  _et al._ [2019] R. T. Co, A. Pierce, Z. Zhang, and Y. Zhao, Dark Photon Dark Matter Produced by Axion Oscillations, Phys. Rev. D 99, 075002 (2019), arXiv:1810.07196 [hep-ph] .
  * Dror  _et al._ [2019] J. A. Dror, K. Harigaya, and V. Narayan, Parametric Resonance Production of Ultralight Vector Dark Matter, Phys. Rev. D 99, 035036 (2019), arXiv:1810.07195 [hep-ph] .
  * Kolb and Long [2021] E. W. Kolb and A. J. Long, Completely dark photons from gravitational particle production during the inflationary era, JHEP 03, 283, arXiv:2009.03828 [astro-ph.CO] .
  * Alexander  _et al._ [2021] S. Alexander, L. Jenks, and E. McDonough, Higher spin dark matter, Phys. Lett. B 819, 136436 (2021), arXiv:2010.15125 [hep-ph] .
  * Blinov  _et al._ [2021] N. Blinov, M. J. Dolan, P. Draper, and J. Shelton, Dark Matter Microhalos From Simplified Models, Phys. Rev. D 103, 103514 (2021), arXiv:2102.05070 [astro-ph.CO] .
  * Bouatta  _et al._ [2004] N. Bouatta, G. Compere, and A. Sagnotti, An Introduction to free higher-spin fields, in  _1st Solvay Workshop on Higher Spin Gauge Theories_ (2004) arXiv:hep-th/0409068 .
  * Sorokin [2005] D. Sorokin, Introduction to the classical theory of higher spins, AIP Conf. Proc. 767, 172 (2005), arXiv:hep-th/0405069 .
  * Copeland  _et al._ [1995] E. J. Copeland, M. Gleiser, and H. R. Muller, Oscillons: Resonant configurations during bubble collapse, Phys. Rev. D52, 1920 (1995), arXiv:hep-ph/9503217 [hep-ph] .
  * Zhang  _et al._ [2020] H.-Y. Zhang, M. A. Amin, E. J. Copeland, P. M. Saffin, and K. D. Lozanov, Classical Decay Rates of Oscillons, JCAP 07, 055, arXiv:2004.01202 [hep-th] .
  * Babichev  _et al._ [2016] E. Babichev, L. Marzola, M. Raidal, A. Schmidt-May, F. Urban, H. Veermäe, and M. von Strauss, Heavy spin-2 Dark Matter, JCAP 09, 016, arXiv:1607.03497 [hep-th] .
  * Lee  _et al._ [2016] H. Lee, D. Baumann, and G. L. Pimentel, Non-Gaussianity as a Particle Detector, JHEP 12, 040, arXiv:1607.03735 [hep-th] .