Can the poltergeist mechanism produce observable GWs?

Abstract

The enhancement of induced gravitational waves (GWs) occurs due to a sudden transition from an early matter-dominated era to the radiation-dominated era. We analyze the impact of the transition rate on the scalar potential. We find that the finite transition duration suppresses the oscillation amplitude of the scalar potential, consequently suppressing the amplitude of the GW energy spectrum. By numerically solving the background and perturbation equations, we demonstrate that the physically motivated models, such as the evaporation of primordial black holes and the decay of Q-balls, cannot produce an observable GW signal.

Introduction. Gravitational waves (GWs) serve as a pivotal tool for investigating the early Universe before big bang nucleosynthesis. There are many potential sources, such as inflation [1, 2, 3, 4, 5], first-order phase transitions [6, 7, 8, 9], preheating after inflation [10, 11, 12, 13], topological defects [14, 15, 16], etc. A notable scenario is GWs induced by the first-order scalar potential, which has been a hot topic for a long time [17, 18, 19, 20, 21, 22, 23, 24, 25, 26]. Generally, to have an observable induced GW signal, there are two scenarios: (1) Primordial curvature perturbations are significantly enhanced on small scales due to some mechanisms [27, 28, 29, 30, 31, 32, 33, 34, 35, 36]. (2) There is a long-lived early matter-dominated (eMD) era followed by a fast transition to the radiation-dominated (RD) era (called the poltergeist mechanism) [37, 38, 39, 40, 41, 42, 43, 44, 45, 46]. Such an eMD epoch is motivated in a wide variety of contexts, including the decaying dark sector [47, 48, 49, 50, 51, 52], primordial black holes (PBHs) [53, 37, 38, 54, 55, 56], post-inflationary reheating [57, 58, 59] and solitons such as Q-balls [60, 61]. Currently, there are no constraints on this epoch before nucleosynthesis [62].

During the eMD epoch, the gravitational potential does not decay, even on subhorizon scales. In this scenario, the amplification of induced GWs comes from the fast oscillation of the scalar potential after the transition, especially the mode entering the horizon deep in the eMD era. The earlier the mode enters the horizon, the faster the scalar potential oscillates in the RD era, and the stronger induced GWs. Since the density contrast of matter $\delta_{m}\propto a(\tau)$ in the eMD era, there is a cutoff scale $k_{cut}$ that keeps the linear perturbation theory valid, which is also the peak frequency of GWs.

In this letter, we revisit the GW amplification mechanism stemming from a rapid transition from the eMD to RD eras. We find that the evolution of the scalar potential is sensitive to the duration time of the transition. For those modes that reenter the horizon deep in the eMD era, even a short duration time can lead to fast decay of the scalar potential (before they oscillate in the RD era) [42, 43]. So, the behavior of the decay rate around the transition is very important. However, in the previous work, they either use the sudden transition approximation or use phenomenological parameterizations, which loses the evolution detail of the decay rate. We omit the parameterized background and analytical expression of the scalar potential. Instead, for the first time, we fully numerically solve the evolution of the background and perturbations in the context of PBHs and Q-balls, which is a belief that the decay rate is fast enough to create a sudden transition process. Our numerical method resolves transient dynamics, revealing that neither PBH evaporation nor Q-ball decay can produce transitions that meet the required rapidity criterion. The finite transition duration emerges as an irreducible physical parameter, thereby establishing that previous estimates of GW amplification in such scenarios are systematically overestimated.

Scalar potential. In order to show the duration time of transition is important to the evolution of the scalar potential, we begin with a parameterized equation of state

$$\omega(\tau)=\frac{1}{6}\left[\tanh{\frac{\tau-\tau_{eq}}{\Delta\tau}}+1\right],$$ (1)

which characterizes the transition from the eMD ($\omega=0$) era to the RD ($\omega=1/3$) era. Here $\tau_{eq}$ is the equal time of matter and radiation, and $\Delta\tau$ denotes the transition width. Combining the Friedmann equation and continuity equation, we can solve the background in terms of

$$\frac{a^{\prime\prime}}{a}+\frac{3\omega-1}{2}\left(\frac{a^{\prime}}{a}\right)^{2}=0.$$ (2)

The equation of motion for the scalar potential is

$$\phi^{\prime\prime}_{k}+3\mathcal{H}(1+\omega)\phi^{\prime}_{k}+[2\mathcal{H}^{\prime}+(1+3\omega)\mathcal{H}^{2}]\phi_{k}+\omega k^{2}\phi_{k}=0,$$ (3)

where $\mathcal{H}=a^{-1}da/d\tau$ is the conformal Hubble parameter. We assume $\omega=c_{s}^{2}$ for simplicity and ignore the entropy perturbation $\delta S$. We choose two sets of $\Delta\tau$ that denote different transition speeds to solve the background numerically, and we calculate the mode $k=100/\tau_{eq}$ as an example, which is an enhanced mode that reenters the horizon in the eMD era. The results are shown in Fig. 1.

Figure 1: Evolution of the equation of state with two sets of parameters denote the slow transition (dashed blue) and fast transition (dashed red). The solid lines represents the evolution of $\phi_{k}$ with $k=100/\tau_{eq}$.

We can see that when the transition rate is slower ($\Delta\tau=0.1$), the mode oscillates during the transition, and their amplitude is smaller compared to the quick transition case ($\Delta\tau=0.01$). We also show the impact of transition speed on the source term of GWs, which is defined by

$$\begin{split}&f(k_{1},k_{2},\tau)=2\phi_{k_{1}}(\tau)\phi_{k_{2}}(\tau)+\frac{4}{1+3\omega(\tau)}\\ &\times\left(\phi_{k_{1}}(\tau)+\frac{\phi^{\prime}_{k_{1}}(\tau)}{\mathcal{H}(\tau)}\right)\left(\phi_{k_{2}}(\tau)+\frac{\phi^{\prime}_{k_{2}}(\tau)}{\mathcal{H}(\tau)}\right).\end{split}$$ (4)

Figure 2: Evolution of the source term in the slow transition (blue line) and fast transition (red line).

In Fig. 2, we choose a specific configuration $k_{1}=k_{2}=100/\tau_{eq}$ and we can see that the reduction is more significant. The result infers that the duration time of the transition is very important to the prediction of the amplitude of GWs, so we need to specifically investigate the transition rate around the transition time.

Decay rates. PBHs and Q-balls are well-motivated candidates to trigger eMD epochs in the Universe. They behave as decaying nonrelativistic fluids. In the early stage, they are stable and their energy density may dominate the Universe if their stable time is long enough. However, when the decay rate becomes stronger, they will decay to radiation and heat the Universe. As mentioned above, the decay rates of these objects around the transition time critically influence the transition speed from matter to radiation domination, which in turn determines the resonant enhancement of GWs induced by primordial curvature perturbations. In general, the decay rate can be expressed as

$$\Gamma\equiv-\frac{1}{M}\frac{dM}{dt}=\frac{n}{t_{eva}-t},$$ (5)

where $M$ is the mass of the matter component, $t_{eva}$ is the evaporation time, when most of their mass become radiation. The decay rate coefficients $n=1/3,3,1,3/5$ for PBHs, thin-wall, thick-wall, and delayed Q-balls [43, 63, 64, 65], respectively. PBHs decay through Hawking radiation, with their mass evolving as $M_{PBH}\propto(t_{eva}-t)^{\frac{1}{3}}$. Q-balls are stable scalar field configurations stabilized by a conserved charge. Their decay mechanisms depend on their type. By comparison of the decay rate coefficients, PBHs have the smallest $n$ thus the fastest transition rate, which is supposed to have a stronger GW signal. So we take PBHs as an example to show how to fully numerically solve the background and perturbations. The discussion is easily extended to the case of Q-Balls.

Evolution of background and perturbations. In order to solve the background in conformal time, we rewrite the decay rate [5] as differential equations. Background evolution can be expressed as

$$\mathcal{H}^{2} = a^{2}(\rho_{r}+\rho_{m})/3,$$ (6)

$$\rho_{m}^{\prime} = -(3\mathcal{H}+\tilde{\Gamma})\rho_{m},$$ (7)

$$\rho_{r}^{\prime} = -(4\mathcal{H}-\tilde{\Gamma})\rho_{r},$$ (8)

$$\tilde{\Gamma}^{\prime} = \mathcal{H}\tilde{\Gamma}+n\tilde{\Gamma}^{2},$$ (9)

where $\tilde{\Gamma}\equiv a\Gamma$ is the conformal decay rate,, $\rho_{m}$ and $\rho_{r}$ are energy density of matter and background radiation. Here $n=1/3$ since we focus on the PBH-dominated case.

We start our numerical calculation at the formation time of PBHs $\tau_{i}$ in the eRD era, where the abundance of PBHs is characterized by $\beta_{i}=\rho_{m}(\tau_{i})/\rho_{r}(\tau_{i})$. We assume $\beta_{i}\ll 1$ throughout this letter. We can ignore the evaporation of PBHs at initial time and the scale factor can be expressed as

$$\frac{a(\tau)}{a_{i}}=\frac{1}{\beta_{i}}\left[(2+\beta_{i}-2\sqrt{1+\beta_{i}})\left(\frac{\tau}{\tau_{i}}\right)^{2}+2(-1+\sqrt{1+\beta_{i}})\frac{\tau}{\tau_{i}}\right].$$ (10)

We rescale the conformal time and scale factor at $\tau_{i}$, where $\tau/\tau_{i}\rightarrow\tilde{\tau},a(\tau)/a_{i}\rightarrow\tilde{a}(\tilde{\tau})$. The energy density can be dimensionless as $\rho_{m}\tau_{i}^{2}/M_{pl}^{2}\rightarrow\tilde{\rho}_{m},\rho_{r}\tau_{i}^{2}/M_{pl}^{2}\rightarrow\tilde{\rho}_{r}$. From now on, we omit the symbol “$\tilde{\quad}$”, just keep in mind that all variables are dimensionless and rescaled at $\tau_{i}$. Then the initial conditions can be expressed as

$$\begin{split}&a(\tau_{i})=1,\quad a^{\prime}(\tau_{i})=\frac{2(1+\beta_{i}-\sqrt{1+\beta_{i}})}{\beta_{i}},\\ &\rho_{r}(\tau_{i})=\frac{27\beta_{i}^{2}}{4(-1+\beta_{i}+\sqrt{1+\beta_{i}})^{2}},\quad\rho_{m}(\tau_{i})=\beta_{i}\rho_{r}(\tau_{i}).\end{split}$$ (11)

The initial condition for $\tilde{\Gamma}$ is more tricky. At formation time, $t_{i}\ll t_{eva}$ and $\Gamma(t_{i})\approx 1/(3t_{eva})$, which is determined by the initial mass of PBHs. However, $t_{eva}$ cannot be determined since we need to integrate the scale factor in conformal time, which is unclear before solving the background equations. In practice, we can pre-set the numerical quantity of $t_{eva}$, then solve the background equations until the end of the transition $\tau_{eva}$, where $\rho_{m}(\tau_{eva})\ll\rho_{r}(\tau_{eva})$ and $\omega(\tau_{eva})\approx 1/3$. Then the actual value of $t_{eva}$ can be determined by integrating the conformal time from $0$ to $\tau_{eva}$. Until now, we know exactly the initial mass of PBHs. In general, only two parameters $\beta_{i}$ and $M_{PBH,i}$ can determine whether there will be the eMD era or not, as well as the depth of the eMD era, if any.

We then calculate perturbations in the conformal Newtonian gauge [66, 38]

$$\delta_{m}^{\prime} = -\theta_{m}+3\phi^{\prime}-\tilde{\Gamma}\phi,$$ (12)

$$\theta_{m}^{\prime} = -\mathcal{H}\theta_{m}+k^{2}\phi,$$ (13)

$$\delta_{r}^{\prime} = -\frac{3}{4}(\theta_{r}-3\phi^{\prime})+\tilde{\Gamma}\frac{\rho_{m}}{\rho_{r}}(\delta_{m}-\delta_{r}+\phi),$$ (14)

$$\theta_{r}^{\prime} = \frac{k^{2}}{4}\delta_{r}+k^{2}\phi-\tilde{\Gamma}\frac{3\rho_{m}}{4\rho_{r}}(\frac{4}{3}\theta_{r}-\theta_{m}),$$ (15)

$$\phi^{\prime} = -\frac{k^{2}\phi}{3\mathcal{H}}-\mathcal{H}\phi-\frac{\mathcal{H}}{2}(\frac{\rho_{m}}{\rho_{tot}}\delta_{m}+\frac{\rho_{r}}{\rho_{tot}}\delta_{r}),$$ (16)

where we neglect anisotropic stress, so $\phi$ is the only scalar perturbation of metric. Here $\rho_{tot}=\rho_{m}+\rho_{r}$ is total energy density. Solve the above equations in superhorion limit ($\tilde{\Gamma}=0$), we get the initial condition in the eRD era

$$\delta_{m,i}=-\frac{3}{2}\phi_{i},\quad\delta_{r,i}=-2\phi_{i},\quad\theta_{m,i}=\theta_{r,i}=\frac{k^{2}\tau}{2}\phi_{i},$$ (17)

where $\phi_{i}$ is the primordial perturbation generated from inflation. It is worth mentioning that our choice is the adiabatic initial conditions.

Numerical method to calculate induced GWs. Given that both the background and perturbations are computed numerically, induced GWs must also be derived through numerical integration. The energy density of induced GWs is given by

$$\Omega_{GW}(k,\tau_{0})h^{2}=\Omega_{r,0}h^{2}\frac{1}{24}(\frac{k}{\mathcal{H}_{c}})^{2}\overline{\mathcal{P}_{h}(k,\tau_{c})},$$ (18)

where $\Omega_{r,0}h^{2}\approx 4.2\times 10^{-5}$ is the current radiation energy parameter, $\mathcal{P}_{h}$ is the power spectrum of GWs. The subscript “c” represents the time when the mode of interest becomes well inside the horizon and the energy density of GWs becomes almost constant. $\overline{\mathcal{P}}_{h}$ is the time average per period $T$ around $\tau$, explicitly in $\bar{x}(\tau)=\frac{1}{T}\int_{\tau-T}^{\tau}x(\tilde{\tau})d\tilde{\tau}$. The power spectrum can be calculated by

$$\begin{split}\mathcal{P}_{h}(k,\tau)&=\frac{4}{81a^{2}(\tau)}\int_{\begin{subarray}{c}|k_{1}-k_{2}|\leq k\leq k_{1}+k_{2}\end{subarray}}d\ln{k_{1}}\,d\ln{k_{2}}I^{2}(k,k_{1},k_{2},\tau)\\ &\times\frac{(k_{1}^{2}-(k^{2}-k_{2}^{2}+k_{1}^{2})^{2}/(4k^{2}))^{2}}{k_{1}k_{2}k^{2}}\mathcal{P}_{\mathcal{R}}(k_{1})\mathcal{P}_{\mathcal{R}}(k_{2}),\end{split}$$ (19)

where $\mathcal{P}_{\mathcal{R}}(k)=\Theta(k_{max}-k)A_{s}(k/k_{*})^{n_{s}-1}$ is the power spectrum of primordial curvature perturbations with $A_{s}\approx 2.1\times 10^{-9}$ being the amplitude at the pivot scale, $n_{s}\approx 0.97$ the spectral tilt and $k_{*}=0.05\mathrm{Mpc^{-1}}$ the pivot scale [67]. $I(k,k_{1},k_{2},\tau)$ is the kernel function defined by

$$\begin{split}&I(k,k_{1},k_{2},\tau)=4k^{2}\int_{0}^{\tau}d\tilde{\tau}a(\tilde{\tau})G_{k}(\tau,\tilde{\tau})\Bigg{[}2T_{k_{1}}(\tilde{\tau})T_{k_{2}}(\tilde{\tau})\Bigg{.}\\ &\Bigg{.}+\frac{4}{3(1+\omega(\tilde{\tau}))}\left(T_{k_{1}}(\tilde{\tau})+\frac{T_{k_{1}}^{\prime}(\tilde{\tau})}{\mathcal{H}(\tilde{\tau})}\right)\left(T_{k_{2}}(\tilde{\tau})+\frac{T_{k_{2}}^{\prime}(\tilde{\tau})}{\mathcal{H}(\tilde{\tau})}\right)\Bigg{]},\end{split}$$ (20)

where $G_{k}(\tau,\tilde{\tau})$ is the Green function for tensor perturbation. $T_{k}(\tau)$ is the transfer function of the scalar potential, where $\phi_{k}(\tau)=T_{k}(\tau)\phi_{i}$. $\omega(\tau)=\rho_{r}(\tau)/3(\rho_{m}(\tau)+\rho_{r}(\tau))$ is the equation of the state parameter.

We can combine two independent homogeneous solution to obtain the Green function

$$G_{k}(\tau,\tilde{\tau})=\frac{v_{1k}(\tau)v_{2k}(\tilde{\tau})-v_{1k}(\tilde{\tau})v_{2k}(\tau)}{v_{1k}^{\prime}(\tilde{\tau})v_{2k}(\tilde{\tau})-v_{1k}(\tilde{\tau})v_{2k}^{\prime}(\tilde{\tau})}\Theta(\tau-\tilde{\tau}),$$ (21)

where the two homogeneous solutions $v_{ik}=ah_{k}$ can be obtained by equation of motion of tensor perturbation

$$\left(\partial_{\tau}^{2}+k^{2}-\frac{1-3\omega(\tau)}{2}\mathcal{H}^{2}\right)v_{ik}=0,$$ (22)

The average square of the kernel can be obtained as [68, 69]

$$\begin{split}\overline{I^{2}(k,k_{1},k_{2},\tau)}\approx I_{2}^{2}(k,k_{1},k_{2},\tau)\overline{v_{1k}^{2}(\tau)}\\ +I_{1}^{2}(k,k_{1},k_{2},\tau)\overline{v_{2k}^{2}(\tau)}\\ -2I_{1}(k,k_{1},k_{2},\tau)I_{2}(k,k_{1},k_{2},\tau)\overline{v_{1k}(\tau)v_{2k}(\tau)},\end{split}$$ (23)

where $I_{n}(k,k_{1},k_{2},\tau)$ is obtained by splitting the Green function into $g_{nk}$ contributions.

We conclude the steps to calculate GWs as follows.

• Use Eqs. (6)-(9) to solve the background equations from $\tau_{i}$ to $\tau_{eva}$ with the initial condition Eq. (11). Continuous background derived analytically during the late RD epoch, valid for $\tau>\tau_{eva}$.

• Find the cut-off scale $k_{cut}$ by solving the perturbation evolution Eqs. (12)-(16), where the cut-off scale is defined as $\delta_{m,k_{cut}}\sim 1$ at the end of transition.

• Choose the range of mode $[k_{min},k_{max}]$ in which we are interested. $k_{min}\tau_{rh}=1/100$ is the mode that enters the horizon in the late RD era. $k_{max}=k_{cut}$ is the cutoff scale.

• For a given mode $k\in[k_{min},k_{max}]$, we obtain the numerical solution of $v_{1k},v_{2k}$ using the analytical solution as the initial condition when they are well outside the horizon during the eRD epoch.

• For each wave number $k$ computed in the previous step, we discretize the $k_{*}\in[10^{-2},10^{2}]k$ interval into 500 logarithmically spaced bins, numerically solve the perturbation equations Eqs. (12)-(16) for each mode, and systematically archive the resultant data.

• Combining $k_{1}$ and $k_{2}$ in the integration domain defined in Eq. (19), where we evaluate the associated kernel functions specified in Eq. (20). Numerical integration is implemented through a trapezoidal quadrature scheme with a bidimensional grid of approximately $500\times 500$ sampling points.

• Calculate the time average of the GW power spectrum.

Following the above steps, we can fully numerically get the power spectrum of GWs.

Numerical results. We choose a typical set of parameters: $\beta_{i}=10^{-6},M_{PBH,i}=10^{3}\mathrm{g}$, in which PBHs can dominate the Universe and evaporate at $9.7\times 10^{5}\mathrm{GeV}$. The background evolution is shown in Fig. 3 and Fig. 4. We label several scales for which we are interested: (1) Matter and radiation equal time $\tau_{eq1}$ (red dashed) and $\tau_{eq2}$ (green dashed). (2) The reentry time of the cutoff scale $\tau_{cut}$, where the modes $k\leq k_{cut}=14.88\mathrm{Hz}$ are safe in linear theory. (3) The time $\tau_{eva}$ when the Universe goes back to the RD epoch.

Figure 3: Evolution of the equation of state and energy density of matter and radiation. The vertical lines represents former matter radiation equal time $\tau_{eq1}$ (red dashed), the time that cutoff scale reenter the horizon $\tau_{cut}$ which defined as $k_{cut}=\mathcal{H}(\tau_{cut})$ (blue), the latter matter radiation equal time $\tau_{eq2}$ (green dashed) and transition finish time $\tau_{eva}$, respectively.

We can see that $\tau_{eq2}$ is very close to $\tau_{rh}$, which infers that the comoving decay rate $\tilde{\Gamma}$ increases rapidly around $\tau_{eq2}$. We can clearly see this rapid increase in Fig. 4. Although $\tau_{eq2}$ and $\tau_{eva}$ seem very closed, we cannot state that it is the justification for the assumption of a sudden transition. As we mentioned above, the modes that reenter the horizon deep inside the eMD era are also sensitive to the duration time of the transition. So we should dive into the evolution details of $\tilde{\Gamma}$ around $\tau_{eq2}$. When $\tilde{\Gamma}\ll\mathcal{H}$, $\tilde{\Gamma}^{\prime}\propto\mathcal{H}\tilde{\Gamma}$ infers that $\tilde{\Gamma}\propto\tau$ in the early stage. On the other hand, when $\tilde{\Gamma}\gg\mathcal{H}$, $\tilde{\Gamma}\propto(1/3)\tilde{\Gamma}^{2}$, until then $\tilde{\Gamma}$ start to grow rapidly. We can see that, at time $\tau_{eq2}$, when $\omega(\tau_{eq2})=1/6$ and $\tilde{\Gamma}(\tau_{eq2})=\mathcal{H}(\tau_{eq2})/2$, $\tilde{\Gamma}$ is still approximately proportional to $\tau$. This indicates that most of the PBH mass had already evaporated by $\tau_{eq2}$.

We also show the evolution of a mode($k=90/\tau_{eq2}$) that reenters the horizon deep inside the eMD era in Fig. 5 and compared to the analytical solution which uses a sudden transition assumption. The numerical results reveal that during the transition phase, the scalar potential experiences rapid decay and evolves proportionally to $\tau$ in the late RD era, showing a substantially reduced amplitude relative to the analytical prediction. It should be noted that the decrease is even greater in the power spectrum since $\mathcal{P}_{h}(k)\propto\langle\phi^{4}(k,\tau)\rangle$. So, the power spectrum of GWs is much smaller than expected, as shown in Fig. 6. The enhancement of the power spectrum is only two orders of magnitude increase. In this case, the peak amplitude $\Omega_{GW}(k,\tau_{0})h^{2}\sim 10^{-21}$ of GWs cannot be detectable in the current observatory.

Figure 4: Evolution of comoving decay rate and comoving Hubble parameter. The zoomed figure shows the details around $\tau_{eq2}$.

Figure 5: Evolution of the scalar potential $\phi_{k}(\tau)$. Blue line is the numerical solution of mode $k=90/\tau_{eq2}$. Orange line is the analytical solution of the same mode which use sudden transition assumption. Green dashed line labels $\tau_{eq2}$.

Figure 6: Energy spectrum of GWs at the time $\tau_{c}$ where the mode of interest becomes well inside the horizon and the energy density of GWs becomes almost constant. The choice of parameters are $\beta_{i}=10^{-6},M_{PBH,i}=10^{3}\mathrm{g}$.

Conclusions and discussions. We have revisited the production of induced GWs during the transition from the eMD era to the RD era. By numerically solving the coupled equations governing background evolution and scalar perturbations, we show that the finite duration of the transition phase plays a decisive role in suppressing GW amplitudes. Unlike the sudden transition approximation, which predicts the resonant enhancement of GWs, realistic PBH evaporation and Q-ball decay leads to rapid decay of the scalar potential during reheating, drastically reducing the source term for GW generation. Our analysis identifies the decay rate of PBHs as a key parameter controlling the transition speed, with slower transitions further diminishing the GW signal. The results underscore the necessity of incorporating realistic transition timescales in modeling early universal GWs, as idealized assumptions systematically overpredict observable signals.

Our result proves that the poltergeist mechanism originating from the PBH evaporation cannot significantly enhance the GW to the observable level in the linear region. In addition, it has been proven that the mode in the non-linear region can also significantly enhance GWs [44]. It infers that if the idealistic poltergeist mechanism (a sudden transition) occurs in the Universe, the nonlinear mode will be even stronger and can easily break the big bang nucleosynthesis constraint of GWs. Our results ease the tension.

In this letter, we have performed numerical calculation in the context of PBHs, especially. This is because the transition rate at late time is much higher compared to other physically motivated candidates such as Q-balls. If PBHs cannot establish an idealistic poltergeist mechanism, neither can Q-balls.

We mainly focus on the adiabatic initial condition throughout this letter. However, PBHs themselves may serve as isocurvature perturbations. In this case, the existence of PBHs can produce relative entropy perturbations and source the curvature perturbation, leading to the enhancement of the power spectrum. Our result just shows that it is hard to enhance induced GWs to an observable level without enhancing the power spectrum of the curvature perturbation.

Acknowledgements. This work is supported in part by the National Key Research and Development Program of China Grant No. 2020YFC2201501, in part by the National Natural Science Foundation of China under Grant No. 12475067 and No. 12235019.