This online article is a supplement to our paper arXiv:1412.7152.

The object of study is a self-gravitating Bose-Einstein condensate, which hypothetically may form a dark matter star or a dark matter core inside a normal star.

The dimensionless GPE-Poisson system

The simulations presented on this page use a dimensionless form of the GPE, together with Poisson's equation for gravity:

\begin{align}i\frac{\partial\psi}{\partial t}&=-\frac{1}{2}\nabla^2\psi+V\psi+c|\psi|^2\psi-\mu\psi,\\ \nabla^2 V&=4\pi G|\psi|^2.\end{align}

This form of the GPE-Poisson system is possible if we choose units such that the reduced Planck constant is $\hbar=1$ and the particle mass is $m=1$.

These choices fix two out of three units of measurement (e.g., the units of mass and action), leaving the choice of the third unit (e.g., length) arbitrary.

We normalize the wavefunction such that $\int |\Psi|^2~dV=N$ is the total particle number. That is, after restoring units, $M=Nm$ is the total condensate mass.

Typically, for small $m$ and large $M$, $N$ becomes very large, which is computationally inconvenient. However, we note that the GPE is invariant under the following set of transformations:

\begin{align}t\rightarrow\lambda^2 t,\qquad r\rightarrow\lambda r,&\qquad V\rightarrow\lambda^{-2}V,\qquad\mu\rightarrow\lambda^{-2}\mu,\\ c\rightarrow\lambda^{-2}\kappa^{-2}c,\qquad G&\rightarrow\lambda^{-4}\kappa^{-2} G,\qquad\psi\rightarrow\kappa\psi.\end{align}

As a result of these rescalings the GPP system will read

\begin{align}\frac{\kappa}{\lambda^2}i\frac{\partial\psi}{\partial t}&=-\frac{1}{2}\frac{\kappa}{\lambda^2}\nabla^2\psi+\frac{\kappa}{\lambda^2}V\psi+\frac{\kappa}{\lambda^2}c|\psi|^2\psi-\frac{\kappa}{\lambda^2}\mu\psi,\\ \lambda^{-4}\nabla^2V&=\lambda^{-4}4\pi G|\psi|^2. \end{align}

In particular, choosing a value of $\kappa\ne 1$ amounts to rescaling $N\rightarrow\kappa^2N$ and $G\rightarrow\kappa^{-2}G$. We shall use this later on to normalize the values of $G$ and $N$ for computationa l convenience.

A Bose star

We wish to simulate a BEC star. Its mass has to be roughly the same as one solar mass, or about $2\times 10^{30}$ kg. The Schwarzschild radius corresponding to this mass is roughly 3 km. Our BEC star must be significantly larger than the Schwarzschild radius in order for the non-relativistic GPE to work. We choose $R=50$ km.

The simulation volume must be significantly larger than the simulated star, to ensure that the star does not interact too much with the boundaries of the volume. After experimentation, we choose a simulation box with coordinates that extend from $-440$ to $+440$ km.

Our simulation code divides space into cells. Too many cells slow down the simulation unacceptably. As a practical limit, we cannot use much more than 100 cells per one side, as the computation time increases as the cube of this number. For practical simulations, we use 80 cells, so the cell size is $11\times 11\times 11$ km.

If we take $11~{\rm km}=\frac{1}{2}\lambda$ to be the Compton half-wavelength of a particle, the corresponding mass is $m=h/c\lambda\simeq 10^{-46}~\rm{kg}\simeq 5.6\times 10^{-11}$ eV. This is the maximum Bose particle mass that we can simulate using a spatial resolution of 11 km.

The time that corresponds to a Compton wavelength of 22 km is approximately 70 μs; this determines the maximum allowable simulation timestep. In most simulations, we shall use τ = 10 μs as our simul ation timestep.

With these considerations in mind, we choose the following units for the purpose of simulating a BEC stellar remnant or stellar core:

1. We set our unit of mass, $[{\rm M}]=10^{-46}$ kg, the mass of a $5.6\times 10^{-11}$ eV Bose-particle.
2. We set $\hbar=1$. Therefore, $10^{-34}~{\rm m}^2{\rm kg}/{\rm s}=1~[{\rm L}]^2[{\rm M}]/[{\rm L}]=10^{-46}~[{\rm L}]^2{\rm kg}/[{\rm T}]$, or $1~{\rm m}^2/{\rm s}=10^{-12}[{\rm L}]^2/[{\rm T}]$.
3. We choose our unit of length: $[{\rm L}]=1~{\rm km}=1000~{\rm m}$. So $1~{\rm m}^2/{\rm s}=10^{-12}[1000~{\rm m}]^2/[{\rm T}]=10^{-6}{\rm m}^2/[{\rm T}]$. Therefore, $[{\rm T}]=10^{-6}~{\rm s}$.

So we measure masses, lengths and time in the following units:

\begin{align} 1~[{\rm M}]&=10^{-46}~{\rm kg}\simeq 5.61\times 10^{-11}~{\rm eV},\\ 1~[{\rm L}]&=1~{\rm km}=1000~{\rm m},\\ 1~[{\rm T}]&=10^{-6}~{\rm s}=1~\mu{\rm s}.\end{align}

In particular, the size of our simulation volume is measured in kilometers; the simulation step is measured in microseconds; and $N=2\times 10^{76}$ simulated particles constitute one solar mass.

In these units, $G=6.67\times 10^{-78}$. Furthermore, a Bose star with one solar mass, $M=2\times 10^{30}$ kg, contains $N=2\times 10^{76}$ particles.

We now choose $\kappa=10^{-38}$, allowing us to rescale $G\rightarrow 0.0667$ and $N\rightarrow 2$. It must be emphasized that the physical conclusions are unaffected by this rescaling. (However, this rescali ng not only ensures that intermediate calculations do not exceed the double-precision range, it allows all calculations to stay within range of single-precision floating point, making it possible in the future t o improve performance by a GPU-based single-precision reimplementation.)

Angular velocity is measured in radians per microsecond. For an $R=50~{\rm km}$ sphere, the equatorial circumference is $2\pi R\simeq 314~{\rm km}$. A tangential velocity of 300,000 km/s (speed of light) corresponds to 955 revolutions per second, or 6,000 radians per second. That is, 0.006 radians per microsecond (our unit of time). In order to remain within the approximately nonrelativistic regime, the refore, angular velocities must be less than $10^{-3}~[{\rm T}]^{-1}$.

Numerical stability

As we perform a simulation using a finite simulation volume, finite spatial grid size, and finite time step, it is important to characterize our understanding of the limitations of our simulation tool. Specifically, we need to understand how nonphysical parameters of the simulation, specifically the finite size of the simulation volume element, the finite timestep, and the finite extent of the simulation volume, affect the stability of the simulation. For this purpose, we ran four very long simulations, varying these nonphysical parameters only.

The parameters in these simulations were identical (the physical parameters described a compact one solar mass Bose star with a 50 km radius, simulated in a cubical volume with 880 km long edges), except for the nonphysical values of the time step and spatial resolution, which are summarized in the table below:

Below is a set of brief videos that show the density and phase evolution of the four cases.

High resolution simulations and rotation

To investigate the physics of a Bose-star, we ran four cases. In three cases, we investigated the effects of varying the rotation rate Ω; in the fourth case, we ran a simulation with a modest rotation rate, but much higher spatial resolution (at the cost of reducing the overall simulation volume.) All four cases ran for 1 second of simulated time, using a $60\times 60\times 60$ spatial grid.

The table below summarizes the four cases that were investigated in this manner:

Videos for these four simulation runs can be viewed below:

Finally, here is a video of of the phase evolution of 36 Bose stars, as described in our paper (arXiv:1412.7152). These simulations were performed using a simulation volume of 480×480×480 km and a spatial grid of 60×60×60, while varying $N$, $\Omega$, $R$ and $a$: