Galaxy Formation with BECDM: I. Turbulence and relaxation of idealised haloes

TL;DR

This paper investigates the properties of Bose-Einstein condensate dark matter haloes, revealing their turbulent behavior, core-halo relations, and energy distributions through theoretical analysis and simulations.

Contribution

It provides a detailed analysis of BECDM halo structure, turbulence, and core-halo mass relations, advancing understanding of ultralight bosonic dark matter.

Findings

Haloes have a soliton core and NFW-like outer profile.

Core-to-halo mass relation: M_c/M ≈ 2.6 Ξ^{1/3}.

Turbulence exhibits a k^{-1.1} velocity power spectrum.

Abstract

We present a theoretical analysis of some unexplored aspects of relaxed Bose-Einstein condensate dark matter (BECDM) haloes. This type of ultralight bosonic scalar field dark matter is a viable alternative to the standard cold dark matter (CDM) paradigm, as it makes the same large-scale predictions as CDM and potentially overcomes CDM's small-scale problems via a galaxy-scale de Broglie wavelength. We simulate BECDM halo formation through mergers, evolved under the Schr\"odinger-Poisson equations. The formed haloes consist of a soliton core supported against gravitational collapse by the quantum pressure tensor and an asymptotic $r^{-3}$ NFW-like profile. We find a fundamental relation of the core=to-halo mass with the dimensionless invariant $\Xi \equiv \lvert E \rvert/M^3/(Gm/\hbar)^2$ or $M_{\rm c}/M \simeq 2.6 \Xi^{1/3}$, linking the soliton to global halo properties. For $r \geq 3.5 \,r_{\rm c}$ core radii, we find equipartition between potential, classical kinetic, and quantum gradient energies. The haloes also exhibit a conspicuous turbulent behavior driven by the continuous reconnection of vortex lines due to wave interference. We analyse the turbulence 1D velocity power spectrum and find a $k^{-1.1}$ power-law. This suggests the vorticity in BECDM haloes is homogeneous, similar to thermally-driven counterflow BEC systems from condensed matter physics, in contrast to a $k^{-5/3}$ Kolmogorov power-law seen in mechanically-driven quantum systems. The mode where the power spectrum peaks is approximately the soliton width, implying the soliton-sized granules carry most of the turbulent energy in BECDM haloes.

Figures

Equations

$$i\hslash \frac{\partial \psi }{\partial t}=-\frac{{\hslash }^2}{2m}{\nabla }^2\psi +mV\psi$$

$${\nabla }^{2}V=4\pi G(\rho -\overline{\rho })$$

$$\psi =\sqrt{\rho }{\mathrm{e}}^{iS/\hslash },\mathbf{v}=\nabla S/m,$$

$$\frac{\partial \rho }{\partial t}+\nabla \cdot (\rho \mathbf{v})=0$$

$$\frac{\partial \mathbf{v}}{\partial t}+\mathbf{v}\cdot \nabla \mathbf{u}=-\frac{1}{m}\nabla (Q+V)$$

$$Q=-\frac{{\hslash }^2}{2m}\frac{{\nabla }^2\sqrt{\rho }}{\sqrt{\rho }}.$$

$$-\frac{1}{m}\nabla Q\equiv \frac{1}{\rho }\nabla \cdot {\mathsf{p}}_{\mathrm{Q}}$$

$${\mathsf{p}}_{\mathrm{Q}}=-{\left(\frac{\hslash }{2m}\right)}^2\rho \nabla \otimes \nabla \ln \rho$$

$$M=\int \rho d^{3}x,$$

$$E$$

$$=$$

$$\mathbf{L}=\int \mathbf{r}\times \rho \mathbf{v}{d}^{3}x$$

$${\rho }_{\mathrm{soliton}}(r)\simeq {\rho }_0{\left[1+0.091\times {\left(\frac{r}{r_{\mathrm{c}}}\right)}^2\right]}^{-8}$$

$${\rho }_0\simeq 3.1\times 10^{15}{\left(\frac{2.5\times 10^{-22}\mathrm{eV}}{m}\right)}^2{\left(\frac{\mathrm{kpc}}{r_{\mathrm{c}}}\right)}^4\frac{M_{\odot }}{{\mathrm{Mpc}}^3}$$

$$\left\{t,x,V,\psi ,\rho \right\}\to \left\{{\lambda }^{-2}\hat{t},{\lambda }^{-1}\hat{x},{\lambda }^2\hat{V},{\lambda }^2\hat{\psi },{\lambda }^4\hat{\rho }\right\}$$

$$\left\{M,E,L\right\}\to \left\{\lambda \hat{M},{\lambda }^3\hat{E},\lambda \hat{L}\right\}$$

$$\left\{t,x,V,\psi ,\rho \right\}\to \left\{\alpha \hat{t},\hat{x},{\alpha }^{-2}\hat{V},{\alpha }^{-1}\hat{\psi },{\alpha }^{-2}\hat{\rho }\right\}$$

$$\left\{M,E,L\right\}\to \left\{{\alpha }^{-2}\hat{M},{\alpha }^{-4}\hat{E},{\alpha }^{-3}\hat{L}\right\}.$$

$$\Xi \equiv |E|/M^3/(Gm/\hslash {)}^2.$$

$$V=\mathrm{ifft}\left[-\mathrm{fft}\left[4\pi G(\rho -\overline{\rho })\right]/k^2\right]$$

$$\psi \leftarrow \exp \left[-i(\Delta t/2)(m/\hslash )V\right]\psi$$

$$\hat{\psi }=\mathrm{fft}\left[\psi \right]$$

$$\hat{\psi }\leftarrow \exp \left[-i\Delta t(\hslash /m)k^2/2\right]\hat{\psi }$$

$$\psi \leftarrow \mathrm{ifft}\left[\hat{\psi }\right]$$

$$\Delta t\le \max \left[\frac{m}{6\hslash }(\Delta x{)}^2,\frac{h}{m|V{|}_{\max }}\right]$$

$$M_{\mathrm{c}}(r_{\mathrm{c}})=3.59\times 10^7{\left(\frac{2.5\times 10^{-22}\mathrm{eV}}{m}\right)}^2{\left(\frac{r_{\mathrm{c}}}{\mathrm{kpc}}\right)}^{-1}{M}_{\odot }$$

$${\psi }_{\mathrm{soliton}}(r,t)=e^{-i\hat{\gamma }{\lambda }_0^{2}t/(\hslash m{c}^2)}\sqrt{{\rho }_{\mathrm{soliton}}(r)}$$

$${\lambda }_0={\left(\frac{{\rho }_0}{m^2{c}^4/(4\pi G{\hslash }^2)}\right)}^{1/4}$$

$$T$$

$$=$$

$$M_{\rm c}/M\simeq 2.6\Xi^{1/3}=2.6\biggr{(}\frac{\lvert E\rvert}{M^{3}(Gm/\hbar)^{2}}\biggl{)}^{1/3},$$

Full text

††pagerange: Galaxy Formation with BECDM: I. Turbulence and relaxation of idealised haloes–Galaxy Formation with BECDM: I. Turbulence and relaxation of idealised haloes††pubyear: 2017

1 Introduction

The $\Lambda$ cold dark matter ($\Lambda$CDM) model has been very successful at describing the large scale structure of our Universe, including the statistical properties of the cosmic microwave background (CMB), and the cosmic web of galaxies across the ages (Springel et al., 2005). State-of-the-art $\Lambda$CDM simulations – e.g., Illustris (Vogelsberger et al., 2014a, b), Eagle (Schaye et al., 2015), Horizon-AGN (Dubois et al., 2014) – include complex modelling of stellar and gas (baryonic) physics that give rise to the galactic population and provide quantitative predictions over cosmological volumes in the non-linear density regime of density contrast for essentially the entire range of mass scales relevant for galaxy formation.

The $\Lambda$CDM model has a nearly scale-invariant spectrum of density fluctuations with substantial power on small mass scales. While the behavior of the power spectrum on extremely small scales depends on the specific physics of the dark matter particles, generic models based on weakly-interacting massive particles have power spectra that extend without suppression all the way to Earth masses (Green, Hofmann & Schwarz, 2004). This feature of CDM models – significant power at small scales – is the source of a number of enduring inconsistencies with galaxy population statistics observations, including the deficit of dwarf galaxies (the missing satellites problem; Klypin et al. 1999; Moore et al. 1999) and the problem with the abundance of isolated dwarfs (Zavala et al., 2009; Papastergis et al., 2011; Klypin et al., 2015), as well as the too-big-to-fail problem (Boylan-Kolchin, Bullock & Kaplinghat, 2011, 2012) and the cusp-core problem (Moore, 1994; Flores & Primack, 1994; Gentile et al., 2004; Donato et al., 2009; de Blok, 2010).

One widely-considered way to resolve these problems is baryonic feedback, which can both alter dark matter gravitational potential wells via supernovae or black hole feedback and prevent the formation of galaxies in low-mass dark matter haloes (Brooks et al., 2013; Oñorbe et al., 2015). Feedback therefore provides a natural mechanism for breaking the scale-free nature of the underlying dark matter density perturbation spectrum. Given the null detection of dark matter particles, however, we must also consider the possibility that the nature of dark matter is different from that of CDM on small scales, meaning that non-baryonic physics (perhaps in combination with baryonic effects) breaks the scale-free nature of dark matter.

The introduction of dark matter self-interactions (Spergel & Steinhardt, 2000), for example, can alleviate these small-scale problems (Vogelsberger, Zavala & Loeb, 2012; Zavala, Vogelsberger & Walker, 2013; Rocha et al., 2013; Vogelsberger & Zavala, 2013; Vogelsberger et al., 2014c; Elbert et al., 2015; Vogelsberger et al., 2016). Alternatively, a cutoff in the primordial power spectrum can be caused by free-streaming in warm dark matter (WDM) models (Bond, Szalay & Turner, 1982; Bode, Ostriker & Turok, 2001), potentially matching observations better than a pure CDM model (e.g., see Schneider et al. (2016) on the missing dwarfs problem and Lovell et al. (2017) on WDM and too-big-to-fail problem). A cutoff in the power spectrum can also be caused by collisional Silk-like damping if there are significant interactions between dark matter and relativistic particles in the early Universe (e.g. Boehm et al. 2002; Buckley et al. 2014; Foot & Vagnozzi 2016). These type of models offer non-baryonic solutions to the small-scale problems of CDM (e.g. Boehm et al. 2014; Vogelsberger et al. 2016).

An ultralight bosonic scalar field is a completely different and intriguing alternative to the CDM paradigm. A bosonic fluid with a particle mass of $m\sim 10^{-22}~{}{\rm eV}$ suppresses small-scale structure in the Universe owing to macroscopic quantum properties (the uncertainty principle) exhibited by the fluid (Peebles, 2000; Hu, Barkana & Gruzinov, 2000; Lee & Koh, 1996). The typical de Broglie wavelength for such a particle – the largest scale at which quantum mechanical effects will appear – is $\lambda_{\rm dB}\sim 1~{}{\rm kpc}$ (Chavanis, 2011; Suárez, Robles & Matos, 2014), similar to the observed sizes of the stellar distributions in dwarf galaxies (Marsh & Silk, 2014; Bozek et al., 2015). This type of fluid has a high temperature of condensation (critical temperature is $\sim$ TeV due to the high number density of axion particles), which becomes much larger than the mean field temperature as the Universe expands and cools. The fluid will thus form a cosmological Bose-Einstein condensate (BEC; Lundgren et al. 2010; Robles & Matos 2013; Matos, Vázquez-González & Magaña 2009, i.e., a superfluid) in the early Universe. On large scales, however, the scalar field behaves just like a collisionless self-gravitating fluid, identical to CDM, and is therefore consistent with modern large-scale cosmological constraints (Matos, Vázquez-González & Magaña, 2009; Suárez & Chavanis, 2017; Li, Rindler-Daller & Shapiro, 2014). We note here that BECDM is often referred to by a variety of other names, including scalar field dark matter, axion dark matter, fuzzy dark matter, quantum dark matter, and $\psi$DM (for recent reviews on BECDM and its astrophysical and cosmological implications, see Suárez, Robles & Matos 2014; Hui et al. 2017; Marsh 2016).

Aside from the possible astrophysical relevance of BECDM, theoretical physics offers motivation for the existence of ultralight scalar-field particles. The axion was first postulated in Peccei-Quinn theory (Peccei & Quinn, 1977; Weinberg, 1978), as an ultralight scalar particle that resolves the strong CP problem in QCD. String-theory compactifications also predict a class of ultralight axions (Arvanitaki et al., 2010). Such particles are candidates for BECDM.

Fully cosmological simulations of BECDM are in their infancy, in part owing to the demanding constraint on spatial resolution required to evolve kpc scale quantum fluctuations throughout the domain. Woo & Chiueh (2009) simulated axion dark matter in 1 Mpc boxes, and Schive, Chiueh & Broadhurst (2014) presented cosmological simulations of BECDM ($2$ Mpc periodic box, dark matter only) with an adaptive mesh and sufficient resolution to characterise the universal soliton-cored haloes that form. Schwabe, Niemeyer & Engels (2016) performed a detailed investigation of idealised merger simulations, documenting a comprehensive parameter study of two-soliton interactions in a non-periodic box. Although there are analytical and numerical solutions with spherical symmetry that suggest BECDM can solve the CDM’s small-scale problems by forming haloes with central cores (Sin 1994; Ji & Sin 1994; Ureña López 2002; Guzmán & Ureña López 2004; Robles & Matos 2012; Martinez-Medina, Robles & Matos 2015; Matos & Robles 2016, but see Slepian & Goodman (2012) for discussion of the impact of finite-temperature effects), numerical simulations with no assumed symmetries are required to confirm these results in a more realistic scenario.

One major consequence of the macroscopic quantum mechanical effects in a BECDM superfluid is that the fluid admits stable, minimum-energy soliton configurations known to form at the centers of self-gravitating haloes (Gleiser, 1988; Seidel & Suen, 1994; Balakrishna, Seidel & Suen, 1998; Guzmán & Ureña López, 2004, 2006). These kpc-scale soliton cores offer one possible solution to the well-known “cusp-core problem” of CDM. Moreover, these solitons are attractor solutions, i.e., they are solutions to which initially unstable configurations that undergo perturbations will converge (Seidel & Suen, 1990, 1994; Lee, 1989; Hawley & Choptuik, 2000; Gleiser & Watkins, 1989; Chavanis, 2011, 2016). Given the robust stability of the ground state solitons, it is expected that they survive even after merging with other solitons, which was also noticed in Schive, Chiueh & Broadhurst (2014).

Given the difficulties associated with numerical simulations of BECDM, most work on the subject thus far has relied heavily on analytic theory. Recently, Gonzáles-Morales et al. (2016) carried out an analysis to constrain the boson mass $m$ by fitting the luminosity-averaged velocity dispersion of dwarf spheroidal galaxies (dSphs) using an analytic soliton core dark matter profile. Assuming radial symmetry and that the halo has not been modified by baryonic physics, they placed an upper limit on the ultralight boson mass: $m<0.4\times 10^{-22}~{}{\rm eV}$ at the $97.5$ per cent confidence level. This is a tighter bound than the earlier study of dSph constraints (Marsh & Pop, 2015) based on the “mass-anisotropy degeneracy” in the Jeans equations that leads to biased bounds on the boson mass in galaxies with unknown dark matter halo profiles. Recently Ureña-López, Robles & Matos (2017) found that dwarf galaxies in the Milky-Way and Andromeda are consistent with $m\sim 10^{-21}\rm eV$. The boson mass constraint from cosmological structure formation requires $m>10^{-23}~{}{\rm eV}$ to create a relevant cut-off in the power spectrum on small scales and remain consistent with large-scale observations (Bozek et al., 2015), such a constraint comes from the Hubble Ultra Deep Field UV-luminosity function and the optical depth to reionisation as measured from CMB polarisation. Hlozek et al. (2015) establishes a constraint of $m>10^{-24}~{}{\rm eV}$ that comes from CMB temperature anisotropies, which is a more robust lower bound as it depends only on linear physics and less modelling of, e.g., star formation rate and halo mass functions. Lyman-$\alpha$ constraints suggest $m>3.3\times 10^{-22}~{}{\rm eV}$ at the $2\sigma$ level, assuming analogies between BECDM and WDM (Hui et al., 2017). This is in moderate tension with the most recent dSph results of Gonzáles-Morales et al. (2016) but consistent with Ureña-López, Robles & Matos (2017). Numerical simulations of BECDM systems are vital to make important progress, including testing the assumptions behind the analytic estimates and possibly ruling out or confirming the BECDM model. A first attempt at modelling the Lyman-$\alpha$ flux power spectrum cutoff with hydrodynamical simulations (Iršič et al., 2017) suggests a constraint of $m<10^{-22}~{}{\rm eV}$; however, these simulations do not include the full quantum effects of BECDM, only its effects on the initial power spectrum.

This paper is the first of a series aimed at quantifying the small-scale effects of BECDM in a cosmological context. We use idealised numerical simulations to analyse previously unexplored properties of relaxed/virialized BECDM haloes that are not covered by analytic work. We also present a numerical algorithm to simulate BECDM, which will be coupled with baryonic physics in fully cosmological simulations in future work. We compare the profiles of our simulated haloes to the analytic soliton relation and also study in more detail the granular quantum fluctuations that are present in the simulations but absent in analytical modeling. We show that the source of these fluctuations is quantum turbulence, an effect that has been seen in non-self-gravitating BEC systems (Kobayashi & Tsubota, 2005; Baggaley et al., 2012; Baggaley, Laurie & Barenghi, 2012; Tsatsos et al., 2016). The turbulence in our haloes, arising from reconnections of quantum vortex lines that form during the merging of haloes, is similar to what is found in dissipationless quantum superfluids with isotropic turbulence.

The manuscript is organised as follows. In Section 2 we discuss the theoretical background and formulation of BECDM haloes. Our code and simulation setup is described in Section 3. Section 4 considers the properties of radially-averaged profiles of virialized structures. Section 5 presents the scaling of the resulting soliton core mass and radius with other fundamental parameters of the system. Section 6 explores the turbulent properties of the virialized haloes, with a focus on the velocity power spectrum. Concluding remarks are given in Section 7.

2 Theoretical background

An ultra light scalar field of spin-0 at zero temperature is described in the non-relativistic limit by the Schrödinger-Poisson (SP) equations (Seidel & Suen, 1990; Sin, 1994; Lee & Koh, 1996; Hu, Barkana & Gruzinov, 2000; Suárez, Robles & Matos, 2014):

[TABLE]

[TABLE]

where the density of the fluid is defined as $\rho=\lvert\psi\rvert^{2}$ and $\overline{\rho}$ is the mean density. Here $m$ is the mass of the boson, $\psi$ is the wave-function of the particles normalised so its square norm is the density, and $V$ is the gravitational potential. In such a system, all particles share a common wave function $\psi$, hence the physical density of the fluid traces the probability density distribution $\lvert\psi\rvert^{2}$.

The physical system can be studied in the field or fluid representation (Suárez & Matos, 2011) via the Madelung transformation (Madelung, 1927; Chavanis, 2011):

[TABLE]

which yields the fluid representation of the SP equations:

[TABLE]

[TABLE]

where $Q$ is the quantum potential:

[TABLE]

Equivalently, one can equate

[TABLE]

in order to define a non-local quantum pressure tensor

[TABLE]

This quantum pressure tensor offers support against collapse from self-gravity. The support is non-local, as it depends on the gradient of the density.

The system has conserved quantities, including the total mass of the system

[TABLE]

and total energy of the system

[TABLE]

The total (quantum) kinetic energy is $K=K_{\rho}+K_{v}$, where $K_{v}$ is the classical contribution and $K_{\rho}$ is the gradient energy due to the quantum pressure tensor. $W$ is the potential energy. Quantum systems, like their classical counterparts, also obey a (well-known) virial (Ehrenfest) theorem: $0=2\langle K\rangle+\langle W\rangle$. In addition, the total angular momentum

[TABLE]

is also a conserved quantity. The systems that we consider in our simulations have no net angular momentum ($\mathbf{L}=\mathbf{0}$). Note that in most previous analytical works the systems studied are assumed to be spherically symmetric and in equilibrium, with a harmonic time-dependent phase ($S=S(t)$). In general, $S=S(\mathbf{x},t)$, then spatial wave interference may lead to regions of large variations in velocity, which seeds turbulence in the fluid. These nonlinear effects can only be studied with numerical simulations such as the ones presented here.

The SP equations admit stable soliton solutions, in which the Heisenberg uncertainty principle/quantum pressure tensor essentially supports the core against collapse under self-gravity. The soliton core solutions are well-approximated by

[TABLE]

(Schive et al., 2014), where $r_{\rm c}$ is the core radius and $\rho_{0}$ is the central density given by:

[TABLE]

The analytic profile fit to the soliton has a flat slope at the center (a ‘core’), and approaches a slope of $r^{-16}$ at the outskirts. This can be compared with the NFW (Navarro, Frenk & White, 1996) profile for CDM, which has an $r^{-1}$ cuspy center and a fall-off going as $r^{-3}$ at large radii.

The SP equations admit a scaling relation, with scaling parameter $\lambda$ (Ji & Sin, 1994):

[TABLE]

Owing to the scaling, the total mass, energy, and angular momentum of the system transform as

[TABLE]

Furthermore, the system may also be scaled by boson mass $m\to\alpha m$ as:

[TABLE]

[TABLE]

Since we assume no net angular momentum in our simulations, the system is then primarily characterised by a single invariant, $\lvert E\rvert/M^{3}$ (Schwabe, Niemeyer & Engels, 2016), which is unchanged under the scaling symmetry. To make it dimensionless, we define the invariant quantity as:

[TABLE]

Note that our definition for $\Xi$ is invariant not just under the $\lambda$ scaling of the SP equations but also with the boson mass (Eqn. 18), which makes our results even more general and scalable to any boson mass.

BEC superfluid systems (such as superfluid liquid helium) are known to exhibit turbulent behavior in a number of regimes (Kobayashi & Tsubota, 2005; Baggaley et al., 2012; Baggaley, Laurie & Barenghi, 2012; Tsatsos et al., 2016). Turbulence in BEC systems is a young and developing field (for a recent review see Tsatsos et al. 2016). Many such systems are described by the Gross-Pitaevskii equations, which are the Schrödinger equations with a non-linear self-interaction term, and no self-gravity (instead, a static potential ‘trap’ is often assumed). Turbulence is possible due to the advective term $\mathbf{u}\cdot\nabla\mathbf{u}$ in the fluid formulation of the governing equations, which is also the case in classical fluid dynamics.

Quantum turbulence is different from its classical manifestation, however (Baggaley et al., 2012; Baggaley, Laurie & Barenghi, 2012; Tsatsos et al., 2016). A direct consequence of the definition of the fluid velocity $\mathbf{v}=\nabla S/m$ is that the flow is irrotational: $\nabla\times\mathbf{v}=0$. However, there is an exception if $\psi$ is not continuous or does not have first or second derivatives. Vorticity in a quantum fluid is thus restricted to degenerate vortex lines or cores, where the velocity may diverge to infinity but the density is zero thus the solution remains physical. These filamentary vortex line structures naturally reconnect and create Kelvin waves that mediate the cascade of energy to smaller scales (Tsatsos et al., 2016). Vortex generation is possible from configurations initially smooth and without vortices (Galati et al., 2013), but a Kelvin’s conservation of circulation theorem applies to the system, placing topological constraints on the vortex lines that are allowed to form.

Turbulence in Gross-Pitaevskii fluids is known to develop Kolmogorov-like $k^{-5/3}$ velocity power spectra when the fluid is mechanically driven on the largest scale of $k$ (Baggaley et al., 2012) (mechanically driven here refers to driving the BEC fluid by grids or propellers): in this case turbulence arising from large spatial-scale modes cascades to smaller scales. In contrast, the spectrum of thermally-driven turbulence in a BEC fluid (i.e., driven by a small heat flux which introduces a counterflow velocity in the superfluid. There is no stirring length-scale introduced into the problem), which lacks energy on the largest scales, exhibits a ‘bump’ in the velocity power spectrum at intermediate scales (the inter-vortex lengthscale) and has been shown to scale as $k^{-1}$ at large wave numbers (Baggaley et al., 2012; Baggaley, Laurie & Barenghi, 2012; Tsatsos et al., 2016). Chiueh et al. (2011) also find $k^{-1}$ steady-state power-spectrum in the linear Schrödinger equations on small scales due to vortex reconnection. The connection between classical and quantum turbulence is an interesting and growing field (Tsatsos et al., 2016). Although understanding the various sources of turbulence is a challenging task, it has been noticed that when the turbulence arises from small to large scale the velocity spectrum is $k^{-1}$ which is characteristic of fluid with isotropic turbulence.

As some of the behavior of superfluids can be quite complicated to capture analytically, especially turbulence, we rely on numerical simulations to study BECDM haloes as described in the next section.

3 Simulations

We simulate $100$ scenarios of a group of soliton cores that merge to form a final virialized halo to study the general properties of virialized BECDM structures that are statistically meaningful. The simulation setup is based on idealised simulations found in Schive et al. (2014); Schwabe, Niemeyer & Engels (2016), which suggest that the initial conditions to form virialized haloes are largely unimportant for the final outcome. In our simulations we mainly create a virialized halo by merging soliton cores of different initial sizes and central densities, but we have also verified, using additional simulations, that we produce consistent types of cored-rather-than-cuspy final BECDM virialized haloes if we merge initially cuspy NFW profiles. Therefore, the simulations are a very useful tool to study the final product of the relaxation of BECDM haloes. The simulations are evolved with a pseudo-spectral method for solving the SP equations in 3-D. Details of the numerical method are provided in the following subsection.

3.1 Pseudo-Spectral Method

We developed a second-order pseudo-spectral solver for the SP equations, which we have also added (Mocz et al. in prep.) into the arepo code (Springel, 2010). Our approach employs a ‘kick-drift-kick’ techniques, akin to symplectic leapfrog $N$-body solvers (Springel, 2005). The wavefunction $\psi$ is evolved with unitary ‘kick’ and ‘drift’ operators.

The variables $\psi(x)$, $\rho(x)$, $V(x)$, are discretized onto a grid of dimension $N^{3}$. In this subsection, the variables will represent the discretized grid versions. Given the density field $\rho$, the potential $V$ can be calculated by transforming to Fourier-space and back:

[TABLE]

where ${\rm fft}\left[\cdot\right]$ and ${\rm ifft}\left[\cdot\right]$ are Fourier transform and inverse Fourier transform operators, respectively, and $k$ are the wave numbers at the corresponding grid locations.

First, the wavefunction is given a ‘kick’ by half a timestep, due to the potential:

[TABLE]

This is followed by a full ‘drift’ (kinetic) step in Fourier-space:

[TABLE]

[TABLE]

[TABLE]

The timestep is completed with another ‘kick’ step (Eqn. 22), and the system is thus evolved from time $t$ to time $t+\Delta t$.

The valid timestep criterion for stability and accuracy of our method, essentially a Courant-Friedrichs-Lewy (CFL) like condition, is that the unitary operators in Equations 22 and 24 do no change the phase by more than $2\pi$ in each timestep. The timestep criterion of Schwabe, Niemeyer & Engels (2016):

[TABLE]

enforces this property ($\lvert V\rvert_{{\rm max}}$ is the maximum of the absolute value of the gravitational potential). Note the timestep scales as $(\Delta x)^{2}$ rather than $\Delta x$ for gravity and Eulerian fluid solvers, which adds computational cost to the simulations.

We briefly compare the differences in existing codes that have been used to solve the SP equations to simulate BECDM. The advantages of our method include: simplicity, use of unitary operators, a ‘kick-drift-kick’ formulation which makes the method readily integrable into a number of existing cosmological codes, and machine precision control of the total kinetic energy during the drift step. This pseudo-spectral method achieves spectral (exponential) convergence in space and second-order convergence in time. The main limitation of pseudo-spectral methods in general is their restriction to a regular grid. Because of this, Schive, Chiueh & Broadhurst (2014) solve the SP equations on an adaptively refined mesh to achieve high dynamic range in cosmological simulations. The implementation requires Taylor expansion of the unitary operators with modified coefficients in order to minimise the small-scale numerical damping. Schwabe, Niemeyer & Engels (2016) use a 4th-order Runge-Kutta finite-difference solver on a grid to solve the SP equations. In general, our code yields consistent results with those obtained in alternative codes in the comparable regime.

3.2 Simulation Setup

We simulate bound systems with various total energy $E$ and total mass $M$, both of which are conserved in the total system which consists of a cubic box of $1~{}{\rm Mpc}$ on a side. We enforce periodic boundary conditions so there is no loss of energy or mass in the simulation. The periodic box allows us to account for incoming waves from all directions, which is closer to the cosmological case (Schive, Chiueh & Broadhurst, 2014) where waves from other haloes at larger distances would also interfere with a given host halo. A similar setup was used in Schive et al. (2014) where they studied idealised merger simulations with a smaller suite. Notably, we find a different fundamental scaling between the core mass and global quantities of the system, and will give some possible explanations for this variation below. Additionally, our study is different and complementary to Schwabe, Niemeyer & Engels (2016), where they analysed mergers of binary solitons in a finite volume where no wave reflection at the boundaries was allowed, they use the same boundary conditions as analytical studies making their results more comparable to those expected for isolated systems.

The primary variable that defines our systems is the invariant (under $\lambda$) $\lvert E\rvert/M^{3}$. For the initial condition at $t=0$, we randomly place between $4$ and $32$ cores with randomly selected soliton radii $r_{\rm c}\in[8,50]~{}{\rm kpc}$, allowing for multiple mergers at any time. We assume no phase offsets in the wavefunction between the cores.

The simulation uses internal units of $[{\rm L}]={\rm Mpc}$, $[{\rm M}]=M_{\odot}$, $[{\rm v}]={\rm km}~{}{\rm s}^{-1}$. Our highest resolution simulations have a resolution of $512^{3}$ cells and are run until a time $t_{\rm end}=10$ (internal units; the physical units are scalable and we rescale all the simulations presented in the paper to a Hubble time), which gives more than sufficient time for the dark matter structure to virialize over many dynamical timescales. The smallest final cores found in our simulations are resolved by at least 4 cells per linear dimension. In our results and analysis, we rescale the simulations with the scaling parameter $\lambda$ (Eqn. 16), chosen so that the total simulation time is the Hubble time $t_{\rm H}$, hence $\lambda=\sqrt{t_{\rm H}/t_{\rm end}}=26$. We stress that the scaling symmetry is a very important feature of the SP equations and our results can be rescaled to other halo masses in a straightforward way. By the scaling symmetry, lower mass haloes take a longer time to virialize. In the internal units of our code, the typical halo masses of our simulations fall into the range of few $\times 10^{7}~{}M_{\odot}$ to few $\times 10^{9}~{}M_{\odot}$, although we can scale our solutions by using $\lambda$. However, for generality, most of our results are presented in a way in which the scaling is factored out and are thus directly applicable to any mass halo. The core properties in the subsequent sections are analysed at the final time of one Hubble time. The simulations were scaled to such a time to demonstrate that halos of such masses become virialized well before the Hubble time.

For the simulations, we used a boson mass of $m=2.5\times 10^{-22}~{}{\rm eV}$, the same as the fiducial value in Schwabe, Niemeyer & Engels (2016). We stress once again, that the results may be scaled with mass $m\to\alpha m$ according to Eqns. 18 and 19.

Fig. 1 shows a volume rendering of the density field of one of our simulations at $4$ different times ($t=0,0.1t_{\rm H},0.2t_{\rm H},t_{\rm H}$), plotted with yt (Turk et al., 2011). As the solitons merge, they interfere quantum mechanically and create waves and interference patterns in the fluid. Shown also are the energy components of the fluid in the box as a function of time. The total energy $E$ is conserved to within a few percent by our method. The quantum gradient energy dominates over the classical kinetic energy, and provides the support against gravitational collapse on small scales. The system is in approximate virial equilibrium: $0\simeq 2\langle K_{\rho}+K_{v}\rangle+\langle W\rangle$.

4 Halo profiles

The primordial dark matter solitons in our merger simulations all lead to the formation of final haloes with a central soliton core, which is well-described by the analytic form found in previous BECDM simulations in the literature (Eqn. 14). These inner solitons are consistent with stability studies where it was shown to be an attractor solution under perturbations. The soliton core has a flat central density followed by a sharp drop in density (as steep as $r^{-16}$). Importantly, there is only a single parameter that characterises the cores: the core radius, or, equivalently, the core mass as the two are related by

[TABLE]

Once the core radius is chosen, the normalisation of the soliton core is determined by the balance between gravitational collapse and the quantum mechanical support: there is no freedom to choose the normalisation. More massive soliton cores are smaller in radius.

Interestingly, the outer profile of the simulated virialized BECDM haloes is found to follow a $r^{-3}$ power-law, similar to the outer profile of an NFW halo in CDM; the break occurs universally at about the soliton size $r_{\rm soliton}\simeq 3.5r_{\rm c}$. Fig. 2 shows the profiles of each of the $100$ haloes, along with the soliton core fits. After normalising the density profiles to their central density, we observe a unique soliton profile with some scatter beyond $r_{\rm soliton}$ attributed to the turbulent behavior. It is important to note that because we assume periodic boundary conditions and since mass is conserved, the system is continuously being perturbed by the reflecting waves that do not attenuate at the edges of the box, the latter precludes reaching the equilibrium configurations that are found analytically, where perturbations from the outer regions cease in a finite time and the systems may reach a dominant mode or configurations with multiple excited states (Matos & Ureña-López, 2007; Robles & Matos, 2013; Martinez-Medina, Robles & Matos, 2015; Bernal, Robles & Matos, 2017; Bernal et al., 2010; Ruffini & Bonazzola, 1969). Given our assumptions and in virtue of the stability arguments for BECDM haloes, it is then expected that the central ground state soliton will be the only mode that remains after a long evolution and the system has reached approximate virial equilibrium.

The radial energy densities of virialized haloes (potential, classical kinetic, quantum gradient) are plotted in Fig. 3 to highlight general properties. The soliton core is supported against gravitational collapse by the quantum pressure tensor, which is expected from the analytic description of steady-state soliton cores. The classical kinetic energy is sub-dominant in the core. The core is stable and protected against disruption from turbulent perturbations, which start appearing at $r\sim 3.5r_{\rm c}$, just where the quantum pressure and kinetic energies become comparable. Equipartition in the three energy densities is seen, however, in the outer parts of the profile, in the $\rho\propto r^{-3}$ region (i.e., the three energy densities follow the same radial profile up to a constant factor); for larger radii the potential energy becomes subdominant (due to boundary effects). Equipartition can be a characteristic feature of turbulence, as is the case here for this dissipationless fluid. The breakdown of the energy densities is found to be universal across all our simulations.

We note that in the region of the stable soliton core, the wave function actually has a time dependent phase (Guzmán & Ureña-López, 2004):

[TABLE]

where $\hat{\gamma}=-0.69223$, and

[TABLE]

This corresponds to a characteristic period of:

[TABLE]

In our simulations we do see oscillations of the core mass/radius of order a few percent with this characteristic frequency (also seen as oscillations in the energy components in Fig. 1). This may result from the intrinsic phase of the stable soliton being slightly perturbed by interfering constructively and destructively with the turbulent/chaotic uncorrelated phases of the surrounding turbulent medium (see Section 6). We note that the soliton core stays smooth and free of substructure at all times in this oscillation.

5 Soliton Core Mass Scaling

We study how the resulting soliton core mass scales with the other fundamental parameters of the system, namely the total mass $M$ and energy $E$. More precisely, $M$ and $E$ here refer to the total mass and energy in the 1 Mpc box, as they are scale invariants in the SP system. Schive, Chiueh & Broadhurst (2014); Schive et al. (2014) claim a relation $M_{\rm c}\propto(\lvert E\rvert/M)^{1/2}$. This is an interesting relation, because the left-hand side scales as the inverse of the soliton radius $r_{\rm c}^{-1}$ and the right hand side is proportional to the halo velocity dispersion $\sigma_{\rm h}$, which gives the relation $r_{\rm c}\sigma_{\rm h}\sim 1$, a nontrivial type of non-local uncertainty principle. However, Schwabe, Niemeyer & Engels (2016) point out that this relation may not be fundamental to the system, but may in fact be biased, due to the scaling symmetry of the fluid. Both sides of $M_{\rm c}\propto(\lvert E\rvert/M)^{1/2}$ scale as $\lambda$, the scaling parameter, so a linear relation would simply be found between these two parameters due to the scaling parameter alone, with no fundamental origin. Schwabe, Niemeyer & Engels (2016) recommended to look for fundamental relations by looking for relations between scale-free invariants, such as $\lvert E\rvert/M^{3}$. The authors do not find a single scaling that fitted their sample of mostly two-body collision simulations, possibly due to the small range of energies sampled by their suite and the strong dependence on the mass ratio of the two merging solitons.

In our 100 simulations of virialized multi-body mergers, essentially characterised by a single parameter $\Xi\equiv\lvert E\rvert/M^{3}/(Gm/\hbar)^{2}$ set by the initial mass and energy (we have assumed no net angular momentum), we do find a fundamental relation between core mass $M_{\rm c}$ and $\Xi$.

[TABLE]

which reproduces our simulations spanning two orders of magnitude in $E$, as shown in Fig. 4. More precisely, a numerical fit to the data yields $M_{\rm c}/M\simeq 2.89^{\pm 0.26}\Xi^{0.346^{\pm 0.013}}$, or $M_{\rm c}/M\simeq 2.63^{\pm 0.04}\Xi^{1/3}$ if the slope is fixed to $1/3$, where $1$-$\sigma$ errors are reported. Schwabe, Niemeyer & Engels (2016) do not find this result in their mainly two-body merger simulations, because the two-body results are sensitive to the mass ratio and the total angular momentum. More importantly, Schwabe, Niemeyer & Engels (2016) used sponge boundary conditions, in which a quasi steady-state solution is not reached; rather, mass and kinetic energy escape at the boundaries, which intrinsically change the total mass $M$ in Eqn. (32). These boundary conditions are closer to those used in analytical solutions of isolated haloes, while our simulations have closer resemblance to a cosmological scenario. Different boundary conditions and varying angular momentum in the simulations could potentially be a source of the differences in the results, but such detailed comparison is outside the scope of this work.

The relation $M_{\rm c}/M\propto\Xi^{1/3}$ simplifies to $M_{\rm c}\propto\lvert E\rvert^{1/3}$ implying that the soliton core traces the total energy of the system. Equivalently, the relation may be written also as: $r_{\rm c}\propto\lvert E\rvert^{-1/3}$ or $\lvert E_{\rm core}\rvert\propto\lvert E\rvert$. Essentially, the relationship tells us there is tight coupling between core and global properties as seen in Fig. 4. These relations suggest that global quantities of virialized haloes such as the total mass or energy are enough to estimate the expected core size, this will be important in cosmological simulations of BECDM where structures appear at several halo masses.

We could also get an estimated relation for the fundamental parameters in the following way: suppose from potential theory, that the gravitational potential at the center of the soliton is proportional to the gravitational potential at the center of the halo. We derive what this means for the scaling relation of the form $M_{\rm c}/M\propto\Xi^{\eta}$. For the soliton core, $M_{\rm c}r_{\rm c}\sim\hbar^{2}/(Gm^{2})$, so the velocity dispersion is: $v_{\rm c}^{2}\sim GM_{\rm c}/r_{\rm c}\sim G^{2}M_{\rm c}^{2}m^{2}/\hbar^{2}$. For the halo, $v_{\rm h}^{2}\sim GM/R_{\rm h}$ and $\lvert E\rvert\sim GM^{2}/R_{\rm h}$. Assuming a constant ratio $v_{\rm c}/v_{\rm h}\equiv\mu$, one may deduce:

[TABLE]

The power $\eta=1/2$ signifies a constant core to halo velocity dispersion ratio. In our simulations, however, we observe $\eta=1/3$, meaning that there is weak mass, energy dependence in the velocity dispersion ratio, namely $\mu\propto\Xi^{-1/6}$.

We note that a power of $\eta=1/2$ is derived analytically from self-similarity of the potential of the core and halo, which basically is the assumption: $M_{\rm c}/r_{\rm c}\sim M/R_{\rm h}$. On the other hand, $\eta=1/3$ can be derived analytically from self-similarity of the energy of the core and halo: $M_{\rm c}^{2}/r_{\rm c}\sim M^{2}/R_{\rm h}$, which may be a better assumption for the strongly coupled turbulent system found in our simulation. This assumption leads to the observed weak mass, energy dependence of $\mu$.

It has also been suggested, in the literature, that $M_{\rm c}/M$ essentially follows from the mass loss in subsequent binary mergers (Du et al., 2017). That work suggests $M_{\rm c}\propto M^{1.44(\beta-1)}$ where $\beta\leq 1$ is the descendant-to-originator mass-fraction. Assuming the halo radius scales approximately as $R_{\rm h}\propto M^{1/3}$, then Eqn. 32 implies $M_{\rm c}\propto M^{5/9}$, corresponding to $\beta=0.69$. This is consistent with the result $\beta\sim 0.7$ of Schwabe, Niemeyer & Engels (2016); Du et al. (2017).

6 Velocity Power Spectra

We compute the 1D radial superfluid energy spectrum $E_{v^{2}}(k)$ defined by Baggaley, Laurie & Barenghi (2012):

[TABLE]

BEC systems without self-gravity but with the non-linear self-interaction term (Gross-Pitaevskii equations) are known to reproduce a classical Kolmogorov scaling $E_{v^{2}}(k)\propto k^{-5/3}$ if mechanically driven on the box scale, and a shallower $E_{v^{2}}(k)\propto k^{-1}$ turbulent cascade in the counterflow regime for large $k$ beyond the inter-vortex length-scale ‘bump’ and no power on the largest scales (Baggaley et al., 2012).

Fig. 5 shows the energy power spectrum $E_{v^{2}}(k)$ for our simulations. The simulations themselves show sustained chaotic motions (stable kinetic energy with time, equipartition) and a homogeneous filamentary distribution of $v$ (which traces out the vortex lines in the fluid; Fig. 6). No turbulence appears inside the soliton. There is no power on the largest spatial scales, as expected, due to lack of large-scale driving, limited by the simulation box size. The simulations all show a well converged power-law relation of $E_{v^{2}}(k)\propto k^{-1.1}$, closer to the thermally-driven counterflow analog Gross-Pitaevskii system from condensed matter physics than the mechanically driven one, and we also observed a maximum mode that carries most of the energy in the turbulent medium. This energy power spectrum is characteristic of isotropic turbulence where small modes dominate the turbulence. As seen in Fig. 5, we find that with our highest resolution we can capture the scale where the spectrum peaks, very low resolution could lead to missing the peak mode $k_{\rm peak}$ and result in a lack of homogeneous turbulence in the simulation.

In the bottom panel of Fig. 5, we show the relation between the scale $d_{\rm peak}=2\pi/k_{\rm peak}$ where $E_{v^{2}}(k)$ peaks for each halo and the core size of the corresponding soliton. We find $d_{\rm peak}\sim 7.5r_{\rm c}$, this corresponds to the region where the potential energy density contributes equally to the kinetic energy density (Fig. 3). Notice that this wavelength $d_{\rm peak}\sim 2r_{\rm soliton}$, which is the total width of the soliton and it is the typical scale where most of the interference is happening, since $r_{\rm soliton}$ is much smaller than the box, this explains the appearance of the homogeneous turbulence throughout the box.

Evidence for the existence of turbulence everywhere in the domain comes from the identification of vortex lines in our simulations (Fig. 6). These filamentary structures are a source of turbulence in a quantum fluid. Vortex lines are degenerate locations in the fluid that have a discontinuity in the differentiability of $\psi$ and have $\lvert\psi\rvert\to 0$. They contain all the vorticity in the fluid as the velocity field must be curl-free elsewhere. Turbulence persist at all times in the box since the system is closed and the total angular momentum ($\textbf{L}=\textbf{0}$) is conserved. The existence of vortex lines is a necessary condition for quantum turbulence, and their reconnection creates Kelvin waves that drive the turbulent motions. Fig. 6 shows a slice of the magnitude and phase of $\lvert\psi\rvert$. We see clear evidence of $\lvert\psi\rvert\to 0$ filamentary structures which correspond with discontinuities in the phase of $\psi$. The figure also shows a zoom-in on the network of reconnection events homogeneous throughout the domain (except inside the soliton code). The reconnection events (hence reconnection rate) are numerically converged, as the figure shows the same events are identified for simulation resolutions $512^{3}$ and $256^{3}$ at the Hubble time. The excellent convergence is due to the exponential spatial convergence properties of the spectral method.

The turbulent structure is seen everywhere in the domain, except, of course, inside the soliton core, which is protected from velocity fluctuations because it is a stable soliton solution. Both the soliton width and the turbulent structures (peak of the velocity power spectrum) are $\approx 2r_{\rm soliton}$, determined by the de Broglie length-scale of the system, the only length-scale in the problem.

7 Concluding Remarks

We have carried out $100$ simulations of virialized BECDM halo cores with periodic boundary conditions to characterise their properties, which are largely universal and independent of initial condition details. Merging multiple haloes with initially cuspy or cored profiles both lead to the formation of stable soliton cores at the centers of BECDM haloes. Our simulation setup provides a useful numerical laboratory for statistically studying the final product of the relaxation process of BECDM haloes with a wide range of total energies. The structure of the resulting dark matter haloes depends primarily on the total mass and energy of the system. The haloes form a stable soliton core with a turbulent $r^{-3}$ NFW-like outer profile and (in the absence of angular momentum) are characterised by a single dimensionless invariant: $\Xi\equiv\lvert E\rvert/M^{3}/(Gm/\hbar)^{2}$. Contrary to previous works, we find that for all of our haloes the core mass of the soliton scales with this quantity $\Xi$ as $M_{\rm c}/M\propto\Xi^{1/3}$, which implies $M_{\rm c}\propto\lvert E\rvert^{1/3}\Rightarrow|E_{\rm core}\rvert\propto\lvert E\rvert$. Properties of the soliton at the centers of haloes are therefore tightly linked to the global halo properties.

Soliton core profiles are described by a single parameter (the core mass, or equivalently the core radius, as the two are related by $M_{\rm c}\propto r_{\rm c}^{-1}$). The size/mass of the core that forms traces the total energy of the entire virialized dark matter halo. We found that the typical soliton size is $~{}3.5r_{\rm c}$, beyond this radius in the BECDM profile the haloes are found to be turbulent, exhibiting a filamentary distribution of vortex lines that form during merging events and are sustained and reconnect to drive turbulence. No turbulence is seen inside the soliton. Equipartition between the potential energy density, classical kinetic energy density, and quantum gradient energy density is seen in the outer core, maintaining a continuously perturbed medium. The turbulence is characterised by a $k^{-1.1}$ power-law in the velocity power spectrum, characteristic of an isotropic turbulence in the fluid. We find this is because the dominant mode in the 1D superfluid velocity spectrum peaks at a scale twice the soliton radius, which is several times smaller than the total length of the system. We find that the suppression of turbulence inside the soliton and the existence of a maximum mode in the velocity power spectra with a scale equal to the soliton width, could explain why the typical scale for the granules in the density field of BECDM simulations is preferentially the soliton size.

The cuspy halo profile universally found in $\Lambda$CDM simulations (Navarro, Frenk & White, 1996) has seen tension with observations that suggest core-like potentials at halo centers. A recent example of such an observation is of SPT-CLJ2011-5228 (Collett et al., 2017), a $z=2.39$ system gravitationally lensed by a $z=1.06$ cluster along the line-of-sight. The inner density profile falls with radius to the power $-0.38\pm 0.04$ ($1\sigma$ errors) out to a radius of $270^{+48}_{-76}$ kpc. This shallow inner profile is in strong tension with our understanding of relaxed cold dark matter haloes, where NFW predicts a $r^{-1}$ profile, and perhaps the flat slope is suggestive of a central soliton.

An interesting application of the systems studied here may be in neutron star glitch statistics, as neutron star interiors may be modelled as a superfluid by the Gross-Pitaevskii equations (Warszawski & Melatos, 2011, 2013). Glitches may originate from the turbulent nature of the fluid, along with the possible intermittent nature of turbulence.

BECDM has been largely studied in the past and a number of independent constraints exist on the boson mass that would make up this type of dark matter, as outlined in the introduction. However, the analytical arguments ultimately need to be validated by numerical simulations to ensure that the analytic assumptions made in the derivations are valid, as was the case historically for the standard CDM scenario. Some of the analyses have suggested moderate tension in the boson mass; however, only through full BECDM cosmological simulations (ultimately involving full baryonic physics) we will be able to confirm these assumptions and place tighter constraints on the boson mass. For example, Lyman-$\alpha$ constraints have not included quantum density fluctuations, which are present in the BECDM model and seen in our simulations

7.1 Context and outlook

The next step in our work will be to simulate the BECDM model coupled to baryons in a fully cosmological setting (Mocz et al. in prep.) to address the impact of baryons on the predictions of BECDM only simulations. We have implemented the numerical method presented in this paper into the arepo code, so we will be able to run the axion dark matter simulations fully-coupled to baryonic components with feedback in upcoming work.

Fig. 7 shows an example of a BECDM $1~{}{\rm Mpc}$ cosmological simulation ran with the arepo code, at redshift $z=4$, with resolution of $1024^{3}$ cells for a boson mass of $m=2.5\times 10^{-22}~{}{\rm eV}$. The resolution of the simulation is $\sim 1$ kpc, enough to capture turbulence and the soliton cores. We have highlighted in the figure a virialized halo as a result of cosmological mergers. Turbulence is also seen in the cosmological box at the intersections of cosmic web filaments where halos have merged, and our current simulations are able to provide a high resolution characterisation of the phenomenon (i.e., the turbulent cascade is resolved over a broader range). The BECDM cosmological simulations, to be described in detail in Paper II, are created with a realistic axion power spectrum at $z\sim 100$ and we aim at using them to explore the effect of varying the boson mass, and to study the coupling with baryonic physics.

Acknowledgments

This material is based upon work supported by the National Science Foundation (NSF) Graduate Research Fellowship under grant no. DGE-1144152 (PM). PM is supported in part by the NASA Earth and Space Science Fellowship (NNX12AB23C). VHR is supported by a UC MEXUS-CONACYT fellowship. PM would like to thank Jerry Ostriker, Alma Gonzalez, Luis Urena-Lopez, and Mikhail Medvedev for valuable discussions, and Sauro Succi for reading of an earlier version of this manuscript. The computations in this paper were run on the Odyssey cluster supported by the FAS Division of Science, Research Computing Group at Harvard University. MV acknowledges support through an MIT RSC award, the support of the Alfred P. Sloan Foundation, and support by NASA ATP grant NNX17AG29G. MBK acknowledges support from the NSF (grant AST-1517226) and from NASA through ATP grant NNX17AG29G and HST theory grants (programs AR-12836, AR-13888, AR-13896, AR-14282, and AR-14554) awarded by the Space Telescope Science Institute (STScI), which is operated by the Association of Universities for Research in Astronomy (AURA), Inc., under NASA contract NAS5-26555. JZ acknowledges support by a Grant of Excellence from the Icelandic Research Fund (grant number 173929-051).