# A Stable Finite-Volume Method for Scalar-Field Dark Matter

###### Abstract

We describe and test a new numerical method to solve the Schrödinger equation in self-gravitating systems, e.g. Bose-Einstein condensates or “fuzzy”/ultra-light scalar field dark matter. The method is a finite-volume Godunov scheme with stable, higher-order accurate gradient estimation, based on a generalization of recent mesh-free finite-mass Godunov methods. It couples easily to particle-based N-body gravity solvers (with or without other fluids, e.g. baryons), manifestly conserves momentum and energy, is numerically stable, and computationally efficient. We consider a variety of test problems and demonstrate that it can accurately recover exact solutions and remains stable even in noisy, poorly-resolved systems, with dramatically reduced noise compared to other proposed implementations. This is non-trivial because the “quantum pressure” is neither isotropic nor positive-definite and depends on higher-order gradients of the density field. We implement and test the method in the code GIZMO.

###### keywords:

methods: numerical — hydrodynamics — cosmology: theory — elementary particles — dark matter## 1 Introduction

Recently there has been a surge of astrophysical interest in numerical methods to solve the Schrödinger-Poisson equation for scalar fields with high occupation number, driven primarily by recent interest in axionic, “Bose-Einstein condensate” (BEC), “ultra-light” scalar-field or “fuzzy” dark matter models. In these models the dark matter (DM) is a boson with mass eV which forms a “heavy” high occupation-number “condensate” (so acts like a massive “cold” DM candidate on large scales), but exhibits coherent quantum effects on smaller scales. At scales eV, the de Broglie wavelength of the DM reaches parsec to kpc scales, leading potentially to interesting astrophysical consequences (for reviews, see Suárez et al., 2014; Hui et al., 2017). Therefore it is particularly interesting to couple quantum dynamics to traditional N-body gravity methods used for cosmological simulations.

Many numerical methods to date are variations of that in Mocz &
Succi (2015) (although see e.g. Schive
et al. 2014; Schwabe
et al. 2016; Kopp
et al. 2017): they noted that one can write the relevant equations in the Madelung form (e.g. Spiegel, 1980), which amounts to coupling the “usual” collisionless Euler equations with self-gravity (solved by -body methods) to an additional “quantum pressure tensor” (QPT) (i.e. force ). However the QPT is neither isotropic nor positive-definite and depends on first and second-derivatives of the density field (so the force depends on third-derivatives); as such it is difficult to devise accurate and numerically stable schemes. For example, the original Mocz &
Succi (2015) scheme and most derivatives are based on an extension of the smoothed-particle hydrodynamics (SPH) equation-of-motion with evaluated at particle locations using standard SPH gradient estimators, but these prove to be noisy and numerically unstable, which is especially problematic when one wishes to follow non-linear collapse of cosmological halos. Meanwhile spectral methods or those which directly evolve the wavefunction (Schive
et al., 2014; Schwabe
et al., 2016) do not conserve mass or couple easily to standard -body codes. In this paper, we therefore introduce a finite-volume Godunov formulation using higher-order accurate matrix-based gradient estimators. We consider a series of test problems to demonstrate the stability and accuracy of the method, in the code GIZMO.^{1}^{1}1The public version of GIZMO is available at: http://www.tapir.caltech.edu/~phopkins/Site/GIZMO.html

## 2 Numerical Method

The dynamics of a self-gravitating Bose-Einstein condensate/scalar field at high occupation number are well-described by the Schrödinger-Poisson equation, , with particle mass , external potential , gravitational potential obeying (including both self-gravity and other contributions), and wavefunction . The Madelung transformation identifies the density (the particle mass times the quantum probability density of the wavefunction) and velocity (the group velocity of a wavepacket) to re-write this as a pair of equations, the Euler equations with the QPT : i.e. and

(1) | ||||

(2) |

For a de Broglie wavelength much smaller than the horizon there is no effect on the global expansion of the Universe and the correction terms to Eq. 1 in an expanding Universe are identical to those for a normal fluid (they do not effect the QPT; see Hui et al., 2017).

### 2.1 Volume Decomposition: Finite-Volume Formulation

We solve Eq. 1 by modifying the Lagrangian, mesh-free finite-volume Godunov (specifically the “meshless finite mass” or MFM) method from Hopkins (2015) – developed for fluid dynamics – to treat the QPT analogous to other anisotropic pressure terms already well-studied (e.g. in kinetic MHD; Hopkins, 2017). This retains the advantages of Lagrangian methods and equal element masses for N-body solvers, while allowing for higher-order solutions to the normal fluid equations in a manifestly-conservative manner.

Briefly, each traditional -body DM “super-particle” now becomes a mesh-generating point at coordinate used for volume decomposition: there is a well-defined volume associated with each , and conserved quantities (e.g. mass, momentum) carried by the particle represent the integral over all DM particles in the domain (rather than the mass/velocity of a Monte-Carlo “super-particle”). Gravity and external potential terms are operator-split and solved in the usual fashion (see below). Writing Eq. 1 in conservative form, the usual finite-volume transformation applies:

(3) | ||||

where is the (oriented) face area between and neighboring volume , and is the interface value of . The problem is now defined by an exchange of a conserved quantity (momentum) between neighbors across a face, so is manifestly conservative.

In MFM, the second-order accurate effective face is given by:

(4) | ||||

(5) | ||||

(6) |

where is the kernel-density defined local mesh-generating point number density, , and is the kernel function with scale-length set adaptively to match the local kernel-averaged inter-element separation ().^{2}^{2}2In this paper we will use a cubic spline for , with maximum radius of compact support ; other choices have small effects on our results. The volume of domain is , so the density within the domain is . This is a “smoothed” Voronoi decomposition.

### 2.2 Gravity & Force-Softening

As noted above, gravity is operator-split. This can be solved by the usual N-body methods, with one important improvement: since we have already de-composed the volume into a well-defined continuous density field, we must solve the gravity equations for the same density field. This amounts to using the fully-adaptive gravitational softening method (where the softening description is equivalent to the solution of the Poisson equation for the continuous density field) as described in detail in Hopkins (2015); Hopkins et al. (2018) for the DM. This already includes the corrections for the non-linear self-gravity terms as the domain described by a single element expands/contracts, so is manifestly conservative of momentum and energy (see Price & Monaghan, 2007).

### 2.3 Gradient Estimation: Calculating the QPT

Eq. 1 depends on higher-derivatives of , so it is critical to choose accurate, consistent, and stable gradient estimators. We therefore determine the first-derivative via the second-order accurate and consistent matrix gradient estimators

(7) |

This provides significant advantages in stability and accuracy over other numerical gradient estimation techniques (for extensive discussion, see Maron et al., 2012; Luo et al., 2008; Lanson & Vila, 2008a, b; Hopkins, 2015; Mocz et al., 2014; Pakmor et al., 2016). For example, the gradient is robust to noise, as compared to face-centered gradient estimators, and its consistency does not depend on the details of the particle arrangement around the point . As shown in Gaburov & Nitadori (2011), this is also the internally-consistent gradient operator for the MFM face operators. We therefore use this to calculate , then iteratively re-apply to each gradient component to calculate the “gradients of gradients” needed for . This provides the required second-order accuracy, stability, and automatically expands the stencil by iteration (for examples and tests, see Muñoz et al., 2014) – in other words, since already involves a neighbor loop for each neighbor , the neighbor stencil used for the second-derivatives is automatically (implicitly) larger, as needed to give robust results. Thus

(8) |

### 2.4 The Riemann Problem & Numerical Stability

What remains is to define the interface flux . Following the MFM method, this is given by the solution of the appropriate Riemann problem at the face, assuming the face is moving with the contact wave (i.e. defined by the mesh motion such that the net mass flux across the face vanishes, as the mesh moves with the flow). We define the left and right states at the face with a piecewise-constant reconstruction, e.g. , (while higher-order reconstructions are in principle straightforward, the already high-order derivatives in make this very noisy and difficult to stabilize). Following Hopkins (2015) where the difference between exact Riemann solvers and the HLLC solver was negligible, we approximate using the HLLC solver for an anisotropic pressure tensor (well studied in elastic dynamics). If the elements and are receding with a speed in excess of , the maximum signal speed in the left or right state (), then the vacuum solution applies (no contact wave exists). Otherwise, in the frame moving with the contact wave (), the interface pressure is

(9) |

where , , the initial wavespeeds are , (Gaburov & Nitadori, 2011), is the identity matrix, and is a limiter function. For , note that the linearized system Eq. 1 has one wave, a longitudinal mode with dispersion relation , so .

The limiter precedes the numerical diffusivity that arises from up-wind mixing in the Riemann problem. Some diffusivity is required for numerical stability, especially important here because can be negative and mesh-generating points move (which can source well-studied numerical instabilities such as the “tensile instability” from elasto-dynamics; Swegle et al. 1995). Taking and equal to its maximum possible value (for a grid-scale particle separation , ) ensures unconditional stability, but produces a very diffusive scheme. This is especially un-desirable in the regime of interest when QPT forces are sub-dominant to gravity (one does not wish to artificially dissipate gravitational motions). By analogy to standard treatments of the tensile instability in negative-pressure tensors in elasto-dynamics (e.g. Monaghan, 2000) and previous studies of anisotropic viscous stress tensors in the MFM method (Hopkins, 2017), we define

(10) |

where and represent the interface value without numerical dissipation and the dissipation term, respectively; denotes the Frobenius norm, and is a weight ( here, which is sufficient for stability without being too diffusive). We also estimate with ), (evaluated with the interface gradient values). Here , where , is the dimensionless (pair-averaged) kernel function, and is the value of this function evaluated at the mean inter-element separation. The term is discussed in Appendix A, and exists to ensure energy conservation. We stress that the detailed form of and have weak effects on our conclusions (entering only in the numerical dissipation, which becomes vanishingly small in converged solutions). But it respects three important properties: (1) it vanishes if the QPT does, (2) it vanishes for mesh-generating points which are well-separated/weakly interacting, or diverging/at rest (these only feel the QPT term), and (3) it rises rapidly (allowing the “full” numerical dissipation) for mesh-generating points which are approaching at small radii within the kernel.

### 2.5 Integration & Timestepping

We integrate Eq. 1 using the standard explicit, leapfrog scheme in GIZMO using adaptive hierarchical timesteps (see Hopkins 2015 for details). The usual timestep critieria apply: e.g. acceleration-based criteria for N-body methods ( with ; Power et al., 2003), and the Courant criterion needed for adaptive gravitational softening methods (); Hopkins et al., 2018). But the QPT imposes its own timestep criterion: because depends on higher numerical derivatives, and the system admits whistler-type waves with , the fastest wavespeed is and numerical stability under explicit integration require a timestep criterion of the form:

(11) |

where here. We emphasize this is critical for stability unless implicit integration is used, something which has not been widely recognized in this area. Quadratic timestep criteria (for e.g. whistler waves, diffusion/conduction problems, etc.) can be prohibitive at very high resolution (small ), motivating fully-implicit integration schemes, but for the DM problems of interest here the requirement is not particularly costly.

## 3 Numerical Tests

We now consider various test problems, all run using the identical code described above. Unless otherwise specified, we do not include gravity or other forces and work in dimensionless units with . All tests are run in 3D – this is important for methods without a regular (e.g. rectilinear Cartesian) grid, where the error properties can be very different even in problems where the analytic solution depends on only one dimension, if particle positions vary in another dimension. We adopt equal particle masses, as (a) density variations require particle-position variations (e.g. an irregular mesh), which makes accurate gradient estimation more challenging, and (b) this is almost always the choice in cosmological simulations.

### 3.1 Instantaneous (Accuracy/Consistency) Tests

First consider a non-dynamical test of the algorithm, to verify that it recovers (numerically) the correct accelerations. We initialize a density profile in a 3D box of arbitrarily large size, and save the numerically-calculated acceleration for each particle.

Hyperbolic-Tangent Profile: Following Nori & Baldi (2018), consider a profile varying along the -axis with analytic form:

(12) |

where . So the density varies from some at to at , over a width of in . Fig. 1 plots the recovered discrete at various resolutions. Even this simple density profile produces a complicated which is not symmetric and has multiple features. With a few resolution elements “across” the density jump, our method recovers the main qualitative features; we see rapid convergence with decreasing .

For comparison, we show the results using the “SPH-like” formulation in Mocz &
Succi (2015); Veltmaat &
Niemeyer (2016); Zhang et al. (2018); Nori &
Baldi (2018).^{3}^{3}3In Nori &
Baldi (2018),
(13)
where refers to the gradient of the kernel function evaluated at , and is the usual SPH correction for variable smoothing lengths (Springel &
Hernquist, 2002). As discussed in Nori &
Baldi (2018), a variety of SPH formulations are possible (see e.g. Price 2012; Hopkins 2013 for general discussion of these degrees-of-freedom), but of those they consider this was the most accurate so we compare it here. We have specifically also tested the variant formulations in Mocz &
Succi (2015); Veltmaat &
Niemeyer (2016); Zhang et al. (2018) and find they all exhibit comparably poor gradient recovery and are all numerically unstable in dynamical tests.
Even at high resolution, the SPH-like result is noise-dominated and systematically biased (as Nori &
Baldi 2018 also showed). This owes to well-known problems: the SPH gradient estimators and equation-of-motion () cannot be made zeroth-or-first order consistent/accurate without being unstable (Price, 2012), and higher-order SPH gradients (especially the form) are notoriously noisy.

3D Gaussian: We have also compared the 3D profile:

(14) |

where and . Our results for this test are qualitatively identical to the profile so we do not show both.

### 3.2 Static/Steady-State Tests

Now add full time-evolution, but consider steady-state solutions.

Ground State of a Damped Harmonic Oscillator: Particles are laid randomly (from a uniform PDF) in a periodic cube with side-length and mean density , with uniform random velocities from to in each direction. We add to the acceleration a term with , corresponding to a 1D harmonic oscillator potential () with frictional damping () to force decay to the ground state. Fig. 1 plots the density distribution well after all velocities have decayed, and the exact ground state .

Advection of a (1D) Density Mode: We initialize in an arbitrarily large box, with uniform velocity ; this has so should advect as . Comparing at shows good agreement; because the profile is smooth (a second-order polynomial), any numerical inaccuracy in the initial is small (much smaller than our test), so is vanishingly small (any initial noise from the finite-sampling of is resolved within a few timesteps by the diffusion in the Riemann problem). Since our Lagrangian method trivially solves advection and is manifestly Galilean-invariant, this is not a challenging test.

### 3.3 Linear/Wave Dynamics

Oblique Traveling (Whistler-Type) Waves: Consider Eq. 1 with , where , are homogeneous, and linearize in the perturbed quantities then make the usual Fourier ansatz: . The linearized equations feature one non-trivial eigenmode: with . So initialize an eigenmode: with , and , , in a 3D periodic box with side-length unity, , and the particles laid on an initially Cartesian grid rotated by the arbitrarily-chosen Euler angles before being perturbed to generate the density profile, so that both and are oblique to the “grid” and each other. We evolve the system for 40 wave periods ().

Our implementation produces convergence and numerical stability. We have repeated this with and different oblique angles or and obtain essentially identical results. As usual, the numerical dissipation in the Riemann problem produces some damping, but this decreases at higher resolution (as desired). Without dissipation ( in Eq. 9), the low-resolution solutions are closer to exact at early times (since they are undamped) but we see numerical instability set in later and destroy the solution. We have also verified that the solutions become numerically unstable if we remove the quadratic timestep condition (§ 2.5), as expected. With the HLLC diffusion term but no limiter (, ), the wave is always damped, even at high resolution (because the grid-scale wavespeed increases at high resolution as , one cannot rely on resolution alone without a limiter here).

Simple Harmonic Oscillator: We initialize , , and add to the acceleration equation the term (corresponding to ). This has exact solution . Although fully non-linear, we find this is a much “easier” test than the whistler-type waves above, as (1) the dispersion relation is simpler (acoustic-like, rather than whistler-like, making it much easier to stabilize), (2) there is less “structure” in the wave, (3) the external potential “corrects” the system (as opposed to being self-generating), and (4) the large (non-linear) amplitudes of the gradients make it somewhat less sensitive to numerical noise. Therefore we do not show the results but simply note the convergence is substantially faster than for the whistler-type wave.

### 3.4 Strongly Non-Linear Tests

Two-Stream Intersection: Following Hui et al. (2017) consider the non-linear evolution of a converging flow with a Gaussian super-position of plane waves as a convenient (and exactly analytically solveable) proxy for the intersection of oppositely-moving streams. In a large periodic, cubic box, and with is an exact solution at all times . We initialize this at , corresponding to a converging flow on the origin, and run it until , at which point it should be exactly mirrored as a diverging flow (representing the plane waves collisionlessly “passing through” one another, or equivalently a “bounce” as the quantum pressure diverges when the converging flow collapses). The scheme is able to recover the diverging solution, which demonstrates it can capture intersecting streams and quasi-collisionless phenomena, and also that the numerical dissipation is not so large to damp the motion (preventing the “bounce” and simply merging into a non-moving soliton at the origin).

Conservation Test: To “stress-test” conservation and stability in the code, we initialize a periodic box with side-length, mass, gravitational constant, and equal to unity, laying particles down randomly according to a Poisson distribution with random initial velocities drawn from a 3D Gaussian with dispersion (with zero net momentum). We evolve this, including self-gravity, until time . We verify (1) the code remains numerically stable, (2) timesteps and densities remain positive-definite, (3) mass is (trivially) conserved, (4) linear momentum is manifestly conserved to machine accuracy if we adopt a single timestep for all particles, (5) with adaptive (individual) timesteps, variations in timestep between neighbors mean the simple implementation of our scheme is not machine-accurate momentum-conserving, but we find the linear momentum error is less than at all times.

### 3.5 Dark Matter Halos & Cosmological Simulations

Isolated Halo: We initialize a spherically-symmetric, Hernquist (1990) mass profile (), self-gravitating () sphere with an isotropic velocity distribution function, using particles. For truly collisionless particles (no QPT), this is an exact equilibrium state. We then enable the QPT, with . The potential at is in these units in the collisionless case (so velocity dispersion ), and so the ground state of the halo with the QPT, while it has no closed-form solution, should approximately be a similar mass profile at large radii, with a constant-density soliton at the center of size . We consider cases with and without an additional damping : without damping the central density oscillates owing to its initial highly out-of-equilibrium central density, with damping the system quickly converges to a stable equilibrium “ground state” with these properties.

Cosmological Simulation: We consider a cosmological, DM-only “zoom-in” simulation of a DM halo (at ) in a Mpc cubic volume, initialized at with a standard CDM cosmology, and zoom-in (high-resolution) Lagrangian region including all particles within virial radii of the halo at , with particle mass ( particles in the halo). The specific halo is m10q with all details in Hopkins et al. (2018). We use physical units with eV, which should produce a kpc soliton in the halo center. We run to to verify that our method allows us to stably evolve scalar-field DM through non-linear cosmological history.

## 4 Discussion

We have developed and tested a method to solve the Schrödinger-Poisson equation, relevant at high occupation number for simulations of e.g. Bose-Einstein condensates or axionic/ultra-light scalar-field (“fuzzy”) DM. The method is based on the finite-volume, meshless finite-mass hydrodynamic scheme in Hopkins (2015). It couples in a straightforward manner to particle-based N-body gravity solvers, maintains equal particle masses, manifestly conserves mass and momentum, is numerically stable, can trivially be solved with additional fluids present, and adds only a fixed computational overhead for simulations which already compute gravitational forces with adaptive gravitational softenings. We implement the method in the code GIZMO.

Unlike some previous implementations discussed in the literature, we show that the methods here recover gradients accurately and retain stability even given highly dis-ordered particle configurations. These are critical in chaotic situations including N-body gravity and cosmological simulations, but even in simple test problems (e.g. a traveling oblique wave) they are challenging, because the “quantum pressure” force () depends on first, second, and third derivatives of the density field, and the “quantum pressure tensor” () is neither isotropic nor positive definite. We present a series of test problems and show it is viable for high-resolution cosmological simulations of scalar-field DM. In future work, we intend to investigate these DM models in greater detail, in cosmological hydrodynamic simulations.

## Acknowledgments

We thank Victor Robles for a number of helpful discussions. Support for PFH was provided by an Alfred P. Sloan Research Fellowship, NSF Collaborative Research Grant #1715847 and CAREER grant #1455342. Numerical calculations were run on the Caltech compute cluster “Wheeler,” allocations from XSEDE TG-AST130039 and PRAC NSF.1713353 supported by the NSF, and NASA HEC SMD-16-7592.

## References

- Gaburov & Nitadori (2011) Gaburov E., Nitadori K., 2011, \mnras, 414, 129
- Hernquist (1990) Hernquist L., 1990, \apj, 356, 359
- Hopkins (2013) Hopkins P. F., 2013, \mnras, 428, 2840
- Hopkins (2015) Hopkins P. F., 2015, \mnras, 450, 53
- Hopkins (2017) Hopkins P. F., 2017, \mnras, 466, 3387
- Hopkins et al. (2018) Hopkins P. F., et al., 2018, \mnras, 480, 800
- Hui et al. (2017) Hui L., Ostriker J. P., Tremaine S., Witten E., 2017, \prd, 95, 043541
- Kopp et al. (2017) Kopp M., Vattis K., Skordis C., 2017, \prd, 96, 123532
- Lanson & Vila (2008a) Lanson N., Vila J.-P., 2008a, SIAM J. Numer. Anal., 46, 1912
- Lanson & Vila (2008b) Lanson N., Vila J.-P., 2008b, SIAM J. Numer. Anal., 46, 1935
- Luo et al. (2008) Luo H., Baum J. D., Löhner R., 2008, Journal of Computational Physics, 227, 8875
- Maron et al. (2012) Maron J. L., McNally C. P., Mac Low M.-M., 2012, \apjs, 200, 6
- Mocz & Succi (2015) Mocz P., Succi S., 2015, \pre, 91, 053304
- Mocz et al. (2014) Mocz P., Vogelsberger M., Sijacki D., Pakmor R., Hernquist L., 2014, \mnras, 437, 397
- Monaghan (2000) Monaghan J. J., 2000, Journal of Computational Physics, 159, 290
- Muñoz et al. (2014) Muñoz D. J., Kratter K., Springel V., Hernquist L., 2014, \mnras, 445, 3475
- Nori & Baldi (2018) Nori M., Baldi M., 2018, \mnras, in press, arXiv:1801.08144,
- Pakmor et al. (2016) Pakmor R., Springel V., Bauer A., Mocz P., Munoz D. J., Ohlmann S. T., Schaal K., Zhu C., 2016, \mnras, 455, 1134
- Power et al. (2003) Power C., Navarro J. F., Jenkins A., Frenk C. S., White S. D. M., Springel V., Stadel J., Quinn T., 2003, \mnras, 338, 14
- Price (2012) Price D. J., 2012, Journal of Computational Physics, 231, 759
- Price & Monaghan (2007) Price D. J., Monaghan J. J., 2007, \mnras, 374, 1347
- Schive et al. (2014) Schive H.-Y., Chiueh T., Broadhurst T., 2014, Nature Physics, 10, 496
- Schwabe et al. (2016) Schwabe B., Niemeyer J. C., Engels J. F., 2016, \prd, 94, 043513
- Spiegel (1980) Spiegel E. A., 1980, Physica D Nonlinear Phenomena, 1, 236
- Springel & Hernquist (2002) Springel V., Hernquist L., 2002, \mnras, 333, 649
- Suárez et al. (2014) Suárez A., Robles V. H., Matos T., 2014, in Moreno González C., Madriz Aguilar J. E., Reyes Barrera L. M., eds, Astrophysics and Space Science Proceedings Vol. 38, Accelerated Cosmic Expansion. Springer International Publishing, Switzerland, p. 107 (arXiv:1302.0903), doi:10.1007/978-3-319-02063-1_9
- Swegle et al. (1995) Swegle J. W., Hicks D. L., Attaway S. W., 1995, Journal of Computational Physics, 116, 123
- Veltmaat & Niemeyer (2016) Veltmaat J., Niemeyer J. C., 2016, \prd, 94, 123523
- Zhang et al. (2018) Zhang J., Sming Tsai Y.-L., Kuo J.-L., Cheung K., Chu M.-C., 2018, \apj, 853, 51

## Appendix A Ensuring Energy Conservation With Un-Resolved Density Structure (Un-Resolved Quantum Potentials)

A general problem arises when we wish to solve Eq. 1 at finite resolution. The Schrödinger-Poisson equation conserves energy, and this is manifest in the Madelung form (Eq. 1) if we define the total energy , i.e. the sum of the usual kinetic and gravitational potential terms, plus the “quantum potential” . But since depends only on , which depends only on the mass configuration, i.e. on the initial conditions and subsequent , then if there is any numerical error in integration, or non-vanishing numerical dissipation terms , these cannot conserve energy (they change the velocities but not or, therefore, ). In usual hydrodynamics, this is not a problem, because there is an additional thermodynamic variable (e.g. temperature, entropy), so although the numerical terms like may change the kinetic energy, this can always be exactly offset by an appropriate exchange with thermal energy (so the terms are diffusive, but manifestly conservative).

In order to restore manifest conservation, we must introduce a “thermodynamic-like” variable into which the energy lost to (inescapable) numerical dissipation can be “stored.” This is the origin of the term in Eq. 9 – it is (numerically) analogous to a standard thermodynamic pressure. Each element carries a scalar (initialized to ), which evolves and defines according to:

(15) | ||||

(16) |

Here is a “source term” – any work done by the numerical diffusion term is added to the “energy” (because of the form of , it is straightforward to show this is positive-definite). As builds up, this generates an isotropic “pressure” , which appears in the equation of motion (Eq. 9), and therefore does work which must be added/removed from the energy. The coefficient , as it is trivial to show that the “effective equation of state” for (Eq. 1) is under isotropic compression/expansion. Note that this also imposes an additional timestep criterion, where , but this is almost always irrelevant compared to the stricter Eq. 11.

This has a very direct physical interpretation. Imagine a field with a single Fourier-mode perturbation that is compressed by rapid inflow (e.g. two streams colliding, in § 3.4), ignoring self-gravity. This should reach some maximum compression, where , so the energy is entirely “stored” in , then it will “bounce” back. At maximum compression, the specific energy , i.e. if there was a large initial kinetic energy, there must be small-scale (high-) modes in . Since depends only on , if the initial kinetic energy of compression (or e.g. gravitational force driving the compression) is sufficiently large, then eventually the required must exceed , the numerical resolution. Regardless of the details of the numerical dissipation, any numerical method with finite resolution will fail at some point to capture the small-scale oscillations/gradients of the density field ( cannot exceed ), but since the flow must still reach a point where , this means, naively, that energy will be lost somewhere and the “return” or “bounce” will be less energetic. The term stores the energy of “unresolved oscillations” in the density field – i.e. the un-resolved high- modes that contribute to , and therefore should contribute to the quantum pressure. In this example, when two parcels with velocities and intersect within a resolution element (so ), the kinetic energy which is no longer present (and should go into sub-grid-scale oscillations of ) is correctly retained in , which then, correctly, “pushes” alongside the resolved quantum pressure from the numerically-calculated .

It is easy to show that if the maximum (smallest oscillation scale) is resolved, then the terms are completely negligible (as they should be). The field therefore also serves a practical numerical purpose: by comparing to the total resolved energy (kinetic plus gravitational plus computed from the resolved density gradients), one immediately obtains an estimate of convergence (specifically, the fraction of energy which has gone into un-resolved modes). For all the test problems in § 3.1-3.4, is trivially satisfied (the mode structure is resolved); so removing the terms makes no visible difference whatsoever to the results shown in Fig. 1. However, these are simple test problems with well-defined characteristic wavenumbers. For highly non-linear cosmological simulations, where external forces (gravity) are often strongly dominant over the quantum pressure, and especially when one wishes to consider larger boson masses (smaller de Broglie wavelengths), it becomes important to account for the unresolved terms. Otherwise, we find that the (numerically) dissipated energy leads to excessive growth of the central solitons at halo centers, which only converges away at extremely high (often impractical) resolution.

Of course, information is still necessarily lost when the oscillations become sub-resolution. Not only does the contribution of different modes become un-resolved, but also, in the simple approach here, the fact that we only record a scalar means that we have also lost direction information. In 3D, if we compress the field to an unresolved peak in one direction, the restoring force will act isotropically. In future work, we will explore whether we can improve on this by replacing with an appropriate tensor that records the un-resolved components of the compression along each direction.