SCALAR: an AMR code to simulate axion-like dark matter models

TL;DR

SCALAR is a new adaptive mesh code based on RAMSES that accurately simulates axion-like dark matter models by solving the Schr"odinger-Poisson equations, enabling detailed studies of fuzzy dark matter.

Contribution

It introduces a novel Schr"odinger equation solver integrated with RAMSES for simulating axion-like dark matter models with adaptive refinement.

Findings

Demonstrates accurate numerical implementation through tests

Enables simulation of non-linear Schr"odinger equations with external potentials

Provides a versatile tool for dark matter research

Abstract

We present a new code, SCALAR, based on the high-resolution hydrodynamics and N-body code RAMSES, to solve the Schr\"odinger equation on adaptive refined meshes. The code is intended to be used to simulate axion or fuzzy dark matter models where the evolution of the dark matter component is determined by a coupled Schr\"odinger-Poisson equation, but it can also be used as a standalone solver for both linear and non-linear Schr\"odinger equations with any given external potential. This paper describes the numerical implementation of our solver and presents tests to demonstrate how accurately it operates.

Equations

$$i\hslash \frac{\partial \psi }{\partial t}=\left[-\frac{{\hslash }^2}{2m}{\nabla }^2+g|\psi {|}^2+m{V}_{\mathrm{ext}}\right]\psi ,$$

$$\int |\psi {|}^2\text{leavevmode}{\mathrm{d}}^{3}x=N,$$

$$\psi =\sqrt{\frac{\rho }{m}}\exp \left(\frac{i}{\hslash }\theta \right),$$

$$\rho (\mathbf{x},t)$$

$$\text{@vec}v(\mathbf{x},t)$$

$$\frac{\partial \rho }{\partial t}+\nabla \cdot (\rho \text{@vec}v)=0,$$

$$\frac{\partial \mathbf{v}}{\partial t}+(\text{@vec}v\cdot \nabla )\text{@vec}v=-\nabla \left(V_{\mathrm{ext}}+\frac{g}{m^2}\rho +Q\right),$$

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

$$\text{@vec}\nabla \times \text{@vec}v=\text{@vec}\nabla \times \text{@vec}\nabla \theta =0.$$

$$i\hslash \left(\frac{\partial \psi }{\partial t}+\frac{3}{2}H\psi \right)=\left(-\frac{{\hslash }^2}{2m{a}^2}{\nabla }^2+m_a\Phi \right)\psi ,$$

$${\nabla }^2\Phi =4\pi G{a}^2\left(|\psi {|}^2-|\psi (a){|}^2\right),$$

$$\psi (\mathbf{x},t_1)=\hat{U}(t_1,t_0)\psi (\mathbf{x},t_0),$$

$$\hat{U}(t_1,t_0)=\exp \left(-\frac{i}{\hslash }{\int }_{t_0}^{t_1}\hat{H}(\mathbf{x},t^{\prime })\text{leavevmode}\mathrm{d}{t}^{\prime }\right),$$

$${\int }_{t_0}^{t_1}\hat{H}(\mathbf{x},t^{\prime })\text{leavevmode}\mathrm{d}{t}^{\prime }\approx \hat{H}(\mathbf{x},t_0)\text{leavevmode}\Delta t,$$

$$\hat{U}(t_1,t_0)\approx \exp \left(-\frac{i}{\hslash }\hat{H}(\mathbf{x},t_0)\Delta t\right).$$

$$\hat{K}(\mathbf{x},t)\equiv -\frac{{\hslash }^2}{2m}{\nabla }^2,$$

$$\hat{W}(\mathbf{x},t)\equiv m\left(V(\mathbf{x},t)+\frac{g}{m}|\psi (\mathbf{x},t){|}^2\right).$$

$$\hat{U}(t_1,t_0)\approx \exp \left(-\frac{i}{\hslash }\hat{W}(\mathbf{x},t_0)\Delta t\right)\exp \left(-\frac{i}{\hslash }\hat{K}(\mathbf{x},t_0)\Delta t\right).$$

$$\psi (\mathbf{x},t_1)=$$

$${\bar{\psi }}_{i,j,k}^{n+1}$$

$$=\left[\sum _{N=0}^{\infty }\frac{1}{N!}{\left(-\frac{i}{\hslash }\hat{K}(\mathbf{x},t_0)\Delta t\right)}^N\right]{\psi }_{i,j,k}^n$$

$$=\left[1+\left(\frac{i\hslash \Delta t}{2m}{\nabla }^2\right)+\frac{1}{2}{\left(\frac{i\hslash \Delta t}{2m}{\nabla }^2\right)}^2+\dots \right]{\psi }_{i,j,k}^n,$$

$${\nabla }^2{\psi }_{i,j,k}^n=$$

$$\frac{{\psi }_{i,j+1,k}^n+{\psi }_{i,j-1,k}^n-2{\psi }_{i,j,k}^n}{\Delta x^2}+$$

$$\frac{{\psi }_{i,j,k+1}^n+{\psi }_{i,j,k-1}^n-2{\psi }_{i,j,k}^n}{\Delta x^2}.$$

$${\psi }_{i,j,k}^{n+1}=\exp \left(-\frac{i}{\hslash }{\hat{W}}_{i,j,k}^n\Delta t\right){\bar{\psi }}_{i,j,k}.$$

$$\displaystyle\rho_{i,j,k}^{n+1}\big{|}_{S}=m\left|\psi_{i,j,k}^{n+1}\right|^{2}.$$

$$\frac{\partial \rho }{\partial t}+\nabla \cdot \mathbf{j}=0,$$

$$\mathbf{j}\equiv -i\frac{\hslash }{2m}\left({\psi }^{\ast }\nabla \psi -\psi \nabla {\psi }^{\ast }\right).$$

$$\mathbf{j}=\frac{\hslash }{m}\left(\Re [\psi ]\nabla \Im [\psi ]-\Im [\psi ]\nabla \Re [\psi ]\right).$$

$$\frac{{\rho }_{i,j,k}^{n+1}-{\rho }_{i,j,k}^n}{\Delta t}$$

$$+\frac{\left(j_{i,j+\frac{1}{2},k}^{n+\frac{1}{2}}-j_{i,j-\frac{1}{2},k}^{n+\frac{1}{2}}\right)}{\Delta x}$$

$$+\frac{\left(j_{i,j,k+\frac{1}{2}}^{n+\frac{1}{2}}-j_{i,j,k-\frac{1}{2}}^{n+\frac{1}{2}}\right)}{\Delta x}=0,$$

$${\psi }_{i+\frac{1}{2},j,k}^{n+\frac{1}{2}}$$

$${\psi }_{i,j+\frac{1}{2},k}^{n+\frac{1}{2}}$$

$${\psi }_{i,j,k+\frac{1}{2}}^{n+\frac{1}{2}}$$

Full text

1 Introduction

The true nature of dark matter is not known. Weakly interacting massive particles (WIMP) are still considered one of the most likely candidates for cold dark matter (CDM), and several experiments are ongoing to try to detect such particles. These are, however, closing in on the neutrino floor where any signal would be drowned in the solar neutrino background (Monroe & Fisher 2007).

A promising alternative to WIMPs are ultra-light axions, fuzzy dark matter (Marsh 2016; Hui et al. 2017; Turner 1983; Press et al. 1990; Goodman 2000; Peebles 2000; Calabrese & Spergel 2016; Guzmán & Ureña-López 2003; Kim & Marsh 2016; Amendola & Barbieri 2006; Hu et al. 2000), and superfluid dark matter (Berezhiani & Khoury 2015; Khoury 2016, 2015; Sharma et al. 2019). These models have distinct and observable signatures on the small scales of structure formation, they are also are able to solve some of the discrepancies observed in CDM simulations, such as the missing satellites problem (Moore et al. 1999; Klypin et al. 1999), the cusp-core problem (de Blok 2010; Bar et al. 2018, 2019) and the too-big-to-fail problem (Boylan-Kolchin et al. 2011). These disparities could, however, have a solution within baryonic processes, which are usually not included in standard CDM simulations as shown in Macciò et al. (2012), Brooks et al. (2013), Governato et al. (2010), Madau et al. (2014), Teyssier et al. (2013), Governato et al. (2012), Pontzen & Governato (2012), Sawala et al. (2016), and Brooks et al. (2017). To understand this better, one should ideally perform simulations including both of these components.

In order to quantify the effects of axion-like dark matter models, one needs to either solve a Schrödinger-Poisson system or use the Mandelung formulation. The latter consists of a set of traditional hydrodynamics equations, with an additional pressure term which can be solved by using methods such as smoothed particle hydrodynamics (SPH), as proposed by Marsh (2015). However, it is known that numerical methods based on the Madelung formulation of quantum mechanics are troublesome in regions around voids. This formulation, indeed, breaks down where the density approaches zero and at interference nodes, as the quantum pressure term can easily become singular (Marsh 2015; Li et al. 2019).

The Schrödinger-Poisson system has several applications in cosmology. For instance, the six dimensional Vlasov equation describing collisionless self-gravitating matter is approximated by a Schrödinger-Poisson system for a complex wave-function in three dimensions. This was proposed as an alternative way for simulating CDM in Widrow & Kaiser (1993). It was later shown, by solving the Schrödinger-Poisson system and comparing it to the Vlasov solver ColDICE (Sousbie & Colombi 2016) in two dimensions, that one has excellent qualitative and quantitative agreement in the solution (Kopp et al. 2017). A similar study is given in Mocz et al. (2018), where the system was solved using a spectral method, demonstrating that one recovers the classical behaviour in the limit where $\hbar\to 0$.

Unfortunately, the methods employed in the above mentioned papers, despite being very accurate, are too expensive to perform high-resolution simulations in three dimensions. The first cosmological, high-resolution, simulation of fuzzy dark matter in three dimensions was performed in Schive et al. (2014) using the code GAMER (Schive et al. 2010, 2018). There, an explicit method, similar to the one we present in this paper, was used. In Veltmaat et al. (2018), they used the classical wave-function to perform zoom-in simulations to study the formation and evolution of ultralight bosonic dark matter halos from realistic cosmological initial conditions and in Mocz et al. (2017) they studied galaxy formation with Bose-Einstein condensate dark matter using a pseudo spectral method (see also Edwards et al. (2018)). There have also been a handful of papers that have performed simulations more in line with the hydrodynamical formulation. In Zhang et al. (2018), a new technique to discretise the quantum pressure is proposed and shown to reproduce the expected density profile of dark matter halos. In Nori & Baldi (2018), a module AX-GADGET for cosmological simulations using SPH inside the P-GADGET3 code is presented. These methods do not solve for the wave-function, but they have the advantage of being much less expensive to run than a full wave-function solver like ours. There have also been simulations performed by using other numerical techniques in either two (Kopp et al. 2017) or three spatial dimensions (Woo & Chiueh 2009).

In this paper we present SCALAR (Simulation Code for ultrA Light Axions in RAMSES): a general adaptive mesh refinement (AMR) code to solve the Schrödinger-Poisson system. Our implementation is in the hydrodynamics and N-body code RAMSES (Teyssier 2002). The structure of the paper is as follows: in Section 2 we present the equations we are to solve, in Section 3 we present the numerical implementation, in Section 4 we present tests of the code and in Section 5 we discuss possible cosmological applications before concluding in Section 6.

2 Theoretical model

A Bose-Einstein condensate (BEC) is a system of identical bosons, where a large fraction of particles occupies the lowest quantum energy state, or the ground state. This phenomenon typically takes place in gases, at very low temperatures or very high densities and it was observed for the first time in Anderson et al. (1995) and Davis et al. (1995). In the condensate regime, these quantum systems behave as a macroscopic fluid and their peculiar features are a macroscopic manifestation of quantum effects.

In general, when Bose-Einstein condensation occurs, thermal de-Broglie wavelengths of particles start to overlap, as they become grater than the mean inter-particle distance. At this point, a coherent state develops and the system behaves as a macroscopic fluid, where only binary collisions are relevant. The dynamics of BECs is complicated, due to the difficulty in modelling particle self-interactions.

However, in the Hartree mean-field theory and in the limit of $T\to 0$, binary collisions are modelled via an effective potential and the whole quantum system can be described by a single-particle wave-function $\psi(\mathbf{x},t)$ obeying the non-linear Schödinger equation:

[TABLE]

where $m$ is the mass of the boson and $g$ is the self-interaction coupling constant. Often, the trapping potential $V_{\rm ext}(\mathbf{x},t)$ is introduced by hand in order to model the presence of a trap, which is responsible for keeping particles confined.

The single-particle wave-function is normalised such that:

[TABLE]

where $N$ is the total number of particles present in the system. As a consequence, the quantity $|\psi(\mathbf{x},t)|^{2}$ represents the number density of particles.

An alternative description of the macroscopic fluid is provided by the so-called Madelung formulation of the Schrödinger equation. In this case, by expressing the single-particle wave-function in polar coordinates:

[TABLE]

the dynamics of the system is described in terms of mass density and velocity, which are macroscopic physical quantities and they are respectively defined as:

[TABLE]

Thus, the Schrödinger equation can be cast into the following system of equations:

[TABLE]

which are known as the Madelung or quantum Euler equations. We recognise Eq. (6) as a continuity equation which expresses conservation of mass. Although the second Madelung equation, Eq. (7), expresses conservation of momentum, it is not the same as the classical momentum equation, as it contains an additional term $Q$, which is called quantum pressure and it is defined as:

[TABLE]

The quantum pressure is a macroscopic manifestation of quantum effects and it is characteristic of Bose-Einstein condensates.

In this formulation, by defining the velocity as in Eq. (5), we are intrinsically assuming that the fluid is irrotational, since:

[TABLE]

However, during the evolution of the wave-function, the phase can develop discontinuities of multiples of $2\pi\hbar$ and its gradient can subsequently generate vorticity in the field, as shown in Uhlemann et al. (2014, 2019).

In cosmology, these kinds of systems can be used to model the dark matter contribution to the energy budget of the Universe. In particular, in the last few decades, models where dark matter is a light boson, such as ultra-light axions or fuzzy dark matter, have received a lot of attention. Due to the small mass of these bosons, macroscopic quantum effects manifest at astronomically relevant scales. In these alternative dark matter models, new signatures are expected within the structure formation process at highly non-linear scales and, therefore, numerical simulations are required in order to explore these scenarios.

Here, the dynamics of dark matter is also described by a system of identical bosons gravitationally bounded. Therefore, the governing equation is a non-linear Schrödinger equation, Eq. (1), where the external potential is replaced by the gravitational potential. In this class of alternative dark matter models, self-interactions between bosons are often neglected, as the coupling constant $g$ is usually parametrically small. The resulting system of equations describing the dynamics of the dark matter fluid is called Schrödinger-Poisson system and, for an expanding Universe, it reads:

[TABLE]

where $a$ is the scale-factor of the Universe, $H\equiv d\log(a)/dt$ is the Hubble rate of expansion, and $\Phi$ is the gravitational potential. With a change of variables $\psi\to a^{3/2}\psi$, the non-linear Schrödinger equation above takes on the form of Eq. (1).

3 Numerical methods

In this section we provide a brief overview of the code RAMSES and the AMR technique. Then, we discuss in details the numerical aspects of the algorithm we implemented in order to solve the non-linear Schrödinger equation. Throughout this section, the dimensionality of the problem is denoted by $\rm{dim}$ and it can be 1, 2 or 3.

3.1 Overview of RAMSES

The RAMSES code was originally designed for cosmological simulations of structure formation and subsequently extended to astrophysical applications. It consists of an N-body particle mesh (PM) code, which solves the gravitational dynamics of a set of macroparticles, sampling the phase space distribution of the dark matter component in the Universe. Through PM algorithms, the mass of each macroparticle is interpolated on a grid and the Poisson equation is solved in order to compute the gravitational potential. Thus, the gravitational force acting on each element of the system and the new phase space position of each macroparticle are computed by solving the corresponding N-body equation with a leapfrog scheme. In addition, RAMSES can solve the dynamics of the baryonic component present in the Universe. In this case, the grid is also used to sample gas parcels and the evolution of the system is described by the equations of hydrodynamics, which are solved by means of a Godunov scheme. For this purpose, Riemann solvers can be used for computing fluxes of conserved physical quantities among cells.

The RAMSES code implements an AMR strategy, where a hierarchy of nested grids is created in order to increase the local resolution according to a set of refinement criteria. In this way, RAMSES can solve accurately gas dynamics and gravitational potential only where more precision is actually needed. This approach reduces consistently the amount of memory needed in cosmological and hydrodynamical simulations, compared to the case where a uniform high-resolution grid is used.

In SCALAR, we rely on the efficient AMR implementation of RAMSES. In order to solve the dynamics of our theoretical model, the single-particle wave-function is sampled by using the original grid allocated by RAMSES for the Poisson and hydrodynamics equations. Also in this case, the AMR approach provides the opportunity to solve the Schrödinger equation with higher resolution only where features of the wave-function are more demanding.

3.2 Adaptive mesh refinement

The basic unit of the AMR hierarchy is an oct, which consists of a set of $2^{\rm dim}$ cells. At each level in the AMR hierarchy, a grid is a collection of octs with the same resolution. The grid with the coarsest resolution is called domain grid and it covers the whole computational domain. During the evolution of the physical system, when the solution starts to develop features and its tracking requires higher resolution, any cell at a given level can be split into a child oct, with double the resolution of the parent cell.

At each time-step, the AMR structure is modified according to a set of refinement criteria. First, for a generic level of refinement $\ell$, a refinement map is created by marking all those cells satisfying at least one refinement criterion. Also cells violating the strict refinement rule are marked for refinement, in order to guarantee that each child oct at level $\ell+1$ is surrounded by, at least, $3^{\rm dim}-1$ neighbours at the coarser level. However, if a given cell at level $\ell$ does not satisfy any refinement criteria anymore, it is marked for de-refinement and subsequently its child octs are destroyed. Then, a new child oct is created at level $\ell+1$ for each marked cell and all the relevant physical quantities are interpolated from level $\ell$. Coarse-fine data interpolation, in an AMR context, is often called prolongation and it can be done by using any of the interpolation schemes which are described in the section below.

When computing the refinement map, physical quantities can fluctuate around the refinement threshold in subsequent time-steps. This means that, some cells at the boundary of a fine resolution region can be refined and de-refined many subsequent times. In this case, the refinement map tends to be quite noisy, since each interpolation operation where the solution is characterised by strong gradients introduces numerical noise in the solution of the non-linear Schrödinger equation. For this reason, once the refinement map is computed according to a given set of refinement criteria, a mesh smoothing operator is applied. For this purpose, a cubic buffer of $n_{\rm expand}$ cells around the computed map is additionally marked for refinement. In this way, even if cells octs are created and destroyed at coarse-fine boundaries, the interpolation procedure is applied in a regions where the wave-function is smoother and, thus, it introduces a lower level of numerical noise.

3.3 The Schrödinger equation

The SCALAR code evolves the solution of the non-linear Schrödinger equation by using a Taylor’s method, similar to the one designed in GAMER. Given the wave-function $\psi(\mathbf{x},t_{0})$, the formal solution of the non-linear Schrödinger equation, Eq. (1), at time $t_{1}=t_{0}+\Delta t$ reads:

[TABLE]

where $\hat{U}(t_{1},t_{0})$ is the time evolution operator and it maps the solution of the Schrödinger equation at two different times. In the general case, the time evolution operator is defined as:

[TABLE]

where $\hat{H}(\mathbf{x},t)$ denotes the Hamiltonian of the system. The operator $\hat{U}(t_{1},t_{0})$ has the main following properties:

• •

$\begin{aligned} \hat{U}(t,t)=1,\end{aligned}$

• •

$\begin{aligned} \hat{U}(t_{1},t_{2})\leavevmode\nobreak\ \hat{U}(t_{2},t_{3})=\hat{U}(t_{1},t_{3}),\end{aligned}$

• •

$\begin{aligned} \hat{U}(t_{1},t_{2})=\hat{U}^{{\dagger}}(t_{2},t_{1})=\hat{U}^{-1}(t_{2},t_{1}).\end{aligned}$

In the limit of $\Delta t\ll 1$, then the following approximation holds:

[TABLE]

and, therefore, Eq. (13) can be approximated as:

[TABLE]

In the general case, the Hamiltonian $\hat{H}(\mathbf{x},t)$ contains different contributions to the total energy of the system. In particular, we can express $\hat{H}(\mathbf{x},t)$ as a sum of contributions describing kinetic and potential energies. Here, we denote these two operators respectively as $\hat{K}(\mathbf{x},t)$ and $\hat{W}(\mathbf{x},t)$, and they are defined as:

[TABLE]

By means of the Lie-Trotter formula (Trotter 1959), the time evolution operator can be split as well:

[TABLE]

As a consequence, the formal solution of the Schrödinger equation can be written as:

[TABLE]

In SCALAR , the two contributions to the time evolution operator are applied separately. First, the 'drift' due to the kinetic part of the Hamiltonian is approximated via Taylor expansion (here for ${{\rm dim}}=3$):

[TABLE]

where, for a generic operator $\hat{O}$, the notation $\hat{O}^{N}$ denotes $N$ consecutive applications of the same operator. In the SCALAR code, the Taylor expansion is performed up to $\mathcal{O}(\Delta t^{3})$, which is the minimum required by the stability analysis of the numerical scheme. Furthermore, the laplacian operator is discretised by a standard second-order finite difference formula:

[TABLE]

Then, the 'kick' due to the potential is computed and the wave-function at the new time-step reads:

[TABLE]

Here, the advantage of the Lie-Trotter splitting is clear: while the kinetic contribution to the time evolution operator needs Taylor expansion in order to be applied, the potential contribution only provides a phase rotation of the wave-function and it can be computed exactly.

Once the new wave-function is computed, the new mass density is computed according to:

[TABLE]

3.4 The continuity equation

In quantum mechanics, the time evolution operator is unitary, as expressed by its properties. This means that the mass density carried by the wave-function is conserved. This is true also if we consider separately the two contributions to the Hamiltonian. However, the Taylor expansion, Eq. (3.3), breaks the unitarity of the time evolution operator. Therefore, in order to improve the conservation properties of our main numerical scheme, we implement a secondary solver for the continuity equation associated to the non-linear Schrödinger equation.

Eq. (6) can be written in its conservative form:

[TABLE]

where $\rho(\mathbf{x},t)$ represents the mass density. Here, the quantity $\mathbf{j}(\mathbf{x},t)$ is the associated density current, or flux, and it is defined as:

[TABLE]

By explicitly considering real and imaginary part of the wave-function, the density current can also be expresses as:

[TABLE]

In SCALAR, Eq. (24) is discretised by using a first-order Godunov scheme:

[TABLE]

where the time-centered fluxes are computed at cell interfaces. In order to compute the time-centered fluxes at cell interfaces, the wave-function is first computed at half time-step, by advancing the solution of $0.5\Delta t$. Then, the wave-function at cell interfaces is estimated in each dimension by linear interpolation:

[TABLE]

Its gradient, instead, is computed in each dimension by means of the first-order finite difference formula and it reads:

[TABLE]

Thus, the time-centered density current at cell interfaces is computed by means of Eq. (26).

However, this solver is not used to explicitly advance in time the mass density, but only to enforce the conservation of mass. Indeed, by denoting \rho_{i,j,k}^{n+1}\big{|}_{S} and \rho_{i,j,k}^{n+1}\big{|}_{C} the new mass densities computed by the main and the secondary solvers respectively, a correcting factor is computed as follows:

[TABLE]

which is used to rescale the wave-function.

Although by solving the continuity equation on top of the Schrödinger equation sensibly improves the conservation properties of the algorithm, this process does not ensure perfect conservation of mass. Since we truncate at the third order in time the Taylor expansion of the kinetic operator, Eq. (3.3), the kinetic solver introduces truncation errors in the solution of the Schrödinger equation. However, by rescaling the wave-function with the correct amplitude computed from the continuity equation, we are able to significantly reduce truncation errors, as shown by the improvement of several orders of magnitude on the conservation properties of SCALAR in the tests we present in Section 4. Unfortunately, the rescaling procedure is subject to accumulation of round-off errors, leading to an evolution of the error in the conservation of mass itself. However, the round-off errors are of the order of machine precision, in contrast to the much higher amplitude of the truncation errors introduced by the kinetic solver. This solution was already adopted in the GAMER code (Schive et al. 2014).

3.5 The solver

SCALAR solves the Schrödinger equation from the coarser to the finer level in the AMR hierarchy. For a generic refinement level $\ell$, the optimal time-step is chosen as:

[TABLE]

where $|V_{\rm max}|$ denotes the maximum absolute value of the effective potential $V+\frac{g}{m}|\psi|^{2}$. Here, $C_{K}$ and $C_{W}$ are Courant factors which are required to be smaller than one. The first term in the square brackets is determined by the Von Neumann stability analysis of the kinetic part of the solver. The second term, instead, requires that the phase of the wave-function does not rotate by a bigger angle than $2\pi C_{W}$ within a time-step. In general, the Courant factors $C_{K}$ and $C_{W}$ are chosen empirically, depending on the characteristics of the physical system that one aims to model. In our case, we set $C_{K}=C_{W}=0.2$ since it provides a good accuracy on the solution of the non-linear Schrödinger equation, without sacrificing too much computation time. In Appendix A we provide a detailed discussion of the Von Neumann stability analysis of the numerical scheme.

In the original RAMSES code two different options are available regarding the choice of the time-step: a single or an adaptive time-step. While the former consists in using the same time-step for all refinement levels and it is determined by the finest level in the AMR hierarchy, the latter allows to use smaller time-steps for finer refinement levels. However, in case of adaptive time-step, for each coarse time-step at level $\ell$ it is possible to perform only two fine time-steps at level $\ell+1$. In SCALAR, an additional option is available: a flexible time-step, where for each coarse time-step at level $\ell$, the number of fine steps at level $\ell+1$ is flexible and it is determined by a level dependent Courant-Friedrichs-Lewy (CFL) condition. From Eq. (31), when the optimal time-step is chosen by the kinetic CFL condition, the time-step scales with the grid size as $\Delta t\propto\Delta x^{2}$, which represents a stricter condition than the usual case of hydrodynamics equations. Therefore, a flexible time-step can reduce significantly the total amount of coarse time-step in a simulation.

Within a generic level of refinement $\ell$, SCALAR solves the non-linear Schrödinger equation for each oct separately. Thus, in order to advance the solution over a time-step, the solver proceeds as follows:

For a given oct, a cubic buffer of neighbours cells is collected. The equations are actually solved only for the central oct, while the buffer cells are used to compute laplacians at each order of the Taylor expansion. If the central oct lies next to the coarse-fine boundary, the wave-function is interpolated into ghost cells from level $\ell-1$. 2. 2.

The kinetic solver evolves the wave-function at the new time-step by means of Eq. (3.3). First, by advancing the solution by $0.5\Delta t$, the half time-step solution $\bar{\psi}_{i,j,k}^{n+1/2}$ is estimated from $\psi_{i,j,k}^{n}$ and it is used later in order to compute the mass density currents. Then, the full-step solution $\bar{\psi}_{i,j,k}^{n+1}$ and the new mass density \rho_{i,j,k}^{n+1}\big{|}_{S} are computed. 3. 3.

The wave-function at half time-step $\bar{\psi}_{i,j,k}^{n+1/2}$ is interpolated at cell interfaces by using Eqs. (3.4)–(3.4) and the time-centered density currents are computed by means of Eq. (26). 4. 4.

The continuity equation is solved and the new mass density \rho_{i,j,k}^{n+1}\big{|}_{C} is computed via Eq. (3.4). 5. 5.

The rescaling factor given by Eq. (30) is computed and $\bar{\psi}_{i,j,k}^{n+1}$ is rescaled in order to preserve mass conservation. 6. 6.

If the given oct lies next to the coarse-fine boundary, the estimated flux is stored for the subsequent reflux operation. 7. 7.

Finally, the phase rotation due to the potential is computed by applying Eq. (22) to $\bar{\psi}_{i,j,k}^{n+1}$, and $\psi_{i,j,k}^{n+1}$ is evaluated.

All steps, except the first and the last ones, are performed separately for each physical dimension. This procedure is called dimensional splitting and it reduces a $N$-dimensional problem into a system of $N$ one dimensional problems. It has the advantage of relaxing the CFL condition of the solver and, therefore, it allows bigger time-steps. A flowchart of our solver is shown in Fig. (1).

3.6 Refinement strategy

In SCALAR, we implement the same 'quasi-Lagrangian' approach as RAMSES uses for hydrodynamics: when the total mass of a given cell exceeds a given threshold, the cell is marked for refinement. The level dependent density threshold is defined as:

[TABLE]

where $M_{c}$ corresponds to the maximum mass allowed per cell.

In addition, following Schive et al. (2014) and the FLASH code (Fryxell et al. 2000), we implement support for the invariant version of the Löhner error estimator. It is based on the second derivative of a given physical quantity, normalised by the first derivative. Considering a generic physical quantity $f$, the error estimator $E_{L}$ reads:

[TABLE]

where the indices $i$,$j$ run over each physical dimensions. Small fluctuations of the physical quantity $f$ are filtered out due to the presence of the second term at denominator. The quantity $|\bar{f}_{ij}|$ is an average of $f$ over dimensions $i$,$j$ and $\xi$ is a small constant. This error estimator is dimensionless and therefore it can be applied to any physical quantity. Furthermore, in Eq. (33), $E_{L}$ is bounded in the interval $[0,1]$. In SCALAR, we apply the Löhner error estimator separately to ${\Re[\psi]}$ and ${\Im[\psi]}$. Then, the final estimation of the error on the wave-function is given by:

[TABLE]

and if it exceeds a user-defined threshold, the cell is marked for refinement. This threshold can be chosen empirically, depending on the features one needs to resolve in the solution of the governing equations. Although it is currently implemented in SCALAR, we did not employed the Löhner error estimator in the test we present here, in Section 4. In general, while testing separately the implementation of the Löhner error estimator, we find that a value of $E_{L}=0.7$ provides a good balance between computational cost and accuracy for the solution of the non-linear Schrödinger equation.

3.7 Spatial and temporal interpolation.

In SCALAR, interpolation is required when a generic level $\ell$ in the AMR hierarchy is refined and new child octs are created at level $\ell+1$, or when, during the solving process, boundary conditions need to be specified for a fine-grid patch and ghost cells are created. In both cases, the wave-function at the coarse level $\ell$ is interpolated down to level $\ell+1$. In order to solve the equations of motion, when the laplacian operator is applied, any discontinuity in the second derivative of the wave-function introduces an error, which propagates into the solution of the non-linear Schrödinger equation and it can destroy the wave-function. Therefore, high-order interpolation schemes are implemented in order to keep the wave-function as smooth as possible.

In particular, in SCALAR, we implement two high-order interpolation schemes. In both cases, the interpolating function is a fourth-order polynomial, but the coefficients of the polynomials are chosen in different ways. In one case, Lagrange basis polynomials are computed in order to set the coefficients, resulting in a fourth-order Lagrange interpolation scheme. In the second case, fourth-order conservative interpolation is performed and the coefficients of the interpolation are set by requiring that cell averages of the interpolated quantities are preserved. In case of adaptive time integration, linear temporal interpolation can also be applied when computing boundary conditions for a fine level patch, since coarse-grid and fine-grid wave-functions can be discretised at different times.

Furthermore, the interpolation can be performed on two different sets of variables: the original set of variables ${\Re[\psi]}$ and ${\Im[\psi]}$, or derived variables $m|\psi|^{2}$ and $\rm{Arg}[\psi]$, corresponding to mass density and phase of the wave-function. The interpolation schemes and the set of variables used for the interpolation process can be specified by the user in the parameter file.

Further details on the interpolation schemes can be found in Appendix B.

3.8 Artificial viscosity

In the tests shown in the upcoming sections, when they were done at the domain level only, the solution of the non-linear Schrödinger equation remains stable for as long as we could run SCALAR. However, when refinements were included, the solver had the tendency to develop spurious high-frequency waves at coarse-fine boundary, even after improving the order of accuracy of interpolation schemes. In order to artificially dump spurious oscillations, we introduced an empirical viscosity term in the non-linear Schrödinger equation. Thus, incorporating the viscosity term, the non-linear Schrödinger equation, Eq. (1), is replaced by:

[TABLE]

where the constant $\epsilon>0$ quantifies the strength of the damping term. For example, if we consider a single plane-wave:

[TABLE]

the viscosity term acts in a similar way as a Gaussian filter, by dumping the wave-function by a factor of:

[TABLE]

This means that, including the artificial viscosity, the wave-function is simply replaced by:

[TABLE]

In this way, the filter leaves untouched physical low frequency modes in the wave-function, while smoothing the spurious numerical oscillations. In general, an artificial viscosity term would affect mass conservation. However, by solving the continuity equation on top of the non-linear Schrödinger equation, mass conservation is enforced and the artificial viscosity simply acts as a viscous force. Indeed, if we consider the Madelung formulation of quantum mechanics instead, Eq. (6) and Eq. (7), the artificial viscosity term would be placed, together with the quantum force, in the momentum equation. Thus, Eq. (7) would read:

[TABLE]

The new viscous force term:

[TABLE]

helps preventing high-frequency waves to build up in time. With the addition of such a viscosity term, we are able to evolve without any issues the wave-function over hundreds of oscillation periods in our tests and, at the same time, preserving mass, energy and agreement with analytical solutions. There is no unique prescription for solving these issues we encountered and, in general, it is possible to design more elaborate artificial viscosity terms.

In an AMR context, when choosing the value of the dumping constant $\epsilon$, one should keep in mind that spurious high-frequency oscillations appear on scales of the local grid resolution, which is not fixed, but it changes according to designed refinement criteria. Thus, the strength of the dumping term should be decided such that it does not over-suppress the wave-function in high-resolution regions, but rather slightly under-suppress spurious oscillations in low-resolution regions. Such a limitation of the artificial viscosity term we implemented motivates further investigations to develop a more accurate scheme for dumping high-frequency spurious oscillations. Empirically, we find that a value in the range $0.2<\epsilon<1$ ensures stability over a long time in all the test cases we present (except the soliton test case, were we set $\epsilon=0$ and we do not use any artificial viscosity), by preventing the growth of spurious oscillations in the solution of the non-linear Schödinger equation.

3.9 Code units

We adopt the set of 'super-comoving coordinates' introduced in Martel & Shapiro (1998), and already used in RAMSES. Thus, the following change of variables is performed:

[TABLE]

where $H_{0}$ is the Hubble constant, $L$ is the box size, and $\overline{\psi}$ is chosen to ensure that $\int|\tilde{\psi}|^{2}\leavevmode\nobreak\ {\rm d}^{\rm dim}\tilde{x}=1$. As a consequence, the resulting non-linear Schrödinger equation reads:

[TABLE]

This set of coordinates was specifically designed for cosmological applications. However, it can be used for any application by setting the scale factor $a$ to unity and replacing $H_{0}$ by a general inverse time scale $T^{-1}$. In the remainder of the paper, all equations are in these code units.

For the particular case of axion dark matter in a cosmological setting (see Eq. (10)) we have $\tilde{g}=0$, $\overline{\psi}\propto a^{-3/2}$, and the potential is determined via the Poisson equation:

[TABLE]

where $\Omega_{i}$ is the fraction of the energy budget of our Universe that is in matter component $i$ (axions, baryons, CDM, etc.) and the mean value of $\tilde{\rho}_{i}$ over the box is set to unity.

4 Tests of the code

In this section we present the numerical experiments we performed in order to test the main features of SCALAR. When testing the accuracy of our numerical schemes, we rely on three main tests: conservation of mass, energy, and reproduction of analytical solutions. Given the total mass $M\left(t\right)$ in the simulation box and the total energy $E\left(t\right)$, the corresponding errors are respectively defined as:

[TABLE]

Instead, by denoting the analytical solution as $\psi_{\rm{a}}\left(x,t\right)$ and the numerical solution as $\psi_{\rm{n}}\left(x,t\right)$, at a given time, the error with respect to the analytical solution is computed according to the following formula:

[TABLE]

where $\left<\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \right>_{x}$ denotes the mean over the box. For the tests we present here, we always assume that each cell in the AMR grid has the same size in all dimensions: $\Delta x=\Delta y=\Delta z$. Furthermore, in all the tests, we only employ the refinement criterion based on density, as we found no need to use the Löhner error estimator in order to achieve an accurate and stable solution over time.

Accuracy and performances.

While running the tests we discuss in this section, we measured the overall accuracy and performances of the numerical algorithms implemented in SCALAR.

The accuracy of a numerical method is often measured by comparing the analytical solution with the numerical solution. For this purpose, we compute the error with respect to the analytical solution as given above. The overall accuracy of SCALAR is obtained by measuring how the global error, computed at a given time, scales with the grid spacing. In Fig. 2, we show the sample of errors we obtained for one of the numerical experiments we performed, the soliton test. The data points are fitted by the formula:

[TABLE]

where $x$ and $y$ represents the grid size and the error with respect to the analytical solution, respectively. The parameters of the fit $\alpha$ and $\beta$ are determined by means of the Linear Least Square (LLS) method and, in this case, we obtain $\alpha\sim 1.8$ and $\beta\sim 6.7$. While $\beta$ is not relevant when measuring the convergence of a numerical scheme, the parameter $\alpha$ represents the slope of the fitting polynomial and it corresponds to the global accuracy of the solver. Thus, for SCALAR, we achieve an overall second-order accuracy.

We measured general performances of SCALAR on a small cluster. Every node of the cluster is equipped with two Intel E5-2670 processors, with 8 cores each, a clock frequency of $f_{\rm CPU}=2.6$ GHz and a total memory of $128$ GB. Nodes within the cluster are interconnected through a Mellanox 56 Gb FDR Inifniband. In application with gravity enabled, SCALAR can evolve the solution of the non-linear Schrödinger equation by updating a single cell in the AMR hierarchy in $\sim 11$ $\mu$s. This means that, our code is able to reach a value of $\sim 10^{5}$ cell updated per second.

Conservation of mass.

The non-linear Schrödinger equation has the conserved quantity:

[TABLE]

which in the Mandelung formulation is just the total mass of the fluid. Mass and energy are not manifestly conserved by the main solver, therefore monitoring them is a useful test.

In simulations with no refinements and without enforcing mass conservation, we typically find the error on the conservation of mass of the order of (the prefactor is here for $\Delta x=2^{-6}$):

[TABLE]

where $T$ is the oscillation period. When we allow for refinements, the situation is typically worse and it is not good enough for cosmological simulations. However, by solving the continuity equation on top of the Schrödinger equation, we observe an improvement on the conservation of mass up to:

[TABLE]

This does not change when we allow refinements and, even though the error grows linearly in time, it is good enough in order to perform cosmological simulations. This is shown in Fig. (3), where we perform the sine wave test on the domain grid only, with a resolution of $N_{\rm cell}=2^{6}$ cells, corresponding to $\ell=6$ and $\Delta x=2^{-6}$, in one dimension. Details regarding the sine wave test are described in Section 4.1.

Please, refer to Section 3 for an explanation of the temporal evolution of the error in the conservation of mass.

Conservation of energy.

Since we enforce mass conservation by solving the continuity equation, energy conservation is a better accuracy test for our code. By defining kinetic and potential energy as:

[TABLE]

the temporal change in the total energy of the system is expressed by:

[TABLE]

As we can see, in the case where the effective potential has no explicit time-derivatives, the energy $E=K+W$ is conserved under the evolution.

In a cosmological setting, the potential $V$ depends on time via the scale-factor and this leads to a Lazyer-Irvine equation (Kopp et al. 2017):

[TABLE]

which can be monitored by integrating it up while performing the simulation.

Comparison to analytical solutions.

The most stringent test we can perform is to directly compare the numerical solution with an analytical solution. However, the discretised version of the governing equations is a different problem than the theoretical continuous limit and, thus, it admits a different solution. Usually, the main difference between the solutions of the continuous and the discretised non-linear Schrödinger equations leads to the wave-function evolving with slightly different temporal phases. For this reason, we estimate the phase difference between the two solutions and, starting from the analytical solution of the continuous non-linear Schrödinger equation, we compute the theoretical solution of the discretised equation. Then, we compare the numerical solution with the proper solution of the discretised non-linear Schrödinger equation. This represents a way of comparing the numerical solution with the one of the physical problem we are actually modelling. In the limit of $\Delta x\to 0$, the solutions of the continuous and the discretised non-linear Schrödinger equations converge to the same solution. In this way, we show that the numerical solution is in excellent agreement with the theoretical solution of the discretised problem, even in the case when the resolution is not extremely high, where computation time would significantly increase.

In the plots below we show, for different choices of initial conditions and potential, the error on the conservation of mass, the error on the conservation of energy, and the error with respect to the analytical solution.

4.1 Sine wave

The sine wave test evolves a static one dimensional density profile, where the initial wave-function is set as:

[TABLE]

It evolves in a constant potential, which is defined as:

[TABLE]

where $m$ is the mass carried by the wave-function and the period of oscillation is given by:

[TABLE]

The full analytical solution reads:

[TABLE]

This numerical experiment is designed to test the creation of ghost cells when computing fine levels boundary conditions. Here we only refine according to the mass criterion, Eq. (32) and, since the density profile does not evolve in time, there is no dynamical creation or destruction of grids: once the refinement map is computed at the beginning, it does not change. We evolve the solution of the Schrödinger equation over $100$ periods of oscillation. It is possible to show that the solution of the discretised equation - the one we are solving - is the same as Eq. (52), up to second order in space, but with a slightly different period of oscillation. Therefore, to factor out the dependence of the period with resolution (which we test seperately), we correct the analytical solution by replacing $T$ with $T_{\rm discrete}$, where:

[TABLE]

This test was performed for ${\rm dim}=1$. The non-linear Schrödinger equation is solved together with the continuity equation, in order to enforce conservation of mass. The domain grid resolution is $N_{\rm cell}=64$, corresponding to $\ell=6$ and $\Delta x=2^{-6}$, and the maximum refinement level is set to $\ell_{\rm max}=8$. When boundary conditions for fine levels are needed, phase and density are interpolated in ghost cells by means of fourth-order conservative interpolation. We used artificial viscosity with $\epsilon=0.2$. The results from this test are shown in Fig. (4).

We have also performed a similar test using a quadratic potential leading to a Gaussian profile:

[TABLE]

with very similar results.

4.2 Travelling wave

This test simulates a one dimensional wave-packet travelling through a periodic box. Here, we test dynamical creation and destruction of grids, since the AMR hierarchy follows the density profile moving towards the direction of the wave. In this case, we have no potential and the initial conditions are defined as:

[TABLE]

where $k_{1}=2\pi n_{1}$, $k_{2}=2\pi n_{2}$ with $n_{1}\not=n_{2}\in\mathbb{N}$. The oscillation frequency of a single mode is:

[TABLE]

and the analytical solution of the Schrödinger equation reads:

[TABLE]

As a consequence, the density is given by:

[TABLE]

where the oscillation period is defined as:

[TABLE]

The wave-function is evolved in time over $100$ oscillation periods. The non-linear Schrödinger equation is solved together with the continuity equation, in order to enforce conservation of mass. Also in this case, the coarse-fine data interpolation is performed by means of fourth-order conservative interpolation. However, while density and phase are interpolated in ghost cells, new refinements are made by interpolating real and imaginary parts of the wave-function. The domain grid has the same resolution as in the previous test and refinements are allowed up to $\ell_{\rm max}=8$. We used artificial viscosity with $\epsilon=0.2$. The results from this test are shown in Fig. (5).

4.3 Soliton

In a cosmological context, SCALAR can be used to simulate the structure formation process with fuzzy dark matter. In this case, the density profiles of the dark matter halos differs from the case of the standard CDM. We can find a stationary solution which can be tested by taking:

[TABLE]

and solving the resulting ODE for $\chi(x)$. A numerical fit to the density profile of a soliton in three dimensions was first suggested in Schive et al. (2014) and then in Marsh & Pop (2015). Despite it is an approximated solution, it is useful when coupling the Poisson equation to the non-linear Schrödinger equation. In this works, the density profile of the soliton was found to be on the form:

[TABLE]

where $r_{\rm core}$ can be chosen as a free parameter, see Appendix C for more details.

We set this density profile analytically and evolve the system. The density profile remains approximately stationary while the wave-function oscillates as:

[TABLE]

This test was performed for $\rm{dim}=3$ and the non-linear Schrödinger equation is solved together with the continuity equation. The domain grid contains $N_{\rm cell}=64^{3}$ cells, corresponding to $\ell=6$ and $\Delta x=2^{-6}$, and the maximum refinement level allowed is $\ell=8$. In this case, both refinement and ghost cells are made by fourth-order conservative interpolation on density and phase. The artificial viscosity term is set to $\epsilon=0$. We found that self-gravity is able to stabilise the wave-function against spurious numerical oscillations. We tested the same physical case with different values of artificial viscosity, but we did not find any improvement in terms of accuracy. In the other tests, where we set $\epsilon=0.2$, the artificial viscosity term was introduced just in order to ensure stability and accuracy of the solution to the non-linear Schrödinger equation over a very large number of oscillation periods. The results for this test are shown in Fig. (6).

5 Cosmological applications

The SCALAR code was originally developed in order to perform numerical simulations of structure formation with fuzzy dark matter. To demonstrate the capabilities of our code, in this section we present a test involving a full 3D cosmological setup.

The simulation box models a $B_{0}=1\leavevmode\nobreak\ \rm{Mpc}/h$ portion of the Universe expanding over time, according to the Einstein-de Sitter model. We take $\Omega_{\rm{\Lambda}}=0.7$ and $\Omega_{\rm{M}}=0.3$, and the Hubble constant is set to $H_{0}=100\leavevmode\nobreak\ h\leavevmode\nobreak\ \rm{km}\leavevmode\nobreak\ \rm{s}^{-1}\leavevmode\nobreak\ \rm{Mpc}^{-1}$, where $h=0.67$. The base resolution, which defines the domain level, is set to $N_{\rm cell}=256^{3}$ and up to one level of refinement is allowed. This means that, in this test, we achieve a maximum resolution of $\Delta x\sim 4\leavevmode\nobreak\ h^{-1}\rm{kpc}$. The mass of the boson is set to $m=10^{-21}\leavevmode\nobreak\ \rm{eV}$ and conservation of mass is enforced by solving the continuity equation on top of the non-linear Schrödinger equation. In order to better appreciate the differences between CDM and fuzzy dark matter, the initial conditions are computed for the case of CDM. The initial density and velocity fields are computed by using the Zel’dovich approximation. Then, following Kopp et al. (2017), we convert them into an initial wave-function.

The evolution starts at redshift $z=1000$ and we run the simulation as long as we are able to resolve the quantum force with at least two cells, meaning until redshift $z=30$. We did not run this simple test further because, due to lack of resolution, we cannot resolve the inner part of collapsed objects. This will be the subject of an upcoming paper.

In Fig. 9, Fig. 10 and Fig. 11 we show the dark matter density field, the real and the imaginary parts of the wave-function for a selection of redshifts, $z=200,100,50,30$. As the field clusters under the effect of gravity, the wave-function develops the wave patterns which are characteristic of this class of models.

By taking the density contrast $\delta\left(\mathbf{x}\right)$ inside the simulation box, we expand it in Fourier modes as follows:

[TABLE]

The matter power spectrum is defined by means of the autocorrelation function, which can be expressed as:

[TABLE]

In Fig. 7, we plot the dimensionless power spectrum, defined as:

[TABLE]

for $z=200,100,50,30$. The results here show the same quantitative behaviour as seen in Woo & Chiueh (2009, see Fig. 2), that performed the same kind of simulation as we do here. As the field gravitationally collapses, the quantum pressure leaves immediately its imprints on $\Delta\left(k\right)$ by producing the characteristic suppression of power at small scales, above $k\sim 300\leavevmode\nobreak\ h^{-1}\text{Mpc}$. Given the mass of the boson, the suppression scale is expected to be around the redshift dependent Jeans wave-number, which can be defined as:

[TABLE]

where $\Omega_{a}$ is the dimensionless density parameter of axions and $h$ is the dimensionless Hubble constant. Between $k=10\leavevmode\nobreak\ h^{-1}\text{Mpc}$ and $k=100\leavevmode\nobreak\ h^{-1}\text{Mpc}$, the power spectrum describes modes in the density field still in the linear regime and, therefore, $\Delta^{2}\left(k\right)$ evolves with redshift according to linear theory:

[TABLE]

in a similar way to the CDM case.

In Fig. 8, the evolution of errors in conservation of mass and energy are plotted against redshift. The error in the conservation of mass slowly evolves in time, as described in Section 4. Furthermore, we track the evolution of the error in the conservation of energy, by integrating the Lazyer-Irvine equation along the simulation and checking at which level of accuracy Eq. (48) is satisfied. As shown in the same figure, in the conservation of energy does not grow significantly, thus ensuring that no energy is numerically dissipated by the solver or the artificial viscosity term.

6 Conclusions

In SCALAR, we implemented a set of numerical algorithms developed in order to solve the non-linear Schrödinger equation in an AMR framework. Eq. (1) can be used to describe the dynamics of a Bose-Einstein condensate, a system of identical bosons in the ground energy state, by means of a single-particle wave-function in a mean field approach. Here, the non-linearity arises from an effective potential, which can contain both a self-interaction term and a generic external potential. Bose-Einstein condensates find their application in several fields. As an example, alternative dark matter models involving Bose-Einstein condensates have been recently developed, such as ultra-light axion dark matter, fuzzy dark matter and superfluid dark matter. The Schrödinger equation is solved with a Taylor method, similar to the algorithm developed in GAMER. In order to improve the conservation properties of the numerical scheme, the continuity equation is solved on top of the non-linear Schrödinger equation and mass conservation is enforced by construction. Empirically, by running several tests, we found that our numerical method is second-order accurate.

In order to test the main components of SCALAR, a test suite was designed. In particular, we tested the performances of the solver with and without solving the continuity equation on top of the non-linear Schrödinger equation, the creation of ghost cells when boundary conditions need to be specified for levels with fine resolution and the dynamical creation and destruction of grids during the process of mesh refinement. For this purpose, we tracked mass and energy conservation properties of our numerical schemes during the evolution of the system. The result is that both mass and energy are well conserved. While the latter remains roughly constant in all the different cases, the former evolves in time. Indeed, the error on the conservation of mass grows according to:

[TABLE]

for all our test-cases as we advance the solution in time. However, by solving the continuity equation on top of the non-linear Schrödinger equation, the error on the conservation of mass improves of several orders of magnitude, growing as:

[TABLE]

and it remains significantly small even for cosmological simulations. Furthermore, we compared the numerical solutions found by SCALAR with the analytical solutions of all the test cases. We show that the numerical solution tracks very well the analytical one over a long evolution time. In this case also, the error with respect to the analytical solution remains roughly constant over time. We showed that high-frequency spurious oscillations created at coarse-fine boundaries by interpolation schemes are efficiently dumped by an artificial viscosity term. However, the long term evolution of the single-particle wave-function still represents a challenge in case the artificial viscosity term is not included. We also run a small cosmological simulation where we show that SCALAR is able to capture the relevant features of models like fuzzy dark matter on cosmological scales.

Future work will aim at developing new high-order interpolation schemes which will not require the inclusion of an artificial viscosity term. Furthermore, we plan to compare the performance and accuracy of SCALAR with similar codes. In a following paper, we plan to exploit SCALAR to explore the non-linear regime of the structure formation with alternative dark matter models. In particular, we want to run a set of high-resolution cosmological simulation in order to verify and provide further predictions of fuzzy dark matter.

The code will be soon publicly available through our GitHub repository111http://github.com/mattiamina.

Acknowledgements.

We thank the Research Council of Norway for their support. Computations were performed on resources provided by UNINETT Sigma2 – the National Infrastructure for High Performance Computing and Data Storage in Norway. HAW was supported by the European Research Council through 646702 (CosTesGrav). We also thank the anonymous referee for several suggestions that helped improve the paper.

Appendix A Stability analysis

The stability condition of a generic PDE solver can be easily found by means of the Von Neumann stability analysis (Charney et al. 1950). For this purpose, the numerical error is decomposed in Fourier modes and a condition on time-step is computed in order to propagate each mode accurately.

We start by considering a generic mode in the Fourier decomposition of the numerical error:

[TABLE]

As mentioned in the previous sections, in SCALAR we use a second-order finite difference formula in order to approximate the laplacian of the wave-function. Thus, for a generic quantity, in one physical dimension we have:

[TABLE]

As a consequence, the second derivative of the error can be written as:

[TABLE]

and the amplification factor can be computed as:

[TABLE]

where $\beta$ corresponds to:

[TABLE]

In order to avoid exponential growth, we require that $|\xi|<1$. Therefore, the stability condition reads:

[TABLE]

where $\cos_{n}$ and $\sin_{n}$ denote to $n^{\rm th}$ order Taylor polynomials of $\cos(x)$ and $\sin(x)$, respectively. Furthermore, $n$ corresponds to the order of the Taylor expansion of the kinetic contribution to the time evolution operator, Eq.(3.3).

In particular, we find that for $n<3$ the numerical scheme is unconditionally unstable. For $n=3$, instead, the stability condition is satisfied as long as:

[TABLE]

The generalisation to $D$ can be done by replacing $(\Delta x)^{2}\to(\Delta x)^{2}/D$ in the formula above.

We also require that the phase angle does not rotate more than $2\pi C_{W}$ within a time-step. Thus, for the kinetic term we require that:

[TABLE]

while for the potential term:

[TABLE]

Combining the three conditions above, the optimal time-step is chosen as:

[TABLE]

where we require the safety factors to be $C_{W},C_{K}<1$.

Appendix B Prolongation operators

The details related to the interpolation schemes we implement in SCALAR are given for the one dimensional case. In case of multidimensional interpolation, the same formulas derived in this section are applied sequentially in each direction.

For simplicity, in order to derive the interpolation formulas we use in SCALAR, we assume that there is an odd number $2N+1$ of coarse data points $\{x_{i},y_{i}\}$ with $-N\leq i\leq N$ and the interpolation is always done for children of the central cell, as shown in Fig. (12).

In this way, we can construct an interpolating polynomial of order $2N$ starting from its general definition:

[TABLE]

In order to find an explicit expression for $P(x)$, the coefficients $\{a_{n}\}$ have to be computed and, depending on the constraints imposed on the interpolating polynomial, different schemes can be implemented.

Lagrange polynomial interpolation

For this interpolation scheme, instead of explicitly computing the coefficients $a_{n}$, the interpolating polynomial is expressed as a linear combination of Lagrange basis functions $l_{j}(x)$, with $j=-N,-N+1,\dots,N-1,N$, which are defined as:

[TABLE]

Thus, coarse data points $y_{i}$ weight the Lagrange basis functions and the interpolating polynomial is constructed as follows:

[TABLE]

In this way, the interpolating polynomial is forced to pass through the data points in the sample. In order to find an explicit formula for a fourth-order Lagrange interpolating polynomial, we directly compute the Lagrange basis functions by means of Eq. (76). Here, we express the differences between cell positions in Eq. (76) as a function of $\Delta x$. Then, the interpolating polynomial is computed at children cell positions, leading to the following interpolation formulas:

[TABLE]

Conservative polynomial interpolation

In this case, constraints on the coefficients of the interpolating polynomial are set by imposing that the mean of the interpolated data on fine cells is equal to the data stored in the coarse cell.

By using Eq. (77), the coarse data points are redefined by means of their cell average:

[TABLE]

and the corresponding linear system is solved in order to derive an explicit expression for the coefficients $\{a_{n}\}$ of the interpolating polynomial. Thus, assuming that the parent cell is split into two children cells, the fine data is obtained by solving the following integrals:

[TABLE]

where $\tilde{y}_{m}$ and $\tilde{y}_{p}$ denote the left and right childern cells respectively.

For a fourth-order polynomial, the solution of the linear system reads:

[TABLE]

Thus, the corresponding interpolation formulas are:

[TABLE]

It is trivial to check the arithmetic average of $y$ over the children cells corresponds exactly to the value stored in the parent cell.

Appendix C Soliton solutions

Self-gravitating bosonic fields can support stable and localised field configurations, where the density profile is static. Such configurations, called solitons, are ubiquitous in models of axion dark matter and exist for $\rm dim>1$. Starting with the Schrödinger-Poisson system (in code units):

[TABLE]

we take the spherically symmetric ansatz:

[TABLE]

The ODE describing the static configuration of the field can be derived by replacing Eq. (92) into the Schrödinger-Poisson system:

[TABLE]

where the laplacian is now expressed in polar coordinates:

[TABLE]

and $\kappa$ represents a free parameter (in our code units $\kappa$ equals $\frac{3}{2}\Omega_{m}a$ for a cosmological simulation with axion-like dark matter). Thus, the system can be solved by considering the initial conditions:

[TABLE]

where $\chi^{\prime\prime}(0)$ is a free parameter, which is set by requiring asymptotic vanishing solution $\chi(\pm\infty)=0$. The oscillation period $T$ can then be computed from the resulting solution, for more details see Schive et al. (2014) and Marsh & Pop (2015).

The soliton profile, normalised such that the total mass is unity, is given by:

[TABLE]

where $r$ is the distance from the center of the box and:

[TABLE]

with $\alpha=0.230$ and $\gamma=-0.692$. In our tests, even if this derivation assumes an infinite box, we use either fixed or periodic boundary conditions. As long as the core radius of the soliton is small enough compared to the simulation box, this solution represents a good approximation.