Tests of Gravity Theories Using Supermassive Black Holes

TL;DR

This paper investigates how supermassive black holes can be used to test theories of gravity that violate the strong equivalence principle, using observational signatures like black hole offsets in galaxy clusters.

Contribution

It proposes a method to detect strong equivalence principle violations via black hole displacements and provides new constraints on galileon gravity parameters from existing observations.

Findings

Black hole offsets can indicate violations of the strong equivalence principle.

Current observations of M87's black hole constrain galileon gravity parameters.

Black holes may escape galaxies within a billion years under certain conditions.

Abstract

Scalar-tensor theories of gravity generally violate the strong equivalence principle, namely compact objects have a suppressed coupling to the scalar force, causing them to fall slower. A black hole is the extreme example where such a coupling vanishes, i.e. black hole has no scalar hair. Following earlier work, we explore observational scenarios for detecting strong equivalence principle violation, focusing on galileon gravity as an example. For galaxies in-falling towards galaxy clusters, the supermassive black hole can be offset from the galaxy center away from the direction of the cluster. Hence, well resolved images of galaxies around nearby clusters can be used to identify the displaced black hole via the star cluster bound to it. We show that this signal is accessible with imaging surveys, both ongoing ones such as the Dark Energy Survey, and future ground and space based surveys. Already, the observation of the central black hole in M~87 places new constraints on the galileon parameters, which we present here. $\mathcal{O}(1)$ matter couplings are disfavored for a large region of the parameter space. We also find a novel phenomenon whereby the black hole can escape the galaxy completely in less than one billion years.

Equations

$${\nabla }^2\varphi +\frac{r_c^2}{3}\left[{\left({\nabla }^2\varphi \right)}^2-{\nabla }_i{\nabla }_j\varphi {\nabla }^i{\nabla }^j\varphi \right]=8\pi \alpha G\rho ,$$

$$\frac{\mathrm{d}\Phi }{\mathrm{d}r}=\frac{GM}{r^2}\left[1+2{\alpha }^2{\left(\frac{r}{r_{\mathrm{v}}}\right)}^{\frac{3}{2}}\right]$$

$$\ddot{r}=-\frac{4\pi }{3}G{\rho }_{0}r+a_{\mathrm{BH}}$$

$$\bar{r}=1\text{ kpc }\left(\frac{a_{\mathrm{BH}}}{2000(\text{km/s}{)}^2/\text{kpc}}\right){\left(\frac{{\rho }_0}{0.1M_{\odot }{\text{pc}}^{-3}}\right)}^{-1}.$$

$$T=10^7\text{ yr}{\left(\frac{R}{\text{kpc}}\right)}^{\frac{1}{2}}{\left(\frac{a_{\mathrm{BH}}}{2000(\text{km/s}{)}^2/\text{kpc}}\right)}^{-\frac{1}{2}}.$$

Full text

1 Introduction

The study of infrared modifications of general relativity (GR) has undergone a renaissance in the last decade, driven partly by the dark energy mystery and partly by recent developments in the construction of healthy higher-derivative scalar-tensor theories (see Clifton et al. (2012); Joyce et al. (2014); Koyama (2015); Berti et al. (2015); Bull et al. (2016); Joyce et al. (2016); Burrage & Sakstein (2016) for recent reviews), which have opened up a new realm of possibilities for driving the acceleration of the cosmic expansion. On small scales, consistency with tests of GR is achieved by utilising screening mechanisms, which use non-linear effects to hide the modifications of gravity in the solar system. One particularly interesting and well-studied paragon for these theories is the galileon (Nicolis et al., 2009), which self-accelerates cosmologically but screens locally using the Vainshtein mechanism (Vainshtein, 1972). On large scales, the galileon can be tested by the modified rate of structure growth (Barreira et al., 2012, 2013); on small scales, the Vainshtein mechanism renders its predictions largely indistinguishable from those of GR.

One interesting exception is that of Hui & Nicolis (2012), who have pointed out that the no hair theorem (Hui & Nicolis, 2013) implies that black holes do not couple to the galileon, and therefore there is a violation of the strong equivalence principle (SEP) whereby they fall at different rates to non-relativistic matter in external fields. In particular, if part of a galaxy’s motion is due to an external galileon field then the supermassive black hole (SMBH) that lies at its center does not feel this, and therefore the SMBH lags behind as the galaxy moves. Hui & Nicolis (2012) have argued that, for constant density cores, the restoring force due to the dark matter and baryons in the galaxy will eventually compensate for the lack of the galileon force, giving rise to an offset in the opposite direction to the galaxy’s acceleration. Hui & Nicolis (2012) considered galaxies moving in the cosmic field, and predicted a small offset, which is difficult to look for observationally.

In this letter, we look for novel scenarios for testing galileon gravity using the predicted SEP violation that are more amenable to observational testing. We consider two new effects:

1. (1)

Galaxy clusters: Satellite galaxies in-falling into galaxy clusters pass through a partially unscreened regime where the galileon field generated by the cluster pulls on the galaxies but not their resident black holes. This can result in an offset which is accessible to ongoing and planned imaging surveys. In this work we use the analysis of M 87 presented in (Asvathaman et al., 2015) to place new constraints. 2. (2)

Escaping black holes: Realistic dark matter haloes whose density falls with distance from the center have a maximum restoring force. Galileon forces larger than this drive the SMBH to escape the galaxy in an observable time-scale (less than one billion years).

2 Galileon Gravity Theories

Galileons are scalar-tensor extensions of GR that are invariant under the Galilean transformation $\phi(x^{\mu})\rightarrow\phi(x^{\mu})+b_{\mu}x^{\mu}+c$ (Nicolis et al., 2009). This symmetry restricts the form of their action to four unique terms (in four dimensions), the so-called quadratic, cubic, quartic, and quintic galileons Nicolis et al. (2009), and ensures that the equations of motion are second-order so that the Ostrogradski ghost instability is absent. In what follows, it will be sufficient to restrict our attention to the simplest model, the cubic galileon, although we note that the results presented here apply equally to all galileon theories including their curved space generalisations (Deffayet et al., 2009, 2011; Gleyzes et al., 2014a, b). In more general cases, there will be additional parameters and the tests we discuss later will constrain combinations of these. The equation of motion for the cubic galileon field sourced by a static non-relativistic object is

[TABLE]

where $i,\,j=1,\,2,\,3$. There are two constant parameters111If one considers the cosmological evolution then these become weak functions of redshift but this dependence is sub-leading provided that one considers the most general galileon theories (as one should), which include an $\mathcal{O}(1)$ constant contribution to $\alpha$ coming from a Brans-Dickie-like coupling. The tests we will discuss here are only sensitive to the sum of the constant and time-dependent pieces, and we take them to be dominated by the constant piece for simplicity. : a dimensionless coupling constant $\alpha\sim\mathcal{O}(1)$, which parameterizes the non-minimal coupling to matter, and the crossover scale $r_{c}$, which determines the size of the non-linear galileon self-interactions. The coupling scale $\Lambda_{3}=(6{M_{\rm{pl}}}/r_{c}^{2})^{1/3}$ is often used instead of $r_{c}$. For a point-mass $M$, one can define the Vainshtein radius $r_{\rm v}^{3}\equiv\frac{4}{3}\alpha GMr_{c}^{2}$, which is the scale where the non-linear galileon terms become important. The total (Newtonian plus galileon) force profile inside the Vainshtein radius is

[TABLE]

so that deviations from GR are suppressed by a factor $(r/r_{\rm v})^{3/2}$. Outside the Vainshtein radius the galileon force becomes unsuppressed, leading to a total force profile that matches the Newtonian one, except for an overall enhancement by a factor of $(1+2\alpha^{2})$. A solar mass object has $r_{\rm v}\sim\mathcal{O}(100\textrm{ pc})$ (Khoury, 2013) and so the solar system lies well within the Sun’s Vainshtein radius. For extended objects, such as the dark matter haloes we discuss in this work, the transition between the screened and unscreened regimes is a rather gradual one (see Schmidt (2010) and fig. 1 below). This is one reason why the outskirts of galaxy clusters are a prime testing ground for galileon theories.

In this theory, the equation of motion for an object of mass $m$ moving in an external Newtonian and scalar potential, $\Phi_{\rm N}^{\rm ext}$ and $\phi^{\rm ext}$ respectively, sourced by some other object is $m\ddot{\vec{x}}=-m\nabla\Phi_{\rm N}^{\rm ext}-\alpha Q\nabla\phi^{\rm ext}$, where $Q$ is the scalar charge that characterises the strength of the object’s coupling to the galileon field. One can show (Hui et al., 2009) that $Q=\int\,{\rm d}^{3}x\,T^{0}_{0}$ i.e. the scalar charge is equal to the baryonic mass222We use the term “baryons” as a proxy for anything that contributes to the matter energy-momentum tensor (as opposed to the pseudo-energy-momentum tensor which includes the contribution from gravitational binding energy). This could include actual baryons as well as dark matter. of the object. This implies that $Q=m$ for non-relativistic objects but that $Q<m$ for compact objects, whose mass receives a significant contribution from the gravitational binding energy. A black hole is the extreme example where $Q=0$. Hui & Nicolis (2013) showed that a black hole cannot source galileon hair, and therefore also does not couple to a galileon external field. Since our primary interest concerns a black hole residing in a galaxy, one might worry that the Vainshtein mechanism is enough to greatly diminish the galileon force on the stars (and dark matter) of the host galaxy, leading to essentially no difference between the falling motion of the stars and that of the black hole. A key feature of the galileon comes to the rescue: the galileon equation of motion is invariant under the galileon symmetry and so, given any solution, one can always add a component with a constant field gradient to obtain a second solution. For instance, for a galaxy located within a cluster, the galileon field sourced by the cluster behaves as a constant-gradient field, since its scale/wavelength is much longer than the size of the galaxy (Hui et al., 2009). The resident stars and dark matter of the galaxy thus respond to this cluster-sourced galileon field, while the black hole does not, setting up an astronomical version of the Eötvös experiment. Note that in this set-up, it is important to account for the effect of Vainshtein mechanism on the cluster-sourced galileon field itself.

3 Galaxy Clusters

A cluster carries sufficient mass that, despite some Vainshtein suppression, the galileon force sourced by it is large enough to give interesting observational effects at distances D\;\mbox{\raisebox{-2.15277pt}{\stackrel{{\scriptstyle>}}{{\scriptstyle{\sim}}}} }\;0.1R_{200} (Schmidt, 2010). As an example, consider a model of the Virgo cluster with mass $M_{200}=10^{15}M_{\odot}$. We model the mass distribution of the cluster with a Navarro-Frenk-White (NFW) profile (Navarro et al., 1996) (with concentration $c=5$) inside $R_{200}$; outside it we model the “2-halo” mass distribution using the fits to N-body simulations of Diemer & Kravtsov (2014). We define $R_{200}$ using the critical density, which corresponds to $R_{200c}$ in Diemer & Kravtsov (2014); we take $H_{0}=72$ km/s/Mpc. Note that we consider haloes at $z=0$ since the tests described presently are intended to be applied to low redshift ($z<0.05$) clusters. Any potential time-evolution of $\alpha$ over this range will be negligible, although we note that more detailed modelling may be required if a sample of useful clusters at redshifts $z\gg 0.1$ were to be found. Previous studies have consistently found that the properties of halos are largely unmodified in galileon cosmologies (Barreira et al., 2014b, 2015), especially inside the Vainshtein radius. Outside, at distances larger than $\sim 2$ Mpc, more detailed modelling may be required, although M 87, which we consider below, lies well within the Vainshtein radius of the Virgo cluster so this is not a caveat to the results we obtain here.

The Newtonian, galileon, and cosmological root mean square (RMS) galileon forces (the scenario considered by Hui & Nicolis (2013)) are plotted in Fig. 1 for models with $r_{c}=500$ and $6000$ Mpc333The latter corresponds to self-accelerating cosmological galileons while the former is often used as a paradigm for Vainshtein screening (Schmidt, 2010). These models are not ruled out by local observational constraints (Dvali et al., 2003; Khoury, 2013). The constraints of Schmidt (2009); Barreira et al. (2014a) coming from cosmological probes do not directly apply to our model since they are derived for a sub-class of models where the constant contribution to $\alpha$ is absent. Future efforts to constrain more general galileon models, in particular those with a linear coupling to matter, will be able to use the results we obtain below as a consistency check.. The galileon force was calculated by solving the field equations for the cubic galileon exactly using spherical symmetry. One can see that there is a large galileon force for r\;\mbox{\raisebox{-2.15277pt}{\stackrel{{\scriptstyle>}}{{\scriptstyle{\sim}}}} }\;0.1R_{200}\sim 0.3 Mpc. Satellite galaxies are typically in-falling towards the cluster within \;\mbox{\raisebox{-2.15277pt}{\stackrel{{\scriptstyle<}}{{\scriptstyle{\sim}}}} }\;4R_{200}\sim 13 Mpc (Diemer & Kravtsov, 2014). These galaxies would feel a large galileon force that their central SMBHs would not, thus allowing the effects of the SEP violation to manifest. Note that the galileon force is suppressed relative to the Newtonian force when r\;\mbox{\raisebox{-2.15277pt}{\stackrel{{\scriptstyle<}}{{\scriptstyle{\sim}}}} }\;R_{200} but it is the amplitude of the galileon force that is important for SEP violations (not this ratio) and so large effects are still expected at distances above $0.1R_{200}$. The cluster forces are at least an order of magnitude larger than RMS cosmological force. One therefore expects offsets of $\mathcal{O}(\textrm{ kpc})$, which are accessible to imaging surveys of nearby galaxy clusters.

The dynamics of this situation are easily exemplified by considering a galaxy where the central density $\rho_{0}$ is approximately constant. In this case, the black hole is subject to the galileon acceleration $a_{\rm BH}$ away from the center and a linear restoring force given by $4/3\pi Gr\rho_{0}$, in which case the equation of motion for the black hole (in the rest frame of the galaxy) is

[TABLE]

so that the black hole undergoes oscillations about the equilibrium point $\bar{r}=3a_{\rm BH}/4\pi G\rho_{0}$ with a time-period $T_{0}\sim(G\rho_{0})^{-1/2}$. In practice, the galileon force turns on slowly over a time $\sim\dot{a}_{\rm BH}/a_{\rm BH}\sim v/d\sim\mathcal{O}(10^{10}\textrm{ years})$ (where $v\sim 100$ km/s is the typical in-fall velocity and $d$ is the distance from the cluster center) whereas the oscillation period is $\mathcal{O}(10^{7}\textrm{ years})$ so that the amplitude of these oscillations is small and the black hole tracks the equilibrium position adiabatically. One then expects the black hole to be offset from the center by a distance

[TABLE]

Clearly the size of the offset depends on the central density $\rho_{0}$ and the external force; we have used a fiducial external force that is typical for the outskirts of a galaxy cluster similar to the one modelled in fig. 1. As the galaxy in-falls, the galileon force increases and so does the offset of the black hole (up to around $R_{200}$); the direction of the offset would be in the opposite direction to the galaxy’s direction of acceleration i.e. away from the center of the cluster.

An example of this (taking $\alpha=1$) is shown in Fig. 2; one can see that observable offsets of $\mathcal{O}(\textrm{kpc})$ can develop at a variety of distances from the cluster’s centre. They are at least an order of magnitude larger than one would expect for galaxies moving in the cosmological field, and so this scenario is well suited to testing galileon gravity. As is evident form the figure, as the galaxy in-falls, the offset becomes larger due to an increased galileon acceleration. This reaches a maximum and begins to fall off as the galaxy approaches the cluster’s center, where the galileon modifications are more screened.

Asvathaman et al. (2015) considered the effect of the sub-clump of the Virgo cluster centered on M 84 and M 86 on the position of the black hole within M 87. Using the same assumptions as Asvathaman et al. (2015), in particular the dynamical model of Walsh et al. (2013) and the photometry of Ferrarese et al. (2006), we find that the strength of the galileon field in the plane of the sky is less than about $700~{}(\mathrm{km/s})^{2}/\mathrm{kpc}$ (or $1000~{}(\textrm{km/s})^{2}/\textrm{kpc}$ in three dimensions). We find the SMBH to be located at the centre of the light to within 0.03 arcseconds (0.03 is the upper limit at 1$\sigma$). We are only sensitive to the motion of the SMBH across the sky because we are measuring its position on the sky relative to the center of the galaxy; consequently, only the component of the galileon field in the plane of the sky is measurable. The distribution of light and mass-to-light ratio in the context of a Hernquist model for the central region of the galaxy was used to obtain an estimate of the restoring force.

Combining this with a model for the Virgo cluster as depicted in Fig. 1 i.e. NFW + 2 halo modelling with $c=5$ and $M_{200}=10^{15}M_{\odot}$ 444masses in the range $M_{200}=1$–$2\times 10^{15}M_{\odot}$ have been reported (Fouque et al., 2001; Peirani & de Freitas Pacheco, 2006), which is less uncertain than effects of modelling the complex 3D structure (Mei et al., 2007) with an NFW profile so we adopt the conservative value of $10^{15}M_{\odot}$. we obtain constraints on the values of $\alpha$ and $r_{c}$ as depicted in Fig. 3. These were found by scanning the parameter space to find the region where the galileon force at $R_{200}$ is smaller than $1000~{}(\textrm{km/s})^{2}/\textrm{kpc}$. The constraints therefore assume that M 87 is located within $R_{\rm 200}$ of the Virgo cluster, but is otherwise insensitive to its precise location because the galileon force is fairly constant for $r<R_{\rm 200}$. We have checked that varying the concentration over the range $c=5$–$10$ produces a minimal change in the constraints. Finer constraints could be obtained by measuring the displacements of many SMBHs relative to the centers of their galaxies and the local density field and averaging over an ensemble. Such an analysis could probe theories with $r_{c}>1500$ Mpc, however, the observations of M 87 already yields tighter constraints on this model than solar system measurements. There are several caveats to our analysis, including the use of spherical symmetry in the predictions (Asvathaman et al., 2015). We leave a detailed analysis of the observational strategy and implementation for future work.

The constant density approximation works well at small distances but the offsets predicted in Fig. 2 lie well in the regime where the fall-off of the density profile with distance becomes important. The density near the center of galaxies tends to fall off with some power $\rho\propto r^{-\beta}$ with $\beta<1$ (Trujillo et al., 2004; Ferrarese et al., 2006) that steepens at larger radii. This results in a rising restoring force that reaches some maximum and begins to decrease when the profile steepens555Profiles with $\beta\geq 1$ result in maximum or divergent forces at the centre but, in practice, the SMBH will drag a disc of stars with it so that the restoring force vanishes at the centre (Asvathaman et al., 2015). The situation described above is therefore generic. . Small galileon forces, such as those generated by the cosmological large scale structure, are not sufficient to overcome the maximum restoring force and hence lead to an offset. Larger forces, such as those that can be generated by clusters can overcome the maximum restoring force and can therefore cause the black hole to escape the galaxy completely. Even for the case that the restoring force is larger at the center, a black hole that is displaced by other mechanisms can be driven out of the galaxy by the galileon force. Such a situation might arise as a result of a merger event or by recoil from gravitational wave “kicks”.

In practice, if the in-falling galaxy is initially at a position where the galileon force is smaller than the restoring force then the black hole will be offset but, as the galaxy moves closer to the cluster, the galileon force will increase and may exceed the restoring force so that there is nothing halting the SMBH’s motion. One can estimate the time-scale for the black hole to escape the galaxy entirely by neglecting the restoring force. Galaxies falling into clusters experience an increasing galileon force whilst the restoring force remains the same and so this quickly becomes a good approximation. Treating the galileon acceleration as constant, the time-scale for the black hole to move a distance $R$ is

[TABLE]

The black hole can then escape the galaxy in less than a billion years. Note that the typical velocity of the black hole is $\mathcal{O}(\textrm{km/s})$, far less than the typical velocity dispersion of the galaxy so that dynamical friction can safely be ignored (Binney & Tremaine, 2008). Precisely which galaxies would allow escape depends on the environment and the galaxy central density. We leave a systematic investigation for the future.

4 Observational Tests

The novel features identified above present new avenues for testing galileon gravity and constraining the parameters $\alpha$ and $r_{c}$ (recall that $\alpha$ parameterises the strength of the galileon coupling to matter and $r_{c}$ parameterises the galileon’s self-interactions).

Displaced SMBHs in nearby galaxies: The gravitational acceleration is several times larger at the outskirts of galaxy clusters than in the field, and is directed towards the cluster center. The predicted offset is detectable via the star cluster bound to the SMBH, which would be offset from the centroid of the stellar light at the center of the galaxy. The correlation with the expected displacement direction from modified gravity, i.e. the opposite direction to the galaxy cluster, is critical in such an exercise as it narrows the search zone and can distinguish it from other sources of fluctuations.

The two key observational parameters are the size of the offset relative to the point spread function (PSF), which determines whether it is resolved, and the change in flux due to the star cluster bound to the SMBH. At $z\simeq 0.05$, an offset of 1 kpc corresponds to 1 arcsecond on the sky, which is larger than the (full width half maximum) PSF of the best ground based telescopes. Out to this distance, hundreds of clusters can be found for which the offsets in in-falling galaxies are resolved. A large fraction of these are already imaged by the SDSS and DES optical surveys and by X-Ray surveys. By identifying suitable galaxies at the outskirts of each cluster, a large sample of galaxies can be assembled. These galaxies could be followed up with high resolution imaging, ideally from the Hubble Space Telescope, and multi-wavelength observations. The sample size may be essential in handling the second observational challenge: detecting the star cluster given the typical flux variations across a galaxy, and other sources of error. For the central parts of elliptical galaxies, these are at the percent level (e.g. Bernardi et al. (2017)), which in many cases is smaller than the displaced SMBH’s star cluster (e.g. Asvathaman et al. (2015)’s study of M87). Another challenge discussed by Asvathaman et al. (2015) is centroiding the galaxy light using the outer isophotes. Finally, we note that the predicted offset is sensitive to the central density profile of the host galaxy, which is challenging to determine. A range of offset values must therefore be considered, and for lower values only the most nearby clusters may be suitable for our test.

Perturbations to the galaxy light profile: while we have not investigated morphological features in this study, a SMBH that is displaced or drifting through the stellar disk will produce characteristic distortions that can be measured by analyzing the images of a large sample of galaxies. Even for galaxies above $z\sim 0.1$, where the offset is no larger than the PSF, the model fitting approach described above can be attempted. Weak lensing studies that measure the low order moments of the surface brightness of galaxies may also be well suited to extracting the skewness or “flexion” in the light distribution, and correlating it with the direction of the external force vector on the galaxy.

Missing SMBHs: for larger galileon forces, or if initially displaced, the SMBH would leave the visible galaxy in less than a billion years, the timescale over which galaxies move with coherent velocities. Hence, a fraction of galaxies would not have central SMBHs in such a scenario. During galaxy mergers, SMBHs are displaced from their center and then occupy steeper parts of the density profile until they lose energy to dynamical friction and sink to the center. The timescales may be sufficient for the modified gravity effect to act on the black hole. The observational finding that the overwhelming majority of galaxies above a certain mass have central SMBHs may provide new limits on $\alpha$ and $r_{c}$ for typical values of the density profile.

Note that for many of these observable tests, we can estimate the direction of the effect from observations of the galaxies surrounding the cluster of interest by estimating the direction of the local gravity acceleration vector. The direction vector is more reliable for galaxies in-falling into galaxy clusters. This is crucial for distinguishing modified gravity effects from other astrophysical processes such as acceleration due to asymmetric jets, recoil from gravitational wave emission, Brownian motion, gravitational slingshot due to mergers, and perturbations due to massive objects such as globular clusters. All of these can displace the central black hole for varying periods of time, but in a direction uncorrelated with the galaxy’s acceleration. These processes are important to study even in the absence of modified gravity and may shed light on the dynamics of galaxy mergers (Merritt & Milosavljevic, 2005).

5 Conclusions

To summarize, we have identified novel methods for testing gravity using the equivalence principle violations for black holes first noted by Hui & Nicolis (2012). Galaxies in-falling into massive clusters may host SMBHs that are offset from the center by $\mathcal{O}(\textrm{kpc})$ and may even be absent altogether. The offset is at least an order of magnitude larger than for galaxies moving in the cosmological field, and is in the opposite direction to the galaxy’s acceleration, which can help distinguish this effect from other astrophysical displacement mechanisms. We have discussed how a sample of nearby galaxy clusters can be used to obtain bounds on the model parameters $\alpha$ and $r_{c}$, which would have strong implications for models of modified gravity that try to explain cosmic acceleration. Indeed, using observations of the central black hole in M 87 we have been able to place new constraints on the parameter space that push into the interesting region for self-accelerating models.

*Acknowledgements: * We are grateful to Eric Baxter, Mariangela Bernardi, Gary Bernstein, Chihway Chang, Neal Dalal, Marla Geha, Jenny Green, Mike Jarvis, Gordon Richards, Ed Moran, Nadia Zakamska, Justin Alsing and colleagues at the Center for Computational Astrophysics. JS is supported by funds provided to the Center for Particle Cosmology by the University of Pennsylvania. BJ is supported in part by the US Department of Energy grant DE-SC0007901. JH is supported by a Discovery Grant from the Natural Sciences and Engineering Research Council of Canada. LH is supported in part by the NASA grant NXX16AB27G and the DOE grant DE-SC0011941.