Pulsar Timing Probes of Primordial Black Holes and Subhalos

TL;DR

Pulsar timing arrays can detect or constrain primordial black holes and subhalos across a wide mass range, offering a new method to probe dark matter subcomponents and their properties.

Contribution

This study demonstrates the sensitivity of pulsar timing arrays, especially the Square Kilometer Array, to primordial black holes and subhalos, providing constraints on their abundance and density.

Findings

SKA can constrain primordial black holes as a dark matter component.

PTAs can detect less dense objects than lensing techniques.

Constraints on subhalo concentration parameters are achievable.

Abstract

Pulsars act as accurate clocks, sensitive to gravitational redshift and acceleration induced by transiting clumps of matter. We study the sensitivity of pulsar timing arrays (PTAs) to single transiting compact objects, focusing on primordial black holes and compact subhalos in the mass range from $10^{-12} M _{\odot}$ to well above $100~M_\odot$. We find that the Square Kilometer Array can constrain such objects to be a subdominant component of the dark matter over this entire mass range, with sensitivity to a dark matter sub-component reaching the sub-percent level over significant parts of this range. We also find that PTAs offer an opportunity to probe substantially less dense objects than lensing because of the large effective radius over which such objects can be observed, and we quantify the subhalo concentration parameters which can be constrained.

Figures

Tables

Table 1. Table 1 : Summary of timing parameters that characterize pulsar timing capabilities. SKA projections are taken from  Rosado  et al. ( 2015 ) . Current constraints are compiled from various sources  Desvignes  et al. ( 2016 ); Verbiest  et al. ( 2016 ); Arzoumanian  et al. ( 2018 ); Reardon  et al. ( 2016 ) as described in the text.

	$T$ [yr]	$t_{\mathrm{RMS}}$ [ns]	$\Delta t$ [wk]	$d$ [kpc]	$N_P$
Current	$5-30$	$50-{10}^4$	$1-4$	$0.5-5$	73
SKA	20	50	2	5	200
Optimistic	20	25	1	10	1000

Equations

$$\varphi (t)={\varphi }_0+\nu t+\frac{1}{2}\dot{\nu }{t}^2+\frac{1}{6}\ddot{\nu }{t}^3+...$$

$$\delta \varphi \equiv \int dt\delta \nu (t)$$

$$t_{\mathrm{RMS}}\equiv \sqrt{\frac{1}{N}\sum _n(t_n^{\text{data}}-t_n^{\text{fit}}{)}^2},$$

$${\left(\frac{\delta \nu }{\nu }\right)}_D$$

$${\left(\frac{\delta \nu }{\nu }\right)}_S$$

$${\left(\frac{\delta \nu }{\nu }\right)}_D$$

$${\left(\frac{\delta \nu }{\nu }\right)}_S$$

$${\left(\frac{\delta \nu }{\nu }\right)}_D$$

$${\left(\frac{\delta \nu }{\nu }\right)}_S$$

$$\text{Pr}(\text{SNR}<x{)}^{N_P}=\text{Erf}{\left(\frac{x}{\sqrt{2}}\right)}^{N_P}=0.95,$$

$${\sigma }_{\ddot{\nu }/\nu }=6\sqrt{\frac{2800\Delta t}{T}}\frac{t_{\text{RMS}}}{T^3}.$$

$$\mathrm{SNR}\equiv \frac{|\ddot{\nu }/\nu |}{{\sigma }_{\ddot{\nu }/\nu }}>4.$$

$$\frac{\ddot{\nu }}{\nu }$$

$$=\frac{GMv}{r_0^3}\left(\hat{\mathbf{v}}+3\frac{v{t}_0}{r_0}{\hat{\mathbf{r}}}_0\right)\cdot \hat{\mathbf{d}},$$

$$r_{\text{min}}$$

$$\sim \frac{10^{-2}\text{ pc}}{{\left(N_{P}f\right)}^{\frac{1}{3}}}{\left(\frac{M}{10^{-6}{M}_{\odot }}\right)}^{\frac{1}{3}}.$$

$$\frac{\ddot{\nu }}{\nu }\simeq \frac{2GMv}{r_{\text{min}}^3}\sim 3\times 10^{-32}\left(\frac{N_{P}f}{200}\right){\text{ Hz}}^2.$$

$$\frac{\ddot{\nu }}{\nu }=\frac{8GM{v}_{\times }^3}{r_{\times }^3}{\left(\frac{t_0{v}_{\times }}{r_{\times }}\right)}^3\left(1-\frac{3{\tau }^2}{t_0^2}\right).$$

$$r_{\times ,\text{min}}$$

$$\sim \frac{0.2\text{ pc}}{\sqrt{N_{P}f}}{\left(\frac{M}{M_{\odot }}\right)}^{\frac{1}{2}}{\left(\frac{\text{kpc}}{d}\right)}^{\frac{1}{2}}.$$

$$\frac{\ddot{\nu }}{\nu }$$

$$\sim 8\times 10^{-33}{\left(\frac{N_{P}f}{200}\right)}^{\frac{3}{2}}{\left(\frac{M_{\odot }}{M}\right)}^{\frac{1}{2}}{\left(\frac{d}{\text{kpc}}\right)}^{\frac{3}{2}}{\text{ Hz}}^2.$$

$$⟨\delta t(t_1)\delta t(t_2)⟩$$

$$⟨\widetilde{\delta t}(f)\widetilde{\delta t}(f^{\prime })⟩$$

$$⟨\widetilde{\dot{\delta t}}(f)\widetilde{\dot{\delta t}}(f^{\prime })⟩$$

$${\text{SNR}}^2=4{\int }_0^{\infty }df\frac{|\tilde{h}(f){|}^2}{S_{\dot{\delta t}}(f)},$$

$$\mathrm{SNR}\simeq \left(\frac{GM}{\tau v^2}\sqrt{\frac{T^3}{12t_{\mathrm{RMS}}^2\Delta t}}\right)\hat{\mathbf{b}}\cdot \hat{\mathbf{d}}.$$

$${\tau }_{\min }$$

$$\sim \frac{20\mathrm{yr}}{\sqrt{N_{P}f}}{\left(\frac{M}{10^{-9}{M}_{\odot }}\right)}^{\frac{1}{2}}{\left(\frac{20\text{ yr}}{T}\right)}^{\frac{1}{2}}.$$

$$M\lesssim 4\times 10^{-8}{M}_{\odot }\left(\frac{N_{P}f}{200}\right){\left(\frac{T}{20\mathrm{yr}}\right)}^3$$

$$\mathrm{SNR}$$

$${\tau }_{\text{min}}$$

$$\sim \frac{20\mathrm{yr}}{N_{P}f}\left(\frac{M}{10^{-4}{M}_{\odot }}\right)\left(\frac{20\text{ yr}}{T}\right)\left(\frac{\text{kpc}}{d}\right).$$

$$M\lesssim 10^{-2}{M}_{\odot }\left(\frac{N_{P}f}{200}\right){\left(\frac{T}{20\mathrm{yr}}\right)}^2\left(\frac{d}{\mathrm{kpc}}\right).$$

$$f_{\text{D, dyn}}^L$$

Full text

I Introduction

Gravity, as the only known coupling of the Dark Matter (DM) with ordinary matter, currently provides the sole window on the nature of the DM. It is via gravitational probes that we have inferred the DM abundance and its behavior as a cold, collisionless fluid. While these macroscopic properties of dark matter have, so far, given us relatively limited information on the theory of DM, it is still possible, despite the weakness of the gravitational couplings, to unveil much more about its properties through those interactions. Every dark matter model predicts small scale structure inside galaxies, and the type of structure, including their density profile and radius, gives information on the DM’s cosmological history, couplings to the visible sector, and to itself. The most popular DM candidate, the Weakly Interacting Massive Particle (WIMP), predicts a power spectrum with weak clustering on small scales, which is difficult to observe experimentally. On the other hand, there are many models that predict an abundance of gravitationally bound structure (which we shall refer to as compact objects) on small scales, and the details of those structures can point to specific models.

There are many models with dynamics in the early universe that give rise to enhanced matter power on small scales. For example, phase transitions (as for axion (or scalar) miniclusters Hogan and Rees (1988); Kolb and Tkachev (1993, 1994a, 1994b, 1996); Zurek et al. (2007); Hardy (2017); Fairbairn et al. (2018); Enander et al. (2017)), early matter domination Erickcek and Sigurdson (2011); Barenboim and Rasero (2014); Fan et al. (2014); Dror et al. (2018), and vector bosons produced during inflation Graham et al. (2016), can all generate enhanced small scale structure. Perhaps the simplest DM substructure consists of primordial black holes (PBH) Hawking (1971); Carr and Hawking (1974); Carr (1975), which is (largely) characterized by a single parameter, the PBH mass, $M$. They arise in many theories with inflationary dynamics producing large density perturbations on small scales, e.g. Carr et al. (1994); Kawasaki et al. (1998); Clesse and García-Bellido (2015); Orlofsky et al. (2017); Inomata et al. (2017). They are one extreme of the generic density profile and hence serve as a benchmark for gravitational probes of small scale structure. The most direct astrophysical probes of DM substructure are currently derived from a wide variety of lensing experiments. These experiments together constrain monochromatic compact object abundance in the mass range $10^{-10}-10~{}M_{\odot}$ to be a subdominant component of DM. Current searches include lensing of the Large Magellanic Cloud (e.g., MACHO Allsman et al. (2001), EROS Tisserand et al. (2007), OGLE Wyrzykowski et al. (2011)), Andromeda (e.g., SUBARU Niikura et al. (2017)) and stars in the local neighborhood (e.g., from KEPLER Griest et al. (2011, 2014), and Gaia Van Tilburg et al. (2018)). Lensing of distant supernovae Zumalacarregui and Seljak (2017) and quasars Wilkinson et al. (2001) have also been considered.

There are also a variety of constraints specific to PBHs. From $10^{-13}~{}M_{\odot}$ down to $10^{-15}~{}M_{\odot}$ (below $10^{-15}~{}M_{\odot}$ PBHs evaporate in a time shorter than the age of the Universe), the existence of white dwarfs Graham et al. (2015) currently constrains PBHs, and femtolensing may do so in the future Katz et al. (2018) (see Gould (1992); Nemiroff and Gould (1995) for earlier work). Above $100~{}M_{\odot}$ a variety of constraints from structure formation and Planck fairly severely constrain the PBH abundance, see  Carr et al. (2010) for a review. There remains a controversial window between $\sim 1-100~{}M_{\odot}$ where PBHs Bird et al. (2016); Clesse and García-Bellido (2016); Sasaki et al. (2018), and not astrophysical black holes or neutron stars, could give rise to the events in LIGO Abbott . The event rate due to PBHs seems consistent with the LIGO event rate, though it has been pointed out that a myriad of other constraints apply to PBHs in this mass window, including disruption of compact stellar systems such as Eriadnus II Brandt (2016) and Segue I Koushiappas and Loeb (2017).

From this discussion we see that there are three regimes where current constraints on compact objects are lacking or limited, and hence where any prospect for setting constraints is particularly intriguing. First, in the mass window between $\sim 1-100~{}M_{\odot}$ relevant for LIGO signals, where constraints exist but there are substantial astrophysical uncertainties. Second, in the mass window below $10^{-11}~{}M_{\odot}$, where lensing constraints currently do not apply. Finally, in the regime where lensing constraints are significant for monochromatic PBHs ($10^{-11}-10~{}M_{\odot}$), but also suffer from astrophysical uncertainties. The constraints in this regime will not have reach to more diffuse subhalos, owing to the requirement that the radius of the object be smaller than the Einstein radius.

In this paper, we explore an astrophysically clean measurement of DM compact objects via pulsar timing measurements across the entire mass window $10^{-12}-100~{}M_{\odot}$, by combining several different gravitational redshift and timing effects in measurements of pulsar periods. Pulsars with millisecond periods, observed over time scales of decades, are known to be remarkably stable clocks. While their periods fluctuate over short times, these fluctuations do not substantially accumulate. In practice one can define a pulse phase of the signal,

[TABLE]

where $\nu$ is the frequency and $\dot{\nu}$, $\ddot{\nu}$ are its first and second derivatives. The most stable pulsars have frequencies of ${\cal O}({\rm kHz})$ and a spin-down rate of the pulsar, $\dot{\nu}/\nu$, ranging from roughly $10^{-23}-10^{-20}\text{ Hz}$, both of which can be fit from the data. Empirically, it is found that $\ddot{\nu}/\nu$ can be below $10^{-31}\text{ Hz}^{2}$ Liu et al. (2018) and is typically not included in fits to the data, allowing one to place upper bounds on processes that would produce a non-negligible $\ddot{\nu}$. Furthermore, any process which induces a modification of the phase,

[TABLE]

can be constrained using pulsar timing measurements.

The quality of pulsar timing data is determined by three parameters. The first parameter is the root-mean-square (RMS) timing residual, $t_{\text{RMS}}$. This is determined after finding the frequency, $\nu^{\text{fit}}$, and its derivative, $\dot{\nu}^{\text{fit}}$, which minimizes the residual between the timing data, $t_{n}^{\text{data}}$, and the timing model, $t_{n}$, where $t_{n}$ is found via the relation $2\pi n=\phi(t_{n})$ from Eq. (1). This gives

[TABLE]

where $N$ is the number of data points, and $t_{n}^{\text{fit}}$ is $t_{n}$ with $\nu=\nu^{\text{fit}}$, $\dot{\nu}=\dot{\nu}^{\text{fit}}$ and all higher order terms dropped. The minimized residual is typically $t_{\text{RMS}}\sim\mu\text{sec}$. The other two parameters are the observation time of the pulsar, $T\sim 10$ years, and the time between measurements, $\Delta t\sim 2$ weeks (also known as the cadence). Clearly the pulsars with the most power to constrain substructure are those with smaller RMS noise, longer observation times, and shorter cadence.

Pulsar timing data can probe DM compact objects since a transit near the timing system will give rise to a change in the observed frequency of the pulsar. We consider changes in the observed frequency of the pulsar due to two effects. First, there can be a gravitational time delay due to a changing gravitational potential affecting the photon geodesic as it moves along the line of sight – this is known as a Shapiro time delay, and was proposed as a probe of dark matter in Siegel et al. (2007). Second, the presence of compact objects can lead to an acceleration of the Earth or pulsar, also changing the observed pulsar period – this is the Doppler effect, and was proposed as a signal of dark matter in Seto and Cooray (2007). These accelerations are optimal for studying smaller masses and are typically more sensitive than Shapiro delays, though in some parameter space, as we will explore in detail, Shapiro delays can be more sensitive due to the long baseline.

The signal from a transiting compact object will look different depending on the relevant timescale, $\tau$, associated with the motion of the compact objects (here we use this variable schematically but give it an explicit, mass-dependent meaning in later sections). If we denote the observation time of a pulsar as $T$, then dynamic signals correspond to $\tau\ll T$, and will appear as blips in the pulsar timing data (analogous to glitches which have been observed in millisecond pulsar data Cognard and Backer (2004); McKee et al. (2016)). Static signals, with $\tau\gg T$, will not be observable as blips but instead as a non-negligible contribution to the second derivative of the frequency, $\ddot{\nu}$.

The idea of using pulsar timing to probe dark matter substructure has a long history. The static contribution of the Shapiro time delay was suggested as a probe of PBHs in Clark et al. (2016); Schutz and Liu (2017), while searches for dynamic signals were considered for single events in Siegel et al. (2007); Seto and Cooray (2007); Kashiyama and Seto (2012), and multiple events in Baghram et al. (2011). None of these analyses, however, considered how the signals were related to each other in the relevant regime of validity. Our results extend, and differ from, previous results as follows. First, we carry out the first analysis to correctly consider all forms of timing signatures, in the dynamic and static limit, and for both Doppler and Shapiro effects. We comment on the interplay between these four signals and their complementary sensitivity in different mass ranges. The comparative analysis has important implications for signals; for example, in contrast to previous work, we find that the Doppler signal dominates in the static limit, substantially modifying the derived constraint. Second, we perform the first study of the single event ‘blip’ signal shapes and compute these shapes in three dimensions; this extends and improves on the previous limits derived in Seto and Cooray (2007); Kashiyama and Seto (2012); Kashiyama and Oguri (2018). Third, we perform projections for current and future pulsar timing experiments in all of the signal regimes, correctly incorporating the impact of the measurement cadence on the constraint for the first time. Lastly, we study the impact of the size of compact objects, parameterized in terms of the profile, on the constraints derived. Note that we do not consider a multi-event (or statistical) signal, as studied in Baghram et al. (2011). While we expect that such an analysis will extend the reach at the low mass end (below ${\cal O}(10^{-9}\,M_{\odot})$ for Doppler signals and below ${\cal O}(10^{-4}\,M_{\odot})$ for Shapiro signals), due to the more complicated nature of the signal, we reserve study for future work Dror et al. (2019).

The outline of this paper is as follows. In Sec. II we describe static and dynamic signatures of transiting compact objects, for both Doppler and Shapiro effects, being careful to delineate the dividing line between the regimes. Next, in Sec. III, we detail the size of the signals expected in the dynamic and static regimes for both Doppler and Shapiro signals. Then we present the analytic and numerical results in Sec. IV, projecting constraints on the fraction of DM in PBHs (or PBH-like subhalos) which can be probed using pulsar timing. These results are extended to more diffuse subhalos in Sec. V, where we show that PTAs have sensitivity to much more extended objects than lensing searches. Finally, in Sec. VI, we summarize our results and suggest ways in which the analysis can be extended.

II Pulsar Timing Signatures from Doppler and Shapiro Effects

Transiting compact objects give rise to two different effects in the time of arrival of pulses from pulsars. The first, the Doppler effect, arises from an acceleration of the Earth or the pulsar. The Shapiro effect, on the other hand, is a gravitational redshift effect along the photon geodesic. Both of these effects cause the photon arrival time to be shifted from the unperturbed propagation value. The constant terms inside of these time shifts are unobservable as they can be absorbed by a redefinition of the unperturbed travel time. We thus consider time-dependent changes which generate a shift in the pulsar frequency, $\delta\nu$. For the Doppler and Shapiro signals, we have, 111Here we assume a weak field approximation, $\Phi\ll 1$, a slowly varying potential during the interaction time scale ($\Phi(r+vr)\simeq\Phi(r)$), where $r$ is the distance of closest approach, and large orbit eccentricity.

[TABLE]

where $\Phi$ is the gravitational potential due to the compact object and $\mathbf{v}$ is its velocity, while $\mathbf{\hat{d}}$ is the direction from the Earth to the pulsar and $z$ parameterizes the path the light takes from the pulsar to the Earth. These expressions can be simplified by assuming the compact object is a PBH of mass $M$,

[TABLE]

where $\mathbf{r}$ is the position of the compact object relative to the pulsar and $\times$ subscript denotes crossing with $\mathbf{\hat{d}}$, ${\bf r}_{\times}\equiv{\bf r}\times{\bf\hat{d}}$. Physically, the Doppler delay derives from integrating over the gravitational field from the compact object and taking the component of the pulsar (Earth) acceleration towards the Earth (pulsar), while the Shapiro delay depends only on components of the position and velocity of the compact object in the direction perpendicular to $\mathbf{\hat{d}}$, as only this gives a time dependent shift to the metric affecting the photons.

As shown in Appendix A these expressions can be further simplified to

[TABLE]

where we have taken the motion of the transiting object as ${\bf r}={\bf r}_{0}+{\bf v}t$. We define $x_{D}\equiv(t-t_{D,0})/\tau_{D}$, $x_{S}\equiv(t-t_{S,0})/\tau_{S}$ as normalized time variables. Here, the width of each signal is given by $\tau_{D}\equiv\left|{\bf r}_{0}\times{\mathbf{v}}\right|/v^{2}$ and $\tau_{S}\equiv\left|\mathbf{v}_{\times}\times\mathbf{r}_{\times}\right|/v_{\times}^{2}$. The times for the passing object to reach its point of closest approach are given by $t_{D,0}\equiv-{\bf r}_{0}\cdot{\mathbf{v}}/v^{2}$, $t_{S,0}\equiv-\left(\mathbf{v}_{\times}\cdot\mathbf{r}_{\times}\right)/v_{\times}^{2}$. For the Doppler delay, the vector pointing from the pulsar to the point of closest approach is given by ${\mathbf{b}}_{D}\equiv{\bf r}_{0}+{\mathbf{v}}t_{D,0}$. For the Shapiro delay the relevant vector points from the line of sight to the point of closest approach, and is given by $\mathbf{b}_{S}=\mathbf{\hat{d}}\times(\mathbf{r}_{\times}+\mathbf{v}_{\times}t_{S,0})$. From here on we will drop the $D,S$ subscripts which will be apparent by context.

The signal shapes are shown in Fig. 1. The Doppler signal has two components depending on the orientation of the incoming object, a transient signal ($\propto\hat{{\mathbf{v}}}\cdot\hat{{\mathbf{d}}}$) and a non-transient signal ($\propto\hat{{\mathbf{b}}}\cdot\hat{{\mathbf{d}}}$). The Shapiro signal is always transient regardless of orientation.

Note that one may be tempted to conclude immediately that a Shapiro signal is always subdominant to the Doppler signal, as it is suppressed by $v^{2}$. However, the Shapiro signal is amplified by the long baseline ($\sim$ kpc) resulting in a much shorter typical timescale, and is able to probe a complementary mass window. We consider this in detail in the next sections.

III Signal Analysis

Both the Doppler and Shapiro signals have a characteristic time scale, $\tau$, corresponding to the time for the compact object to pass the line of sight. We note that the signal width ($\tau$) and time to the center of the blip ($t_{0}$) are parametrically the same scale, with the differences being due to the objects’ orientation.

If $\tau\gtrsim T$, where $T$ is the observing time, which we define as the static limit, we will observe only a small section of the signal, which will have a power series expansion in small $t/\tau$. As we discussed in Sec. I, in the static limit, the first two terms in the expansion are unobservable as they are degenerate with the frequency and its derivative. However, a traversing compact object can still be detected by its higher order contributions to Eq. (1) corresponding to the coefficient of the ${\cal O}(t^{2})$ term (this is degenerate with a measurement of pulsar frequency second derivative, which is known to be small).

If, on the other hand, $\tau\lesssim T$ (the dynamic limit), the whole signal shape is seen and therefore a power series expansion will no longer hold. Since the typical compact object spacing is smaller for lower masses, the distance to the pulsar timing system is also smaller making the dynamic signal dominantly present at lower masses. In this limit, one can look for the entire signal shape in pulsar timing data, analogous to searches for gravitational waves and stellar microlensing events. Note that deep in this same regime multiple events will typically transit the line of sight over the observation period, where a statistical approach is relevant as proposed in Ref. Baghram et al. (2011); we leave analysis of such a multi-event signal for future work Dror et al. (2019), where we expect improved reach at lower masses.

We now compute the observables in pulsar timing experiments for the four different searches (Doppler dynamic and static, Shapiro dynamic and static) and their corresponding signal to noise ratio (SNR). The SNR in each case can be estimated analytically by assuming that the constraints are dominated by the object that comes closest to the pulsar (for the Doppler delay) or the line of sight (for the Shapiro delay). We call this the closest-object approximation and it holds in most of the parameter space for all four searches. We also derive an analytic estimate for the split between the dynamic and static limits for the different searches, highlighting their corresponding sensitivity regions.

In the following sections, we apply the following statistical procedure to determine the reach. In order to ensure no false positives among the entire pulsar timing array we need to set a threshold for the SNR. In the absence of a signal the SNR at each pulsar is a one-sided Gaussian random variable and therefore we can calculate the $95\%$ confidence threshold value, $x$, by

[TABLE]

which gives $x=3.66$ with $N_{P}=200$. We fix the threshold value to be four for simplicity, which ensures no false positives with greater than $95\%$ confidence. In our Monte Carlo simulations, we also require that a signal manifests in $90\%$ of randomly generated universes. Correspondingly, in our analytic estimates we use the $90$th percentile relevant length scale, denoted by a ‘min’ subscript, derived in Appendix B.

III.1 Static Limit

In the static regime the constraint is derived from requiring that $\ddot{\nu}/\nu$ is small enough to be consistent with the fit shown in Eq. (1). In setting constraints, we assume a dedicated analysis where $\nu$, $\dot{\nu}$ and $\ddot{\nu}$ are fit simultaneously (as opposed to the usual procedure which only fits $\nu$ and $\dot{\nu}$, and one assumes $\ddot{\nu}$ is small). This is necessary since otherwise the fits for $\nu$ and $\dot{\nu}$ will absorb part of $\ddot{\nu}$, diminishing the signal.

Assuming the data can be characterized by white noise and a signal of the form of Eq. (1), the RMS noise can be taken as the uncertainty in each measurement, and the total expected uncertainty obtainable by a least squares fit on the second derivative is found to be Liu et al. (2018)222We disagree with the limit setting procedure employed in Schutz and Liu (2017) which requires the cubic term in the timing residuals to be below $t_{{\rm RMS}}$ (corresponding to the condition, $\ddot{\nu}/\nu\lesssim 6t_{{\rm RMS}}/T^{3}$) since this does not account for the “sampling factor” of $\sqrt{2800\Delta t/T}$. Coincidentally for pulsar timing array data this is an ${\cal O}(1)$ correction for most pulsars since $10^{3}\text{ wk}\sim 10\text{ yr}$.

[TABLE]

Current pulsar data have an uncertainty of ${\cal O}(10^{-31}~{}{\rm Hz}^{2})$ while sensitivities are projected to reach ${\cal O}(10^{-33}~{}{\rm Hz}^{2})$ for a single pulsar. This allows us to define a suitable SNR,

[TABLE]

There are several observational challenges in implementing this analysis. First, in addition to DM compact objects there are other sources which produce a contribution to $\ddot{\nu}$, such as the existence of dark planets near to the pulsar, as well as a genuine spin down of the pulsar Liu et al. (2018)333Objects in our solar system are not an important background since they experience a yearly modulation and are fit for in the analysis.. Given this, a static search presents challenges as a discovery method, though it can be reliably used to set constraints on the existence of compact objects. Interestingly, for some mass ranges, compact objects predict a static Doppler signal in conjunction with a dynamic Shapiro signal, discussed later in this section, which would increase the confidence in the measurement and potentially provide some information about the object size (see Sec. V).

We now calculate the expressions for the static contributions of Doppler and Shapiro signal shapes (when the transiting objects are PBHs), and subsequently estimate their contribution to $\ddot{\nu}$. In the Dynamic subsection we will discuss the division between the static and dynamic signals.

III.1.1 Doppler

To obtain the size of the Doppler signal shape, we take the second time derivative of Eq. (8) evaluated at $t=0$, and use the relation between $\tau$ and $t_{0}$, $\tau^{2}+t_{0}^{2}=r_{0}^{2}/v^{2}$, to find

[TABLE]

where in the final step we used the relation $\left|{\mathbf{b}}\right|=\tau v$.

Note that a static signal can never be isolated and will always be due to a collection of compact objects. However, to understand the sensitivity analytically, we can still make progress by employing the closest-object approximation. A concern in this approach is whether the impact of far away objects is truly small since the number of compact objects at a given distance grows with distance. However, this is usually a small effect for two reasons. First, the signal size has a steep power of $r_{0}$ in the denominator. Second, the contribution to $\ddot{\nu}$ does not grow coherently with number of objects and the contribution from a single object can be positive or negative (depending on the object’s trajectory). Nevertheless, we note that this approach breaks down when the signal is not deep within the static regime, and where the contribution from multiple objects is needed to adequately estimate the signal size.

Therefore to estimate the signal size from Eq. 14, we calculate the minimum typical distance of a DM compact object. This is derived from the minimum distance of randomly distributed points, around $N_{P}$ pulsars, to the closest pulsar, and we assume each compact object is of mass $M$. The result is derived in Appendix B and we quote the result here,

[TABLE]

Roughly speaking, the static signal condition can be taken as $r_{{\rm min}}\gtrsim vT$. However, as discussed above, this is the condition that the static analytic estimates are valid, not that a static search cannot be performed, as the static analytic estimates make use of the closest-object approximation; but when many objects are near the static limit boundary, many objects make similar contributions to $\ddot{\nu}$.

With this expression for $r_{{\rm min}}$ we can now estimate the size of this signal (dropping angular factors),

[TABLE]

Notice that the mass dependence has dropped out due to scaling of the minimum distance as $M^{\frac{1}{3}}$. We note that single pulsar measurements have $\sigma_{\ddot{\nu}/\nu}\gtrsim 10^{-31}~{}{\rm Hz}^{2}$ and so do not have sufficient sensitivity to see this static signal. However, as we will see future arrays profit significantly from an increase in both observation time (as the uncertainty drops as $\propto T^{-7/2}$), and from many more pulsars in the, and so will be capable of measuring such tiny deviations.

In general compact objects will generate a signal if they gravitationally interact with either the Earth or the pulsar, which we label as the ‘Earth term’ and ‘pulsar term’, respectively. Using correlations between pulsars one can reduce the noise that affects only the signals due to a compact object interacting with the Earth. The constraints on a Doppler signal are, however, always dominated by the pulsar term, as opposed to the Earth term. This can be understood in the limit where all pulsars are identical since in that case for the Earth term, $\sigma_{\ddot{\nu}/\nu}\propto 1/\sqrt{N_{P}}$, while the signal is constant, such that the SNR scales as $\sqrt{N_{P}}$. On other hand, for the pulsar term the noise is independent of the number of pulsars, while the signal size grows linearly with $N_{P}$, such that the SNR scales as $N_{P}$. Therefore, for large $N_{P}$, the pulsar term dominates.

III.1.2 Shapiro

The computation of the static Shapiro signal is analogous to the static Doppler signal. Taking the second time derivative of Eq. (9), evaluating at $t=0$, and simplifying with $t_{0}^{2}+\tau^{2}=r_{\times,0}^{2}/v_{\times}^{2}$, gives

[TABLE]

As expected, $\ddot{\nu}$ appears parametrically suppressed compared to the Doppler contribution, though (as commented previously) this is deceiving due to the different, and typically smaller, distance scale in the denominator. We also find the same power of $r_{0}$ in the denominator, suggesting that we can again use a closest-object approximation up to the boundary between the static and dynamic regimes where multiple objects are crucial for obtaining the correct signal size.

In a manner similar the Doppler case, we compute the smallest expected distance of a DM compact object to the line of sight toward some pulsar, $r_{\times,\text{min}}$. This is done in Appendix B, with the final result,

[TABLE]

As before, we are now able to estimate the size of the cubic $\ddot{\nu}/\nu$ in the closest-object approximation. Omitting angular factors, this is

[TABLE]

Note that the mass dependence does not drop out, as it did in the Doppler static case. This can be traced to the geometry involving the distance to the line of sight, which results in the minimum distance scaling as the square root (rather than 1/3 power) of the number density of DM compact objects. This mass scaling agrees with the analysis of Ref. Schutz and Liu (2017) (though as we commented previously the constraint on $f$ we obtain does not agree because of a difference in the limit setting procedure). Because the Shapiro static signature is small and subdominant to the Doppler static signal, current pulsar data are unable to constrain a static signal from DM compact objects. Nevertheless in Sec. IV we show that future arrays may be able to observe such tiny contributions.

III.2 Dynamic limit

In the dynamic limit a compact object is close enough to the line of sight that it crosses in a time smaller than the observation time. This means that pulsar timing experiments see the entire signal shape, rendering the expansions in the static limit invalid. Specifically, we take the dynamic constraint to be $\tau<t_{0}<T-\tau$ and we note that this implies $\tau<T/2$. To extract small signals out of a noisy background we use the prescription employed in gravitational wave searches known as the Matched Filter procedure Moore et al. (2015a, b). The idea is to take the time-of-arrival data, apply a filtering procedure (namely, we convolute the data with the Weiner filter), and extract an optimal signal to noise ratio (SNR). For simplicity, we work in the limit $\Delta t\ll\tau,\left|t_{0}\right|\ll T$ such that the measurement is unaffected by cadence or finite width effects (adding these is straightforward but complicates the expressions). Furthermore, we assume that the timing residual noise is white with a variance given by $t_{\text{RMS}}$, i.e. Moore et al. (2015a)

[TABLE]

The signal we consider here is $\delta\nu/\nu=\dot{\delta t}$, whose power spectrum is given by

[TABLE]

Using a one-sided power spectral density for the noise we identify, $S_{\dot{\delta t}}(f)\equiv 8\pi^{2}t_{\text{RMS}}^{2}\Delta t\,f^{2}$, giving a SNR Moore et al. (2015a),

[TABLE]

where $\widetilde{h}(f)$ is the Fourier transform of the $\delta\nu/\nu$ signal. We now compute this for Doppler and Shapiro signals given in Eqs. (8) and (9).

Unlike the static signal, we expect the backgrounds in the dynamic case to be less worrisome. Most importantly, the characteristic signal shape is unlikely to have significant overlap with other sources of noise. Perhaps the most prominent candidate to mimic a dynamic signal are pulsar glitches, which have recently been observed in millisecond pulsars Cognard and Backer (2004); McKee et al. (2016). However, pulsar glitches are well parameterized by an instantaneous peak in the phase with a subsequent falling exponential. Thus they have a different frequency structure than the signals of interest here. A more troubling background is dark baryonic objects. The baryonic mass distribution is, however, peaked near a solar mass, whereas the objects we consider in the dynamic limit have masses $M\lesssim 10^{-2}\,M_{\odot}$.

III.2.1 Doppler

To find the SNR of the Doppler signal we insert the signal shape of Eq. (8) into the expression for the SNR in Eq. (23). This is valid because, in contrast to the static case, the fitting procedure for $\nu,~{}\dot{\nu}$ is not degenerate with the signal. This gives an SNR of

[TABLE]

Note that we have dropped the term proportional to $\hat{{\mathbf{v}}}\cdot\hat{{\mathbf{d}}}$ in Eq. 8, as it is parametrically suppressed by $\tau/T$, which is small in the dynamic limit.

The signal size, as well as the transition between the dynamic and static regimes, can be understood by employing the closest-object approximation, as discussed previously. Inspection of Eq. (24) suggests that a closest object approximation should hold due to the $1/\tau\propto 1/r_{0}$ in the denominator (and the small spread in the velocity distribution), and we have checked this using a Monte Carlo simulation, which we discuss in Sec. IV. In order to obtain an estimate of the SNR, we compute an estimate of the minimum $\tau=\tau_{\text{min}}$, generated by a random set of points. Note that $\tau$ also corresponds to the minimum impact parameter since, $\left|{\mathbf{b}}\right|=\tau/v$. We calculate this explicitly in Appendix B and quote the result here,

[TABLE]

Combining Eqs. (24) and (25) gives a good estimate of the largest SNR for a given mass and DM density. We can further estimate the condition for the nearest object to be in the dynamic limit, meaning $\tau_{{\rm min}}\lesssim T/2$,

[TABLE]

For such small masses, these signals are sufficiently quick to appear as a transient in pulsar timing experiments, while for larger masses the Doppler signal can only appear in a static search which was detailed previously.

Again one can compare contributions from the pulsar and Earth terms assuming the pulsars are identical. For dynamic signals, the Earth term scales as $S_{n}\propto 1/N_{P}$, while the signal is constant. However, for the pulsar term the noise is independent of the number of pulsars, while the signal size grows linearly with $N_{P}$. Therefore, the SNR scales identically with $N_{P}$ for the pulsar and Earth term.444This is in contrast to the results presented in Kashiyama and Seto (2012), which claims to achieve more powerful constraints with the Earth term at lower masses. The discrepancy can be traced to parametrically different estimates of the impact parameter, $b$. Ref. Kashiyama and Seto (2012) assumes $b\sim vT$, based on dimensional analysis, whereas we derive more specific estimates. Therefore, deep in the dynamic limit their constraints should be comparable. However, for the pulsar term the dynamic condition, $\tau<T/2$, is easier to satisfy since $\tau\propto 1/\sqrt{N_{P}}$, whereas for the Earth term $\tau$ is independent of $N_{P}$. Thus for larger masses we expect the pulsar term to become more sensitive. For simplicity we only use the pulsar term but note that one could achieve improved sensitivity deep in the dynamic limit by studying both of these contributions.

III.2.2 Shapiro

Finally we arrive at the dynamic Shapiro signal. To find the SNR of the dynamic Shapiro signal we insert the signal shape of Eq. (9) into the expression for the SNR in Eq. (23) as described in Appendix A. This gives an SNR of

[TABLE]

The SNR is independent of $\tau$ in the $\tau\ll T$ limit. This can be understood intuitively. For a transient signal, the SNR scales linearly with the width of the signal however the signal size scales as the $1/r$-potential and hence $1/\tau$, conspiring to produce an SNR independent of $\tau$ at leading order. Nevertheless, keeping higher order terms in the SNR leads to an expression that eventually decreases as $\tau$ approaches $T$ from below.

Since far objects can produce a signal, the closest object approximation breaks down more quickly than for the dynamic Doppler signal, and one should account for the multiple events in the SNR, reserved for future work Dror et al. (2019) with the methods proposed in Baghram et al. (2011). Accounting for multiple events necessarily involves a random signal shape and therefore the matched filter procedure used here will not be applicable. However, even with a different signal analysis technique an SNR accounting for the multiple events is expected to be larger than this SNR, such that the constraints quoted here are conservative. With a single blip analysis this independence results in a minimum mass at which no signal will be seen for any given pulsar, which will later result in a hard cutoff in the projected constraints.

As before, the minimum $\tau$ is related to the minimum impact parameter of a set of randomly generated points near an infinite line, $\tau=\left|{\mathbf{b}}\right|/v$. Thus we obtain a minimum signal width,

[TABLE]

We can now use this to estimate the maximum mass which will generate a dynamic Shapiro signal, $\tau_{{\rm min}}\lesssim T/2$:

[TABLE]

IV Constraints on Primordial Black Holes

We are now prepared to compute the sensitivity of PTAs to PBHs using the signatures we have discussed. Assuming a null result, we study the capability of the searches to set constraints on DM compact objects in the $(M,f)$ plane, where $f\equiv\Omega/\Omega_{{\rm DM}}$ denotes the fraction of dark matter contained in PBHs of mass $M$.

Before setting constraints let us briefly comment on current and future pulsar timing capabilities. To estimate the capabilities of current pulsars we compiled data from PPTA Reardon et al. (2016), EPTA Desvignes et al. (2016), Nanograv Arzoumanian et al. (2018), as well as the combined international collaboration, IPTA Verbiest et al. (2016), culminating in 73 unique pulsars (for 13 of these pulsars no distance was quoted, so we assume a distance typical of the rest of the set, 1 kpc from the Earth). For the 2015 data releases we assume an additional three years’ observing time in order to derive constraints corresponding to current data. In setting our limits with current data we use parameters from this set without any additional approximations (in particular, we do not resort to approximating the current data as an identical set of pulsars with some chosen parameters). Since pulsar timing precision improves quickly with observation time, continuing to observe these pulsars results in a rapid improvement in the ability of PTAs to discover DM compact objects. Nevertheless, with the upcoming construction of the Square Kilometer Array (SKA) Keane et al. (2015) (and its already running precursor, MeerKAT Bailes et al. (2018)), the number of high precision millisecond pulsars is expected to dramatically increase with the potential of uncovering every millisecond pulsar beaming toward Earth in the entire Milky Way. The particular capabilities of the future millisecond pulsar set are not well-known, due to uncertainties both in final capabilities of the array and the number of detectable pulsars in our galaxy. In forming our projections we use the Phase II numbers in Ref. Rosado et al. (2015) which correspond to 200 millisecond pulsars with 50 ns timing, and two week cadence. Furthermore, we take the typical distance of a pulsar from the Earth to be five kpc. We also present results for a more optimistic case, where SKA finds 1000 millisecond pulsars, with 25 ns timing, which can be observed with weekly cadence and have a typical distance from the Earth of ten kpc. The assumed experimental parameters are summarized in Table 1.

We now present our constraints on $f$ using the analytic formulae derived in the previous sections. Note that the analytic formulae drop angular factors, assume a velocity, $v$, of $250$ km$/$s, and use the SKA PTA parameters given in Table 1. Four of the subsequent constraint equations arise from equating the relevant SNR to four. The other two are simply reformulations of Eqs. (26), (29), which indicate the transition between the dynamic and static regions. Our results are summarized in Fig. 2; these results are generated with a numerical simulation (detailed further below), but are consistent, to ${\cal O}(1)$ numbers, with the analytic results quoted in detail next. The comparison between the analytic results and numerical simulation is discussed further (with a plot detailing differences) in Appendix C.

To begin we consider the Doppler search in the dynamic limit. In this case dropping the angular factors in the SNR equation, Eq. (24), equating the SNR to four, and substituting $\tau=\tau_{\text{min}}$ from Eq. (25), constrains $f$ to

[TABLE]

The $L$ superscript denotes that this analytic constraint corresponds to the left-hand side of the triangular “Doppler-dyn” wedge in Fig. 2, labeled $\propto M^{-1}$ at low masses. This behavior does not continue indefinitely, but is cut off when the closest object no longer satisfies our dynamic condition, $\tau_{\text{min}}\lesssim T/2$, where $\tau_{\text{min}}$ is given by Eq. (25). This is equivalent to Eq. 26 and constrains $f$ to

[TABLE]

where the $R$ superscript indicates that this analytic constraint corresponds to the right hand side of the same triangular wedge in Fig. 2, labeled $\propto M$.

We now repeat the above arguments for the other searches. The Shapiro constraint in the dynamic limit is obtained from equating the SNR, Eq. (27), to four, which sets a lower bound on the masses reachable with a single event,

[TABLE]

corresponding to the left-hand side of the wedge labeled “Shapiro-dyn” in Fig. 2.

Furthermore, the dynamic condition, $\tau_{\text{min}}<T/2$, with $\tau_{\text{min}}$ again given by Eq. (28) can be written, similar to Eq. 29, as a condition on $f$ as,

[TABLE]

corresponding to the right-hand side of the Shapiro dynamic wedge in Fig. 2 labeled $\propto M$.

Similarly analytic scalings in the static limit can also be derived. Equating Eq. (16) to four and substituting $r_{0}=r_{\text{min}}$ from Eq. (15) yields a constraint on $f$ from the static Doppler search,

[TABLE]

This corresponds to the curve labeled “Doppler-stat” and $\propto M^{0}$ in Fig. 2. Likewise equating Eq. (19) to four and substituting $r_{\times}=r_{\times,\text{min}}$ from Eq. (18) sets the constraint on $f$ from the static Shapiro search,

[TABLE]

which corresponds to the “Shapiro-stat” curve in Fig. 2 labeled with scaling $\propto M^{1/3}$. Note that, as mentioned earlier, the scaling of the static results at their low mass end, where the static approximation is breaking down, is not trivial since it does not follow the closest-object approximation. We do not attempt to study this in detail analytically but note that over most of the interesting parameter space, the static search in the highly static regime ($\tau\gg T$) has the most promising reach.

While the analytic results give the correct scaling and approximately correct magnitude, to set more rigorous constraints, which take into account the dark matter velocity distribution and all angular factors, we employ a Monte Carlo (MC) simulation which takes the PBHs to be monochromatic (of a single mass) and randomly distributed with density $\rho_{{\rm DM}}=0.46~{}{\rm GeV}/{\rm cm}^{3}$ Sivertsson et al. (2017), and with a Maxwell Boltzmann velocity distribution with a mean of $220~{}{\rm km/s}$. The statistics employed are the same as discussed in the beginning of Sec. III. In both the dynamic and static limits we take the signal to be the largest SNR generated in the entire array of pulsars.

We consider four different scenarios composed of the different pulsar sets. In the first scenario (SC1) we compute the constraints for only the current 73 pulsar set, assuming the search is done today. In scenario two (SC2) we assume the same set of pulsars, except the observation time of each is increased by ten years. In the third scenario (SC3) we consider the current pulsar set observed for 30 more years, and the addition of the set of pulsars from SKA measured for 20 years. This effectively assumes SKA will start taking data ten years from now with all current pulsars continuing to be studied. Lastly, scenario four (SC4) is the same as scenario three but with the optimistic parameters. In this case the new pulsar measurements dominate the current pulsar measurements.

We present our full results in Fig. 3. The main peaks at $10^{-9}-10^{-8}\,M_{\odot}$ and $10^{-4}-10^{-3}\,M_{\odot}$ arise from the dynamic Doppler and Shapiro searches. The non-smooth behavoir of the dynamic Shapiro curves is (largely) due to current pulsars each having different noise. While for other searches, the pulsar timing array will have an “effective noise” leading to smooth curves, the single blip approximation in Shapiro leads to a sharp cut in sensitivity for each pulsar in the set at a different mass (as a consequence of the $\tau$-independence of the SNR) leading to the features shown in the figure. For SC3 and SC4, where the constraints dominantly arise from a new (assumed identical) set of pulsars from SKA, these features do not arise. At larger masses, the static searches become important with Doppler static becoming approximately constant at large masses and static Shapiro only becoming relevant for SC4 (where it gains due to its large assumed baseline). For comparison we include constraints from lensing experiments, which, as we will show in the next section, only apply for very compact objects. We see that even with the current set of pulsars, a dedicated search could begin to probe $f<1$ within ten years. With the inclusion of additional pulsars from SKA, pulsar timing can scan a huge mass range, from as low as $\sim 10^{-12}\,M_{\odot}$ and could constrain PBHs to sub-percent fractions of the DM. As already mentioned, in addition to the correct scalings, we also find that the MC and analytic results are in agreement up to $\mathcal{O}(1)$ factors, with the main differences being due to the angular factors, as well as all PBHs sharing a uniform velocity in the analytic approximation.

Lastly, we note that there are several ways that these constraints could be drastically improved besides the addition of pulsars with better timing parameters, as can be seen in Eqs. (30), (32), (34), and (35). For example, if a millisecond pulsar is found much farther away (say $d\sim 10-100\,\text{ kpc}$), the constraints from the Shapiro delay will improve, while a pulsar in a DM dense region (such as near the galactic center or in globular clusters) will yield stronger constraints for all proposed signals.

V Constraints on Subhalos

We now turn to constraining more diffuse DM subhalos, which as we now show, can be detected using the proposed searches for compact objects. In principle, diffuse subhalo signals can be calculated using the same procedure we invoked for compact objects, namely, computing $\delta\nu/\nu$ for the Doppler and Shapiro signals, and finding the SNR using Eqs. (12) or (23), depending on the mass range of interest. The induced strain by a passing object is only dependent on the gradient of the gravitational potential due to the passing object (see Eqs. (4) and (5)). For a subhalo a distance $r$ from the pulsar timing system, Gauss’s law states that the gradient of the potential is given by,

[TABLE]

where the integral runs from the center of the subhalo to the point of interest. Substituting this expression into (4) or (5) gives $\delta\nu/\nu$ for a generic DM subhalo profile. Computing the signal shape and size is, however, rather involved; after all, $M(r)$ is time dependent as the halo moves. On the other hand, one can set conservative bounds by replacing $M$ in Eqs. (14), (17), (24), and (27), by the minimum $M(r)$ over the time of observation, which is calculated at the point of closest approach of the halo to the pulsar (Doppler) or line of sight (Shapiro).

The above lower bound is a particularly good estimate for dynamic signals. In this case the signal shape can be split up into two components, the one from the inner ring encapsulated by the impact parameter $b$ and an additional piece for the outer ring:

[TABLE]

A matched filtering prescription is optimized to look for a specific signal shape and hence would remove (part of) the additional contribution from the ring, leaving behind primarily a signal from the inner circle of radius $b$. While more sophisticated analysis could improve and include these effects, it is beyond the scope of this work.

In the case of static signals, the non-compactness of the subhalo manifests itself as additional contributions in time derivatives of $M(r)$ in both the Doppler and Shapiro signals, as well a deviation in the integral over the line of sight, in the case of Shapiro. We have verified that these corrections contribute $\mathcal{O}(1)$ to the lower bound evaluated simply by taking $M(r)$ as evaluated at the initial position, $M(r_{0})$.

To relate the limits computed earlier for PBHs to subhalos, it is convenient to define a sensitivity distance, which is the typical distance of a compact object to the pulsar (Doppler) or line of sight (Shapiro) in order to induce a particular SNR (e.g. four). Limits on a particular DM fraction at a given mass are set when the minimum distance is smaller than this sensitivity distance.

The sensitivity distance in the case of the dynamic Doppler signal can be computed using Eq. (24), and substituting SKA parameters given in Table 1 while taking $v\sim 10^{-3}$, gives

[TABLE]

The Shapiro dynamic sensitivity distance is a more complicated function of mass, as the SNR is $\tau$ independent in the $\tau\ll T$ limit. The curve shown in Fig. 4 is derived by taking the SNR given in Appendix A and finding the $r,M$ curve which satisfies $\text{SNR}(\tau=r/v,M)=4$.

In the static limit, Eqs. (14) and (17) give the distances,

[TABLE]

We plot the sensitivity distances as a function of mass in Fig. 4. Note that in Fig. 4 the dynamic curves end at $r=vT$ since, if the object passes at larger $r$, it would not satisfy the dynamic condition discussed earlier.

If the subhalo radius is smaller than the sensitivity distance then its effects on pulsar timing searches are identical to a PBH of the same mass. On the other hand, if the DM subhalo has a radius larger than the sensitivity distance, it may still be constrained if $r_{\text{enc}}$, the inversion of $M(r)$, is less than the sensitivity distance for some mass $M$. Physically, this means that the whole subhalo is too diffuse to measure but the core may be compact enough to measure.

To explore this possibility we consider two halo profiles of the generic form,

[TABLE]

where $M_{\text{vir}}$ is the virial mass of the halo, $c\equiv r_{{\rm vir}}/r_{s}$ is the concentration parameter, $r_{\text{vir}}\equiv(3M_{\text{vir}}/800\pi\rho_{c})^{1/3}$ is the virial radius, and $\rho_{s}$ is an overall normalization factor fixed by requiring that the total mass inside of the virial radius is the virial mass. The standard NFW profile corresponds to taking $\alpha=1,~{}\beta=2$, and an ultracompact minihalo Fillmore and Goldreich (1984); Bertschinger (1985) corresponds to $\alpha=9/4,~{}\beta=0$ (though see, e.g., Delos et al. which suggests an $\alpha=3/2$, $\beta=3/2$ profile can be a better fit to numerical simulations for halos produced from gravitational collapse of some primordial power spectra).

In Fig. 4 Left (Right) we plot $r_{\text{enc}}$ for NFW (UCMH-like $\alpha=9/4,~{}\beta=3/4$) halos of virial mass $M_{\text{vir}}=10^{-7}\,M_{\odot}$, $10^{-1}\,M_{\odot}$ and concentration parameters $c=100$, $10^{8}$ and the PBH-limit, $c\rightarrow\infty$. For a given subhalo, if $r_{\text{enc}}$ passes below a particular sensitivity curve at a mass, $M^{\ast}$, then the search is sensitive to an effective subhalo of mass $M^{\ast}$. These $M^{\ast}$ mass subhalos make up $f^{\ast}=f\,M^{\ast}/M_{\rm vir}$ of the DM. Hence the limits of interest, $(M_{\text{vir}},f)$, can be (conservatively) extracted from $(M^{\ast},f^{\ast})$, which is constrained for PBHs. Thus finally we obtain,

[TABLE]

where $f_{{\rm PBH}}(M^{\ast})$ it the limit extracted from Fig. 3. As we can see in Fig. 4, a Doppler search is sensitive to the largest radii, followed by a Shapiro delay search.

A similar procedure can be adopted to translate microlensing constraints for PBHs to subhalos. This time the sensitivity radius corresponds the Einstein radius, $r_{E}$, above which microlensing experiments will not see modulations in the source brightness. The Einstein radius for a source and lens at distance $D_{S}$ and $D_{L}$ respectively is given by,

[TABLE]

The sources considered by microlensing experiments range from nearby stars Griest et al. (2011, 2014), the Large Magellanic Cloud (e.g.,Allsman et al. (2001); Tisserand et al. (2007)), Andromeda (e.g., Niikura et al. (2017)), and even distant supernova Zumalacarregui and Seljak (2017). For stellar sources, the distances are $\mathcal{O}(1-100\text{ kpc})$ while supernovae are sensitive to much larger distances, ${\cal O}(\text{Gpc})$. The typical Einstein radii of these sources are,

[TABLE]

and are illustrated in Fig. 4. For diffuse halos with radii much larger than this distance microlensing is unable to see a significant signal. Note that because the Einstein radii are much smaller than the PTA sensitivity distance, time delays due to the creation of multiple source images (as considered in Ref. Katz et al. (2018)) are subdominant to the effects considered here.

Finally, we use Eq. (44) and the intersections in Fig. 4 to find limits on the dark matter fraction $f$ as a function of the total DM subhalo mass, $M_{{\rm vir}}$, in Fig. 5 for the SKA pulsar set defined in Sec. IV as SC3. The lensing constraints are scaled in the same way as pulsar timing, where we take the sensitivity line to be the characteristic Einstein radius from objects in the Andromeda galaxy for the Subaru constraint and from objects in the Milky Way for the MACHO/Eros/Ogle constraint curve.

We observe in Fig. 5 that in the PBH-limit, in most of the mass range, constraints in SC3 are stronger for lensing compared to pulsar timing; however as the concentration parameter is decreased, lensing drops off in sensitivity relative to PTAs. We find that for NFW halos with $c=10^{8}$ the Subaru search is only sensitive to halos of $M\gtrsim 10^{-9}~{}M_{\odot}$. For a CDM-inspired Sánchez-Conde and Prada (2014) concentration parameter $c=100$, we no longer observe any sensitivity from lensing, while PTAs can constrain a non-negligible $f$. Note, however, that this sensitivity occurs only for very low halo masses (significantly below an Earth mass) where one expects halo disruption. On the other hand, for slightly higher concentration parameters $c=10^{3}$, sensitivity to $f<1$ occurs similarly to $c=100$ for low mass halos, but also (from a static Doppler search) for $M\gtrsim 10^{-3}\,M_{\odot}$.

While our analysis of diffuse halos is schematic and suffers from ${\cal O}(1)$ corrections, it serves to emphasize the complementarity between lensing and pulsar timing probes. Fully exploiting the potential of PTAs to constrain diffuse halos and specific models of structure formation is an intriguing problem which we leave for future work Dror et al. (2019).

VI Conclusions

In this work we considered pulsar timing constraints on DM compact objects, focusing on primordial black holes and subhalos. We studied pulsar timing signatures over the mass range from $10^{-12}\,M_{\odot}$ to $100\,M_{\odot}$ finding that, depending on an objects mass, different searches are required to detect it. We examined four different types of searches, dynamic and static signals, each arising from Doppler and Shapiro time delays. Importantly, we computed the signals in three-dimensions and highlighted their relation to one another. Furthermore, using a Monte Carlo analysis we performed projections for pulsar timing capabilities using current and future pulsar timing experiments and understood their scaling using analytic techniques. With dedicated searches we found that current pulsar timing arrays can, over the next decade, set non-negligible constraints through dynamic searches. Farther into the future, we expect sub-percent level constraints over the entire range with upcoming pulsar timing arrays.

There are two primary ways in which the capabilities of detecting DM compact objects from pulsar timing can be drastically improved. First, we assumed that all DM compact objects are in regions with DM density comparable to our local density. If instead pulsars are discovered in DM-dense regions (e.g., close to the galactic center or within dwarf galaxies), then it is possible to quickly improve the power of Doppler signals. Similarly if pulsars are discovered with a line of sight passing through a DM-dense region, then the capabilities of a Shapiro search will be greatly enhanced. Second, the Shapiro search potential is sensitive to the distance to the pulsar (in the case of uniform density, limits on the fraction of the DM constrained scale linearly with the baseline), such that if pulsars are discovered significantly farther from our local neighborhood (e.g., extra-galactic pulsars), then the Shapiro search will quickly become more powerful.

The constraints studied here apply to substructure which has survived to the present day, and we do not attempt to relate these substructures to specific astrophysical or particle physics models. Relating structure on such small scales to a model is a challenging exercise due to the uncertainties on the survival of low-mass subhalos to the present epoch. Previously, Refs. Baghram et al. (2011); Kashiyama and Oguri (2018) considered evolution of subhalos in a vanilla Cold DM (CDM) paradigm with Stable Clustering and spherical halo models to predict halo substructure, and came to opposite conclusions about the feasibility of detecting DM substructure in the CDM model with PTAs. The difference in the conclusions of these papers is likely partly due to the DM clustering model and partly due to a difference of methods with respect PTA constraints. Utilizing the methods proposed in this paper, we plan to consider PTA constraints on vanilla CDM, and other models such as axion miniclusters, in future work Dror et al. (2019). If theory predictions can be made reliably, pulsar timing will become a powerful tool to probe the nature of DM.

Finally we note that our analysis was entirely focused on single events. For lower masses (below $\sim 10^{-9}\,M_{\odot}$ for the Doppler dynamic search and below $\sim 10^{-4}\,M_{\odot}$ for the Shapiro dynamic search), we expect multi-event signals that are not detectable as single events to become important. Nevertheless, they may leave an imprint in pulsar timing arrays that can detected using a statistical prescription as considered in Baghram et al. (2011). We leave further study of this limit to future work Dror et al. (2019). Taken together, however, a coherent picture is emerging for how compact objects can be constrained across a huge mass range with one observational tool.

Acknowledgments:

We thank Adrienne Erickcek, Andrey Katz, Joachim Kopp, Vikram Ravi, Katelin Schutz, Sergei Sibiryakov, Wei Xue, Xingjiang Zhu, and Miguel Zumalacarregui for useful discussions. We especially thank Stephen Taylor for extensive discussions and collaboration during early parts of this work. JD, HR, TT and KZ are supported in part by the DOE under contract DE-AC02-05CH11231. KZ thanks the CERN theory group for hospitality for the duration of this work. Some of this work was also completed at KITP, supported in part by NSF Grant No. NSF PHY-1748958, and at the Aspen Center for Physics, which is supported by NSF grant PHY-1607611.

Appendix A Simplification of Doppler and Shapiro signals and SNR

This Appendix contains many of the details for the results presented in Sec. II, on the signal shapes for the Doppler and Shapiro delays. In both cases, we work in the limit where the trajectory of the compact object is unaffected by the presence of the pulsar-Earth system (this amounts to assuming highly unbound orbits for the Doppler signal and a large baseline for the Shapiro signal). In either case we define a “time until blip center”, $t_{0}$, and a “signal width”, $\tau$, corresponding to the dot-product and cross-product of the distance-velocity vectors respectively.

A.1 Doppler

A.1.1 Signal

We begin by considering the Doppler effect, where our goal is to compute the velocity of the pulsar as a function of time,

[TABLE]

where ${\bf r}$ is the vector pointing from the pulsar to the compact object and is taken to be unaffected by the presence of the pulsar, ${\bf r}={\bf r}_{0}+{\mathbf{v}}t$. It is convenient to introduce the time variables,

[TABLE]

such that $t_{0}^{2}+\tau^{2}=r_{0}^{2}/v^{2}$. In this case the magnitude of the position is

[TABLE]

where we have defined a normalized time, $x\equiv(t-t_{0})/\tau$. This gives a resultant strain,

[TABLE]

where we defined the impact parameter vector ${\bf b}\equiv{\bf r}_{0}+t_{0}{\mathbf{v}}$, obeying $\left|{\mathbf{b}}\right|=v\tau$. We thus obtain the result in Eq. (8).

A.1.2 SNR

Here we derive Eq. (24) from the signal Eq. 52 in the $\tau,\left|t_{0}\right|\ll T$ limit. The signal must be windowed in order to account for the finite time of the experiment and we do this with a top-hat function in the time domain. The SNR from Eq. (23) in the $\left|t_{0}\right|\ll T$ limit then becomes

[TABLE]

the frequency integral can be regulated with a fictitious “mass” term in the denomator parameterized by $\epsilon$

[TABLE]

where the ellipses denote terms which are odd in $x$ or $y$ and therefore vanish when the $x,y$ integrals are taken (including the divergent piece). Inserting this into the above and carrying out the $x,y$ integrals we find the $\tau\ll T$ limit,

[TABLE]

in agreement with Eq. (24).

A.2 Shapiro

A.2.1 Signal

We now compute the signal shape in the case of a Shapiro time delay starting with (5). The gradient of the PBH potential, $\Phi$, is evaluated along the line of sight (taken in the ${\bf\hat{z}}$ direction), giving

[TABLE]

It is convenient to use coordinates that make the cylindrical symmetry of the problem manifest, ${\bf r}_{\times}\equiv{\bf r}\times{\bf\hat{d}}$ and ${\mathbf{v}}_{\times}\equiv{\mathbf{v}}\times{\bf\hat{d}}$, which results in time variables with the same interpretations as for Doppler:

[TABLE]

Carrying out the integral for signals close to the line of sight we can rewrite the strain in terms of these coordinates as

[TABLE]

As in the Doppler signal we can rewrite the magnitude of the signal shape as

[TABLE]

We have again defined $x\equiv(t-t_{0})/\tau$. Inserting this into (58) gives the signal shape

[TABLE]

A.2.2 SNR

Here we derive Eq. (27) from Eq. 61 in the $\left|t_{0}\right|,\tau\ll T$ limit. Analogous to our approach in deriving the Doppler SNR, in the time domain we window our function for a finite observing time, giving

[TABLE]

As for the Doppler signal we regulate using 54. In the $\tau\ll T$ limit the SNR becomes

[TABLE]

in agreement with Eq. (27).

Appendix B Minimum Distance Statistics

This Appendix is dedicated to deriving the minimum distances and times quoted in Sec. III. Each of the signals considered (Doppler/Shapiro, dynamic/static) depends on the relevant distance of objects described by a random variable, $B$. The dominant signal comes from the object which has the minimal $B$, and therefore to gain an analytic understanding of the dominant signal we calculate the $100\times p$th percentile of $B_{\text{min}}\equiv\text{min}\left(B_{1},B_{2},...,B_{N}\right)$. Every statistic of $B_{\text{min}}$ can be calculated with its cumulative distribution function (CDF), $F_{B_{\text{min}}}$, and because the $B_{i}$’s are independent and identically distributed random variables we can solve for $F_{B_{\text{min}}}$ solely in terms of $F_{B}$:

[TABLE]

The $100\times p$th percentile of $B_{\text{min}}$ is calculated by solving $F_{B_{\text{min}}}(b_{p})=p$. For each signal we take $N$ objects populated inside the relevant volume. At the end of this Appendix we show the comparison between calculating the constraints on $f$ using the distances derived here and our Monte Carlo simulation.

B.1 Doppler

B.1.1 Static

For this signal the relevant distance is simply the distance from the object to a pulsar (we take the pulsar to define the origin). Assuming the objects are uniformly distributed in a sphere of radius one means $F_{B}(b)=b^{3}$. Therefore in the large $N$ limit,

[TABLE]

For a sphere of radius $R$, we multiply $b_{p}$ by $R$ and take $N=4\pi N_{P}nR^{3}/3$, giving

[TABLE]

Substituting $p=0.9$ gives (15).

B.1.2 Dynamic

Here the relevant distance is the impact parameter of the passing object. Being only concerned with an order of magnitude estimate, we simplify the problem of finding the CDF of $B$ by taking the objects to move in the same direction. The impact parameter is then uniformly distributed across the cross sectional area. Therefore inside of a circle of radius one, in the cross section centered at the pulsar, the CDF of $B$ is $F_{B}=b^{2}$ and therefore in the large $N$ limit,

[TABLE]

To find $N$ we calculate how many objects are in the dynamic limit. The dynamic limit is defined by $T-\tau>t_{0}>\tau$. Geometrically this constrains how far an object can be as a function of it’s impact parameter and defines a conical volume. Therefore in time $T$, $N=N_{P}n\pi vTR^{2}/3$. Multiplying (67) by $R$ and substituting in $N$ gives

[TABLE]

Substituting $p=0.9$ and dividing by $v$ gives Eq. (25).

B.2 Shapiro

B.2.1 Static

The relevant distance here is the distance from a PBH to the central axis of a cylinder. Assuming the objects are uniformly distributed in a cylinder of radius one, and length $d$, means $F_{B}(b)=b^{2}$. Therefore in the large $N$ limit,

[TABLE]

Taking $N=N_{P}n\pi dR^{2}$, while multiplying by $R$ to for a cylinder of radius $R$, gives

[TABLE]

Taking $p=0.9$ gives (18).

B.2.2 Dynamic

Analogous to the Doppler dynamic signal, the relevant distance here is the impact parameter defined from the central axis of a cylinder. Again, being only concerned with an order of magnitude estimate, we simplify the problem of finding the CDF of $B$ by taking the objects to move in the same direction. The impact parameter is then uniformly distributed across the cross sectional area. Therefore inside a rectangle of width $d$ and height one, the CDF of $B$ is $F_{B}=2b$ such that in the large $N$ limit,

[TABLE]

The dynamic limit is defined by $T-\tau>t_{0}>\tau$. Geometrically this constrains how far an object can be as a function of its impact parameter, defining a parallelogram. Therefore in time $T$, $N=N_{P}nvTdD/2$. Multiplying (71) by $D$, for a cylinder of diameter $D$, and substituting $N$ gives

[TABLE]

Taking $p=0.9$ and dividing by $v$ gives Eq. (28).

Appendix C Comparison of Analytic and Monte Carlo Constraints

We now make an explicit comparison in Fig. 6 between the analytic and numeric results presented in Sec. IV. Six of the eight analytic curves are described by Eqs. 30, 32, 34, and  35 for the four signal types, and by Eqs. 31 and 33 which are the boundary of the dynamic regimes. The other two (dotted) curves are $r_{\text{min}}=vT$, where $r_{\text{min}}$ is found in Eq. (15), and $r_{\times,\text{min}}=vT$, where $r_{\times,\text{min}}$ is found in Eq. (18). At first glance these lines seem to be in disagreement with the Monte Carlo. However one must remember that the these only indicates when the closest object is in the static limit, and the Monte Carlo is summing the contributions to $\ddot{\nu}$ of every object. Therefore this left hand side tells us where we would should expect the Monte Carlo and analytic methods to start diverging, consistent with what is observed.

As can be seen in Fig. 6, at least deep in the static or dynamic regimes where analytic approximations are expected to hold, the analytic estimates differ from the full Monte Carlo results by $\mathcal{O}(1)$ factors in the scaling regimes. This is due to the two main differences in the analytic and Monte Carlo results. The first difference is that in the full Monte Carlo the velocity of each object is drawn from a Maxwell Boltzmann distribution whereas the analytic results assume $v=250$ km/s. The second difference is because analytic results drop the angular dependence of some expressions, such as $\mathbf{\hat{b}}\cdot\mathbf{\hat{d}}$ in Eq. (24).