DETECTABLE SEISMIC CONSEQUENCES OF THE INTERACTION OF A PRIMORDIAL BLACK HOLE WITH EARTH

The determination of the constituents of dark matter is an outstanding problem that continues to receive great attention (Feng 2010). Constraints derived from the gravitational microlensing survey Expérience de Recherche d'Objets Sombres (Alcock et al. 1998; Tisserand et al. 2007) limit at most 8% of our galaxy's dark matter halo to be comprised of objects with mass M^PBH between 10²⁶ g and 2.4 × 10³⁴ g (for comparison, Earth's mass is approximately 6 × 10²⁷ g). Primordial black holes (PBHs), which could constitute a significant fraction of these massive objects at a range of masses (Zel'Dovich & Novikov 1967; Hawking 1971; Carr & Hawking 1974), are subject to Hawking radiation losses (Hawking 1974), a process by which they emit photons of a characteristic frequency ∝1/M^PBH. PBHs with masses M^PBH < 10¹⁵ g would thus evaporate within the age of the universe, and those with M^PBH < 10¹⁷ g would emit potentially detectable gamma rays (Page & Hawking 1976; Sreekumar et al. 1998; Carr et al. 2010). Based on these criteria, PBHs in the mass range 10¹⁷ g < M^PBH < 10²⁶ g are unconstrained and the entirety of dark matter in our galaxy could in principle be constituted by them. One way to narrow these constraints is to search for special oscillations of the Sun (and potentially other stars) induced by the gravitational interaction with proximal and impacting PBHs. The impacts of PBHs of the size of asteroids (M^PBH = 10²¹ g) with the Sun would be detectable by existing solar observatories, but the rate of occurrence of such events is estimated at a disappointingly low 10⁻⁸ yr⁻¹ (Kesden & Hanasoge 2011).

The prospect of black holes (BHs) striking Earth excites considerable public interest. The PBHs we consider here are of masses small enough that they do not pose an existential threat to our planet, and because they are at the lower end of the allowable mass range, will cause negligible non-seismic effects. However, for masses the size of asteroids and larger, devastation similar in some manner to the impact of an asteroid may ensue. Consider, for example, the Tunguska impact and explosion in Siberia on 1908 June 30. No comet ash has ever been recovered around this region, leading to the speculation that it might have occurred due to an impact with a PBH (Jackson & Ryan 1973). Subsequently, this theory was dismissed because of (among other reasons) the lack of secondary signatures around the exit location of the supposed PBH (Beasley & Tinsley 1974; Burns et al. 1976a). Another example is the putative identification of supersonic waves in recorded seismograms, hypothesized to have been caused by the impact of a "quark nugget" (QN) with Earth (Anderson et al. 2003) on 1993 November 24. This claim was later falsified when it was shown that the signals were actually due to a hitherto unreported earthquake (Selby et al. 2004). QNs, also termed strangelets, are a hypothetical phase of matter comprised of roughly equal parts of up, down, and strange quarks that have lower binding energy than iron, and thus may be the "ultimate" stable state of matter (Witten 1984). Because a QN is comprised of materials at nuclear density and can be of any mass up to a neutron star mass, it will have a very similar gravitational effect as a PBH when passing through Earth, and present a similar seismic signature. However, smaller QNs, of order M = 10¹⁵ g, could still produce observable explosions in the atmosphere, in contrast to PBHs of the same mass. This could be used to distinguish these two classes of compact object collisions with Earth.

The speed of PBHs in unbound orbits (Lisanti & Spergel 2011) is on the order of 200–400 km s⁻¹, indicating a maximum transit time through Earth of less than a minute. During its passage, a PBH would exert gravitational forces on the whole globe and generate seismic waves that propagate through Earth's interior. Due to the fast movement of the PBH compared to typical seismic wave speeds, such an event may be treated effectively as a supersonic source that generates Mach waves. The signatures of Mach waves are different from those of earthquake-induced seismic waves, thereby providing a means of detecting PBHs. Amplitudes of such waves have also been estimated (Khriplovich et al. 2008), albeit in a simplistic scenario where Earth is treated as a homogeneous fluid ball.

Gravitational forces in the near-field of the PBH are very strong, leading to nonlinear effects such as the cracking of rocks and heating. However, in this paper we focus on the far-field, treating the PBH as a moving Newtonian gravitational source to the seismic wave equation (Kesden & Hanasoge 2011), which we solve using the publicly available SPECFEM3D_GLOBE software package. An outline of the rest of the paper is as follows. In Section 2, we introduce the governing equations and describe the numerical solution method. In Sections 3 and 4 we discuss the results, in the latter section focusing on the unique spherical normal modes induced by a PBH. In Section 5, we provide some order-of-magnitude estimates of event rates and non-seismic effects. In Section 6, we suggest a strategy for detection, before concluding in Section 7. We relegate to Appendices A and B convergence tests and a few additional simulation results.

2.1. Governing Equations

In seismology, the displacement field s of particles within Earth's volume Ω is governed by the seismic wave equation

$$\begin{equation} \rho \, \partial ^2_t {\mathbf {s}}- \mbox{\boldmath $\bf \nabla $}\cdot {\mathbf {T}}= {\mathbf {f}} , \end{equation} \tag{ 1 }$$

where ρ is the mass density, T is the symmetric second-order stress tensor, and f is the source. Equation (1) is a linearized version of conservation of linear momentum, which proves to be sufficiently accurate for seismic wave propagation in Earth's interior.

In an elastic material, the stress tensor T is determined in terms of the strain via Hooke's law,

$$\begin{equation} {\mathbf {T}}= {\mathbf {c}}{\bm :} {\bm \varepsilon } , \end{equation} \tag{ 2 }$$

where c is the fourth-order elastic tensor with at most 21 independent components, describing the elastic properties of the material, and ${\bm \varepsilon }$ is the strain tensor

$$\begin{equation} {\bm \varepsilon =\frac{1}{2}\left[\mbox{\boldmath $\bf \nabla $}{\mathbf {s}}+\left(\mbox{\boldmath $\bf \nabla $}{\mathbf {s}}\right)^T\right]} . \end{equation} \tag{ 3 }$$

In an isotropic elastic medium, the fourth-order tensor reduces to two independent parameters—the bulk modulus κ and the shear modulus μ. Furthermore, the shear modulus μ vanishes in an acoustic medium (e.g., the outer core). Therefore, only the parameter κ is relevant in fluids. It is important to note that because of this difference between solids (inner core, crust, and mantle) and liquids (outer core), the governing equation takes slightly different forms, as will be further discussed later. Complications associated with attenuation are readily captured by an absorption band model (Dahlen & Tromp 1998) and are accommodated in numerical simulations based on the introduction of memory variables (e.g., Komatitsch & Tromp 1999).

In global/regional seismology, the source term f represents earthquakes, whereas in exploration seismology, active sources, such as man-made explosions, are usually employed. To study the effects of PBHs, neither source type is relevant. Following Kesden & Hanasoge (2011), we consider a classical BH, for which the source term f may be written as the gradient of the gravitational potential ϕ of the BH:

$$\begin{equation} {\mathbf {f}}= \mbox{} {-}\rho \, \mbox{\boldmath $\bf \nabla $}\phi . \end{equation} \tag{ 4 }$$

The potential ϕ obeys Newton's law of gravitation

$$\begin{equation} \phi ({\mathbf {x}},t)=\mbox{}{-}\frac{G\,M^\mathrm{PBH}}{ \|{\mathbf {x}}-{\mathbf {x}}^\mathrm{PBH}(t)\|} , \end{equation} \tag{ 5 }$$

where G is the universal gravitational constant and M^PBH is the mass of the PBH. Note that the gravitational potential (hence the gravitational force) changes with time as the PBH is fast moving with a trajectory

$$\begin{equation} {\mathbf {x}}^\mathrm{PBH}(t)={\mathbf {x}}^\mathrm{PBH}_0-{\mathbf {v}}^\mathrm{PBH} \, t , \end{equation} \tag{ 6 }$$

where x^PBH₀ is its initial position (when we start the simulation) and v^PBH is the velocity of the PBH. The energy loss of the PBH during its passage through Earth is so small that the velocity of the PBH remains constant. Consequently, the source term f takes the form

$$\begin{equation} {\mathbf {f}}({\mathbf {x}},t) = \mbox{}{-} G\,M^\mathrm{PBH}\,\rho ({\mathbf {x}})\,\frac{{\mathbf {x}}-{\mathbf {x}}^\mathrm{PBH}(t)}{ \|{\mathbf {x}}-{\mathbf {x}}^\mathrm{PBH}(t)\|^3} . \end{equation} \tag{ 7 }$$

Given its mass M^PBH and trajectory (6), the gravitational force that the PBH exerts on particles at any location within Earth may be calculated.

On Earth's surface ∂Ω, traction-free boundary conditions must be satisfied,

$$\begin{equation} \hat{{\mathbf {n}}}\cdot {\mathbf {T}}= \mathbf {0} \quad ({\mathbf {x}}\in \partial \Omega, \forall t) , \end{equation} \tag{ 8 }$$

where $\hat{{\mathbf {n}}}$ denotes the unit outward normal along the boundary, i.e., on the free surface ∂Ω in this case. It is critical to note that Earth consists of several major shells: the crust and mantle that behave as solids on seismological timescales, the outer core that freely flows like a liquid, and the inner core that is a dense solid comprised of mostly iron. Hence, the boundaries between these layers, namely, the core–mantle boundary (CMB) and the inner-core boundary (ICB), involve solid–fluid interactions. The boundary conditions on fluid–solid internal boundaries are continuity of normal displacement and traction,

$$\begin{equation} [\hat{{\mathbf {n}}}\cdot \,\,{\mathbf {s}}]^{+}_{-}=0 \quad ({\mathbf {x}}\in \Sigma, \forall t) , \end{equation} \tag{ 9 }$$

$$\begin{equation} [\hat{{\mathbf {n}}}\cdot {\mathbf {T}}]^{+}_{-}=\mathbf {0} \quad ({\mathbf {x}}\in \Sigma, \forall t) , \end{equation} \tag{ 10 }$$

where [ g ]⁺ ₋ denotes the jump of any function g across an internal boundary Σ (either the CMB or ICB). Together with the initial conditions

$$\begin{equation} {\mathbf {s}}({\mathbf {x}},0)=\mathbf {0} , \quad \partial _t{\mathbf {s}}({\mathbf {x}},0)=\mathbf {0} , \end{equation} \tag{ 11 }$$

the system of governing equations presented in this section may now be solved numerically.

2.2. Weak Form

The wave equation (1) is in differential form, which is generally referred to as the strong formulation. Classical numerical methods, such as finite-difference and pseudo-spectral methods, usually involve a strong formulation. However, to account for solid–fluid boundary conditions properly, a weak formulation or integral form of the wave equation is preferable. The weak formulation is obtained by taking the dot product of Equation (1) with an arbitrary test vector w, and integrating over Earth's volume Ω:

$$\begin{equation} \int _\Omega \big(\rho \partial ^2_t {\mathbf {s}}-\mbox{\boldmath $\bf \nabla $}\cdot {\mathbf {T}}\big) \cdot {\mathbf {w}}{d}^3{\mathbf {x}}= \int _\Omega {\mathbf {f}}\cdot {\mathbf {w}}{d}^3{\mathbf {x}} . \end{equation} \tag{ 12 }$$

2.2.1. Crust and Mantle

The crust and mantle occupy the solid outer parts of Earth. Using the divergence theorem, Equation (12) in the crust and mantle may be simplified as

$$\begin{eqnarray} &&\int _{\Omega _\mathrm{M}} \rho \partial ^2_t {\mathbf {s}}\cdot {\mathbf {w}}{d}^3{\mathbf {x}}+ \int _{\Omega _\mathrm{M}} \mbox{\boldmath $\bf \nabla $}{\mathbf {w}}:{\mathbf {T}}{d}^3{\mathbf {x}}- \int _{\partial \Omega } \hat{{\mathbf {n}}}\cdot {\mathbf {T}}\cdot {\mathbf {w}}{d}^2{\mathbf {x}}\nonumber\\ &&\quad+ \int _{\mathrm{CMB}} \hat{{\mathbf {n}}}\cdot {\mathbf {T}}\cdot {\mathbf {w}}{d}^2{\mathbf {x}}=\int _{\Omega _\mathrm{M}} {\mathbf {f}}\cdot {\mathbf {w}}{d}^3{\mathbf {x}} , \end{eqnarray} \tag{ 13 }$$

where Ω_M denotes the volume of the crust and mantle, CMB the core–mantle boundary, and ∂Ω Earth's free surface.

The traction continuity condition on the CMB, i.e., Equation (10), reduces to

$$\begin{equation} \hat{{\mathbf {n}}}\cdot {\mathbf {T}}=-p \hat{{\mathbf {n}}}\quad \quad ({\mathbf {x}}\in \mathrm{CMB}, \forall t) , \end{equation} \tag{ 14 }$$

where p is the pressure field in the outer core.

Taking into consideration the free surface boundary condition (8) and traction continuity condition on the CMB (14), Equation (13) may be rewritten as

$$\begin{eqnarray} &&\int _{\Omega _\mathrm{M}} \rho \partial ^2_t {\mathbf {s}}\cdot {\mathbf {w}}{d}^3{\mathbf {x}}+ \int _{\Omega _\mathrm{M}} \mbox{\boldmath $\bf \nabla $}{\mathbf {w}}:{\mathbf {T}}{d}^3{\mathbf {x}}- \int _{\mathrm{CMB}} p \hat{{\mathbf {n}}}\cdot {\mathbf {w}}{d}^2{\mathbf {x}}\nonumber\\ &&\quad=\int _{\Omega _\mathrm{M}} {\mathbf {f}}\cdot {\mathbf {w}}{d}^3{\mathbf {x}} . \end{eqnarray} \tag{ 15 }$$

It is the pressure field on the CMB that facilitates the exchange of information from the outer core to the mantle.

2.2.2. Outer Core

The outer core is a liquid in which the governing equation simplifies significantly. The displacement field in the outer core s^fluid may be expressed in terms of a scalar potential χ as (Komatitsch & Tromp 2002a)

$$\begin{equation} {\mathbf {s}}^\mathrm{fluid}=\frac{1}{\rho } \mbox{\boldmath $\bf \nabla $}\chi . \end{equation} \tag{ 16 }$$

Substituting this definition into Equation (1), we may obtain a scalar governing equation for fluids

$$\begin{equation} \frac{1}{\kappa }\partial ^2_t\chi -\mbox{\boldmath $\bf \nabla $}\cdot \left(\frac{1}{\rho }\mbox{\boldmath $\bf \nabla $}\chi \right) =-\frac{1}{\kappa }\phi . \end{equation} \tag{ 17 }$$

The weak form of this equation is

$$\begin{eqnarray} &&\int _{\Omega _\mathrm{O}} \frac{w}{\kappa } \partial ^2_t\chi {d}^3{\mathbf {x}}+\int _{\Omega _\mathrm{O}} \frac{1}{\rho } \mbox{\boldmath $\bf \nabla $}w\cdot \mbox{\boldmath $\bf \nabla $}\chi {d}^3{\mathbf {x}}-\int _{\mathrm{CMB}} \frac{w}{\rho } \mbox{\boldmath $\bf \nabla $}\chi \cdot \hat{{\mathbf {n}}}{d}^2{\mathbf {x}}\nonumber\\ &&\quad+\int _{\mathrm{ICB}} \frac{w}{\rho } \mbox{\boldmath $\bf \nabla $}\chi \cdot \hat{{\mathbf {n}}}{d}^2{\mathbf {x}}=-\int _{\Omega _\mathrm{O}} \frac{w}{\kappa }\phi {d}^3{\mathbf {x}} , \end{eqnarray} \tag{ 18 }$$

where Ω_O denotes the volume of the outer core and w is an arbitrary scalar test function. The continuity condition for the normal component of the displacement field on the CMB and ICB (9) reduces to

$$\begin{equation} \frac{1}{\rho }\mbox{\boldmath $\bf \nabla $}\chi \cdot \hat{{\mathbf {n}}}=\hat{{\mathbf {n}}}\cdot {\mathbf {s}}^\mathrm{solid} \quad \quad ({\mathbf {x}}\in \mathrm{CMB\ or\ ICB}, \forall t) , \end{equation} \tag{ 19 }$$

where s^solid refers to the displacement on the solid side of the boundary. Now Equation (18) may be rewritten as

$$\begin{eqnarray} &&\int _{\Omega _\mathrm{O}} \frac{w}{\kappa } \partial ^2_t\chi {d}^3{\mathbf {x}}+\int _{\Omega _\mathrm{O}} \frac{1}{\rho } \mbox{\boldmath $\bf \nabla $}w\cdot \mbox{\boldmath $\bf \nabla $}\chi {d}^3{\mathbf {x}}-\int _{\mathrm{CMB}} w \hat{{\mathbf {n}}}\cdot {\mathbf {s}}{d}^2{\mathbf {x}}\nonumber\\ &&\quad+\int _{\mathrm{ICB}} w \hat{{\mathbf {n}}}\cdot {\mathbf {s}}{d}^2{\mathbf {x}}=-\int _{\Omega _\mathrm{O}} \frac{w}{\kappa }\phi {d}^3{\mathbf {x}} . \end{eqnarray} \tag{ 20 }$$

It is the normal component of the displacement field at the CMB and ICB that passes information from the mantle and the inner core to the outer core, respectively. By definition, the pressure field in the outer core is

$$\begin{equation} p=\mbox{}{-}\kappa \mbox{\boldmath $\bf \nabla $}\cdot {\mathbf {s}}=\mbox{}{-}\kappa \mbox{\boldmath $\bf \nabla $}\cdot \left(\frac{1}{\rho }\mbox{\boldmath $\bf \nabla $}\chi \right)=\mbox{}{-}\partial ^2_t\chi -\phi , \end{equation} \tag{ 21 }$$

where in the last equality we have used Equation (17).

2.2.3. Inner Core

The weak form in the solid inner core is similar to the crust and mantle, namely,

$$\begin{eqnarray} &&\int _{\Omega _\mathrm{I}} \rho \partial ^2_t {\mathbf {s}}\,{\bm \cdot}\, {\mathbf {w}}{d}^3{\mathbf {x}}\,+ \int _{\Omega _\mathrm{I}} \mbox{\boldmath $\bf \nabla $}{\mathbf {w}}:{\mathbf {T}}{d}^3{\mathbf {x}}\,+ \int _{\mathrm{ICB}} p \hat{{\mathbf {n}}}\,{\bm \cdot}\, {\mathbf {w}}{d}^2{\mathbf {x}}\nonumber\\ &&\qquad=\int _{\Omega _\mathrm{I}} {\mathbf {f}}\,{\bm \cdot}\, {\mathbf {w}}{d}^3{\mathbf {x}} , \end{eqnarray} \tag{ 22 }$$

where Ω_I denotes the volume of the inner core. Equations (15), (20), and (22) only hold under simple circumstances, because they do not take into account rotation or self-gravitation. General but more complex formulations may be found in Komatitsch & Tromp (2002a, 2002b); such complications are accommodated by the spectral-element solver SPECFEM3D_GLOBE.

2.3. Black Hole Singularity

One of the main differences between earthquake and BH simulations lies in the fact that we have to deal with different source terms in Equations (15), (20), and (22), namely,

$$\begin{equation*} \int _{\Omega _\mathrm{O}} \frac{\phi }{\kappa } w {d}^3{\mathbf {x}}\quad \mathrm{and} \quad \int _{\Omega _\mathrm{M}+\Omega _\mathrm{I}} {\mathbf {f}}\cdot {\mathbf {w}}{d}^3{\mathbf {x}} , \nonumber \end{equation*}$$

where ϕ and f are determined by Equations (5) and (7), respectively. Unlike earthquake simulations, where sources are associated with fault planes, here we encounter volumetric body forces. When a point within Earth is close to the location of the BH, the gravitational force and potential become singular. To avoid such unphysical singularities, we implement the source term as

$$\begin{equation} \widetilde{{\mathbf {f}}}({\mathbf {x}},t)\!=\! \left\{ \begin{array}{@{}@{}ll} \!{-}\,{{\it GM}^\mathrm{PBH}\,\rho ({\mathbf {x}})}\,\frac{{\mathbf {x}}-{\mathbf {x}}^\mathrm{PBH}(t)}{ \|{\mathbf {x}}-{\mathbf {x}}^\mathrm{PBH}(t)\|^3} , &\|{\mathbf {x}}-{\mathbf {x}}^\mathrm{PBH}(t)\| \ge R_\mathrm{crit} , \\ \ms \ms \!{-}\,{{\it GM}^\mathrm{PBH}\,\rho ({\mathbf {x}})}\,\frac{{\mathbf {x}}-{\mathbf {x}}^\mathrm{PBH}(t)}{R_\mathrm{crit}^3} , &\|{\mathbf {x}}-{\mathbf {x}}^\mathrm{PBH}(t)\| < R_\mathrm{crit} , \end{array} \right. \end{equation} \tag{ 23 }$$

$$\begin{equation} \widetilde{\phi }({\mathbf {x}},t)=\left\{ \begin{array}{@{}l} {-}\frac{G\,M^\mathrm{PBH}}{ \|{\mathbf {x}}-{\mathbf {x}}^\mathrm{PBH}(t)\|} , \\ \ms \ms \qquad\|{\mathbf {x}}-{\mathbf {x}}^\mathrm{PBH}(t)\| \ge R_\mathrm{crit} , \\ \ms \ms {-}\frac{3 G\,M^\mathrm{PBH}}{2 R_\mathrm{crit}} +\frac{G\,M^\mathrm{PBH}}{2 R_\mathrm{crit}^3} \|{\mathbf {x}}-{\mathbf {x}}^\mathrm{PBH}(t)\|^2 ,\\ \ms \ms \qquad \|{\mathbf {x}}-{\mathbf {x}}^\mathrm{PBH}(t)\|<R_\mathrm{crit} , \end{array}\right. \end{equation} \tag{ 24 }$$

where R_crit is a chosen critical length scale within which the amplitude of the gravitational force decreases linearly. It may be shown mathematically that $\widetilde{\phi }$ and $\widetilde{{\mathbf {f}}}$ capture all features of the original source terms ϕ and f, as long as the critical length is less than the size of the numerical grid. In our simulations, the size of an element is ∼30 km, while R_crit is around ∼6 m, i.e., five thousand times smaller than the element size. Convergence tests using different R_crit illustrate that further decreasing R_crit does not change the simulation results (Appendix A).

The last but equally important issue is the selection of the time step for the numerical simulations. As mentioned previously, the speed of the BH is typically a hundred times larger than the seismic wave speeds, and it takes only a few minutes for the PBH to traverse Earth. Therefore, we need a much smaller time step than usual to resolve the trajectory of the BH. Figure 9 explains why this is important. In practice, we choose a time step that is one thousandth of the original time step when the PBH is close to Earth.

The setup of our simulations is illustrated in Figure 1. Different trajectories as well as two distinct BH speeds, namely, 200 km s⁻¹ and 400 km s⁻¹, are considered. The simulations are carried out on a cubed-sphere mesh (Komatitsch & Tromp 2002a, 2002b) comprising 8.2 million elements. Given the wave-speed structure of Earth, the shortest period the mesh can resolve is around 15 s. We set the BH mass M^PBH = 10¹⁵ g, noting that amplitudes of the resultant seismic waves scale linearly with M^PBH.

Zoom In Zoom Out Reset image size

An instructive case is that of a PBH whose trajectory coincides with the center of Earth, i.e., d = 0 in Figure 1. We use the one-dimensional Preliminary Reference Earth Model (Dziewonski & Anderson 1981), aligning the PBH trajectory with Earth's rotation axis. Four snapshots from this simulation are shown in Figure 2, illustrating generation, propagation, transmission, and reflection of PBH-induced seismic waves. Seismograms from this simulation are displayed in Figure 3. Due to symmetry, east–west ground motions are very weak compared to the other two components, and are thus not shown. The north–south component shows not only distinctive arrivals due to cylindrical waves from the supersonic motion of the PBH, but also classical Rayleigh surface waves and reflected/transmitted core phases. The vertical component, besides containing signatures similar to the north–south component, registers unique waves that arrive at almost the same time everywhere on Earth's surface, distinct from those arising due to earthquakes. These waves are generated at the core–mantle boundary (CMB) and ICB due to free-slip interactions between the solid mantle and inner core and the liquid outer core.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We study four additional PBH trajectories, two of which pass through Earth's interior (d = 0.4 R_⊕ and d = 0.8 R_⊕) while the other two pass nearby (d = 1.2 R_⊕ and d = 1.6 R_⊕). A faster PBH traveling at a speed of 400 km s⁻¹ is also considered. To avoid redundancy, only d = 0.8 R_⊕ for the faster PBH is shown in Figure 4 (more images may be found in Appendix B). On the north–south component, direct arrivals are clearly visible due to the supersonic motion of the source as well as Rayleigh surface waves, which behave in the same manner as those shown in Figure 3. Distinctions from the simple example also emerge. For example, the entrance and exit points of the PBH change from the poles to latitudes around ±40°, as indicated by the focus of the surface waves. Waves arriving after the surface waves are coherent but more difficult to interpret. The vertical component is largely unchanged from the simple example. Besides other signatures found on the north–south component, linear features due to the spherical waves exist as well. For cases where the PBH does not hit our planet, we see no coherent signals on the horizontal components. But spherical waves are still visible on the vertical component, since they are gravitationally induced at the CMB and ICB. Simulations in three-dimensional Earth model S40RTS (Ritsema et al. 2011) preserve these key features, as illustrated in Appendix B.

Zoom In Zoom Out Reset image size

For all cases, ground motions are on the order of 10⁻⁶ m. For comparison, ground motions due to the Bolivia M_w = 8.1 earthquake on 1994 June 9 are on the order of 10⁻³ m. From an energy perspective, seismic energy injected by the 10¹⁵ g PBH is six orders in magnitude lower than that of the Bolivia earthquake, corresponding to an earthquake of magnitude M_w = 4.

Another distinctive feature of PBH-induced seismic waves is the appearance of unusual spheroidal normal modes (which are not excited by earthquakes) in spectra of the time series, as shown in Figure 5.

Zoom In Zoom Out Reset image size

The significance of the power spectra arises from normal-mode analysis (Dahlen & Tromp 1998). Any motions on Earth can be decomposed into a set of periodic oscillations—Earth's free oscillations or normal modes:

$$\begin{equation} {\mathbf {s}}({\mathbf {x}},t)=\mathrm{Re}\sum _k (i\omega _k)^{-1} {\mathbf {s}}_k({\mathbf {x}}) \int _{-\infty }^t C_k(t^{\prime }) \mathrm{exp}[i\omega _k(t-t^{\prime })] {{d}}t^{\prime } , \end{equation} \tag{ 25 }$$

where s_k denotes the kth mode, ω_k its eigenfrequency, and C_k its excitation.

4.1. Earthquake Excitation

For an earthquake, the excitation term C_k is controlled by the location of the earthquake and its mechanism:

$$\begin{equation} C_k(t)=\int _{\Sigma ^t} {\mathbf {m}}({\mathbf {x}},t) {\bm :} {\bm \varepsilon }_k^*({\mathbf {x}}) {{d}}^2{\mathbf {x}} , \end{equation} \tag{ 26 }$$

where Σ^t represents the fault plane at time t, m the moment density tensor, and ${\bm \varepsilon }_k$ the strain tensor associated with the kth mode

$$\begin{equation} {\bm \varepsilon }_k=\frac{1}{2}\left[\mbox{\boldmath $\bf \nabla $}{\mathbf {s}}_k+\left(\mbox{\boldmath $\bf \nabla $}{\mathbf {s}}_k\right)^T\right] . \end{equation} \tag{ 27 }$$

When the point source approximation applies, as is often the case in global seismology, the excitation term reduces to

$$\begin{equation} C_k(t)\approx {\mathbf {M}}(t) {\bm :} {\bm \varepsilon }_k^*({\mathbf {x}}_s) , \end{equation} \tag{ 28 }$$

where the strain tensor ${\bm \varepsilon }_k$ at the location of the point source x_s interacts with the moment tensor

$$\begin{equation} {\mathbf {M}}(t)=\int _{\Sigma ^t} {\mathbf {m}}({\mathbf {x}},t) {{d}}^2{\mathbf {x}} . \end{equation} \tag{ 29 }$$

Equation (28) determines not only which modes can be excited by an earthquake, but also how strongly excited these modes are. It is critical to note that earthquakes only occur in the crust and upper mantle. Therefore, modes that are not sensitive to the shallow parts of Earth cannot be excited by natural earthquakes; examples include the Slichter mode (₁S₁) and Stoneley modes trapped at the CMB and ICB.

4.2. PBH Excitation

A PBH triggers normal modes in a different way than earthquakes. The excitation follows the form of a distributed force

$$\begin{equation} C_k(t)=\int _\Omega {\mathbf {f}}({\mathbf {x}},t)\cdot {\mathbf {s}}_k^*({\mathbf {x}}) {{d}}^3{\mathbf {x}} , \end{equation} \tag{ 30 }$$

where f satisfies

$$\begin{equation} {\mathbf {f}}({\mathbf {x}},t)=G M^\mathrm{PBH} \rho ({\mathbf {x}}) \mbox{\boldmath $\bf \nabla $}\frac{1}{ \|{\mathbf {x}}-{\mathbf {x}}^\mathrm{PBH}(t)\|} . \end{equation} \tag{ 31 }$$

Substituting Equation (31) into (30), we obtain

$$\begin{eqnarray} &&C_k(t) =\int _\Omega G M^\mathrm{PBH} \rho ({\mathbf {x}}) {\mathbf {s}}_k^*({\mathbf {x}})\nonumber\\ &&\quad\cdot \mbox{\boldmath $\bf \nabla $}\left[\frac{1}{r^{\prime }}\sum _{l=0}^{\infty }\left(\frac{4\pi }{2l+1}\right)\left(\frac{r}{r^{\prime }}\right)^l\sum _{m=-l}^{l}Y_{lm}(\theta,\phi) Y_{lm}^*(\theta ^{\prime },\phi ^{\prime })\right] {{d}}^3{\mathbf {x}}\nonumber \\ &&\quad=\int _\Omega G M^\mathrm{PBH} \rho ({\mathbf {x}}) {\mathbf {s}}_k^*({\mathbf {x}})\cdot \left(\hat{{\mathbf {r}}}\partial _{r}+r^{-1}\mbox{\boldmath $\bf \nabla $}_1\right)\nonumber\\ &&\quad\times \left[\frac{1}{r^{\prime }}\sum _{l=0}^{\infty }\left(\frac{4\pi }{2l+1}\right)\left(\frac{r}{r^{\prime }}\right)^l\sum _{m=-l}^{l}Y_{lm}(\theta,\phi) Y_{lm}^*(\theta ^{\prime },\phi ^{\prime })\right] {{d}}^3{\mathbf {x}},\nonumber \\ \end{eqnarray} \tag{ 32 }$$

where r ≡ ||x||, r' ≡ ||x^PBH||, $\mbox{\boldmath $\bf \nabla $}_1\equiv \hat{\bm \theta } \partial _\theta +\hat{\bm \phi } (\sin \theta)^{-1}\partial _\phi$ denotes the surface gradient on the unit sphere, and Y_lm(θ, ϕ) are spherical harmonics (Edmonds 1960). Note that this expansion is valid only for r' > r. When the PBH is within Earth, r may be larger than r' for some points x. In that case, we need to interchange r and r', leading to a more complex result than the one considered here.

Mode s_k may be a spheroidal mode ^Ss_k or a toroidal mode ^Ts_k:

$$\begin{eqnarray} ^S{\mathbf {s}}_k({\mathbf {x}})&=&{}_nU_l(r) Y_{lm}(\theta,\phi)\, \hat{{\mathbf {r}}}+\frac{1}{\sqrt{l(l+1)}}\, {}_nV_l(r) \mbox{\boldmath $\bf \nabla $}_1Y_{lm}(\theta,\phi) ,\nonumber\\ \end{eqnarray} \tag{ 33 }$$

$$\begin{eqnarray} ^T{\mathbf {s}}_k({\mathbf {x}})&=&\frac{\mbox{}-1}{\sqrt{l(l+1)}}\, {}_nW_l(r) (\hat{{\mathbf {r}}}\times \mbox{\boldmath $\bf \nabla $}_1)Y_{lm}(\theta,\phi) , \end{eqnarray} \tag{ 34 }$$

where l, m, and n are the degree, order, and overtone number of the mode, respectively; _nU_l, _nV_l, and _nW_l are the corresponding radial eigenfunctions. Due to the relations

$$\begin{eqnarray} \hat{{\mathbf {r}}}\cdot \mbox{\boldmath $\bf \nabla $}_1&=&0 , \end{eqnarray} \tag{ 35 }$$

$$\begin{eqnarray} \hat{{\mathbf {r}}}\cdot \left(\hat{{\mathbf {r}}}\times \mbox{\boldmath $\bf \nabla $}_1\right)&=&0 , \end{eqnarray} \tag{ 36 }$$

$$\begin{eqnarray} \mbox{\boldmath $\bf \nabla $}_1\cdot \left(\hat{{\mathbf {r}}}\times \mbox{\boldmath $\bf \nabla $}_1\right)&=&0 , \end{eqnarray} \tag{ 37 }$$

it is obvious that no toroidal modes can be generated by a PBH, because the excitation term vanishes. For spheroidal modes, Equation (32) reduces to

$$\begin{eqnarray} C_k(t)&=& {}_nA_l(t) Y^*_{lm}(\theta ^{\prime }(t),\phi ^{\prime }(t)) , \end{eqnarray} \tag{ 38 }$$

$$\begin{eqnarray} {}_nA_l(t)&=&\frac{4\pi G M^\mathrm{PBH}}{(2l+1) \left(r^{\prime }(t)\right)^{l+1}} \int _0^{R_\oplus } \rho (r) r^{l+1}\nonumber\\ &&\times \left[l\, {}_nU_l(r) +\sqrt{l(l+1)}\, {}_nV_l(r)\right] {{d}}r , \end{eqnarray} \tag{ 39 }$$

where we have used the orthogonality of spherical harmonics. Note that the PBH location x^PBH varies with time, and therefore r', θ', and ϕ' are also functions of time. Equations (38) and (39) differ from Equation (28) in that the eigenfunctions everywhere along Earth's radius are involved. The implication is that even modes whose sensitivities are restricted to the deeper parts of Earth, i.e., modes that cannot be triggered by earthquakes, may be excited by a PBH.

Substituting Equations (38) and (39) into Equation (25), we obtain

$$\begin{eqnarray} {\mathbf {s}}({\mathbf {x}},t)&=&\mathrm{Re}\sum _k (i\omega _k)^{-1}{} ^S{\mathbf {s}}_k({\mathbf {x}})\nonumber\\ &&\times \int _{-\infty }^t {}_nA_l(t^{\prime }) Y^*_{lm}(\theta ^{\prime }(t^{\prime }),\phi ^{\prime }(t^{\prime }))\, \mathrm{exp}[i\omega _k(t-t^{\prime })] {{d}}t^{\prime }.\nonumber\\ \end{eqnarray} \tag{ 40 }$$

Note again that the summation over k actually denotes triple summations over l, m, and n. In order to evaluate Equation (40) efficiently, the summation over m is carried out analytically using the spherical-harmonic addition theorem

$$\begin{equation} \sum _{m=-l}^{l} Y^*_{lm}({\theta ^{\prime },\phi ^{\prime }}) Y_{lm}(\theta,\phi)=\left(\frac{2l+1}{4\pi }\right) P_l(\cos \Theta) , \end{equation} \tag{ 41 }$$

where cos Θ = cos θcos θ' + sin θsin θ'cos (ϕ − ϕ'). Finally, Equation (40) reduces to

$$\begin{eqnarray} {\mathbf {s}}({\mathbf {x}},t)&=&\left(\frac{2l+1}{4\pi }\right)\mathrm{Re}\sum _{l,n} (i {}_n\omega _l)^{-1}\nonumber\\ &&\times \bigg[{}_nU_l(r)\, \hat{{\mathbf {r}}}\int _{-\infty }^t {}_nA_l(t^{\prime }) P_{l}(\Theta (t^{\prime })) \mathrm{exp}(-i {}_n\omega _l t^{\prime }) {{d}}t^{\prime } \nonumber \\ &&+ \frac{1}{\sqrt{l(l+1)}}\, {}_nV_l(r) \int _{-\infty }^t {}_nA_l(t^{\prime }) \mbox{\boldmath $\bf \nabla $}_1 P_{l}(\Theta (t^{\prime }))\nonumber\\ &&\times\mathrm{exp}(-i {}_n\omega _l t^{\prime }) {{d}}t^{\prime } \bigg] \mathrm{exp}(i {}_n\omega _l t). \end{eqnarray} \tag{ 42 }$$

Equation (42) is difficult to implement for a moving PBH, however, it is straightforward to consider a stationary and transient PBH effect. In this case, instead of Equation (31), the force term is

$$\begin{equation} {\mathbf {f}}({\mathbf {x}},t)=G M^\mathrm{PBH} \rho ({\mathbf {x}}) \mbox{\boldmath $\bf \nabla $}\frac{1}{ \|{\mathbf {x}}-{\mathbf {x}}_0^\mathrm{PBH} \|} \delta ^{\prime }(t) , \end{equation} \tag{ 43 }$$

where δ'(t)—rather than δ(t)—is used to avoid inducing a permanent displacement. In this simpler case, where r', θ', and ϕ' are time independent, Equation (42) takes the following form:

$$\begin{eqnarray} {\mathbf {s}}({\mathbf {x}},t)&=&\left(-\frac{2l+1}{4\pi }\right)\mathrm{Re}\sum _{l,n} {}_nA_l \bigg[{}_nU_l(r) P_{l}\left(\cos \Theta \right) \hat{{\mathbf {r}}}\nonumber\\ &&+\frac{1}{\sqrt{l(l+1)}}\, {}_nV_l(r) \mbox{\boldmath $\bf \nabla $}_1 P_{l}\left(\cos \Theta \right)\bigg] \mathrm{exp}(i {}_n\omega _l t).\nonumber\\ \end{eqnarray} \tag{ 44 }$$

We shall use this configuration solely for benchmarking purposes, testing the generation of spherical waves by comparing the results of a normal mode calculation with those obtained based on a spectral-element method (SEM) calculation. For earthquake sources and impacts, the SEM has already been benchmarked extensively against normal modes (e.g., Komatitsch & Tromp 2002a; Meschede et al. 2011). Figure 6 shows a comparison between a normal-mode calculation and a spectral-element simulation for two stationary and transient PBHs, which are symmetrically placed with respect to Earth's center. Using two PBHs instead of one ensures that Earth's center remains stationary.

Zoom In Zoom Out Reset image size

In this section, we provide some order-of-magnitude estimates of event rates and non-seismic effects of a BH or QN moving through the Earth. The radius of a BH is

$$\begin{equation} r_\mathrm{BH}\approx 1.5 \times 10^{-15}\ {\rm m}\ \overline{M} , \end{equation} \tag{ 45 }$$

where $\overline{M}\equiv M/10^{15}$ g is its mass in units of 10¹⁵ g. Assuming that a QN is a uniform sphere with density $\rho _\mathrm{QN} = 10^{21}\ \rm {g} \ \rm {m}^{-3} \rho _0$, with ρ₀ of order unity (Witten 1984), it has a radius of

$$\begin{equation} r_\mathrm{QN}\approx 6 \times 10^{-3}\ \rm {m} \left(\,\frac{\overline{M}}{\rho _0}\,\right)^{1/3} . \end{equation} \tag{ 46 }$$

Restricting the discussion to masses less than ∼10⁻³ M_⊕, one can see that though the QN is significantly larger than the BH, as gravitational sources both can still be considered point particles to good approximation. Assuming that a fraction χ of the estimated dark matter density in the vicinity of the solar system—see for example Lisanti & Spergel (2011)—is in the form of compact objects of mass $\overline{M}$, the expected event rate is

$$\begin{equation} N\approx 4\times 10^{-7} \,{\rm yr}^{-1} \:\frac{\chi \ \bar{v}\ \bar{d}^{\,2}}{\overline{M}} , \end{equation} \tag{ 47 }$$

where $\bar{v}\equiv v/200$ km s⁻¹ is the velocity in units of 200 km s⁻¹ and $\bar{d}=b_{\max }/R_\oplus$ is the maximum distance, in units of R_⊕, out to which the passage of a compact object with impact parameter b_max will create the kind of impulsive gravitational response that gives rise to the unique seismic signatures discussed here. A rough estimate suggests $\bar{d}\approx 10\hbox{--}20$. These events are thus quite rare for masses that could produce observable seismic signatures. However, as we will show below, on the lower mass end of BHs that are seismically observable, the corresponding non-seismic signatures will be negligible, and seismology may offer a unique way to observe them. This is not the case for QNs, which, because of their larger radii, could produce strong atmospheric explosions.

The Bondi radius R_B for accretion of material onto a compact object is

$$\begin{equation} R_B \approx 3\times 10^{-9}\ {\rm m} \: \frac{\overline{M}}{\bar{v}^2} . \end{equation} \tag{ 48 }$$

For the masses considered here, this is negligibly small for BHs and is inside the radius of a QN, thus for QNs the relevant radius is simply r_QN. Following Burns et al. (1976b), the radius R_m out to which the gravitational impulse of the passing compact object could melt surface rock is

$$\begin{equation} R_m \approx 5\times 10^{-7}\ {\rm m}\: \frac{\overline{M}}{\bar{v}} , \end{equation} \tag{ 49 }$$

and following Labun et al. (2011), the transverse distance L out to which rock will be shattered due to the gravitational tidal field is

$$\begin{equation} L \approx 10^{-5}\ {\rm m} \ \overline{M} . \end{equation} \tag{ 50 }$$

Again, both R_m and L are within the radius of the QN and negligible for smaller mass BHs. Applying the results of Ruderman & Spiegel (1971), we may estimate the energy that will be deposited into the atmosphere by the passage of the compact object as it passes through. For BHs, the main effect is due to the impulsive gravitational acceleration of air molecules, giving

$$\begin{equation} \Delta E^\mathrm{atm}_\mathrm{BH} \approx 0.3 \ {\rm J} \ \left(\,\frac{\overline{M}}{\bar{v}}\,\right)^2 . \end{equation} \tag{ 51 }$$

To create an explosion such as the Tunguska event of ∼10¹⁵–10¹⁷ J, for example, would require a BH of mass $10^{8}\hbox{--}10^{9}\, \overline{M}$. QNs cause a similar gravitational impulse, however their larger surfaces experience a drag that can deposit significantly larger amounts of energy compared to BHs for the smaller mass objects:

$$\begin{equation} \Delta E^\mathrm{atm}_\mathrm{QN} \approx 4 \times 10^{10} \ {\rm J} \ \bar{v}^2 \left(\,\frac{\overline{M}}{\rho _0}\,\right)^{2/3} . \end{equation} \tag{ 52 }$$

Again, to create Tunguska-like explosions requires massive QNs, however, in contrast to BHs, smaller QNs can still produce sizable explosions, e.g., an $\overline{M}=1$ QN imparts the energy equivalent of ∼10 tons of TNT.

The event rates above are exceedingly low for massive ($\overline{M}\gg 1$) compact objects, though as discussed in Labun et al. (2011), direct impacts may leave features distinguishable from asteroid impacts, and hence the effective baseline for observation can be extended by looking at geological features on Earth, Moon, Mercury, and Mars. Also, as BH and QN detectors the gas giants have larger collecting areas, and may be useful in constraining compact object dark matter candidates in regimes not readily accessible by other observations.

Based on numerical simulations, we find that PBH transits through Earth generate unique and significant seismic signals. Our theoretical results may be used to develop a strategy for detecting PBHs. Unlike earlier suggestions (Anderson et al. 2003), we propose to use spherical rather than cylindrical waves as the key detection mechanism. There are several advantages to this approach. First, detection is more reliable because the simultaneous global arrivals induced by slip at the CMB and ICB are unique and cannot be produced by earthquakes. In contrast, cylindrical-wave arrivals may resemble certain earthquakes, especially when the number of available stations is small. Such a misinterpretation, involving an unreported earthquake, occurred previously (Selby et al. 2004). Second, in order to detect cylindrical waves, travel time measurements are involved, which presents a by-no-means trivial task. When the trajectory is unknown, it is difficult to align seismic stations, which further increases the challenge of picking the relevant arrivals. In contrast, the linear travel-time feature associated with simultaneous global arrivals can be simply revealed by comparing seismograms side by side in a so-called record section. Because the signals have the same arrival time, how we align the related seismograms is less critical. Third, due to reverberations of the CMB- and ICB-induced spherical waves, unique free-oscillation spectra are expected. These spectra involve a frequency-spacing between modes excited by a PBH governed by the pole-to-pole compressional-wave travel time, and the excitation of modes not affected by earthquakes, which are restricted to Earth's crust and upper mantle. Stacking of global seismic spectra may be required to bring out these free-oscillation characteristics. Finally, the detector we propose applies equally to scenarios where the PBH does not collide with Earth, but passes sufficiently close to generate detectable spherical waves. Of course cylindrical waves are still important, because fitting the arrival times of Mach waves provides constraints on the PBH trajectory which cannot be obtained by considering spherical waves only.

Seismograms recorded by the global seismographic network are routinely monitored to detect landslides and glacial earthquakes (Ekstrom et al. 2003), and it may be feasible to include detection criteria for PBHs in this monitoring process.

We study the effects of PBHs passing through or nearby Earth using numerical simulations based on a spectral-element method. The fast motion of the PBH causes rapidly varying classical gravitational forces that are carefully numerically accommodated. We also satisfactorily resolve numerical issues due to the singularity in the gravitational near-field of the PBH. We predict three characteristic features. The first feature involves cylindrical waves due to a supersonic source, generating Mach waves which are readily understood and have been studied previously (Anderson et al. 2003). The second feature involves spherical waves generated by free-slip at solid–fluid boundaries, namely, the core–mantle and inner-core boundaries. These spherical waves present a unique signature at global seismographic stations that cannot be produced by earthquakes. Third, free-oscillation spectra are predicted to exhibit evenly spaced peaks, with a spacing determined by the pole-to-pole compressional-wave travel time. Modes not excited by earthquakes, such as ₂S₂, should be detectable in spectra induced by PBHs. Probably the most likely scenario is the passage of a PBH nearby Earth. In such a near-miss scenario, the cylindrical-wave features vanish, but the spherical waves remain. Therefore, we may still use global seismograms and their spectra to detect a PBH, although it would be difficult to determine its trajectory. We also note that spherical waves could in principle be excited by asteroids along trajectories proximal to Earth, potentially containing new seismic information about the interior.

In this article, we only consider gravitational interactions between a PBH and Earth. Potential PBH impact phenomena, which are still debated, are not considered. Related seismic effects may be readily taken into consideration, as exemplified by the meteorite impact study of Meschede et al. (2011). Such effects are of course irrelevant for near misses.

As discussed in Section 2.3, we introduce a critical length scale R_crit in order to eliminate unphysical singularities. This length scale determines the maximum force we allow in the numerical simulations, as well as the behavior of the source time function. On one hand, amplitudes of simulated seismograms will be too small if R_crit is too large—the forces the PBH exerts on Earth would be severely underestimated. On the other hand, the smaller R_crit, the more the effective source time function looks like a delta function, which contains high-frequency information that cannot be resolved by our finite-size elements. The second issue is less problematic, since high-frequency noise can be safely removed by applying low-pass or band-pass filters in post processing. For these reasons, R_crit should be chosen to be as small as possible, as long as the maximum force does not exceed the limits of computer representation. However, we will show that when the critical length scale is far less than the size of the elements, differences in the weak form implementation vanish.

Figure 7 illustrates an element which contains the PBH at an arbitrary moment.

Zoom In Zoom Out Reset image size

The source term involves a volume integral of the gravitational force over element Ω_e, that is,

$$\begin{eqnarray*} \int _{\Omega _e} {\mathbf {f}}\cdot {\mathbf {w}}{{d}}^3{\mathbf {x}}&=&\int _{\Omega _e} \widetilde{{\mathbf {f}}}\cdot {\mathbf {w}}{{d}}^3{\mathbf {x}}+ \int _{\Omega _\mathrm{PBH}}({\mathbf {f}}-\widetilde{{\mathbf {f}}})\cdot {\mathbf {w}}{{d}}^3{\mathbf {x}}\nonumber \\ &=& \int _{\Omega _e} \widetilde{{\mathbf {f}}}\cdot {\mathbf {w}}{{d}}^3{\mathbf {x}}+\int _{\Omega _\mathrm{PBH}} G M^\mathrm{PBH} \rho ({\mathbf {x}})\nonumber\\ &&\times \left[-\frac{{\mathbf {x}}-{\mathbf {x}}^\mathrm{PBH}(t)}{ \|{\mathbf {x}}-{\mathbf {x}}^\mathrm{PBH}(t)\|^3}+\frac{{\mathbf {x}}-{\mathbf {x}}^\mathrm{PBH}(t)}{R_\mathrm{crit}^3}\right]\nonumber\\ && \cdot {\mathbf {w}}({\mathbf {x}}) {{d}}^3{\mathbf {x}}\nonumber \\ &\approx & \int _{\Omega _e} \widetilde{{\mathbf {f}}}\cdot {\mathbf {w}}{{d}}^3{\mathbf {x}}+\int _{\Omega _\mathrm{PBH}} G M^\mathrm{PBH} \rho ({\mathbf {x}}^\mathrm{PBH}(t))\nonumber\\ &&\times \left[-\frac{{\mathbf {x}}-{\mathbf {x}}^\mathrm{PBH}(t)}{ \|{\mathbf {x}}-{\mathbf {x}}^\mathrm{PBH}(t)\|^3}+\frac{{\mathbf {x}}-{\mathbf {x}}^\mathrm{PBH}(t)}{R_\mathrm{crit}^3}\right] \nonumber\\ &&\cdot {\mathbf {w}}({\mathbf {x}}^\mathrm{PBH}(t)) {{d}}^3{\mathbf {x}}\nonumber \\ &=& \int _{\Omega _e} \widetilde{{\mathbf {f}}}\cdot {\mathbf {w}}{{d}}^3{\mathbf {x}}+G M^\mathrm{PBH} \rho ({\mathbf {x}}^\mathrm{PBH}(t)) {\mathbf {w}}({\mathbf {x}}^\mathrm{PBH}(t))\nonumber\\ &&\cdot \int _{\Omega _\mathrm{PBH}} \left[-\frac{{\mathbf {x}}-{\mathbf {x}}^\mathrm{PBH}(t)}{ \|{\mathbf {x}}-{\mathbf {x}}^\mathrm{PBH}(t)\|^3}+\frac{{\mathbf {x}}-{\mathbf {x}}^\mathrm{PBH}(t)}{R_\mathrm{crit}^3}\right] {{d}}^3{\mathbf {x}}, \nonumber \\ \end{eqnarray*}$$

where we have used the fact that when R_crit ≪ Δh, Ω_PBH can be treated as a single point compared to element Ω_e, i.e., ρ(x) ≈ ρ(x^PBH(t)) and w(x) ≈ w(x^PBH(t)) for any x ∈ Ω_PBH at time t.

Define r ≡ x − x^PBH(t), and hence r ≡ ||r|| = ||x − x^PBH(t)||. Then the volume integral reduces to

$$\begin{eqnarray} \int _{\Omega _e} {\mathbf {f}}\cdot {\mathbf {w}}{{d}}^3{\mathbf {x}}&=& \int _{\Omega _e} \widetilde{{\mathbf {f}}}\cdot {\mathbf {w}}{{d}}^3{\mathbf {x}}+G M^\mathrm{PBH} \rho ({\mathbf {x}}^\mathrm{PBH}(t)) {\mathbf {w}}({\mathbf {x}}^\mathrm{PBH}(t))\nonumber\\ &&\cdot \int _{\Omega _\mathrm{PBH}} \left(-\frac{{\mathbf {r}}}{r^3}+\frac{{\mathbf {r}}}{R_\mathrm{crit}^3}\right) {{d}}^3{\mathbf {x}}\nonumber \\ &=& \int _{\Omega _e} \widetilde{{\mathbf {f}}}\cdot {\mathbf {w}}{{d}}^3{\mathbf {x}}, \end{eqnarray} \tag{ A1 }$$

because the integral over Ω_PBH is non-singular and simply vanishes. This can be easily understood from a physical perspective: within a small spherical volume around the PBH, the net gravitational force from the PBH is zero because gravity is isotropic in the tiny sphere. Equation (A1) is important, because it provides a way to calculate the source term involving a singular force f using a non-singular force $\widetilde{{\mathbf {f}}}$, which is easily implemented in the spectral-element method.

We design an experiment to test this result. Three complementary simulations are carried out, each with a different R_crit: the first one with R_crit = 10⁻⁶ R_⊕, the second one with R_crit = 10⁻⁷ R_⊕, and the last one with R_crit = 10⁻⁸ R_⊕, where R_⊕ = 6371 km is the radius of Earth. One seismogram from each of the three simulations is shown in Figure 8. The three simulations give the same results, thereby demonstrating that the critical length scale R_crit = 10⁻⁶ R_⊕ used in our simulations is sufficiently small for our theoretical derivation to be valid. Using $\widetilde{{\mathbf {f}}}$ instead of f does not only removes singularities, but also gives reliable seismograms.

Zoom In Zoom Out Reset image size

Another issue involves capturing the widely differing timescales associated with the motion of the PBH and the periods of the induced seismic waves. Figure 9 explains why this is important.

Zoom In Zoom Out Reset image size

A few additional simulation results are shown in Figures 10–12.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size