Formation, Gravitational Clustering and Interactions of Non-relativistic Solitons in an Expanding Universe

TL;DR

This paper explores how non-relativistic scalar field solitons form, cluster, and interact under gravity in an expanding universe, revealing complex dynamics like mergers and binary formation through simulations and analytic estimates.

Contribution

It provides new insights into the formation, stability, and gravitational clustering of solitons in an expanding universe, combining numerical and analytical methods.

Findings

Rapid soliton formation driven by self-interactions

Gravitational clustering leads to mergers and binary formation

Analytic estimates of instability scales and soliton profiles

Abstract

We investigate the formation, gravitational clustering and interactions of solitons in a self-interacting, non-relativistic scalar field in an expanding universe. Rapid formation of large number of solitons is driven by attractive self-interactions of the field, whereas the slower clustering of solitons is driven by gravitational forces. Driven closer together by gravity, we see a rich plethora of dynamics in the soliton "gas" including mergers, scatterings and formation of soliton binaries. The numerical simulations are complemented by analytic calculations and estimates of (i) the relevant instability length and time scales, (ii) individual soliton profiles and their stability, (iii) number density of produced solitons, and (iv) the two point correlation function of soliton positions as evidence for gravitational clustering.

Equations

$$\left[i\left({\partial }_t+\frac{3}{2}H\right)+\frac{1}{2a^2}{\nabla }^2-U_{\mathrm{nl}}^{\prime }(|\psi {|}^2)-\Phi \right]\psi =0,$$

$$\frac{{\nabla }^2}{a^2}\Phi =\frac{{\beta }^2}{2}\left[|\psi {|}^2+\frac{1}{2a^2}|\nabla \psi {|}^2+U_{\mathrm{nl}}(|\psi {|}^2)\right]-\frac{3}{2}{H}^2,$$

$$H^2=\frac{{\beta }^2}{3}\overline{\left[|\psi {|}^2+\frac{1}{2a^2}|\nabla \psi {|}^2+U_{\mathrm{nl}}(|\psi {|}^2)\right]},$$

$$\beta \equiv \frac{M}{m_{\mathrm{pl}}}\ll 1.$$

$$U_{\mathrm{nl}}(|\psi {|}^2)=-\frac{|\psi {|}^2}{2}\frac{|\psi {|}^2}{1+|\psi {|}^2}.$$

$$\psi (t_{\mathrm{in}},\mathbf{x})=\bar{\psi }(t_{\mathrm{in}})+\frac{1}{\sqrt{L^3}}\left(\frac{m}{M}\right)\sum _{\mathbf{k}}\delta {\psi }_{\mathbf{k}}{e}^{i\mathbf{k}\cdot \mathbf{x}},$$

$$\mathrm{with}\bar{\psi }(t_{\mathrm{in}})=1\mathrm{and}⟨|\delta {\psi }_{\mathbf{k}}{|}^2⟩\sim \frac{1}{2}.$$

$$\psi (t,\mathbf{x})=\bar{\psi }(t)\left[1+\epsilon \frac{\delta {\psi }_{\mathbf{k}}(t)}{\bar{\psi }(t)}{e}^{i\mathbf{k}\cdot \mathbf{x}}\right],$$

$$\left({\partial }_t^2+\frac{k^2}{4}\left[k^2+4|\bar{\psi }{|}^2{U}_{\mathrm{nl}}^{\prime \prime }(|\bar{\psi }{|}^2)\right]\right)\frac{\delta {\psi }_{\mathbf{k}}}{\bar{\psi }}=0.$$

$$k^2<-4|\bar{\psi }{|}^2{U}_{\mathrm{nl}}^{\prime \prime }(|\bar{\psi }{|}^2),$$

$$\displaystyle\mu_{k}=\Bigg{|}i\frac{k}{2}\sqrt{k^{2}+4|\bar{\psi}|^{2}U_{\rm nl}^{\prime\prime}(|\bar{\psi}|^{2})}\Bigg{|}\,.$$

$$k_{\star }$$

$$\text{with}{\mu }_{k_{\star }}$$

$$\frac{{\mu }_k}{H}\sim \frac{1}{\beta }\frac{1}{a^{3/2}}\gg 1\mathrm{for}\text{rapid growth}.$$

$$k<k_J\approx \sqrt{\sqrt{2}\beta |\bar{\psi }|}\mathrm{with}{\mu }_k=\sqrt{\frac{1}{2}{\beta }^2|\bar{\psi }{|}^2-\frac{k^4}{4}}.$$

$$\frac{k}{a}\lesssim \{\begin{array}{ll}& a^{-3/2}\text{self-interactions},\\ & \sqrt{\beta }{a}^{-3/4}\text{gravity}.\end{array}$$

$$\frac{{\mu }_k}{H}\sim \{\begin{array}{ll}& {\beta }^{-1}{a}^{-3/2}\text{self-interactions},\\ & 1\text{gravity}.\end{array}$$

$$\frac{m}{M}{k}^{3/2}⟨|\delta {\psi }_{\mathbf{k}}{|}^2{⟩}^{1/2}\sim \bar{\psi },$$

$$a^3{n}_{\mathrm{sol}}\sim (k_{\mathrm{nl}}/2\pi {)}^3,$$

$$\psi (t,r)=e^{-i\nu t}\Psi (r).$$

$$\left[\nu +\frac{1}{2r^2}{\partial }_r(r^2{\partial }_r)-U_{\mathrm{nl}}^{\prime }({\Psi }^2)-\Phi \right]\Psi =0,$$

$$\frac{1}{r^2}{\partial }_r(r^2{\partial }_r)\Phi =\frac{{\beta }^2}{2}\left[{\Psi }^2+\frac{1}{2}({\partial }_r\Psi {)}^2+U_{\mathrm{nl}}({\Psi }^2)\right].$$

$$\mathcal{E}$$

$$=\mathcal{O}[10^2]\times {\left(\frac{M}{m}\right)}^{2}m,$$

$$\frac{d\mathcal{N}}{d(-\nu )}>0\text{where}\mathcal{N}\equiv \int d^{3}r{\Psi }^2(r),$$

$${\xi }_{\mathrm{LS}}(r)=\frac{DD}{RR}-2\frac{\mathcal{N}-1}{\mathcal{N}}\frac{DR}{RR}+1,$$

$${\xi }_{\mathrm{LS}}(r)\propto \frac{1}{r^2},$$

$$\frac{d\ln (a^3{n}_{\mathrm{sol}})}{d\ln a}\simeq 0.1,$$

$$S=\int \frac{d^{4}x}{\hslash c^2}\sqrt{-g}\left[\frac{R}{16\pi G}-\frac{1}{2}{g}^{\mu \nu }{\partial }_{\nu }\varphi {\partial }_{\mu }\varphi -V(\varphi )\right],$$

$$V(\varphi )=\frac{m^2{c}^2}{2{\hslash }^2}{\varphi }^2+V_{\mathrm{nl}}(\varphi ),$$

$$d{s}^2=\left(1+2\frac{\Phi }{c^2}\right)(cdt{)}^2-\left(1-2\frac{\Phi }{c^2}\right)d{\mathbf{x}}^2.$$

$$\frac{1}{\sqrt{-g}}{\partial }_{\mu }\left(\sqrt{-g}{g}^{\mu \nu }{\partial }_{\nu }\varphi \right)+{\partial }_{\varphi }V(\varphi )=0,$$

$$\frac{{\partial }^2\varphi }{\partial (ct{)}^2}-\left(1+4\frac{\Phi }{c^2}\right){\nabla }^2\varphi -\frac{4}{c^2}\frac{\partial \Phi }{\partial (ct)}\frac{\partial \varphi }{\partial (ct)}+\left(1+2\frac{\Phi }{c^2}\right){\partial }_{\varphi }V=0.$$

$${\nabla }^2\Phi =\frac{4\pi G}{c^2}{T}^0{}_0,$$

Full text

I Introduction

Solitons are self-localized, persistent configurations in nonlinear field theories which have been studied intensely in a broad range of contexts including cosmology, high energy physics, nonlinear optics and cold-atom physics, condensed matter physics, fluid mechanics and mathematics Lee and Pang (1992); Vilenkin and Shellard (2000); Manton and Sutcliffe (2004); Weinberg (2012); Liebling and Palenzuela (2012); Kivshar and Malomed (1989).

In cosmology, for example, solitons can emerge naturally at the end of inflation and dominate the energy density (e.g. Amin et al. (2012)), or related configurations can form in the axion field that might constitute the entirety or part of the dark matter (e.g. Kolb and Tkachev (1994)). Depending on the context, they can act as new sources of gravitational waves Zhou et al. (2013); Antusch et al. (2018); Helfer et al. (2018); Palenzuela et al. (2017); Lozanov and Amin (2019), potentially lead to the formation of primordial black holes Khlopov et al. (1985); Cotner and Kusenko (2017); Cotner et al. (2018), be involved in baryogenesis Enqvist and McDonald (1998); Lozanov and Amin (2014), and change the approach to radiation domination in the early universe Adshead et al. (2015); Lozanov and Amin (2017, 2018); they can also provide novel insights into the small scale problems in the cold dark matter paradigm Hu et al. (2000); Alcubierre et al. (2002); Marsh and Pop (2015); Schive et al. (2014); Hui et al. (2017).

To explore many of these implications, it is important to consider their formation, and their interactions resulting from gravity and self-couplings of the field. In this paper we explore the gravitational clustering and gravitational as well as non-gravitational interactions of non-relativistic solitons, starting with the formation of solitons from cosmological initial conditions.111Here, by cosmological initial conditions we mean an almost homogeneous field with small perturbations. See Fig. 1 for a visual overview of soliton formation and clustering in an expanding universe.

We focus on nontopological solitons in a non-relativistic scalar field theory. We include strong self-interactions in the theory, while gravity is included under the assumption that it is weak. In our simulations, the rate of expansion of space is determined by the average energy density of the field.

There is a large amount of diverse literature on nontopological solitons in real and complex scalar field theories in a cosmological context; this paragraph is a sample rather than a comprehensive review of the literature. For work on individual solitons, see, for example, Kaup (1968); Coleman (1985); Bogolyubsky and Makhankov (1976); Copeland et al. (1995); Amin and Shirokoff (2010); Kasuya et al. (2003); Amin (2013); Chavanis (2011); Chavanis and Delfini (2011); Chavanis (2012, 2018). For the early universe, soliton formation in relativistic fields in an expanding universe but ignoring gravitational interactions has been considered, for example, in Kusenko and Shaposhnikov (1998); Farhi et al. (2008); Amin et al. (2012); Gleiser et al. (2011). In the late universe context, gravitational interactions are included in the non-relativistic limit, but self-interactions are ignored or typically assumed to be very weak (e.g. Schive et al. (2014); Mocz et al. (2017); Schwabe et al. (2016)). In this non-relativistic, non-interacting limit, halos and solitons within them have been shown to form. Binary soliton collisions or interactions and their implications have also been explored under controlled initial conditions (e.g. Palenzuela et al. (2007); Amin et al. (2014); Schwabe et al. (2016); Helfer et al. (2018)). The fate of a “prepared” collection of relativistic solitons (oscillons) with random velocities was considered in two dimensions and without gravity in Hindmarsh and Salmi (2008). The mergers of a small group of pre-existing non-relativistic solitons, with gravity included but without self-interactions, was explored in Mocz et al. (2017); Schwabe et al. (2016).

In our work we present the following results for the first time: We simulate and analyze the case of soliton formation with strong self-interactions, starting with cosmological initial conditions. Thereafter, a “gas” of solitons emerges in a self-consistently expanding universe, followed by gravitational clustering of solitons and eventual, dynamically rich, close encounters. We provide a quantitative understanding of the formation, gravitational clustering, individual properties and interactions of solitons based on simulations and analytic calculations. We note that quite generally, we can use the results in the present work to understand the formation and gravitational clustering dynamics of the non-relativistic limits of oscillons, Q-balls, and boson stars with strong self-interactions.222The emergent, approximate $U(1)$ symmetry in the non-relativistic limit makes oscillons and Q-balls almost identical Kasuya et al. (2003), or at the very least obtainable from one another.

The present work is somewhat related to (but does not rely on) a recent exploration of gravitational perturbations from oscillons and transients Lozanov and Amin (2019). In Lozanov and Amin (2019), soliton formation in a relativistic Klein-Gordon equation in an expanding background was investigated; however, gravitational perturbations were calculated passively (i.e. gravitational clustering was not present). Here, we focus on non-relativistic fields, but clustering due to gravity is included. While the models and context are not identical, a qualitative comparison between relativistic and non-relativistic models and results is discussed in the Appendix.

The rest of the paper is organized as follows. In Section II we discuss the model for a non-relativistic, self-interacting field in an expanding universe with weak field gravity. In Section III, we briefly discuss the lattice simulation and our numerical algorithm. The initial conditions for the simulations are provided in Section IV. We analyze linear instabilities from self-interactions and gravitational interactions in Section V. The numerically calculated power spectrum for the field perturbations is provided in Section VI. In Section VII we discuss the formation of solitons, followed by a discussion of their individual profiles and stability in Section VIII. The gravitational clustering of solitons is discussed in Section IX, and resulting strong soliton interactions are explored in Section X. Finally, we present our conclusions and future directions in Section XI. In the Appendix we discuss connections to a related relativistic system (at the level of the equations, instabilities, solitons and initial conditions).

II The Model

We use the following equations of motion (and constraint equations) to explore the dynamics of a non-relativistic, self-interacting, self-gravitating scalar field in an approximately homogeneous and isotropic universe:

[TABLE]

where $\overline{[\ldots]}$ indicates a spatial average, $a(t)$ is the scalefactor, $H(t)=\dot{a}(t)/{a}(t)$ is the Hubble rate, $\psi(t,{\bf x})$ is the complex field amplitude, $\Phi(t,{\bf x})$ is the Newtonian potential, and $U_{\rm nl}(|\psi|^{2})$ encodes the self-interactions of the field.333We have checked that qualitatively similar results are obtained even if we set $U_{\rm nl}\rightarrow 0$ in the Poisson and Friedmann equations, but keep $U_{\rm n}^{\prime}(|\psi|^{2})\equiv\partial_{|\psi|^{2}}U_{\rm n}(|\psi|^{2})$ in the nonlinear Schrödinger equation.

All variables and parameters appearing in the above equation are dimensionless. We have expressed time $t$ in units of $\tau_{m}=\hbar/mc^{2}$, lengths in units of $\lambda_{m}=\hbar/mc$, the Newtonian gravitational potential $\Phi$ in units of $c^{2}$ and $|\psi|^{2}$ in units of $m^{2}M^{2}c^{3}/\hbar^{3}$. Note that $m^{2}M^{2}c^{3}/\hbar^{3}$ has dimensions of mass density. We assume that the parameter

[TABLE]

There are three relevant scales in the equations (not easily discernible in the non-dimensional version): $m=$ mass of particles of our field (without self-interactions), $M$ determines the strength of the self-interactions, and ${m_{\rm pl}}$ is the reduced Planck mass which determines the strength of gravity. We work in a parameter regime with: $m\ll M\ll{m_{\rm pl}}$. The fiducial value used in the present paper is $M=0.03\,{m_{\rm pl}}$ (though we have also varied $M$ by a factor of a few). This particular parameter regime can be natural when identifying $\psi$ as the non-relativistic approximation to the inflaton field Lozanov and Amin (2018) (with $m\simeq 2\times 10^{-4}M$). The hierarchy $m\ll M\ll{m_{\rm pl}}$ is also natural for an axion-like field where $M$ plays the role of the decay constant $f$; in this case $m$ can be much smaller (e.g. Arvanitaki et al. (2010); Hlozek et al. (2015)). We note that $m$ is essentially setting units of quantities in our equations, and the behavior we explore will be qualitatively valid for any energetically dominant, cosmological scalar field regardless of the particular value of $m$ (modulo initial conditions).

For the purpose of this paper, we chose $U_{\rm nl}(|\psi|^{2})$ with a saturated nonlinearity :

[TABLE]

The saturated nonlinearity refers to the fact that for $|\psi|\gg 1$, $U_{\rm nl}^{\prime}(|\psi|^{2})\rightarrow\rm const.$ which means that the nonlinearity appearing in the equation of motion for $\psi$ is bounded. This form is not strictly necessary, and different powers of $|\psi|^{2}$ in the denominator of $U_{\rm nl}$ [for example, $(1+|\psi|^{2})^{\alpha}$ or $(1+|\psi|^{2\alpha})$ with $\alpha>0$], are also worth exploring; however, we do not consider these here.

Note that for $|\psi|^{2}\ll 1$, the above choice yields $U^{\prime}(|\psi|^{2})=-|\psi|^{2}$, which makes the first equation in (LABEL:eq:MasterEq) analogous to the usual nonlinear Schrödinger equation with attractive interactions (ignoring gravity). Equation (LABEL:eq:MasterEq) also matches the equations of motion for axions, or symmetric inflationary potentials in this non-relativistic, small amplitude limit.

While not necessary for our present purposes, we explore the connection of our non-relativistic equations to those obtained from a relativistic theory in the Appendix. We also refer the reader to (for example) Namjoo et al. (2018); Eby et al. (2018) for more detailed discussions of the non-relativistic limit of relativistic scalar field systems (typically in the weak interaction limit). At the leading order, the non-relativistic limit of real or complex scalar fields should yield equations similar to ours.

III Lattice Simulations

We solve our Schrödinger-Poisson system in a self-consistently expanding background (see eq. (LABEL:eq:MasterEq)) on an $N=400^{3}$ lattice.444Smaller lattices were also used to check for convergence, and other numerical checks. For example, we halved the resolution, and the locations of solitons did not change. In our highest resolution simulations, solitons contain $\gtrsim 10$ pixels per linear dimension ($\mathcal{O}[10^{3}]$ pixels per soliton volume). The field evolution uses a second-order in time (exponential convergence in space) ‘kick’-‘drift’-‘kick’ spectral method of Mocz et al. (2017). For our numerical method, the total run time scales as $\mathcal{O}[N^{5}]$, which limits $N$ from being too large. The initial box size is $L\sim H^{-1}$, and we run our simulations from $a_{\rm i}=1$ to $a_{\rm f}=20$ (with corresponding $t_{\rm f}-t_{\rm i}\simeq{\rm few}\times 10^{3}\,m^{-1}$). The box size ($L$), resolution ($\Delta x=L/N$), and time duration of the simulations are chosen so that (i) the relevant instability scales (discussed below) are captured in the simulation, (ii) our solitons are resolved ($a_{\rm f}\Delta x\lesssim\mathcal{O}[1]$), and (iii) we have a sufficient number ($\mathcal{O}[10^{2}]$) of solitons in our simulation volume to make statistically significant statements about their properties, interactions and clustering.

We find our solitons in the numerical simulations by locating local maxima in the $|\psi|^{2}$ field (by comparing each pixel to its nearest neighbor in a $3\times 3\times 3$ region), and taking all points with comoving density about some threshold. We look at the radial density profiles about these points and verify that they fall on the central amplitude–radius relation predicted for solitons shown in Fig. 4. In practice, we found that we could distinguish solitons from other local inhomogeneities (which are less dense), with our threshold of $a^{3}|\psi|^{2}>25$. The results are invariant to the particular choice of threshold over a range of values: $10$s–$100$s. A lower threshold would start including extraneous linear fluctuations, and a higher threshold would start excluding solitons.

IV Initial Conditions

We begin with an almost homogeneous field with small spatial perturbations (mimicking zero-point fluctuations) of the form

[TABLE]

where $|\delta\psi_{\bf k}|$ are drawn from a Raleigh distribution, and the phases for $\delta\psi_{\bf k}$ are drawn from a uniform distribution.555Note that $\delta\psi_{\bf k}$ is in units of $m^{1/2}$ (with $\hbar=c=1$). Recall that $\psi$ is measured in units of $mM$ and $L$ in units of $m^{-1}$ which together lead to the appearance of the $m/M\ll 1$ coefficient. To arrive at the above initial conditions, we found it easiest to start from the relativistic case with the relativistic field $\phi=(\sqrt{2}/m)\Re[\psi e^{-imt}]$ (see the Appendix for details). For the initial conditions, we ignore self-interactions, as well as fast time variations and assume $k\lesssim m$. Refinements are possible (such as $|\delta\psi_{\bf k}|^{2}\sim(1/2)\sqrt{k^{2}+m_{\rm eff}^{2}}$), but are not expected to change the results qualitatively. We assume that $a(t_{\rm in})=1$. The gravitational perturbations and $H$ are then obtained self-consistently using eq. (LABEL:eq:MasterEq). We choose $\bar{\psi}(t_{\rm in})=1$ because (as we will see) for $\bar{\psi}(t_{\rm in})\lesssim\beta^{-1}$, instabilities due to self-interactions are ineffective. On the other hand, for $\bar{\psi}(t_{\rm in})\gtrsim\beta^{-1}$ we are forced to introduce a timescale $H^{-1}$ via the Friedmann equation which is comparable to $\tau_{m}$, thus potentially entering a fast timescale, relativistic regime.

To remain consistent with our non-relativistic approximation, we introduce a cutoff in the initial spectrum $\langle|\delta\psi_{\bf k}|^{2}\rangle=0.5e^{-k^{2}}$ which removes relativistic ($k\gg 1$) modes. We have checked that our results are qualitatively insensitive to order unity changes in amplitudes of the initial perturbations as well as the cutoff.

V Linear Instabilities

As seen in Fig. 1, there is a rapid growth in field/density perturbations on a characteristic length scale, which results in the formation of solitons. We calculate and compare this instability with gravitational instability below.

V.1 Self-Interaction Instability

Let us consider small spatial perturbations around a homogeneous solution $\bar{\psi}(t)$:

[TABLE]

where $\varepsilon=(m/M)L^{-3/2}$. Sufficiently long wavelength perturbations of the field are unstable due to self-interactions of the field $U^{\prime}(|\psi|^{2})$. To see this, let us first ignore expansion and gravitational interactions (that is, $a=1$, $H=0$, $\Phi=0$), and substitute eq. (5) into eq. (LABEL:eq:MasterEq). At the background level, we find $\bar{\psi}(t)=\bar{\psi}(0)e^{-i\nu t}$ with $\nu=U^{\prime}(|\bar{\psi}|^{2})<0$. At linear order in the perturbation, we find666To obtain this equation, we found it useful to first derive the first-order equations for the real and imaginary parts of the perturbation $e^{i{\bf k}\cdot{\bf x}}\delta\psi_{\bf k}/\psi$ and then combine them to get the second-order-in-time equations for each part. The real and imaginary parts satisfy the same second-order linear equation; thus, we arrive at eq. (6).

[TABLE]

Note that $U_{\rm nl}^{\prime\prime}(|\bar{\psi}|^{2})<0$ for our case. Thus, we have unstable, exponentially growing perturbations $|\delta\psi_{\bf k}/\bar{\psi}|\propto e^{\mu_{k}t}$ for

[TABLE]

For a given $|\bar{\psi}|$, the mode that grows the fastest has a wavenumber

[TABLE]

The corresponding (approximate) expressions in an expanding universe, are obtained via ${\bf k}\rightarrow{\bf k}/a$. Moreover, in an expanding universe $\bar{\psi}\propto a^{-3/2}$ and $H\sim\beta a^{-3/2}$.

In an expanding universe, this growth rate should be compared to $H$ to ascertain whether the growth of perturbations can compete with expansion related dilution. Using our expressions for $U_{\rm nl}(|\psi|^{2})$ in eq. (LABEL:eq:SI) and $H^{2}$ from the Friedman equation (LABEL:eq:MasterEq), we need

[TABLE]

In the above expression we have assumed that $|\bar{\psi}|\lesssim 1$.

V.2 Gravitational Instability

Spatial perturbations of the field also grow due to gravitational interactions (we ignore self-interactions for the moment). Again, ignoring expansion, usual linear instability analysis of eq. (LABEL:eq:MasterEq) reveals that the unstable perturbations grow exponentially $|\delta\psi_{\bf k}/\bar{\psi}|\sim e^{\mu_{k}t}$ when Hu et al. (2000)

[TABLE]

Heuristically, including expansion means that $|\bar{\psi}|$ and $k$ redshift, and $k$ above should be interpreted as a physical wavenumber $k/a$.777We recognize that including expansion more carefully, the gravitational instability is power-law type rather than exponential, and the fractional overdensity must grow as $\sim a$ for $k<k_{J}$; however, our argument is sufficient to capture the slowness of gravitational instability compared to the self-interaction one Easther et al. (2011).

We end this section by noting that there are two instability scales associated with self-interactions and gravity respectively (see eqs. (LABEL:eq:SI) and (10)). Assuming $|\bar{\psi}|\sim a^{-3/2}\lesssim 1$, the instabilities are active on physical wave-numbers

[TABLE]

The unstable modes have characteristic “growth rates”:

[TABLE]

This simple scaling analysis reveals that for $\beta\ll 1$, the self-interaction instability will dominate at early times.

VI Power Spectrum

The power spectrum of the field perturbations is shown in Fig. 2. The initial spectrum (black) is based on our initial conditions (see eq. (LABEL:eq:IC), including an exponential cutoff which removes $k\gg 1$ modes at this time).

The dashed blue line is the expected power spectrum at $a=1.5$ based on our instability analysis in Section V. This calculated power spectrum is consistent with the numerically evaluated spectrum at the same time which was obtained using the full lattice simulation, both with local gravitational interaction included (solid line) and turned off (dotted line).888We note that there is some power on $k/a\gtrsim 1$ in the power spectrum; part of this is from initial conditions where we were not aggressive in removing all $k/a\gtrsim 1$ modes, and part from rescattering due to nonlinear evolution. However at late times, most of the power is on $k/a\lesssim 1$.

Soon after, the perturbations start becoming nonlinear, and backreaction of the perturbations on the homogeneous evolution of the field becomes significant. The scalefactor when the perturbations become nonlinear can be obtained from the following heuristic criterion which compares the amplitude of field perturbations to the background homogeneous field:

[TABLE]

where the left-hand side is an estimate of the variance of fluctuations on a scale $l\sim k^{-1}$. The above criterion is satisfied by a combination $(a_{\rm nl},k_{\rm nl})$ such that the field perturbations on the comoving scale $k_{\rm nl}$ become nonlinear first. For $\beta=0.03$, we analytically estimate $a_{\rm nl}\simeq 2.1$ and $k_{\rm nl}\simeq 0.7$. Note this scale $k_{\rm nl}/a\simeq 0.35$ in the spectrum in Fig 2 (see the blue curves). A characteristic scale is also visible in the second column ($a=2$) of the snapshots of the field evolution shown in Fig. 1.

VII Soliton Formation

Once the perturbations become nonlinear, the attractive self-interactions lead to the formation of localized, roughly spherical energy density configurations (our solitons) at the peaks of the density perturbations. The comoving number density of such peaks (and hence of solitons) is crudely given by

[TABLE]

at the time of formation (see Amin (2010); Amin et al. (2010)). Using $k_{\rm nl}\simeq 0.7$, we get $a^{3}n_{\rm sol}\sim 10^{-3}$, consistent with our simulations (see Fig. 3 ).

The formation of solitons following the initial linear instability is clearly visible in the snapshots shown in Fig. 1. While we do not show the $a=3$ snapshot, the formation of solitons is complete by this time. The $a=4$ snapshot shows well-formed, and separated solitons, with typical overdensity inside the solitons of $\mathcal{O}[10]$.

In more detail, Fig. 3 shows the comoving number density of solitons as a function of time in our simulations. The initial number density established by the formation of the solitons is independent of self-gravity. However, gravity is strong enough to lead to subsequent mergers or disruptions which leads to a small drop in number density of solitons at late times. In addition, we cannot rule out that gravity is causing some individual solitons to become unstable. The drop in comoving number density is evident in the difference between the dashed (ignoring gravitational interactions) and solid lines ($\lesssim 10\%$ per Hubble time).

We find that a large fraction ($\sim 70\%$) of the energy in a comoving volume of the universe is locked up in solitons. We only count regions with over-densities $\gtrsim 4$ as part of solitons for this estimate. This result is consistent with related earlier simulations using the relativistic nonlinear Klein-Gordon equation in an expanding universe (but ignoring gravitational clustering); see for example Amin et al. (2010, 2012).

VIII Individual Solitons

The first two equations in eq. (LABEL:eq:MasterEq) (ignoring expansion) admit spatially localized, spherically symmetric, solitonic solutions of the form

[TABLE]

We substitute this ansatz into (LABEL:eq:MasterEq) to obtain equations for the profile $\Psi(r)$ and gravitational potential $\Phi(r)$:

[TABLE]

Note that $\nu$ can be absorbed into the definition $\tilde{\Phi}=\Phi-\nu$. We then find smooth, localized, node-free solutions for $\Psi(r)$ for each $\Psi(0)$, by appropriately adjusting $\tilde{\Phi}(0)$.999If needed, we can recover $\nu=\Phi-\tilde{\Phi}$ by insisting that $\Phi(r)\rightarrow 0$ for $r\rightarrow\infty$. In practice, recovering accurate values of $\nu$ is not easy since $\tilde{\Phi}$ falls off as a power law.

We note that by going to the large $r$ limit of the profile equations, $\Psi(r)$ decays in an exponential fashion at large radii (see Schiappacasse and Hertzberg (2018)). This will be relevant when discussing soliton interactions.

In Fig. 4 we plot the $1/e$ width of these soliton profiles as a function of the central amplitude (solid black curve) using the profiles obtained from the above procedure. Note that the width is non-monotonic in the central amplitude. The data points in this plot correspond to solitons extracted from our simulations and are in excellent agreement with the calculated analytic expectation. Note that for early times ($a=2$), not all high density regions are solitons yet; hence, they do not lie on the analytic curve initially.

While we have done the above calculation including gravity, the gravitational potential remains small for most of the solitons: $|\Phi(0)|=\mathcal{O}[10^{-3}]$ for $\beta=\mathcal{O}[10^{-2}]$, and gravity does not significantly affect profiles for central amplitudes $\Psi(0)\lesssim\rm few$. The same is true in our simulations. We also show the gravitational potential at the center of these solitons Fig. 4 (top axis).

The mass (or energy) per soliton is101010Note that ignoring the gradient and potential terms only changes the answer by a factor of few. We briefly restore units with $\hbar=c=1$ to clarify that each soliton contains a large number of $m$ particles.

[TABLE]

for the range of central amplitudes shown in Fig. 4 and seen in simulations. Note that with $m\ll M$, $\mathcal{E}\gg m$. We find that the energy is a non-monotonic function of $\Psi(0)$, with a minimum near $\Psi(0)\simeq 1$.

Stability

From our calculated profiles, we find that for $-\nu\gtrsim 0.05$ (correspondingly, $\Psi(0)\gtrsim 0.9$):

[TABLE]

whereas it is smaller than zero at smaller amplitudes. This Vakhitov-Kolokolov stability criterion Vakhitov and Kolokolov (1973) guarantees stability for solitons with $\Psi(0)\gtrsim 0.9$ against long-wavelength perturbations.

The stability criterion elegantly explains the dearth of solitons with central amplitudes below $\Psi(0)\lesssim 1$ in Fig. 4(see also Eby et al. (2018), where this criterion is argued to hold even when including gravity in the non-relativistic limit).111111A long-wavelength stability analysis for relativistic solitons (oscillons) was carried out in Amin et al. (2010); Amin and Shirokoff (2010) (albeit in a different self-interaction potential, and without gravity), which also showed that the above stability criterion correctly predicted the survival of large amplitude oscillons in simulations. We further note that three-dimensional oscillons in sine-Gordon potentials (for axions, but without gravity) are not stable and have a relatively short lifetime, compared to flattened potentials Salmi and Hindmarsh (2012); Amin et al. (2012). Oscillons in flattened potentials can last longer than $10^{7}m^{-1}$Salmi and Hindmarsh (2012), whereas the duration of our simulations is $t_{\rm f}-t_{\rm i}\sim{\rm few}\times 10^{3}m^{-1}$. See the Appendix for further references on lifetimes in the relativistic case. A more detailed stability analysis including gravity for our saturated potentials would be useful.121212For a related analysis in the case of axions, see Hertzberg (2010); Visinelli et al. (2018).

IX Gravitational Clustering

For $\beta\ll 1$, gravitational clustering is expected to become important at late times (significantly after the solitons have formed, see eq. (12)). At these late times, this universe essentially behaves as a matter dominated universe ($a(t)\propto t^{2/3}$), with solitons becoming our new non-relativistic dust particles on scales much larger than their size. As a result, our zeroth order expectation is that the gravitational clustering of these solitons should proceed in a manner similar to dust in an expanding universe. Moreover, we can ignore non-gravitational forces between the solitons at separations much larger than $2r_{e}$ because we expect them to be Yukawa-like, with the force falling away exponentially with separation.131313This is also reminiscent of the force between solitons as analyzed by Manton (1979).

We construct the two-point correlation function of soliton locations obtained from our simulations to quantitatively investigate the effects of gravitational clustering. In Fig. 5, we show the two-point correlation function of the solitons, calculated with the Landy-Szalay estimator Landy and Szalay (1993); Wall and Jenkins (2012):

[TABLE]

where there are $\mathcal{N}$ solitons (the data $D$), and $\mathcal{N}$ uniform randomly chosen points $R$, and $DD$ is the number of soliton pairs in a given comoving radial separation bin, $RR$ is the mean count for the random points over several realizations $R$, and $DR$ is the cross-correlation statistic.

As seen in Fig. 5, the measured two point correlation function is the same for the cases with and without gravitational interactions at early times soon after soliton formation ($a\lesssim 4$). The distribution is close to Poissonian on large scales: $\xi_{\rm LS}(r\gtrsim 10)\approx 0$. However, the comoving scale $r_{\rm nl}\sim k_{\rm nl}^{-1}$ which is the typical separation of solitons when they first form manifests itself in a negative correlation function on small scales (we find very few solitons with separations less than $k_{\rm nl}^{-1}$).

If we allow for gravitational interactions, solitons begin to cluster. This clustering can be quantified in our simulations at late times as excess power in $\xi_{\rm LS}$ (for $a\gtrsim 10$). Consistent with clustering of point particles in a matter dominated universe starting with uncorrelated positions Saslaw (1980), we find

[TABLE]

where $r$ is a comoving separation.141414We checked that if we replace the solitons by point particles after $a=4$, the correlation function evolves in a qualitatively similar manner Fitting the model $\xi_{\rm LS}\propto a^{\alpha}r^{\beta}$ for our 6 simulations in the range of $a=10$ to $a=20$, we find $\alpha=1.7\pm 0.3$, $\beta=-2.1\pm 0.2$. It would be interesting to explore this clustering further in detail, since it might reveal differences from the point particle case at late times.

X Strong Soliton Interactions

Self-gravity plays the important role of bringing solitons together at late times (i.e., significantly after their formation), and allows them to interact.151515There are interactions at early times when gravity is ignored as well, but this is not so at late times in our simulations. We find that the some solitons have a significant velocity at early times with and without gravity, which will be investigated quantitatively in the future. Figure 6 shows three different types of interactions that are achieved from our cosmological initial conditions.

Solitons “repel or bounce off” each other when the relative phase of the interacting solitons $|\theta_{1}-\theta_{2}|\approx\pi$ where $\psi_{a}(t,{\bf x})=\Psi_{a}({\bf x})e^{-i(\nu_{a}t+\theta_{a})}$ with $a=1,2$. We have verified this phase structure in our simulations during such a repulsive interaction. 2. 2.

A few solitons merge to form more massive solitons (typically when the relative phase is $\sim 0$), resulting in a change in the number density of solitons. Such interactions are typically accompanied by the generation of a burst of scalar waves as the solitons settle into new configurations. 3. 3.

A small fraction of solitons form orbiting binaries, and we even see an occasional three-body interaction. 4. 4.

Only ${\rm few}-10\%$ of the number of solitons in our simulations undergo strong encounters per Hubble time.161616We inspected 6 numerical runs with different initial conditions to get this number. This is consistent with the rate of change in the comoving number density of solitons

[TABLE]

as seen from Fig. 3.

We re-iterate that bouncing, binary formation, and merging of solitons are self-consistently obtained from our cosmological initial conditions. Evidently, the dynamics of these strong interactions are quite rich, and deviate from the expectations of treating these solitons as just point particles. The relative phase of the solitons plays an important role in these close encounters.

We note that at late times ($a_{\rm f}\simeq 20$), we have about 10 pixels per linear dimension of the soliton ($\sim 10^{3}$ pixels per volume of the soliton). As a result, the detailed dynamics (such as post-interaction kicks at late times) of individual strong interactions should be interpreted with some care. While it is not easy to improve the resolution significantly for the entire simulation, zoom-in, higher resolution simulations focusing on soliton interactions using initial conditions from our simulations would be useful. A more detailed investigation of the rich dynamics of close encounters with higher resolution simulations is left for future work. For an early, and detailed investigation of $Q$-ball interactions (relativistic complex field valued analogs of our solitons), but without gravity, see Axenides et al. (2000); Battye and Sutcliffe (2000).

The repulsive and attractive behavior of such solitons as a function of relative phase can be heuristically understood as follows. Consider a probe soliton moving past another stationary soliton (in the absence of gravity). The nonlinearity in the Schrödinger equation ($\propto|\psi|^{2}$ for $|\psi|^{2}\ll 1$) can be thought of as a nonlinear refractive index.171717This is more than an analogy since nonlinear Schrödinger equations are used to model light pulse propagation in nonlinear media Aitchison et al. (1991); we learned of the above heuristic explanation from the same paper. For soliton formation and interactions in yet another context (Bose-Einstein condensates), see for example Nguyen et al. (2017). If the two solitons are in phase, we expect this term to be larger in the region between the solitons than the case when the stationary soliton is absent. It also increases towards the stationary soliton. As a result, this larger refractive index, and its gradient, will cause the core of the probe soliton to bend towards the stationary one; i.e. there will be attraction between the solitons. On the other hand, when our two solitons are out of phase, the $|\psi|^{2}$ between the two solitons will be smaller and have to go to zero in the middle (from symmetry), causing the probe soliton to move away from the stationary one (hence “repulsion”). A more detailed, effective potential based analysis at large separations is provided by Palenzuela et al. (2007); Cotner (2016).

XI Conclusions & Future Directions

We investigated the dynamics of non-relativistic scalar fields in an expanding background. By including self-interactions and gravitational interactions, we demonstrated the formation of solitons driven by self-interactions from cosmologically relevant initial conditions, followed by gravitational clustering of solitons. We showed that this clustering leads to dynamically rich interactions between solitons including scattering, merging and binary formation at late times (which is absent in the case when gravity is not included). The highly nonlinear dynamics were explored by numerically solving the Schrödinger-Poisson system of equations with self-interactions and weak field gravity in a self-consistently expanding universe.

We provided analytic results and estimates for (i) the time scales and length scales associated with soliton formation, (ii) the spatial distribution of solitons, (iii) the number density of solitons, (iv) the individual properties of our three-dimensional solitons, including their stability, and (v) the two-point function related to the gravitational clustering of solitons.

We showed agreement between our analytic calculations and numerical simulations. The estimates and analytic results also provide an understanding of how the results depend on essential physical parameters in our problem, allowing for broader applicability beyond that of the fiducial models considered in this paper. In the Appendix we discuss the connection of our work to the case where the fields satisfy a relativistic Klein-Gordon equation in an expanding universe (in particular, Lozanov and Amin (2019)). A more careful comparison with relativistic simulations, and many associated subtleties and caveats, is left for future work.

Our work points towards a number of new avenues of exploration: (1) What is the end state of a gravitationally and non-gravitationally interacting “soliton gas”? What is the velocity and angular momentum distribution? This investigation is not purely gravitational because of the close encounters of the solitons in an expanding universe, where the phase plays a dominant role (see Schwabe et al. (2016) for the non-interacting case). (2) For our initial conditions, individual solitons seem to be far from forming black holes. However, rare, accidental over-densities or over-densities driven by gravitational clustering and mergers might make it more favorable to form black holes. Numerically intensive, general relativistic simulations of soliton formation from cosmological initial conditions and strong self-interactions have not yet been done Helfer et al. (2017); Widdicombe et al. (2018). (3) The close encounters could be a source of stochastic gravitational waves from solitons in the early universe, in addition to those from formation of the solitons in the early universe. (4) It is possible to consider a different expansion history of the background (for example, radiation domination) and an axion-like potential as well as inhomogeneous initial conditions, which would make parts of our analysis relevant for the formation of quasi-stable axitons Kolb and Tkachev (1994) and axion miniclusters Hogan and Rees (1988) in the early universe.181818Radiation domination makes soliton formation and clustering more difficult starting from approximately homogeneous initial conditions.

Acknowledgments

We thank Kaloian D. Lozanov for help in comparing the results in this paper to relativistic simulations in an expanding background and for helpful comments on the manuscript, Marcos G. Garcia for useful discussions regarding the connection between non-relativistic/relativistic field theories, Rohith Karur for checking the importance of gravitational interactions in soliton profiles, and Eugene E. Lim for helpful comments on the draft. Part of the simulations were carried out on the Odyssey cluster supported by the FAS Division of Science, Research Computing Group at Harvard University. MA is supported by a DOE grant DE-SC0018216. Support (PM) for this work was provided by NASA through Einstein Postdoctoral Fellowship grant number PF7-180164 awarded by the Chandra X-ray Center, which is operated by the Smithsonian Astrophysical Observatory for NASA under contract NAS8-03060. Part of this work was carried out at the Aspen Center for Physics, which is supported by National Science Foundation grant PHY-1607611.

XII Appendix

XII.1 Connection to a Relativistic Model

In the main body of the paper we did not include a detailed analysis of the non-relativistic limit of strongly self-interacting relativistic theories (if it exists). We took certain non-relativistic field equations with strong self-interactions as given, and explored the solutions. As we discuss below, the equations we use can be justified as being obtained by integrating out the fast time variation in the weakly interacting limit. While we believe that some aspects of the relativistic – non-relativistic connection persists at large self-interactions as well, a rigorous mapping is beyond the scope of the present work. We note that even with strong self interactions, the spatio-temporal variations of the solutions of the system under consideration remain non-relativistic, and the gravitational potential remains small, making the exploration in the main body of the paper self-consistent in this respect.

To derive our equations of motion (LABEL:eq:MasterEq) from a relativistic scalar field theory (in a particular limit discussed below), consider a real scalar field $\phi$ within general relativity. Consider a real scalar field minimally coupled to gravity with the action

[TABLE]

where $\phi$ has dimensions of energy, $R$ is the Ricci scalar, $g_{\mu\nu}$ is the metric, $g$ is the determinant of the metric, and $d^{4}x=(cdt)d^{3}x$. We are interested in potentials of the form

[TABLE]

where $V_{\rm nl}(\phi)$ contains the non-quadratic part of the potential, whose shape is controlled by a scale $M$. As a concrete example, we can consider the potential $V(\phi)=(m^{2}M^{2}/2)\tanh^{2}(\phi/M)$ Kallosh and Linde (2013a, b), although the precise form is not necessary for most of the discussion that follows.

XII.1.1 The weak field approximation – non-expanding spacetime

In the weak field limit (i.e. for $\Phi/c^{2}\ll 1$ where $\Phi$ is the Newtonian gravitational potential) and in the absence of expansion, the metric is determined by the line element of the form

[TABLE]

Note that we are ignoring anisotropic stress, as well as vector and tensor perturbations. In the linear regime, anisotropic stress is absent and will be absent away from solitons at the very least. In a time averaged sense, the anisotropic stress will be small inside the solitons.

The equation of motion (nonlinear Klein-Gordon equation in curved spacetime) satisfied by the field $\phi$ is

[TABLE]

which, to leading order in $\Phi/c^{2}$, yields

[TABLE]

In turn, the Einstein equations reduce to the Poisson equation

[TABLE]

with

[TABLE]

XII.1.2 The non-relativistic limit

In order to consider the “non-relativistic” limit, it is convenient to redefine the real scalar $\phi$ in terms of a complex field $\psi$, factoring out the rest energy contribution

[TABLE]

where the normalization constant is chosen so that (27) reduces to the usual non-relativistic Poisson form. Straightforward substitution into (26) and (28) yields

[TABLE]

and

[TABLE]

Let $\tau_{m}=\hbar/mc^{2}$ and $\lambda_{m}=\hbar/mc$. Now let us assume that $|\tau_{m}\partial_{t}\psi|\ll|\psi|$; similarly $|\tau_{m}\partial_{t}\Phi|\ll|\Phi|$. We now average $T^{0}_{0}$ over a period $2\pi\tau_{m}$ assuming that $\Phi$ and $\psi$ do not change appreciably over this period.191919This part is not entirely rigorous, and it deserves to be handled with care. This yields

[TABLE]

It is convenient to define

[TABLE]

Assuming $|\Phi/c^{2}|\ll 1$, and $|\tau_{m}\partial_{t}\psi|\ll|\psi|$, the above expression simplifies to

[TABLE]

To get the Schrödinger-like equation, we multiply eq. (30) by $e^{it/\tau_{m}}$ and average over a period $2\pi\tau_{m}$, again assuming that $\Phi$ and $\psi$ do not change appreciably over this period. This temporal averaging will get rid of the h.c part in (30). Moreover, we divide the resulting equation by $mc^{2}$, to get

[TABLE]

We will show in the next subsection that $\frac{\tau_{m}}{\sqrt{2}}\langle e^{it/\tau_{m}}\partial_{\phi}V\rangle=\psi\partial_{|\psi|^{2}}U_{\rm nl}(|\psi|)^{2}=\psi U_{\rm nl}^{\prime}(|\psi|^{2})$.

Treating $\Phi/c^{2}$ and $\tau_{m}\partial_{t}$ as separate small parameters, and keeping leading order terms in each (but ignoring $\Phi/c^{2}\times\tau_{m}\partial_{t}$), we have

[TABLE]

After some re-arranging

[TABLE]

On the one hand it is clear that $\nabla^{2}\psi\gg(\Phi/c^{2})\nabla^{2}\psi$. However, for large amplitude solitons $\lambda_{m}^{2}\nabla^{2}\psi\sim\psi$; hence, it is not clear that we can drop this term compared to the $\psi$ in the term with the $\Phi$ coefficient (also see Namjoo et al. (2018)). A similar argument holds for $\langle\ldots\rangle$ terms. For the discussion that follows, we will use $\nabla^{2}\psi\gg(\Phi/c^{2})\nabla^{2}\psi$, to arrive at

[TABLE]

Using our result for $\langle T^{0}_{0}\rangle$, we also have the Poisson equation at the lowest order in $\Phi/c^{2}$

[TABLE]

These are our master equations used for time evolution of the field and for determining the metric potential (see eq. (LABEL:eq:MasterEq)). In arriving at eq. (LABEL:eq:MasterEq) in this limit (ignoring expansion for the moment), we assumed weak field gravity and restricted ourselves to scalar metric perturbations without anisotropic stress.

Since we were not interested in reproducing the limit of a particular relativistic theory in the main body of the text, we simply took $U_{\rm nl}$ to be an effective potential for our theory. Nevertheless, by using (33) we can link $U_{\rm nl}$ to $V_{\rm nl}$ at least for small amplitudes. We turn to this task next.

The time averaging procedure in eq. (33) is mathematically well defined for any potential which admits a Taylor expansion and has a quadratic minimum.

[TABLE]

The non-quadratic (nonlinear part) of this potential is

[TABLE]

Using $\phi=\sqrt{2}\Re[\psi e^{-it/\tau_{m}}]$ and taking a time average of this nonlinear part over a period $2\pi\tau_{m}$, we have

[TABLE]

We also need the time average of $\partial_{\phi}V_{\rm nl}$ for the equations of motion:

[TABLE]

It is beneficial to have a fitting function for $\langle V_{\rm nl}\rangle$. For a potential of the form

[TABLE]

we can find that for $|\psi|\leq\pi/(2\sqrt{2})$, an excellent approximation to $U_{\rm nl}(|\psi|^{2})$ is provided by

[TABLE]

Rescaling our field by $mM$, we recover the potential used in the main body of the text. We caution, that the form beyond $|\psi|>\pi/(2\sqrt{2})$ need not be simply connected to the relativistic potential. Moreover, at these large amplitudes, we might benefit by time-averaging over amplitude-dependent frequencies.

To include the effect of background expansion we consider a metric of the form

[TABLE]

where $a(t)$ is the scalefactor. Our complete set of equations then becomes (under the assumption that $H^{-1}\gg\tau_{m}$ and $c/H\gg\lambda_{m}$),

[TABLE]

where $\overline{[\ldots]}$ indicates a spatial average. The third equation is obtained from the Einstein equations for a homogeneous and isotropic universe (the Friedmann equation). This completes our derivation, with caveats, of the master equations (LABEL:eq:MasterEq) that are used in the main body of the paper.

For recent derivations and discussions of the non-relativistic limit, as well as decay rates for solitons with and without weak-field gravity (but in a non-expanding universe), see Namjoo et al. (2018); Eby et al. (2018).

XII.2 Details on Initial Conditions

In this appendix, we derive the vacuum initial conditions for a free non-relativistic field from the appropriate relativistic free field vacuum perturbations. Starting with the definition of the Fourier transform:

[TABLE]

where $\psi^{R,I}_{\bf k}(t)$ are the Fourier transforms of $\psi^{R,I}(t,{\bf x})$. Hence we have

[TABLE]

Similarly we have

[TABLE]

where we have assumed $|\dot{\psi}^{R,I}_{\bf k}(t)|/m\ll|\psi^{I,R}_{\bf k}(t)|$. At an initial time $t=0$, we have

[TABLE]

Now, following the implementation in Defrost Frolov (2008), we can write $\phi_{\bf k}(0)=b_{\bf k}/\sqrt{2\omega_{k}}$ and $\dot{\phi}_{\bf k}(0)=c_{\bf k}\sqrt{\omega_{k}/2}$ where $\langle b_{\bf k}b_{{\bf q}}^{*}\rangle=\delta_{{\bf k}{\bf q}}$ and $\langle c_{\bf k}c_{{\bf q}}^{*}\rangle=\delta_{{\bf k}{\bf q}}$ with $b_{\bf k}$ and $c_{\bf k}$ independent complex numbers (4 identically distributed independent variables). The real and imaginary parts of each are drawn from a zero mean Gaussian distribution, with a variance of $1/2$. Then, we have

[TABLE]

where in the second equality we assumed $k\ll m$ so that we have $\omega_{k}=\sqrt{k^{2}+m^{2}}\approx m$.

We are interested in $\psi_{\bf k}(t)$ which is the Fourier transform of $\psi(t,{\bf x})$. It can be written as

[TABLE]

Hence,

[TABLE]

with amplitude drawn from a Raleigh distribution and the phase drawn from a uniform distribution. This is consistent with the result in the main body of the paper.

XII.3 Comparison of Linear Instability Relativistic and non-relativistic Systems

The instability analysis discussed in the main text is connected to Floquet analysis in the corresponding relativistic theory (see for example, Lozanov and Amin (2018)). However, the instability bands as well as the Floquet exponents can differ from the relativistic case at large amplitudes and relativistic wave-numbers. For the relativistic version (with $a=1$, $H=0$), the perturbation to the homogeneous field satisfies:

[TABLE]

where the field $\phi$ is measured in units of $M$ and spacetime in units of $m^{-1}$. The periodic term in $V^{\prime\prime}_{\rm nl}(\bar{\phi})$ leads to growth of perturbations of the form $\delta\phi_{\bf k}\sim P_{{\bf k}}(t)e^{\Re[\tilde{\mu}_{k}]t}$ where $\tilde{\mu}_{k}$ are the Floquet exponents and $P_{{\bf k}}(t)$ are periodic functions. We find that $\mu_{k}\approx\Re[\tilde{\mu}_{k}]$ for $\bar{\phi},\bar{\psi}\ll 1$ and $k\lesssim 1$. The boundary of the non-relativistic band yields a good approximation to the relativistic case for $\bar{\psi}\lesssim 1$.

XII.4 Non-relativistic Solitons and Oscillons

It is worth making a comparison of our non-relativistic solitons discussed in Section X to the relativistic ones (oscillons). Recall that $\phi=\sqrt{2}\Re[\psi e^{-it}]$. For small amplitude solitons (oscillons), we expect $\phi(t,r)\approx\phi(r)\cos(\omega t)+\ldots$ (with $\omega<1$). The solitons in the nonlinear Schrödinger equation have the form $\psi(t,r)=\Psi(r)e^{-i\nu t}$. Hence

[TABLE]

where $\nu<0$. We caution the reader that this small amplitude analysis should merely be taken as a guide. The actual relativistic solitons can include multiple frequencies, including breathing modes at large amplitudes. We compared the profiles of relativistic solitons obtained from the simulations in Lozanov and Amin (2019), and found good qualitative agreement with Fig. 4, albeit with more scatter around the curve (after appropriate scaling of the parameters).

A numerical study of the lifetime and stability of large amplitude relativistic oscillons (but without gravitational interactions) in flattened potentials like the one we use here has been discussed in Salmi and Hindmarsh (2012). A more detailed connection between non-relativistic solitons and oscillons, as well as analysis of the stability of relativistic cases (typically for small amplitude) can be seen in Fodor et al. (2008); Hertzberg (2010); Mukaida et al. (2017); Saffin et al. (2014).

XII.5 Probability Density Functions of the Density and Gravitational Potential

The probability density function (PDF) of the energy density and the gravitational potential in our simulation is show in Fig. 8. Note that the gravitational potential in the simulation volume remains small $|\Phi|\ll 1$. Moreover the formation of the “shelf” in the density PDF($a\gtrsim 4$) is characteristic of the systems in which soliton formation takes place; the same qualitative behavior was seen when simulating relativistic systems with a related self-interaction potential Lozanov and Amin (2019) (see Fig 3. in that paper; however, note that $\beta\approx 8\times 10^{-3}$ in that figure). Note that $\beta\ll 1$ is required for the instability that generates solitons to be effective in a self-consistently expanding universe. The same $\beta$ also controls the strengths of the gravitational potential. This competition makes it difficult to generate individual solitons with large gravitational potentials via the self-interaction instability.