Crater Morphology of Primordial Black Hole Impacts

One long standing question in astrophysics and cosmology is the nature of dark matter (e.g. Bertone & Hooper, 2018, and reference therein). Candidates for dark matter particles can be divided into two categories: microscopic and macroscopic, where the latter are often referred to as MACHOs (massive compact halo objects). It has been suggested that dark matter is composed (at least in part) of black holes formed from density fluctuations in the early universe (Carr et al., 2016). Such objects are called primordial black holes (PBHs). Since the idea has been proposed, different observations have excluded certain mass ranges for these primordial black holes (Carr et al., 2020; Carr & Kühnel, 2020). Currently, there are several poorly constrained mass windows, notably between $10^{16}\,\rm g$ and $10^{24}\,\rm g$, where the constraints are weak enough that PBHs could be a large fraction of the dark matter, though the exact mass windows are dependent on the PBH mass spectrum. While several constraints have been reported in this mass range some have been contested (Katz et al., 2018), which motivates the present work to develop a new independent constraint in this asteroid-mass to sublunar-mass window.

In this study we consider the possibility that a primordial black hole collided with the moon in the past. As we shall show later, the craters from such collision will have a shape markedly different from an ordinary meteor impact. This feature can readily be detected by surface images and radar scans.

The idea of a collision between a primordial black hole (or another compact object) with planetary bodies has been considered by various authors in the past. Rafelski et al. (2013b) have considered in some detail the physics of compact ultra-dense matter colliding with solar system bodies. They found that relativistically compact matter will generally completely penetrate the body forming both ‘entrance’ and ‘exit’ wounds on a hyperbolic trajectory. Penetration is generally supersonic which may produce a seismic signal (Khriplovich et al., 2008) which has been searched for (Herrin et al., 2006; Banerdt et al., 2007; Horowitz & Widmer-Schnidrig, 2020). Accretion onto a high velocity PBH ‘entering’ the earth has been proposed as an explanation for the Tunguska event (Jackson & Ryan, 1973). ‘Exit’ events have meanwhile been compared to volcanic events, transporting crust material high into the atmosphere (Rafelski et al., 2013b, a). While surface processes on the earth complicate efforts to establish historic events as PBH encounters, the Moon and other solar system bodies without active surface processes hold a record of several billion years of impacts in their surface geology. In this work we develop the theory of cratering dynamics in PBH encounters with solar bodies and describe a possible observational campaign to constrain low mass MACHO dark matter.

While we generally consider PBHs interacting with the Moon in this work, our discussion generalises to many solar system bodies and ultradense dark matter candidates.

The structure of the paper is as follows. In section 2 we present analytic calculations of the shock wave trajectories in the case of a regular cratering event and in the case of a PBH impact. Using the shock trajectories, we calculate how the ejecta will be distributed outside the crater, and show that each one of these cases will yield a unique ejecta distribution. In section 3 we present numeric hydrodynamic simulations that verify our analytic results. Finally, in section 4 we discuss the results and the implications for PBHs as dark matter particles.

In order to verify the analytic results derived in the previous section we performed hydrodynamic simulations using the moving mesh code RICH (Yalinewich et al., 2015). In the first simulation we collide a cold sphere with a cold half space, both filled with ideal gas with the same density. This simulation is supposed to represent a normal cratering event. We chose the adiabatic index to be $\gamma=1.4$, since this value best describes the high pressure shock behaviour in silica (Yalinewich & Schlichting, 2019). Since we want to resolve the impactor and also follow the shock wave to distances much larger than the impactor radius, we arranged 20,000 mesh generating points in a non-uniform way inside the computational domain. In the vicinity of the impactor the resolution is 0.1 of the radius of the impactor, but at larger radii the cells are larger such that the ratio between the cell size and radius is roughly constant around 0.1. With this distribution of points, we were able to extend the boundaries of the computational domain to $10^{6}$ impactor radii. To simplify the simulation, we assume a head on collision, and run the simulation on a 2D grid using cylindrical symmetry. The density of both the impactor and the target is set to 1 while the density in the rest of the computational domain is $10^{-9}$, the velocity with which the impactor travels toward the target is set to 1 and the ambient pressure before impact is $10^{-9}$ everywhere. We ran the simulation to time $10^{7}$, where the unit of time here is the time it takes the impactor to traverse its own radius. The simulation does not include gravity, so the debris expands ballistically from the impact site. In figure 1 we plotted the density of the ejecta as a function of distance from the impact site. Using the results of the previous section we can estimate the slope of this profile analytically. We saw that $v\propto m^{-\beta}$. At distances much larger than the radius of the impactor, the velocity profile is roughly homologous $v\approx d/t$, where $d$ is the distance from the impact site and $t$ is the time since the impact. Combining these relations with the definition of density $m\approx\rho d^{3}$ we find a scaling relation for the density profile $\rho\propto r^{-3-1/\beta}$. For $\beta=2/3$ we obtain $d\ln\rho/d\ln r=-4.5$, which is close to the value obtained in the numeric simulation $d\ln\rho/d\ln r=-4.1$.

The second simulation is intended to describe a head on collision of a primordial black hole with the surface of the moon. The setup for this simulation is the same that described in the previous paragraph, except the density of the impactor is $10^{6}$ times greater than that of the target, and also we increase the resolution along the trajectory of the impactor inside the target, so the total number of points is larger (120,000). The radial density profile from this simulation is shown in figure 2. From this figure we infer the slope $d\ln\rho/d\ln r=-5.9$. Our analytic theory predicts that $\beta=1/3$, and so the slope of the density profile should be $d\ln\rho/d\ln r=-6$. The numeric simulations confirm our analytic prediction, and show that the behaviour of a PBH collisions are markedly different from a normal cratering event. In figure 3 we present sequential log density snapshots from both simulations to highlight the differences between them.

In this work we considered the possibility that a primordial black hole collided with the moon. We show that a crater resulting from such a collision will have a different shape from a regular meteor impact, namely, a compact object impact crater will have steeper ejecta blanket. Most craters on the moon have ejecta blankets whose shape is compatible with a meteor impact. However, there are craters with steeper ejecta blankets whose origin does not involve exotic objects (e.g. Mare Imbrium, see Zhu et al., 2017).

Steep ejecta blankets can therefore be used to flag craters of interest, but are not indisputable evidence of a compact object impact. A stronger evidence for black hole impact is the presence of a high pressure, pyrite-type phase of silica (Kuwayama et al., 2005). Regular impacts on the moon only produce the stishovite and seifertite phases of silica. This is because the new phase of silica appears above 40 TPa, whereas peaks pressure in a typical collision is around 1 TPa (Kayama et al., 2018). We therefore propose that steep ejecta blankets can be used to flag craters of interest, but to confirm a black hole impact origin would require finding the high pressure phase of silica.

If such a crater did exist on the moon, there is a possibility that it will have eroded away, due to meteor impacts or solar wind sputtering. For a one metre crater, the erosion time is over a billion years (McDonnell et al., 1977), and for larger craters take even longer. Hence, we conclude that erosion is not a concern.

The smallest black holes that produce a detectable one metre crater have a mass of about $10^{19}$ g. We estimated that even for these small black holes, the impact rate is so low that, even if all dark matter is composed of such black holes, there’s only about a 10% chance that one collided with the moon throughout its lifetime. With such a low rate we are not able to place constraints on the dark matter distribution. However, the same considerations used here for the moon could apply just as well to any rocky body without an atmosphere. Therefore, with future missions to other rocky bodies in our solar system it will eventually be possible to use the method outlined in this work to constrain the mass distribution of primordial black holes.

AY is supported by the Natural Sciences and Engineering Research Council of Canada (NSERC), funding reference #CITA 490888-16. This work made use of the sympy (Meurer et al., 2017), numpy (Oliphant, 2006) and matplotlib (Hunter, 2007) python packages.

The source code for the numerical simulation and documentation can be found on github at https://github.com/bolverk/huji-rich.