Probing primordial black holes from a first order phase transition through pulsar timing and gravitational wave signals

In this work, we assess the sensitivity reach of pulsar timing array (PTA) measurements to probe pointlike primordial black holes (PBHs), with an extended mass distribution, which originate from collapsed Fermi balls that are formed through the aggregation of asymmetric U(1) dark fermions trapped within false vacuum bubbles during a dark first order phase transition (FOPT). The PBH formation scenario is mainly characterized by the dark asymmetry, strength of the FOPT, rate of FOPT, and the percolation temperature. Meanwhile, for PBH masses of interest lying within $10^{-10} M_\odot - 10^{2}M_\odot$, the relevant signal for PTA measurements is the Doppler phase shift in the timing signal, due to the velocity change induced by transiting PBHs on pulsars. Taking the dark asymmetry parameter to be $10^{-4}$ and $10^{-5}$, we find that percolation temperatures within the 0.1-10 keV range, FOPT rates above $10^3$ times the Hubble parameter at percolation, and FOPT strengths within $10^{-6}-0.1$ can give rise to PBHs that can be probed by an SKA-like PTA observation. On the other hand, the accompanying gravitational wave (GW) signal from the FOPT can be used as a complementary probe, assuming that the peak frequency lies within the $\mathcal{O}(10^{-10})-\mathcal{O}(10^{-7})$ Hz range, and the peak GW abundance is above the peak-integrated sensitivity curves associated with pulsar timing observations that search for stochastic GWs. At the fundamental level, a quartic effective potential for a dark scalar field can trigger the FOPT. By performing a parameter scan, we obtained the class of effective potentials that lead to FOPT scenarios that can be probed by SKA through pulsar timing and GW observations.

Recently, a novel scenario incorporating FOPT and dark sector particles has been proposed [14][15][16][17], in which PBHs are formed from the collapse of macroscopic intermediate states called Fermi balls (FBs).Under this framework, which will be the main focus of this study, PBHs with sub-solar mass and a wide range of abundance can be generated.We briefly describe the physical mechanism behind the formation of PBHs in this scenario.The FOPT is induced by the finite temperature potential of a dark scalar φ, when the dark sector temperature drops below the critical temperature T c .The dark fermions, denoted by χ, need to be trapped in the false vacuum during FOPT, so that the expanding true vacuum bubbles squeeze χ to form FBs. A way to realize the trapping of dark fermions is by introducing a Yukawa interaction of the form g χ φ χχ; when the φ gets a nonzero vacuum expectation value (VEV), v φ , it lifts the mass difference of χ between the true and false vacua.Further requiring g χ v φ T c would trap χ particles in the false vacuum and form FBs. Once the range of the Yukawa interaction, set roughly by the inverse of the thermal mass of φ, becomes comparable with the separation distance between the bound χ particles in the FB, it triggers the collapse of FBs into PBHs.Some constraints on the properties of PBHs can be obtained by considering limits from Hawking radiation and gravitational lensing [5].For M PBH 10 −16 M , Hawking radiation becomes significant, and light particles from PBH evaporation contribute to cosmic ray fluxes in the present epoch.For M PBH 2.5 × 10 −19 M , the lifetime of such PBHs is shorter than the age of the Universe, but they can nevertheless inject entropy into the thermal bath through PBH evaporation, and may introduce distortions in the cosmic microwave background (CMB) and affect Big Bang nucleosynthesis (BBN).On the other hand, for 10 −11 M M PBH 10 M , gravitational microlensing observations, e.g.Subaru-HSC/OGLE-IV/EROS-2 [18,19], are able to set stringent limits on the PBH abundance.
In this work, we consider pulsar timing arrays (PTAs) as a complementary probe of PBHs within the window 10 −10 M M PBH 10 2 M .A PTA consists of a catalog of millisecond pulsars (MSPs), which are astrophysically relevant due to their extremely accurate timing signals.This property of MSPs has been exploited to measure gravitational waves (GWs) in the nanohertz frequency band, by considering a PTA as an interferometer [20,21].However, the timing signal of an individual pulsar in a PTA can be affected by the gravitation of the surrounding PBHs, mainly through Doppler and Shapiro effects [22][23][24][25][26][27].The former happens when a PBH passes through the vicinity of either a pulsar or the Earth, and induces acceleration; as for the latter, the Shapiro effect is due to the time delay in the propagation of the signal when a PBH crosses the line-of-sight between a pulsar and the Earth.To obtain limits on PBH properties, we perform a Monte Carlo simulation to assign masses, positions, and velocities on PBHs, based on the PBH mass and phase space distributions.The amplitude of the phase shift, relative to the intrinsic pulsar timing signal can be used to compute a statistical quantity called the signal-to-noise ratio (SNR).For a single mock simulation, the representative SNR is given by the maximum SNR across all pulsars.The sensitivity of a PTA to a PBH formation scenario is then determined by checking if the maximum SNR is above a threshold value for 90% of all mock simulations performed.
In the case of PBHs formed during a dark FOPT in the early Universe, the PBH mass distribution and abundance can ultimately be derived from the finite-temperature quartic potential and particle content in the dark sector.The scalar potential determines the percolation temperature, FOPT strength, and FOPT rate, which directly impact the PBH mass distribution and abundance.In addition to the Doppler and Shapiro effects induced by PBHs, another useful probe of the FOPT are GWs with frequencies between 10 −10 Hz to 10 −6 Hz, which can potentially be studied in future observations such as SKA [28], THEIA [29], and µAres [30].Similar to the derivation of PTA limits, we calculate the associated SNR of a GW signal from a FOPT over the course of a 20-year observation period.As elaborated in detail in [31], a convenient method to calculate the GW signal SNR is through the use of peak-integrated sensitivity curves (PISC).Under this framework, the relevant parameters are reduced into the peak GW abundance Ω peak h 2 and peak frequency f s ; the PISC for a specific experiment and observation time simply corresponds to an SNR equal to unity.The FOPT parameters, to which PTA is sensitive, provide predictions for Ω peak h 2 and f s ; the SNR can be easily calculated as the ratio of the predicted peak GW abundance to the PISC evaluated at the predicted peak GW frequency.This paper is organized as follows.We revisit the method of obtaining sensitivities on the phase shift of the pulsar timing signals via Doppler and Shapiro in Sec. 2. In Sec. 3 we apply this method, as an illustrative example, to the case of monochromatic PBHs and obtain limits on the PBH fraction and PBH mass.Secs. 4 and 5 cover mainly the PBH formation scenario via dark FOPT, which includes a discussion of the extended PBH mass function and mixed PBH-FB scenario, and the pulsar timing sensitivities projected on the FOPT parameter space.We calculate the limits from GWs produced during a dark FOPT in Sec. 6.We specify the quartic potential in Sec. 7, which realizes the FOPT, and we use it to investigate the correlation between PTA and GW signals.Finally, we summarize the results in Sec. 8.

Pulsar timing as probe of substructure
In this section, we shall provide a detailed discussion, following closely references [25,26], of pulsar timing as a way to determine the properties of a collection of transiting substructures.It is well-known that pulsars emit periodic signals, in the form of electromagnetic waves, that are extremely stable over extended time periods; the stability of pulsar signals is comparable with Earth-based atomic clocks.The pulsar signal arriving on Earth at time t can be characterized by its phase φ(t), which is fitted with a polynomial function, referred to as the timing model, given by where ν is the frequency of the emitted pulsar signal, and ν/ν is the spin-down rate.On the other hand, a transiting object in the vicinity of the pulsar, and/or close to the Earth-pulsar line of sight, can induce a phase shift to the pulsar emission.In Eq. (2.1), the signal from such transiting objects will manifest starting at third-order; the additional contribution to the linear and quadratic terms can be reabsorbed into the fit to the pulsar period and spin-down rate [23].In this study, we shall focus on the Doppler phase shift induced by the passage of compact objects; specifically for this case, we take the dominant signal to come from the velocity shift of the relevant pulsar, due to the acceleration on the pulsar caused by the transiting compact object.Another effect, known as the Shapiro time delay, changes the propagation time of the signal along the line of sight due to the gravitational potential sourced by the transiting object.In the case of transiting PBHs, previous literature, e.g.[23][24][25][26], has shown that the Doppler signal can be used to constrain monochromatic PBH masses as low as ∼ 10 −11 − 10 −10 M , whereas the Shapiro signal loses sensitivity in that mass range.As for the frequency shift induced by the the Doppler and Shapiro effects, it can be shown that they are respectively given by [24] where Φ is the potential generated by the compact object, d is the unit vector pointing from the Earth to the pulsar, v is the velocity of the compact object.We take the assumption of a constant velocity trajectory, r = r 0 + vt: this is a reasonable assumption since altering the velocity of the object requires it to be extremely close to the pulsar, and such strong-field events are rare [23].
For the Shapiro contribution to the frequency shift, the integration in z is performed along the line-of-sight, and the gradient of the potential is understood to be evaluated at each point along the line-of-sight.Given a shift in frequency δν(t), the corresponding phase shift is (2.4) We then sift out the signal amplitude from the phase shift that cannot be reabsorbed into the intrinsic phase signal from the pulsar.This can be done by subtracting away from δφ(t).For a collection of N compact objects, each of which is labelled by an index i = 1, ..., N , the so-called residual signal amplitude h I (t) at the I th pulsar is The signal-to-noise (SNR) ratio, a quantity which we will eventually use to assess the sensitivity of pulsar timing in probing the presence of transiting substructures, is evaluated using h I (t).This is given by (2.7) In the above expression, ∆t is the cadence of the pulsar timing measurement, which refers to the time interval between pulsar timing measurements; T obs is the observation time of the measurement; and t rms is the timing residual of the timing measurement, which refers to the uncertainty in the pulsar timing data.For pulsar timing observations comparable with the reach of SKA, we have ∆t = 2 weeks, t rms = 50 ns, T obs = 20 years [32].The overall SNR is taken to be the maximum SNR across an array of pulsars, so that The SNR, obtained from the residual signal, is completely determined from the positions and velocities of the compact objects, and the fixed positions of the pulsars.
In determining the PTA sensitivity on the distribution of compact objects/substructures in the Galaxy, we shall adopt a frequentist, statistical approach where the initial positions and velocities of the compact objects are drawn from a phase space distribution, which can be factorized into a spatially-dependent part and a velocity-dependent part.We assume that the velocity distribution is a truncated Maxwell-Boltzmann distribution of the form ) where v rms = 325 km/s, v esc = 600 km/s [24], and Meanwhile, the spatial volume distribution of compact objects is assumed to be uniform within some simulation volume; in the Solar neighborhood, this can be taken to be a fraction of the local DM density ρ DM = 0.46 GeV/cm 3 .Ultimately, this makes the SNR a random variable.In the absence of a signal, the SNR follows a one-sided Gaussian [24]; assuming that the SNRs from each pulsar are independent of each other, the joint probability density is just a product of these one-sided Gaussians.Given a threshold p-value, which gives the probability that the alternative hypothesis-in this case, the presence of a residual signal in the phase shift-occurred by random chance, the threshold SNR value, denoted by SNR * , is where N pul is the number of pulsars in the array.Then we have so that For N pul = 200 pulsars, the threshold SNR is SNR * = 3.66 ≈ 4 for p = 0.1.We then construct a sufficiently large number of mock "Universes," taken to be 1000, where we specify the initial positions and velocities of the compact objects, keeping the pulsar positions fixed; if a signal, i.e. when the SNR is at least 4, exists in at least 90% of randomly generated "Universes," then PTA measurements are sensitive in probing the substructure under consideration.
3 Sanity checks: pointlike, monochromatic PBH mass In the limit where the transiting object is pointlike, e.g. in the case of primordial black holes (PBHs), and in the case where the velocity of the transiting object is regarded as a constant, the frequency shifts simplify into (see Appendix A for more details) Here, b is the unit vector pointing from the pulsar to the closest distance approach of the transiting object.In either case, x refers to a rescaled time variable given by The quantities t 0,D and t 0,S refer, respectively, to the time at which the object reaches the closest distance to the pulsar, and the closest distance to the Earth-pulsar line of sight.The time variable τ , which exists for each signal, can be thought of as the characteristic width of the signal.
In the above expressions, we ignore constants of integration, which will simply be reabsorbed by φ 0 in Eq. (2.1).For a given monochromatic PBH mass M P BH , and PBH mass fraction f P BH , one can obtain the local number density of PBHs, which will determine the size of the sampling volume where the locations of the PBHs will be drawn.Note that the shape of the sampling volume will be different for the Doppler signal and for the Shapiro signal.This can be understood by noting that the signal width for the Doppler signal depends on the closest distance to the pulsar, which makes a spherical volume centered at the pulsar a natural shape for the simulation volume.The radius of this spherical volume is 6.099 × 10 −6 kpc N DM 100 The radial position of a PBH within this volume can be specified by taking a random number p i from a uniform distribution on the unit interval, and this is given by Analogously, the Shapiro signal width depends on the closest distance to the Earth-pulsar line of sight; since the object must transit between the Earth and the pulsar, the natural shape of the simulation volume is a cylinder coaxial with the line of sight, with a length equal to the Earth-pulsar separation L. The radius of this cylindrical volume is so that a randomly assigned initial distance from the line of sight can be taken to be where, analogously, p i is taken from a uniform distribution on the unit interval.Fig. 1 shows the result of our Monte Carlo simulation that we implemented using the procedure outlined in Section 2, where we fixed the number of samples within each simulation volume, associated with a single pulsar, to 100.Focusing on the Doppler signal sensitivity curve denoted by the red solid curve, we see that the PBH fraction stays flat for PBH masses of at least ∼ 10 −6 M , while lower PBH masses correspond to weaker sensitivity in the PBH fraction.This trend can be understood by considering the limit in which the PTA sensitivity is determined by a single PBH that yields the smallest signal width τ , where we can ignore coherence/cancellation effects, which occur upon summing over the phase shifts from each PBH, due to the angular coefficients d • b and d • v in Eq. (3.4).As discussed in [24], for low PBH masses, the minimum sample signal width is much less than the observation time scale, which falls under the dynamical Doppler regime.In this case, the limit on the PBH fraction scales as the inverse of the PBH mass, mainly due to the fact that lighter PBHs must be sufficiently close to the pulsar to generate a large SNR.On the other hand, since the required simulation volume becomes larger for heavier PBHs, the typical distance of PBHs from pulsars would be comparable to or larger than the distance traversed by the PBH during the observation period, and this falls under the static Doppler regime.From dimensional analysis, the amplitude of the Doppler phase shift at third order in t is proportional to G N M P BH νT 3 obs /(v 2 cτ 3 D ); taking the smallest τ D to be b min /v, where b min is the smallest impact parameter among PBHs in the sample, one can show that b min scales as the PBH number density proportional to (M P BH /f P BH ) 1/3 , and hence the dependence of the Doppler phase shift on the PBH mass disappears.As a side note, we observe from Fig. 1 that gravitational microlensing provides the most stringent limits for the PBH mass range of interest; on the other hand, pulsar timing provides a complementary probe of these pointlike objects, and the sensitivity of pulsar timing can be improved by, e.g., increasing the observation time or increasing the number of pulsars in the catalog [24].Furthermore, the Doppler signal sensitivity offers better discovery reach than microlensing for PBHs heavier than solar mass.
In the left panel of Fig. 2, we compare the result of a full MC simulation with simply taking the PBH sample, for each pulsar, with the smallest signal width.We expect fairly comparable results between them, based on the idea that the PBH with the smallest τ will provide the dominant contribution to the Doppler phase shift; in the actual plot, there appears to be an O(1) difference in the results because the full MC simulation considers a collection of PBHs, which can introduce some coherent enhancement in the amplitude of the Doppler phase shift.Meanwhile, we investigated the behavior of the fluctuation in the PBH fraction for M P BH 10 −6 M , when we change the number of mock universes in the MC simulation.Taking MC simulations with 500, 1000, 5000, and 10000 randomly generated universes, we found that the mean squared fluctuation σ 2 of f P BH within the considered mass range, decreases with increasing number of mock universes as listed in Table 1.This is an indication that the fluctuations in f P BH are a result of the random nature of the MC simulation itself.The left panel is a comparison of the Doppler PTA sensitivity limits from a full MC simulation (solid), and the result of simply taking the PBH sample that yields the smallest signal width τ (dashed).We both see similar trends in the PBH fraction-PBH mass sensitivity.In the right panel, we show the fluctuation in the limit on the PBH fraction about the mean value, for different numbers of mock universes in the MC simulation.The variance in the PBH fraction decreases with increasing number of mock universes.Table 1.

No. of runs f P BH
Tabulated values of the average PBH fraction and the variance, in the mass range 10 −6 M to 10 1 M , of the PBH fractions plotted in the right panel of Fig. 2.

PBH from a dark first-order phase transition
We now consider a PBH formation mechanism that can eventually give rise to a PBH mass function, which goes beyond the case of a monochromatic PBH mass distribution.Following [33], we consider the onset of a first-order phase transition (FOPT), during the radiation dominated era of the early Universe, within some dark sector.This is triggered by the presence of a dark sector scalar field φ with an associated quartic potential given by V ef f (φ, T ), where T is the temperature of the dark radiation.In general, the visible and dark sector temperatures, which are respectively denoted by T SM and T , can be different; it is then convenient to introduce the temperature ratio parameter which can be time-dependent.Assuming that there is no exchange of energy and entropy between the two sectors, and mass thresholds are irrelevant, we shall take ξ to be roughly constant in the relevant regime of interest.
The VEV of the scalar field φ is the order parameter characterizing the phase transition.At some critical temperature T = T c , the VEV abruptly jumps from the false vacuum (FV) at φ = 0, to the true vacuum (TV) at φ = φ + .In a FOPT, quantum tunneling allows the transition from the false to the true vacuum that are separated by a potential barrier.A bubble of the TV will form within a Hubble volume at the nucleation temperature T = T n < T c , corresponding to a time t = t n such that [34] where Γ(t) is the nucleation rate.The quantity Γ(t) determines the rate at which TV bubbles form, within a comoving volume and within a comoving time.Given the threedimensional instanton action S 3 (T ), the nucleation rate is As time progresses, the dark sector plasma will be populated by TV bubbles, and both phases coexist.For some comoving time t larger than the comoving time t c at which the critical temperature is reached, the comoving radius of a TV bubble, formed at t i , is where v w is the bubble wall velocity.The exact dynamics of bubble wall expansion is beyond the scope of our study (see, e.g.[35][36][37][38]), but we shall assume that the TV bubbles quickly reach their terminal velocity due to the friction from the surrounding plasma.The fraction of the Universe in the FV is then [39] In the last line, we note that the scale factor a(t) goes as t 1/2 in the radiation dominated era, so one can show that the physical size of the TV, formed at Assuming that the phase transition is relatively short, we simply have From the effective potential, we can define useful thermodynamic quantities.The free energy density difference between the FV and TV phases is the corresponding energy and entropy densities are respectively given by The pressure is then Assuming that the scalar field can be described as a perfect fluid, the trace of the energy momentum tensor is Given the trace of the energy momentum tensor, we can introduce which measures the strength of the phase transition; here ρ R is the total radiation energy density, coming from the relativistic species in the dark and SM sectors, i.e.
ρ R (t) = ρ v (t) + ρ d,eq (t), (4.13) We can also define the strength of the phase transition normalized with respect to the dark sector radiation energy density, given by which differs from α tr by a factor In addition to the dark sector scalar field, we introduce a dark sector Dirac fermion χ, possessing some global dark U(1) charge, and this may be regarded as a DM component.One can write down a renormalizable Lagrangian associated with φ and χ, which respects this global U(1) symmetry, and it is given by In this study, we take the fermion mass parameter m χ to be zero.Meanwhile, the trilinear Yukawa term serves a dual purpose of introducing an attractive interaction between dark fermions, and also provides additional mass to fermions within TV regions, equal to g χ φ + (T ), during the onset of the phase transition.We shall assume that during the relevant period of FOPT, there is no change in the effective relativistic degrees of freedom, and in the dark sector, we have If we take the case where g χ φ + (T ) T , the DM particles lying in the FV do not have sufficient kinetic energy to penetrate through the bubble wall; this process is referred to as DM filtering [40].As TV bubbles increase in size and population, the DM particles are trapped in pockets of FV bubbles that eventually shrink, forcing the DM particles to condense into Fermi balls (FBs); similar ideas were explored in, e.g.[41].Once the FBs form, no light degrees of freedom remain in the dark sector.Note that because the process of pair annihilation of dark fermions into scalars is quite efficient in depleting them, it is necessary to introduce an asymmetry in the number of DM particles relative to DM antiparticles; the asymmetry is denoted by the parameter where is the entropy density of the visible sector, which is conserved assuming that no entropy nor energy exchange occurs between the visible and dark sector.The time t = t * > t n refers to the percolation time, where the FV filling fraction is taken to be f F V (t * ) = 0.29 [42].The corresponding number density of FV bubbles with radii lying between R r and R r + dR r is [16,39] where β, having the same units as the Hubble parameter, is roughly the inverse of the time duration of the FOPT, and v w is the bubble wall velocity.In obtaining Eq. (4.21) we assumed that the TV nucleation rate can be written as which implies that the the FV filling fraction, I(t), can be written as In obtaining Eq. (4.23), it is assumed that the saddle point approximation applies in evaluating Eq. (4.6).
As the TV bubbles expand and increase in number during the FOPT, the regions containing the FV shrink and will eventually be completely surrounded by the walls of the TV bubbles.For this reason, the FV bubbles are not exactly spherical: quantitatively, the asphericity of a FV bubble with radius R * is captured by the parameter A, such that the volume of such a FV bubble is A 4πR 3 * /3.This parameter is directly tied to the FV filling fraction through the following normalization condition: The total U(1) charge in a FV bubble of radius R * = R r (t * ) is given by [16,43] Subsequently, the FV bubble shrinks to form a FB containing the same amount of dark U(1) charge.While the individual dark fermions themselves are massless, the mass of the FB is a combination of the Fermi gas kinetic energy, the FV energy, and a negative contribution from the attractive Yukawa interaction.The total energy of a FB is a function of its radius, and the resulting FB configuration, at a given FB temperature, corresponds to the local energy minimum.In the limit where ∆V ef f (T * ) T 4 * , the FB mass is [16,17] where ∆V ef f (T * ) is the energy difference between the false and true vacua, evaluated at the dark sector percolation temperature, denoted by T * ≡ T (t * ).Using Eq. (4.25), Q F B can be traded with R * , so that a relation between the FB mass and R * can be established using Eq.(4.26).
Due to the attractive Yukawa interaction, the FBs can collapse into PBHs.As the FB cools down, the range of the Yukawa interaction increases in the case where the dark scalar mass scales as the square of the temperature; at some collapse temperature T cl , the Yukawa contribution to the FB energy becomes significantly large and the FB configuration ceases to be stable [16,44].Typically, since T cl is extremely close to the percolation temperature, and because of the mild temperature dependence of the FB mass according to Eq. (4.26), the PBHs inherit the mass of the FB progenitor.From Eq. (4.21), one can obtain the corresponding normalized probability distribution of PBH masses produced through this mechanism; this is given by where M is the mass of a PBH and the relationship between R * and M can be obtained from Eqs. (4.25) and (4.26), and identifying M with M F B .To obtain the SM entropy density, we use [45] to evaluate g * s,v (T SM ) ≈ g * ρ,v (T SM ) for a given visible sector temperature T SM .
In implementing the Monte Carlo simulation to calculate the PTA sensitivity curve, we need to construct sample universes by sprinkling each mock universe with PBHs, with the added feature that we have to assign masses to the PBHs in accordance with the distribution provided in Eq. (4.29).Sampling the distribution P (M ) can be done by first taking a random number w from the unit interval [0, 1), and then obtain M such that 1−w = c(x(M )), where c(x(M )) is the complement of the cumulative distribution function (CDF) c(x(M )) corresponding to P (x), such that P (M ) dM = P (x) dx, where The CDF integral is numerically stable when x(M ) I * .For x(M ) I * , the numerical evaluation of the complement function is more stable compared to the numerical integration of the CDF integral.By using the identity For exp(−xI * ) < 10 −7 , cres is at most 1% of c, so we can simply drop cres ; this is demonstrated nicely in the right panel of Fig. 3. Before resorting to a full Monte Carlo simulation, we can estimate the FOPT parameter region that is sensitive to pulsar timing array measurements, by first obtaining the dependence of the average PBH mass and PBH fraction on the FOPT parameters.In most cases, the temperature range at which the phase transition occurs is much smaller than T c , so that we can perform a linear expansion of ∆V ef f around T = T c , i.e.
where c ≡ (T c ) can be interpreted as the latent heat that can be released to the plasma during the phase transition.The strength of the FOPT transition can be estimated as and, in turn, we can express ∆V ef f (T * ) as Furthermore, we can focus our attention to parameter regions which lead to reheat temperatures below the critical point, so that ρ Typically, the lowest temperature T 0 , at which the potential still develops a barrier between the false and true vacua, does not stray away from T c .If we allow T * > T 0 ≥ 0.8T c , then this gives a maximum possible upper bound for α tr of 0.4; in other words, the criterion that we do not reheat the dark plasma close to T c implies that α tr is at most O(10 −1 ).
To simplify our initial analysis, we shall take the following choices for the FOPT parameters: For the temperature mismatch parameter ξ, the above choice can be shown to be consistent with constraints on N ef f ; at the epoch of Big Bang nucleosynthesis (BBN), ∆N ef f ≤ 0.4 [46].Firstly, if the phase transition had occurred before BBN, the absence of light degrees of freedom in the dark sector makes the BBN constraints on N ef f irrelevant; otherwise, the contribution of the dark sector to Here we are assuming that the Chapman-Jouguet condition holds [47], where the bubble walls expand like spherical detonation waves; note that, since the formation of the bubbles occurred in the dark sector plasma, it is appropriate to evaluate the bubble wall velocity at α d , defined in Eq. (4.15), rather than α tr .However, it has been pointed out by [48] that bubble wall expansion in the context of cosmological phase transitions can be characterized by a wider class of detonation scenarios.In [49] they provided contour plots for the bubble wall velocities in terms of the particle physics model parameters that trigger the FOPT.In this work, however, we are only interested in general trends in the limits on the FOPT parameter space, and we shall adopt the simplest approach where we stick with Eq. (4.44).Meanwhile, the above fixed choice for T * /T c implies that α tr 0.14; note that this upper limit on α tr is a matter of choice based on Eq. (4.41).Then the average PBH mass is given by: where At percolation time, the total energy density of PBHs is Relative to the background cosmological cold DM (CDM) density, the total energy density of PBHs at percolation time is Using Eqs.(4.47) and (4.52), we can then obtain the projection on the f P BH − M P BH plane, for certain choices of the FOPT parameters.This is shown in Fig. 1, where we project out contours of constant β/H * = 2.5 × 10 3 , 5.0 × 10 3 , 10 4 , and 2.5 × 10 4 , for DM asymmetry parameters η χ = 10 −5 and 10 −3 .Note that larger asymmetry parameters lead to larger dark U(1) charges in FBs, implying heavier PBHs.However, to maintain the same PBH fraction, an increase in η χ must be compensated by a decrease in the percolation temperature.The dots in Fig. 1 correspond to different values of the temperature at percolation T * ; on a contour with fixed β/H * , the percolation temperature decreases in the direction of increasing average PBH mass, since M P BH scales with the percolation temperature as T −2 * .Meanwhile, smaller values of β/H * lead to larger FV bubble radii, which then correspond to heavier PBH masses; this trend can also be seen in the contours in Fig. 1.The locations of the dots in the contours of Fig. 1, relative to the PTA sensitivity curve obtained by assuming a monochromatic PBH mass distribution, indicate that PTA observations with similar properties as SKA would be sensitive to PBHs formed through a FOPT, with an associated percolation temperature T * lying around O(1) keV, with a FOPT rate around O(10 3 − 10 4 )H * , assuming that the DM asymmetry is 10 −5 ; for η χ = 10 −3 , T * must be around O(10) eV.
Alternatively, we can invert Eqs.(4.47) and (4.52) to determine the FOPT parameters in terms of M P BH and ω P BH, * .We have where c 0 ≡ 1 192 (A 4π/3/f F V (t * )) 2 ∞ I * dx x 3 e −x (1 − e −x ) ln 6 (x/I * ) 70.For fixed percolation temperature T * and DM asymmetry η χ , the estimates of the FOPT parameters α tr and β/H * given in Eqs.(4.54) and (4.56) are useful to roughly pinpoint the parameter region where pulsar timing is expected to be sensitive, even without performing a full Monte Carlo simulation of the PBHs.In Fig. 4 we display the contours of constant percolation temperature T * on the α tr − β/H * plane, obtained from Eqs. (4.54) and (4.56), which correspond to the expected sensitivity reach of an SKA-like PTA observation to measure Doppler signals produced by monochromatic PBHs of mass M and fraction ω P BH * .For a given T * , the FOPT parameter points below the projected sensitivity curve will give rise to PBHs that provide a maximum SNR above threshold for the Doppler signal.The vertical dashed line corresponds to the maximum α tr , for a fixed T * /T c , beyond which the dark sector plasma will be reheated to the critical temperature and the false and true vacua will undergo a prolonged period of phase coexistence [50].Meanwhile, the horizontal line at β/H * = 10 marks a convenient lower bound of phase transition rates that are faster than the Hubble expansion at percolation.To compare the results with a full Carlo simulation, estimates for the viable FOPT parameters from 4, and determine the fraction of mock universes f that produce a maximum SNR above threshold.We expect this fraction to be close to 0.9, so we define a parameter ∆ ≡ f /0.9 − 1 which quantifies how good the estimates are in representing the PTA sensitivity limit in the FOPT parameter space.Looking at right hand side of each panel in Fig. 4, we find that |∆| is at most 5%; this rather satisfactory agreement between the estimates and the full MC result can be attributed to the shape of the PBH mass distribution being peaked around M , so the PBH mass distribution can roughly be characterized by a monochromatic distribution centered at M .Notice also that the projected curves are similar in shape as the corresponding PTA sensitivity curves: the branch where the PBH fraction is fixed is mapped into a line of constant α tr , while the low-mass branch, where the PBH fraction becomes larger, is mapped into a rising branch in the α tr − β/H * plane.We see that the projected curves, for η χ = 10 −5 (10 −4 ), fall on the desired FOPT parameter region, i.e. α tr less than the maximum value to avoid reheating the dark plasma to T = T c , and for β/H * > 10, for percolation temperatures ranging from 1-10 keV (0.1-1 keV).

Mixed PBH-FB scenario
In the previous section, we considered the scenario in which all of the FBs collapse into PBHs.We assume that the intermediate step of FB formation occurs, so we require that, at percolation time, the range of the Yukawa interaction L φ M −1 φ (T * ) is less than the average interparticle distance, which scales as R * /Q 1/3 F B , between χ within the FBs; otherwise, the FV bubbles containing trapped DM particles will immediately collapse into PBHs [16].If FBs are formed, a necessary condition for the collapse to occur is a mechanism that allows FBs to cool down, which can be achieved through the emission of a light particle [16,44], or the emission of neutrinos that can carry away energy from the FB.The latter process, mentioned in [51], may be realized by assuming a portal interaction between the scalar φ, introduced in Eq. (4.17), and the Higgs.If the emitted particles are relativistic, then the energy loss rate scales as T 4 according to the Stefan-Boltzmann law; furthermore, if the cooling rate is faster than the Hubble expansion, this ensures that the FB will inevitably collapse into PBHs [16,52] (see also [53] for a recent discussion on blackbody and evaporative cooling).On the other hand, this criterion is dependent on the model which describes the DM interaction with light particles; if the couplings are sufficiently weak such that, e.g. the cooling time is longer than the age of the Universe, then some of the FBs may not collapse into PBHs.In addition, the mixed PBH-FB scenario may be realized by noting that FBs lying on the low mass tail of the FB distribution will cool down faster than heavier FBs, and form PBHs. Hence, a possible final outcome might involve a situation where a fraction r of substructures from FOPT come in the form of FBs, and the remainder, (1−r), in PBHs.The exact fraction r appears to be dependent on the underlying mechanism or DM microphysics, so we shall treat r as a model independent parameter.The distributions for the FB radii and masses remain unchanged, and can be sampled by picking a random value of x using the distribution in Eq. (4.31); the additional feature in our setup is the probability FB, given by Accounting for the finite size of the FB, the phase signal is multiplied by a form factor as prescribed in Appendix A. For a spherically symmetric, constant volume density FB [54] with radius R F B (M ), the form factor can be obtained using Eq.(A.30), and we find that For a given substructure mass M , the form factor in Eq. ( 5.1) must be evaluated at the instant where the object is closest to the pulsar (in the case of Doppler), or where the object is closest to the line of sight (in the case of Shapiro).Note also that Eq. ( 5.1) tends to unity in the limit where x R F B .Hence, to determine whether a Doppler signal can be used to distinguish between FBs and PBHs in this scenario, we can compare the typical distance between the object and the pulsar, given in Eq. (3.6), versus the average radius of FV bubbles, R , which is an upper bound for the average radius of FBs, R F B .The average FV bubble radius is so that we have (5.3) The above estimate tells us that the typical size of FBs is much smaller than the typical pulsar-object distance, so a FB is effectively indistinguishable from a PBH from the point of view of a pulsar.Similarly, for the Shapiro signal, the relevant length scale is R S given by Eq. (3.8); compared to the typical FB size, we have (5.4) Hence, even the Shapiro signal cannot discriminate between a pointlike object and a finite size FB.

Gravitational wave production
During a FOPT, gravitational waves (GWs) can be produced, and these can serve as complementary signals that we can utilize to probe FOPT scenarios in the early Universe.It has been shown in [37] that the dominant contribution to the GW abundance from a FOPT comes from sound waves.The signal spectrum is therein) Similar to the prescription employed to evaluate the Chapman-Jouguet speed, the efficiency factor κ s , which is generally a function of the bubble wall velocity and phase transition strength, must be evaluated at α d since the FOPT occurred in the dark sector [55,56].
Assuming that the FOPT falls in the nonrunaway regime, we have where One of the usual methods employed to determine the sensitivity of GW experiments to GW signals produced, e.g. through a FOPT, is to check if the predicted GW abundance Ω sig (f )h 2 is above the GW sensitivity curve for a given experiment Ω noise (f )h 2 , for some frequency f within the operating frequency band [f min , f max ] of the experiment.Another method requires that the SNR, defined as is above some threshold value ρ th .Here, τ obs,GW is the expected runtime of the experiment, and n det = 1 for an autocorrelation measurement and n det = 2 for a cross-correlation measurement (e.g.aLIGO, aVirgo, KAGRA, CE, ET, LISA, DECIGO, BBO, PTA).However, reference [31] points out disadvantages of these two methods, and instead advocates the idea of introducing peak-integrated sensitivity curves (PISCs).The main assumption of this framework is the factorizability of the GW abundance spectrum into a peak value-a model dependent quantity-multiplied by a spectral function, whose form is model independent.
Then we can write As an illustration, one can take the case of GW produced through sound waves; the peak GW abundance is given by Eq. (6.2) the spectral Eq. ( 6.3).The ansatz in Eq. (6.8) allows us to write where (6.10) As defined in Eq. (6.10), the PIS only depends on the shape of the noise spectrum intrinsic to the GW experiment or mission, and on the form of the spectral function.Once we have the PIS for an experiment/mission and for a specific source of GW signal, the SNR can immediately be obtained from Eq. (6.9) once we have the model dependent quantities, namely the peak GW signal and the corresponding peak frequency, without repeating the integration Eq. (6.7) for each set of model parameters.
In our scenario where GWs are produced from sound waves during a FOPT in the dark sector, specifying the FOPT parameters will provide a prediction for the peak GW abundance and the corresponding frequency, as provided in Eqs.(6.2) and (6.4), respectively.In Fig. 5, we projected the FOPT parameters from Fig. 4, for T * /T c = 0.9, and η χ = 10 −4 and η χ = 10 −5 , corresponding to PBH formation scenarios that can be probed by PTA, on the Ω peak s h 2 − f s plane.The cross marks indicate the same set of average PBH mass values as in Fig. 4, where M P BH increases with f s .On the other hand, the PISCs corresponding to THEIA, SKA, and µAres are shown as dashed lines.Note that µAres is a proposed space-based GW antenna [30] that will be expected to operate in the µHz frequency band, with better sensitivity compared with LISA.Meanwhile, THEIA [29] is another proposed mission designed to perform high precision astrometry at an accuracy level that is 1000 times better than the ESA Gaia mission.The principle behind detecting stochastic GW signals using astrometry is based on the idea that such GW signals will produce correlated angular displacements of stars, and monitoring of the proper motions of a collection of stars can be used to place limits on the GW abundance [57,58].As for SKA [28], gravitational wave astronomy can be performed in an analogous manner, as laid down in, e.g.[59], by looking at the effect of GWs on the timing signal of a collection of pulsars, and the correlation between the changes in the timing signals of the observed pulsars can be used to infer the strength of the GW signal.For all PISCs, we assume a 20-year observation period; for any other choice for the observation period, the PISC must be rescaled by an appropriate factor that is proportional to τ 1/2 obs,GW .Conversely, a prediction for the GW SNR that is a factor of 10 larger than the SNR threshold for a GW experiment, for a given observation time τ obs , means that the required observation time to reach the SNR threshold is reduced by a factor of 100.
Looking at the PISCs in the top row of Fig. 5, THEIA provides the most stringent constraints on stochastic GW signals in the peak frequency range around 10 −10 − 10 −7 Hz, with SKA limits being competitive in the subnanohertz range.Beyond the microhertz THEIA will be sensitive to any FOPT scenario that yields a GW signal SNR for SKA above unity.For the projected FOPT parameter points, where T * /T c is fixed to 0.9, the corresponding SNR can be directly read off from Fig. 5 using Eq.(6.9).In the case of η χ = 10 −5 , we found that most of the parameter points along the PTA sensitivity curve, corresponding to percolation temperatures of 1-2 keV, can produce GWs that can be probed by SKA, assuming that the threshold SNR is set to 1; for percolation temperatures of at least 5 keV, SKA is not sensitive to the associated GW signal.Meanwhile, for η χ = 10 −4 , SKA is sensitive to FOPT scenarios with T * = 0.1 − 0.2 keV, and the SNR goes below threshold for higher percolation temperatures, as shown in the top right panel of Fig. 5.This trend in the SNR, attributed to the decrease in the peak GW abundance in the frequency range below ∼ 10 −8 Hz, can be understood by looking at Eq. (6.2): the behavior of the peak GW abundance Ω peak s h 2 is mainly driven by the decrease in α tr , which can happen as the result of increasing percolation temperature, according to the α tr estimate in Eq. (4.54).
In light of the recent discovery of low frequency gravitational waves by NANOGrav based on 15-year observational data [60], one can check that that NANOGrav is quite insensitive to the GWs produced in our mechanism, assuming the parameter range of interest.Given the strain sensitivity h c (f ) provided in [61], the corresponding GW abundance can be computed via, e.g.[58,62], We then calculated the PISC and display the resulting curve in Fig. 5, and one can see that the recent NANOGrav result has little to no direct impact on our setup.Focusing on the proposed THEIA experiment, the bottom row panels of Fig. 5 show the corresponding GW SNR for the same set of benchmark FOPT scenarios.We generally find a sharp change in sensitivity for β/H * in the 2-5 keV range for η χ = 10 −5 , and in the 0.2-0.5 keV range for η χ = 10 −4 .In this range of percolation temperatures, the peak frequency where the projected PTA sensitivity curve crosses the THEIA PISC is below ∼nHz where the THEIA PISC is decreasing.Beyond a certain percolation temperature, the peak frequency is above nanohertz where the THEIA PISC is increasing; this can only happen if α tr increases with the percolation temperature.From Eq. (4.54), one can increase the PBH fraction, which corresponds to lighter PBH masses in the PTA sensitivity curve.In turn, lighter PBH masses can be achieved if the FV bubble, which trapped the DM particles that formed the FB progenitor, has a smaller radius: this amounts to a fairly large dimensionless FOPT rate β/H * .This is also the reason why the projected PTA sensitivity curves that lie above the THEIA PISC correspond to the low M P BH tail, for higher percolation temperatures.
It is then worth emphasizing that for benchmark FOPT scenarios with DM asymmetry η χ = 10 −5 (10 −4 ), T * at 2 keV (0.2 keV) or below, and β/H * around the O(10 2 ) − O(10 4 ) range, these produce PBHs that can be probed through an SKA-like PTA observation of Doppler phase shifts; on the other hand, using SKA to look for GW signals from a FOPT will also be sensitive to these FOPT scenarios.This may provide a strong motivation for including the search for Doppler shifts in the pulsar timing signal as one of the goals of SKA, whether or not a stochastic GW signal can be detected by SKA, albeit a lack of a GW signal and a presence of a PTA signal would be inconclusive in validating the PBH formation scenario discussed in this study.It might also be possible that a GW signal will be found, but not a PTA signal: this may occur for FOPT scenarios with high percolation temperatures, which lead to light PBHs beyond the reach of pulsar timing.Certainly, a simultaneous detection of stochastic GW signals and Doppler shift in the phase signal would be an ideal situation in probing this FOPT scenario.

Generic quartic potential
At this point, we were able to identify the FOPT parameters, that produce PBHs that can be probed by PTA, and also generate GWs that can potentially be detected through correlated shifts in positions of a collection of stellar objects, or through a space-based GW antenna.In Sec. 4, we have introduced the effective potential V ef f for the dark sector scalar field φ, from which all relevant quantities that describe the FOPT originate.The goal now is to obtain the class of models which lead to the desired FOPT parameters that can be probed by PTA and stochastic GW observations.We can take a generic form for the effective potential V ef f following [17], so that The parameter T 0 , the destabilization temperature, can be traded with the energy gap between the false and true vacua at zero temperature, which we shall denote as B. One can show that the true vacuum configuration φ 0 at zero temperature satisfies the following equations so that we have where ζ(B) is the positive root to the polynomial equation From the effective potential, we can also extract the critical temperature T c at which the true and false vacua are degenerate in energy; we find that T c is the positive root of the equation Taking the fundamental parameters that characterize the effective potential to be we can identify the physical conditions that guarantee the existence of T 0 and T c .Requiring that we obtain positive roots from Eqs. (7.5) and (7.6), we have From the expression for the effective potential, Eq. (7.1), the VEV, energy gap between false and true vacua, and both their derivatives, are given by ∆V ef f = 1 3 The above expressions can be readily used to calculate the FOPT strength from Eqs. (4.11) and (4.12).Furthermore, obtaining the FOPT rate β requires knowledge of the nucleation rate provided in Eq. (4.3).For a certain class of quartic effective potentials for a single scalar field, [63] obtained a semianalytic expression for the O(3)-symmetric bounce action given by ) where Such quartic effective potentials are assumed to have two local minima, where each minimum can be associated with either the true/false vacuum, separated by a barrier.We require that λ > 0 so that the potential is always bounded from below.Note that 0 < δ < 2, where δ 2 corresponds to the thin wall limit [64], where the two potential minima are nearly degenerate in energy.The FOPT rate β at the percolation temperature T * is then given by The time-temperature relation is obtained from entropy conservation in the SM and the dark sector, assuming that they are only coupled through gravity.Since the duration of the phase transition is rather short relative to the Hubble expansion, and assuming that we are away from mass thresholds, we can assume that the effective number of relativistic degrees of freedom is constant in either sector.Furthermore, the dark and visible sectors do not exchange energy and entropy, so the entropy in each sector is separately conserved.This implies that the temperature ratio in the two sectors is constant, and this condition was invoked in Sec. 4. Noting that one can show that Fixing A = 0.1 and η χ = 10 −5 , we performed the scan over the parameters {B 1/4 , C, D, λ}.We restricted our scan within An initial summary of the results of the scan is shown in Fig. 6, where we plot one dimensional histograms of the physical observables in SKA, namely: the average PBH mass and PBH fraction; and the peak GW abundance and peak GW frequency.Starting from flat priors for the fundamental parameters appearing in the effective potential, the solid black histograms represent the resulting distributions of the physical observables.
We note that we only included parameter points in the regime where |T * − T c | < 0.1T c and α tr is sufficiently small such that the dark plasma will not be reheated back to the critical temperature.The scanned hypervolume leads to a fairly wide range of: average PBH masses, from ∼ 10 −22 M to ∼ 10M ; PBH fractions from ∼ 10 −2 to ∼ 10 2 ; peak GW abundances from ∼ 10 −31 to ∼ 10 −8 ; and peak GW frequencies from ∼ 10 −10 Hz to ∼ 10 −2 Hz.Imposing the criterion that we only include the FOPT scenarios that lead to Doppler phase shifts within reach by SKA (labelled as "SKA-PTA"), we find that the coverage of PBH masses and PBH fractions are significantly reduced as expected, since PTA is sensitive to PBH masses of at least ∼ 10 −10 M and PBH fractions of at least ∼ 0.3.On the other hand, recall from Eq. (4.47) that the average PBH mass goes as the inverse square of the percolation temperature and is weakly dependent on the FOPT strength α tr ; at the same time, a detectable signal in PTA and GW observations will require larger FOPT strengths at low percolation temperatures.Hence, the upper limit on the average PBH mass is a result of imposing a maximum α tr .As for GW observables, the f s histogram peaks at ∼ 10 −8 Hz, and covers a frequency band of ∼ 10 −9 to ∼ 10 −7 Hz, corresponding to the effective frequency of SKA.Focusing on the lower left corner of Fig. 6 where we the histogram plot for Ω peak h 2 , we see that can still generate a wide range of GW abundances even after selecting those FOPT scenarios that produce both PTA signals and PBH fractions below unity; these FOPT scenarios are below the reach of GW searches in SKA.Further applying the filter where we take those points that also give a detectable GW signal in SKA, this narrows down the range of the values of the physical observables; in particular, average PBH masses within ∼ 10 −8 to 10 −5 M can generate signals in PTA and GW observations using SKA.
In Fig. 7 we present the results of our parameter scan as corner plots, where we conveniently show two dimensional slices of the entire scanned hypervolume.The red curves roughly enclose the points which correspond to FOPT scenarios where the resulting PBHs can generate a detectable Doppler shift the pulsar timing signal from an SKA-like PTA Those points that also generate GW signals within the sensitivity reach of SKA are bordered by blue curves, and the subset points where the PBH fraction is below unity are bounded by the green curves.The points outside the red curves correspond either: (i) unphysical points which not yield values for T 0 and/or T c , according to Eq. (7.7), and hence do not lead to a FOPT; or (ii) points which are physical, but are beyond the reach of both pulsar timing probes of PBHs nor searches for stochastic GWs in SKA.We focused solely on filtering our scanned data points based on the sensitivity reach of SKA to both Doppler signals from transiting PBHs and GW signals from FOPT to highlight the potential of SKA to probe FOPT in the early Universe through these two observational signatures.Based on our scan, we can identify the class of effective potentials, which could come from some fundamental physics framework, that yield interesting phenomenological signatures which can be probed through SKA.Analogously, to provide us with insights on how the fundamental parameters the physical observables, we show a similar set of corner plots in Fig. 8.
We can identify a few interesting features from the corner plots of Figs. 7 and 8.In Fig. 7, we observe a lower limit for B 1/4 at ∼ 1 keV, which arises as a result of the physical conditions imposed in Eq. (7.7); similarly the white region in the range 10 −1 ≤ D ≤ 1 in the λ-D panel can be attributed to an unphysical region where no T c can be obtained.As for the PBH masses, PBH fractions, and GW frequency, we see in Fig. 8 that the main parameter that sets the scale for these physical parameters is B 1/4 ; note that we have traded B 1/4 with the destabilization temperature T 0 , so B 1/4 sets the temperature scale at which the FOPT occurs.Furthermore, we observe that the constraint on the PBH fraction below unity, for points that can be probed by SKA through pulsar timing and GW observations, restricts B 1/4 , up to ∼ 20 keV.This can be understood by looking at the upper left corner of Fig. 8, where the PBH fraction appears to be positively correlated with B 1/4 .Intuitively, this makes sense since a larger B 1/4 leads to a larger percolation temperature, which will increase the total energy density of the PBHs, leading to larger PBH fractions.Still in the same panel, the average PBH mass is anticorrelated with B 1/4 ; here, the average PBH mass is expected to decrease for higher percolation temperatures according to Eq. (4.47).Hence, pushing for larger values of B 1/4 will generally bring us to lighter PBHs beyond the reach of PTA.

Conclusions
In this paper, we studied the PBHs produced from the collapse of Fermi balls, which originate from DM particles filtered away from the walls of expanding true vacuum bubbles during a dark FOPT.Earlier work has shown that the generated PBHs follow a mass distribution, dictated by the distribution of the radii of false vacuum bubbles.In turn, the distribution is determined by the: percolation temperature; strength of FOPT; rate of FOPT; bubble wall velocity, assumed to be Chapman-Jouguet; and the DM asymmetry.
We determined the sensitivity reach of pulsar timing arrays to the Doppler phase shift produced by these transiting PBHs, taking the case of monochromatic PBHs as a reference scenario and taking the average PBH mass to be at 10 −10 M .Viable parameter provide values of the SNR, a functional of the phase shift, above a threshold set by the number of pulsars in the PTA. the PTA f P BH − M P BH curve in the case of monochromatic PBHs, and identifying monochromatic PBH mass with average PBH mass M from the PBH mass distribution, we were able to identify the corresponding FOPT parameters that are within reach by an SKA-like PTA measurement.As a cross check, we performed a Monte Carlo simulation of PBH masses sampled from the PBH mass distribution set by the viable FOPT parameters, and found that all of these parameter points lead to a situation where 90% of the mock simulations produce a Doppler signal with an SNR above threshold.This can be explained by the fact that the mass distribution is sharply peaked, and thus most PBHs in the simulation have masses that hover around M .Fixing η χ = 10 −5 and the ratio of the percolation temperature to the critical temperature to be 0.9, the rate of the FOPT ranges from β/H * 10 to 10 4 , covering a wide range of values for the FOPT strength α tr , from ∼ 10 −6 to 10 −1 , within the 1-10 keV temperature range.Since the PBH fraction goes as ∆V 1/4 ef f , and since the PTA sensitivity curve requires a fairly constant PBH fraction for PBH masses above 10 −8 M , a lower percolation temperature implies a lower radiation energy density, requiring a larger FOPT strength.The case of η χ = 10 −4 and percolation temperatures 0.1-1 keV leads to similar sensitivity reach in α tr and an order of magnitude increase in β/H * .For η χ = 10 −5 , and setting the percolation temperature to be fixed around the 1-10 keV range, we require at least β/H * ∼ 10 3 to produce average PBH masses below 10 −8 M .Below this PBH mass, PTA limits on the PBH fraction become weaker, which translates to larger values of the FOPT strength.
The same FOPT that produced the PBHs can also produce GW signals that can be probed through the correlated shifts of a collection of stars via pulsar timing.We concentrated our attention to GWs produced by sound waves since these provide the dominant contribution to the GW signal from FOPT.Their spectral shape is known and is determined by the peak GW abundance and peak frequency.We then obtained the peak-integrated sensitivity curves for THEIA, SKA, and µAres, i.e. the peak GW abundance as a function of the corresponding peak frequency corresponding to a GW signal with SNR = 1 and 20year observation time.We found that the corresponding peak frequency of the GW signal falls within the sub-nHz to sub-µHz range, and that THEIA offers the best sensitivity, while SKA offers sensitivity to our FOPT scenario that is comparable with THEIA.For projected SKA sensitivities to GWs, relevant FOPT rates lie in the range 10 2 β/H * 10 4 provides SNR > 1 for the benchmark cases η χ = 10 −5 , keV-scale percolation temperatures, and T * /T c = 0.9.
Finally we considered a class of generic effective quartic potentials that can realize the FOPT, from which we can directly calculate the FOPT parameters relevant for PBH formation and GW production.Performing a parameter scan over this class of effective potentials, fixing η χ = 10 −5 , we identified points that lead to PBH Doppler signals within reach by an SKA-like PTA measurement, and also those that lead to GW sound wave production that can be seen by SKA.We found that there is a significant region of the parameter space that provides a detectable Doppler phase shift signal in an SKA-like PTA measurement, stochastic GW searches in SKA may also cover.may provide further motivation to utilize the distortions in the pulsar timing signal, beyond the search stochastic GWs, to search for signatures of FOPT in the early Universe, most notably the PBH formation discussed in this work.It is hoped that shifts in the timing signal in PTAs will be used not just to search for stochastic GWs, but also substructures that could be present in the Galactic neighborhood.The results of our scans can also be useful for constraining explicit models with an additional scalar field that can be mapped to our generic quartic potential.
For the Doppler signal, τ D = b/v, where b is the closest distance between the pulsar and the transiting object.Meanwhile, for the Shapiro signal, τ S = b ⊥ /v ⊥ , where b ⊥ is the closest distance to the line of sight, and v ⊥ is the magnitude of the velocity component perpendicular to the line of sight.The phase shifts associated with the Doppler and Shapiro signals are then obtained by integrating the frequency shifts over time, and we get

10 − 10 Figure 1 .
Figure1.Sensitivity curves on the PBH parameter space, corresponding to Doppler signal, assuming SKA-like observational reach, and monochromatic PBH mass.We have included combined microlensing constraints for pointlike substructures from Subaru-HSC, OGLE-IV, and EROS-2 surveys, derived in[18,19].The dot dashed lines correspond to the projections of first order phase transition parameters (discussed in Section 4) on the PBH fraction and the average PBH mass, for DM asymmetry parameters η χ = 10 −5 (left panel) and η χ = 10 −3 (right panel), and bubble wall velocity fixed to the Chapman-Jouguet value of v w (α d ) 1 for α tr = 0.1 and α d ∼ 10 3 .In the left panel, the dot markers correspond to percolation temperature values 0.25, 0.50, 1.00, 2.00, 5.00, and 10.0 keV; in the right panel, we have 2.5, 5.0, 10.0, 20.0, 50.0, and 100.0 eV.Please refer to the text for more details.

10 − 10
Figure 2.The left panel is a comparison of the Doppler PTA sensitivity limits from a full MC simulation (solid), and the result of simply taking the PBH sample that yields the smallest signal width τ (dashed).We both see similar trends in the PBH fraction-PBH mass sensitivity.In the right panel, we show the fluctuation in the limit on the PBH fraction about the mean value, for different numbers of mock universes in the MC simulation.The variance in the PBH fraction decreases with increasing number of mock universes.

Figure 3 .
Figure 3.In the left panel we show the normalized probability distribution of PBH masses, in the case where the temperature of the Universe at percolation time is T * = 200 GeV, g * = 106.75,η χ = 10 −3 , v w = 0.2, β/H * = 60, and [∆V ef f (T * )] 1/4 = 100 GeV.The vertical red, solid line and green, dashed line correspond to the average PBH mass M and PBH mass at the peak of the PBH mass distribution, respectively.Please refer to the main text for more details.The right panel shows the CDF c associated with the PBH mass function, as a function of the universal dimensionless parameter x(M ).Note that this CDF only depends on f F V (t * ), and is unaffected by the choice of FOPT parameters; the FOPT parameters only dictate the form of x(M ).We also show the numerical result for c, c0 , and cres for f F V (t * ) = 0.29; the black vertical line marks the position where exp(−xI * ) = 10 −7 .

Figure 8 .
Figure 8. Similar to Fig. 7, we show the correlations of the fundamental parameters in the effective potential with the physical parameters that can be probed in PBH searches via Doppler shifts in pulsar timing signals and searches.
GeV/cm 3 is the current background cosmological CDM density.If we regard PBHs as a component of CDM, we can assume that the local PBH density relative to the local DM density, at present time, is simply given by