Pulsar Timing Probes of Primordial Black Holes and Subhalos

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.


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 also a variety of constraints specific to PBHs. From 10 −13 M down to 10 −15 M (below 10 −15 M PBHs evaporate in a time shorter than the age of the Universe), the existence of white dwarfs [31] currently constrains PBHs, and femtolensing may do so in the future [32] (see [33,34] for earlier work). Above 100 M a variety of constraints from structure formation and Planck fairly severely constrain the PBH abundance, see [35] for a review. There remains a controversial window between ∼ 10 − 100 M where PBHs [36], and not astrophysical black holes or neutron stars, could give rise to the events in LIGO [37]. 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 [38] and Segue I [39]. Many of these constraints are, however, plagued by astrophysical uncertainties.
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 ∼ 10 − 100 M relevant for LIGO signals, where constraints exist but there are substantial astrophysical uncertainties. Second, in the mass window below 10 −11 M , where lensing constraints currently do not apply. Finally, in the regime where lensing constraints are significant for PBHs (10 −11 − 10 M ) but may 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 , 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, where ν is the frequency andν,ν are its first and second derivatives. The most stable pulsars have frequencies of O(kHz) and a spin-down rate of the pulsar,ν/ν, ranging from roughly 10 −23 − 10 −20 Hz, both of which can be fit from the data. Empirically, it is found thatν/ν can be below 10 −31 Hz 2 [40] and is typically not included in fits to the data, allowing one to place upper bounds on processes that would produce a non-negligibleν. Furthermore, any process which induces a modification of the phase, can be constrained using pulsar timing measurements. The quality of pulsar timing data is determined by three parameters. The first parameter is the root-meansquare (RMS) timing residual, t RMS . This is determined after finding the frequency, ν fit , and its derivative,ν fit , which minimizes the residual between the timing data, t data n , and the timing model, t n , where t n is found via the relation 2πn = φ(t n ) from Eq. (1). This gives where N is the number of data points, and t fit n is t n with ν = ν fit ,ν =ν fit and all higher order terms dropped. The minimized residual is typically t RMS ∼ µsec. The other two parameters are the observation time of the pulsar, T ∼ 10 years, and the time between measurements, ∆t ∼ 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 [41]. 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 [42]. 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, τ , associated with the motion of the compact objects (here we use this variable schematically but give it an explicit, massdependent meaning in later sections). If we denote the observation time of a pulsar as T , then dynamic signals correspond to τ T , and will appear as blips in the pulsar timing data (analogous to glitches which have been observed in millisecond pulsar data [43,44]). Static signals, with τ T , will not be observable as blips but instead as a non-negligible contribution to the second derivative of the frequency,ν.
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 [45,46], while searches for dynamic signals were considered for single events in [41,42,47], and multiple events in [48]. 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 [42,47,49]. 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 multievent (or statistical) signal, as studied in [48]. While we expect that such an analysis will extend the reach at the low mass end (below O(10 −9 M ) for Doppler signals and below O(10 −4 M ) for Shapiro signals), due to the more complicated nature of the signal, we reserve study for future work [50].
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, δν. For the Doppler and Shapiro signals, where Φ is the gravitational potential due to the compact object and v is its velocity, whiled 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 , where r is the position of the compact object relative to the pulsar and × subscript denotes crossing withd, r × ≡ r ×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 tod, as only this gives a time dependent shift to the metric affecting the photons. 1 Here we assume a weak field approximation, Φ 1, a slowly varying potential during the interaction time scale (Φ(r + vr) Φ(r)), where r is the distance of closest approach, and large orbit eccentricity. As shown in Appendix A these expressions can be further simplified to where we have taken the motion of the transiting object as r = r 0 + vt. We define x D ≡ (t − t D,0 )/τ D , x S ≡ (t − t S,0 )/τ S as normalized time variables. Here, the width of each signal is given by τ D ≡ |r 0 × v| /v 2 and τ S ≡ |v × × r × | /v 2 × . The times for the passing object to reach its point of closest approach are given by For the Doppler delay, the vector pointing from the pulsar to the point of closest approach is given by b D ≡ r 0 + vt 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 b S = d × (r × + v × 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 (∝v ·d) and a non-transient signal (∝b ·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 (∼ 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, τ , corresponding to the time for the compact object to pass the line of sight. We note that the signal width (τ ) 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 τ 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/τ . 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 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, τ 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. [48]; we leave analysis of such a multi-event signal for future work [50], 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 closestobject 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 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.

A. Static Limit
In the static regime the constraint is derived from requiring thatν/ν is small enough to be consistent with the fit shown in Eq. (1). In setting constraints, we assume a dedicated analysis where ν,ν andν are fit simultaneously (as opposed to the usual procedure which only fits ν andν, and one assumesν is small). This is necessary since otherwise the fits for ν andν will absorb part ofν, 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 [40] Current pulsar data have an uncertainty of O(10 −31 Hz 2 ) while sensitivities are projected to reach O(10 −33 Hz 2 ) for a single pulsar. This allows us to define a suitable SNR, 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ν, such as the existence of dark planets near to the pulsar, as well as a genuine spin down of the pulsar [40] 3 . 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ν. In the Dynamic subsection we will discuss the division between the static and dynamic signals.

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 τ and t 0 , where in the final step we used the relation, |b| = τ 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ν 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, Roughly speaking, the static signal condition can be taken as r min 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ν. With this expression for r min we can now estimate the size of this signal (dropping angular factors), Notice that the mass dependence has dropped out due to scaling of the minimum distance as M 1 3 . We note that single pulsar measurements have σν /ν 10 −31 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 ∝ 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, σν /ν ∝ 1/ √ N P , while the signal is constant, such that the SNR scales as √ 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.

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 As expected,ν 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 ×,min . This is done in Appendix B, with the final result, As before, we are now able to estimate the size of the cubicν/ν in the closest-object approximation. Omitting angular factors, this is 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. [46] (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.

B. 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 τ < t 0 < T − τ and we note that this implies τ < 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 [51,52]. 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 ∆t τ T such that the measurement is unaffected by cadence or finite width effects (adding these is straightforward but complicates the expressions). Furthermore, we assume the measurement noise is white. In this case the SNR is given by, where S n ≡ ∆t (t RMS /T ) 2 . 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 [43,44]. 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 10 −2 M .

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. (20). This is valid because, in contrast to the static case, the fitting procedure for ν,ν is not degenerate with the signal. This gives an SNR of Note that we have dropped the term proportional tov ·d in Eq. 8, as it is parametrically suppressed by τ /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 the employing a closest-object approximation, as discussed previously. Inspection of Eq. (21) suggests that a closest object approximation should hold due to the 1/τ ∝ 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 τ = τ min , generated by a random set of points. Note that τ also corresponds to the minimum impact parameter since, |b| = τ /v. We calculate this explicitly in Appendix B and quote the result here, Combining Eqs. (21) and (22) 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 τ min T /2, 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 ∝ 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. 4 Therefore, deep in the dynamic limit their constraints should be comparable. However, for the pulsar term the dynamic condition, τ < T /2, is easier to satisfy since τ ∝ 1/ √ N P , whereas for the Earth term τ 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.

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. (20). This gives an SNR of Here the SNR drops off as τ −1/2 and therefore we can again estimate the limits using the closest-object approximation. As before, the minimum τ is related to the minimum impact parameter of a set of randomly generated points near an infinite line, τ = |b| /v. Thus we obtain a minimum signal width, We can now use this to estimate the maximum mass which will generate a dynamic Shapiro signal, τ min T /2: 4 This is in contrast to the results presented in [47], 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. [47] assumes b ∼ vT , based on dimensional analysis, whereas we derive more specific estimates.   [59]. Current constraints are compiled from various sources [53][54][55]60] as described in the text.

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 ≡ Ω/Ω 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 [53], EPTA [54], Nanograv [55], as well as the combined international collaboration, IPTA [56], 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) [57] (and its already running precursor, MeerKAT [58]), 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. [59] 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 I. 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 I. Four of the subsequent constraint equations arise from equating the relevant SNR to four. The other two are simply reformulations of Eqs. (23), (26), 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 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. (21), equating the SNR to four, and substituting τ = τ min from Eq. (22), constrains f to The L superscript denotes that this analytic constraint corresponds to the left-hand side of the triangular "Doppler-dyn" wedge in Fig. 2, labeled ∝ 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, τ min T /2, where τ min is given by Eq. (22). This is equivalent to Eq. 23 and constrains f to where the R superscript indicates that this analytic constraint corresponds to the right hand side of the same triangular wedge in Fig. 2, labeled ∝ 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.
which corresponds to the left-hand side of the wedge labeled "Shapiro-dyn" in Fig. 2 labeled ∝ M −1 . Furthermore, the dynamic condition, τ min < T /2, with τ min again given by Eq. (25) can be written, similar to Eq. 26, as a condition on f as, corresponding to the right-hand side of the Shapiro dynamic wedge in Fig. 2  which corresponds to the "Shapiro-stat" curve in Fig. 2 labeled with scaling ∝ 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 (τ 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 ρ DM = 0.46 GeV/cm 3 [61], and with a Maxwell Boltzmann velocity distribution with a mean of 220 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 and 10 −4 − 10 −3 M arise from the dynamic Doppler and Shapiro searches. 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 ∼ 10 −12 M 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 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. (27), (29), (31), and (32). For example, if a millisecond pulsar is found much farther away (say d ∼ 10 − 100 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 δν/ν for the Doppler and Shapiro signals, and finding the SNR using Eqs. (12) or (20), 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 implies that the gradient of the potential is given by, where the integral runs from the center of the subhalo to the point of interest. Substituting this expression into (4) or (5) gives δν/ν 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), (21), and (24), 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: 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 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 dynamic signals can be computed using Eqs. (21) and (24) with SKA parameters given in Table I, and taking v ∼ 10 −3 , In the static limit, using Eqs. (14) and (17) give the distances, 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 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, where M vir is the virial mass of the halo, c ≡ r vir /r s is the concentration parameter, r vir ≡ (3M vir /800πρ c ) 1/3 is the virial radius, and ρ 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 α = 1, β = 2, and an ultracompact minihalo [62,63] corresponds to α = 9/4, β = 0 (though see, e.g., [64] which suggests an α = 3/2, β = 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 enc for NFW (UCMHlike α = 9/4, β = 3/4) halos of virial mass M vir = 10 −7 M , 10 −1 M and concentration parameters c = 100, 10 8 and the PBH-limit, c → ∞. For a given subhalo, if r enc passes below a particular sensitivity curve at a mass, M * , then the search is sensitive to an effective subhalo of mass M * . These M * mass subhalos make we obtain, where f PBH (M * ) 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, The sources considered by microlensing experiments range from nearby stars [26,27], the Large Magellanic Cloud (e.g., [23,24]), Andromeda (e.g., [25]), and even distant supernova [29]. For stellar sources, the distances are O(1−100 kpc) while supernovae are sensitive to much larger distances, O(Gpc). The typical Einstein radii of these sources are, 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. [32]) are subdominant to the effects considered here.
Finally, we use Eq. (39) 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 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 10 −9 M . For a CDM-inspired [65] 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, e.g. c = 10 3 , sensitivity to f < 1 occurs similarly to c = 100 for low mass halos, but also for M 10 −3 M for a static Doppler search.
While our analysis of diffuse halos is schematic and suffers from 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 [50].

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 to 100 M 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 DMdense 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. [48,49] 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 [50]. 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 ∼ 10 −9 M for the Doppler dynamic search and below ∼ 10 −4 M 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 [48]. We leave further study of this limit to future work [50]. 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.
As in the Doppler signal we can rewrite the magnitude of the signal shape as We have again defined, x ≡ (t − t 0 )/τ . Inserting this into (A7) gives the signal shape,

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 × p th percentile of B min ≡ min (B 1 , B 2 , ..., B N ). Every statistic of B min can be calculated with its cumulative distribution function (CDF), F Bmin , and because the B i 's are independent and identically distributed random variables we can solve for F Bmin solely in terms of F B : The 100 × p th percentile of B min is calculated by solving F Bmin (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.

a. 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, For a sphere of radius R, we multiply b p by R and take N = 4πN P nR 3 /3, giving b p = − 3 4π ln (1 − p) N P n 1 3 .

b. 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, To find N we calculate how many objects are in the dynamic limit. The dynamic limit is defined by T −τ > t 0 > τ . 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πvT R 2 /3. Multiplying (B4) by R and substituting in N gives Substituting p = 0.9 and dividing by v gives Eq. (22).

a. 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, Taking N = N P nπdR 2 , while multiplying by R to for a cylinder of radius R, gives Taking p = 0.9 gives (18).

b. 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, The dynamic limit is defined by T − τ > t 0 > τ . 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 nvT dD/2. Multiplying (B8) by D, for a cylinder of diameter D, and substituting N gives Taking p = 0.9 and dividing by v gives Eq. (25).