Cosmological standard timers from unstable primordial relics

In this article we study a hypothetical possibility of tracking the evolution of our Universe by introducing a series of the so-called standard timers. Any unstable primordial relics generated in the very early Universe may serve as the standard timers, as they can evolve through the whole cosmological background until their end while their certain time-varying properties could be a possible timer by recording the amount of physical time elapsed since the very early moments. Accordingly, if one could observe these quantities at different redshifts, then a redshift-time relation of the cosmic history can be attained. To illustrate such a hypothetical possibility, we consider the primordial black hole bubbles as a concrete example and analyze the mass function inside a redshifted bubble by investigating the inverse problem of Hawking radiation. To complete the analyses theoretically, the mass distribution can serve as a calibration of the standard timers.


I. INTRODUCTION
The cosmology of Λ-cold dark matter (ΛCDM) is widely acknowledged as the most successful theory in depicting the evolution of the Universe, which has been examined in various high-precision observations, namely, the cosmic microwave background (CMB) [1][2][3], large scale structure (LSS) [4], and the acceleration at late times [5,6], etc.However, there remain puzzles in the Universe, such as the recent focuses on the Hubble tension [7][8][9] and σ 8 tension [10], which hint on the existence of unknowns beyond the ΛCDM.
Accompanied with the dramatic developments of observational technologies, various approaches of astronomical measurements have been proposed to study the evolution of the Universe.Type Ia supernovae (SNe) produces consistent peak luminosity after calibration, of which the distances can measure by comparing their absolute and apparent magnitudes, and hence they can serve as a standard candle [11].Baryon acoustic oscillations (BAO) determine the fixed scale that the sound wave can travel before the recombination, through probing the BAO scale at different redshifts, and then the cosmological models can be effectively constrained, and accordingly, BAO can work as a standard ruler [12,13].Gravitational waves (GWs) and electromagnetic (EM) waves from binary compact objects and EM counterpart provide the calibration between luminosity distance and redshift, which makes GWs standard sirens in constraining cosmological parameters [14].While the aforementioned methods can constrain cosmological models, the distinct physics behind them would bring difficulties in excluding the unrecognized systematic uncertainties and may still cause the tension for cosmological parameters.Moreover, most observational signals lie within certain redshift windows.Standard candle via SNe covers redshift window 0.01 < z < 2.3 [15].Standard ruler via BAO by observing galaxies can be applied mostly on 0.106 ≲ z ≲ 0.61 [12,13], the signals from Lyα forest of high-redshift quasars [16] and CMB [10] provide the acoustic scales at z = 2.34 and 1100, respectively.Standard siren observations based on aLIGO are only allowed within z < 1 [17,18].
In contrast, unstable primordial relics can be used to track the whole evolution history of the Universe.If these relics are unstable, the decay of the relics tracks the physical time of the Universe, just as standard timers distributed in the Universe.If we observe distant standard timers, a redshift-time relation is obtained, which can be used to constrain cosmological models.
In this paper, we discuss an explicit model of the standard timer, i.e., primordial stellar bubbles consisting of primordial black holes (PBHs) in tracking the evolution of the Universe.The PBH stellar bubbles can be generated from some primordial physics, such as multi-stream inflation [19,20], which causes a number of PBHs clustering in specific regions.Its primordial origin produces the same PBH mass functions inside PBH stellar bubbles.Since the observed signals from PBH bubbles are determined by the internal mass function evolution and external cosmic expansion, the analysis of the observed signals can construct the one-to-one correspondence between the PBH bubble internal mass functions and the evolution of the Universe.In the pioneer paper [21], we studied the detectability of such PBH bubbles with the current and the near future observations, and presented the detailed analyses of the possible observational signals from such a stellar bubble, including the EM and GW signals, coming from the light PBHs' Hawking radiations and binary mergers of PBHs inside a bubble.Limited by the current resolution of gamma-ray detectors, these PBH bubbles can be regarded as the exotic celestial objects.After having used the point-source differential sensitivity in the 10-year observation of Fermi LAT for a high Galactic latitude (around the north Celestial pole) source [22], we showed that the lowest detectable bubble mass is M bub ∼ 10 32 g with the peak mass 10 15 g of the lognormal distribution of PBH masses inside the bubble.Also, we found that massive PBH binaries with M pk ∼ 10 34 − 10 38 g and M bub ∼ 10 45 − 10 48 g can produce detectable GWs within in the frequency band of LISA and BBO.Hence, EM and GW signals are complementary for light and heavy PBH bubbles.Since PBH bubbles can be regarded as point-sources for the current EM observations (like the Fermi LAT experiment), the redshift information of PBH bubbles is thus imprinted in their EM signals.As we will show below, one can infer the redshift information of PBH bubbles through the inverse process, which is essential to the determination of the time evolution of our Universe.In the following context, we focus on the EM signals coming from PBH bubbles and the role that PBH bubbles can play in , hence we call it the standard timer.
This paper is organized as follows.We present a brief introduction of PBH bubbles in Sec.II; In Sec.III, we study the inverse problem for Hawking radiation to extract the mass function of a PBH bubble from EM observations.Then, we describe the approach of calibrating the redshifts of PBH bubbles from the inverse results and construct the redshift-time calibration through the standard timer.In Sec.IV, we apply the redshift-time calibration from the standard timer to test the standard cosmological model.The conclusion and discussions are given in Sec.V.
In usual cases, PBHs are regarded as individual isolated objects in the Universe, however, there are many scenarios for clustering of PBHs at some certain scales.If the sizes of these regions are small enough, say, smaller than the resolutions of current telescopes, they behave as exotic celestial objects, i.e., the PBH bubbles [20,21].Generally speaking, these stellar bubbles can be generated from some new-physics phenomena that might have occurred in the primordial Universe, such as, quantum tunnelings during or after inflation [56,57], multi-stream inflation [19,58,59], inhomogeneous baryogenesis [60], etc.In these cases, before the bubble-wall tension vanishes, the field values are different between inside and outside of the bubble.Such difference can result in different local physics inside the bubble (for PBH cases, see [20,61] for details), namely the production rate of the exotic species of matter.
In general, the abundance of bubbles is determined by the probability of the tunneling (for phase transition) or bifurcation (for multi-stream inflation).And the size of the bubbles is determined by the comoving scale at which tunneling or bifurcation happened.For example, in the multi-stream inflation model, the radius of the bubble is similarly R b = R 0 exp(−N b ), where R 0 is the radius of the current observable Universe and N b is interpreted as the e-folding number between the beginning of the observable inflation to the bifurcation.Since the bifurcated path eventually merges, the tension of the bubble wall vanishes automatically.The number density of the bubble n b is determined by the shape of the multi-field potential, and the amplitude of the isocurvature fluctuation during inflation.We should keep in mind that the local abundance of PBHs inside a bubble f local ≡ Ω PBH /Ω tot is naturally enhanced by the probability of formation of this bubble, such as the bifurcation probability f bifu < 1, see the Fig. 1 in Ref. [20].We can write f local = f glob f −1 bifu , where f glob < 1 is the normal definition of the abundance of PBHs over the whole observable Universe [26].So, it is expected that the enhanced observable signals coming from these PBH bubbles, like the detectable GW and EM signals with the current observations [20,21].Notably, [21] shows that even for PBH bubbles at quite high reshifts z ∼ 10 3 , their EM signals are possible to be detected by FermiLAT experiment.So that in this sense, PBH bubbles are able to track the history of the Universe over a long period.

III. PBH BUBBLES AS A STANDARD TIMER
The same primordial origins of PBH bubbles are expected to be formed at a nearly single moment with the same initial mass functions [20].Then, PBHs inside a bubble would evaporate through Hawking radiations, which results in the deformation of their mass function insides the bubble.We notice that there are two key quantities encoded in this process.One is the cosmic expansion which causes the redshift of the observed photon energies, the other is the physical evolution time that hides in the deformed mass function.The cosmological redshift depends on the evolution of the Universe, while the physical evolution time is determined by Hawking evaporation and does not rely on the cosmological redshift.Hence, the independent channels of cosmological redshifts and physical times indicate that PBH bubbles can calibrate the redshift-time relation, which works as a standard timer of the Universe.
In order to extract cosmological redshifts and physical times from the EM signals of PBH bubbles, two problems need to be resolved.One is that the redshift cannot be obtained directly from the observed photon energy due to the unknown emission energy of photon.However, the observed photon flux spectrum is a result from the interplay of the Hawking radiations emitted over the PBH mass function and the cosmological redshifts, which enable us to obtain the redshifted PBH mass function from the observed photon flux spectrum and thus the redshift is encoded in it.This method is discussed in this section and the details presented in Appendix A and B. There are some other methods for helping estimate the redshift of PBH bubbles such as the redshifts from the host galaxies or neighbour astrophysical objects [62].The other problem is to determine the physical evolution times of PBHs in lack of their primordial mass function.This issue can be addressed when the redshift is extracted from the observed photon flux spectrum, the physical PBH mass function can be extracted from the redshifted mass function, then physical evolution time between two PBHs bubbles can be obtained by comparing their physical PBH mass functions.In this section, we will introduce the method of extracting the redshifted mass function of PBHs from observed photon flux spectrum, then obtain the redshift from the redshifted mass function, and construct the standard timer array in the Universe.

A. Inverse Problem for Hawking Radiation
The Hawking radiation spectrum R(E) emitted from a PBH bubble follows where H(E, M ) is Hawking radiation kernel and n(M ) is the PBH mass function.In order to study the evolution of PBH mass function n(M ), extracting the PBH mass function from the observed emission rate spectrum is a key step, and we call it the inverse problem for Hawking radiation.The similar problem has already been studied for decades, i.e., determining the temperature distribution on the thermal radiator for a given radiation power spectrum, which is known as the inverse black-body radiation problem [63].There are several inversion methods that have been proposed to solve the inverse black-body radiation problem, such as iteration method and Möbius inversion transformation [63,64] (see Appendix A for more details).
In studying the inverse Hawking radiation problem, the properties of the Hawking radiation kernel H(E, M ) and emission rate spectrum R(E) are crucial in analysis.The Hawking radiation kernel can be divided into two parts, the primary emission H p (E, M ) (the direct Hawking emission) and the secondary emission H s (E, M ) (from the decay of gauge bosons or heavy leptons and the hadrons produced by the fragmentation of primary quarks and gluons [65]), and we write Since the gamma ray observations are more sensitive in the high-energy range and the primary emission dominates such Hawking radiation [66], we focus on the primary emission in the inverse Hawking radiation problem.On the other hand, the secondary emission involves the fragmentation and hadronization of quark and gluon jets, which lacks the precise analytic expressions, especially in the low-energy range [67].Moreover, one can in principle involve the secondary contribution by virtue of the numeric calculations, we leave this to the future study.In practice, the emission rate spectrum R(E) comes from the high-energy gamma ray observation, which means R(E) is in the form of a data array instead of an analytic expression.As a result, we should apply discretization on the Eq. ( 1), which is in form of Here, w j is the weight of the quadrature formula, N is the number of discretization nodes for the mass function n(M ).Then, Eq. ( 1) is transformed to a linear algebraic equation, we can apply the method for the least squares problem (see Ref. [68] for more details) to resolve the mass function vector as follows, Here, R and n are the emission rate vector and mass function vector, respectively.K is the kernel matrix, whose element is Actually, due to the observational error existed in the emission rate vector R, there exist a number of potential vectors n for the Eq. ( 4), and the regularization methods [69] can be applied to determine the physical mass function vector.We also show the formal method for the inverse Hawking radiation problem in Appendix B.
In the observed photon flux from a PBH bubble, the photon energy is redshifted by the cosmic expansion.The photon flux can be expressed as Here, L(E; z) is the intrinsic luminosity emitted from the PBH bubble at redshift z which can be calculated as , where V (z) is the volume of the PBH bubble.Then, Eq. ( 5) can be written as There is 1 + z term appearing in the Hawking kernel, and the unknown redshift causes the uncertainty in determining the PBH mass function n(M ; z) in the inverse problem.In general, we can randomly choose the redshift in the Hawking kernel until the proper mass function returns, however, the error in introducing unknown redshift may increase the error in mass function.Consequently, we can transform 1 + z term from the Hawking kernel to the mass function, which results in the redshift uncertainty in the mass function.In the primary emission of the Hawking radiation H p (E, M ), the E and M terms in H p (E, M ) are symmetry, see Eq. (B2), so we have Then, Eq. ( 6) can be written as Here, M ′ ≡ M (1 + z).Then we have the equation, Following the method for the inverse Hawking radiation problem in Eq. ( 4), we can obtain the redshifted mass function

B. Redshift Calibration
In the redshifted mass function f (M ; z), the redshift z is unknown, which is essential for the standard timer.Considering two PBH bubbles with redshifts z 1 and z 2 , we get the redshifted mass functions Here, a factor 1/(1 + z) appears in the argument of mass function n(M ) which transforms the argument M → M/(1 + z) in n(M ).By comparing the locations of the large-mass tails of these two normalized redshifted mass functions, where the Hawking radiation effect is negligible and the large-mass tails remain nearly unchanged (see the blue and red solid curves in Fig. 1), we obtain the redshift ratio η: For an illustration, we take a normalized lognormal mass function [see Eq. ( 39)] at redshifts z = 1 and z = 2 shown in Fig. 1, which gives the η ≈ 1.5005 with the error 0.03%.For a given redshift z 1 , z 2 = η(1 + z 1 ) − 1, we can reconstruct the normalized PBH mass function at z 1 and z 2 as follows, Here, ñ and f denote the normalized comoving mass functions.Then, the physical evaporation times can be extracted from ñ(M ; z 1 ) and ñ(M ; z 2 ).Suppose that the physical evaporation time between z 2 and z 1 is t m , we write the mass evaporation relation [see Eq. ( 24)] as where δ(t m ) is a term corresponding to the evaporated mass during the period z 2 to z 1 .Then, we construct the relation between ñ(M ; z 1 ) and ñ(M ; z 2 ) as Applying Eq. ( 14) into Eq.( 15), we obtain In the low-mass approximation M z1 ≪ δ(t m ) which holds for the small values of M z1 , since the most of PBH's masses have been evaporated for small M z1 .Thus, we obtain, Here, the term ñ(δ(t m ); z 2 )/δ 2 (t m ) can be extracted from ñ(M z1 ; z 1 ) in low-mass approximation, which gives the δ(t m ) from ñ(M z2 ; z 2 ).The physical evaporation time t m is thus extracted from δ(t m ).Impressively, the low-mass tail on the logarithmic scale possesses the universal slope d log n(M )/d log M ≃ 2, which is independent on the initial PBHs' mass function, arising from the properties of mass evaporation [70], which is also confirmed by our following analysis, see Fig. 5.
In order to determine the redshift z 1 , we need to compare the physical evaporation time t m with the cosmic evolution time t z between z 1 and z 2 .The proper z 1 should be chosen such that the physical evaporation time is same as the cosmic evolution time as follows: where t z is assumed to follow the standard cosmology evolution, which is calculated as Then the redshift of one PBH bubble is calibrated.
The accuracy of the calibrated redshift z depends on the redshift ratio η.A practical method for redshift calibration should be stable, i.e., the error in calibrated redshift can not be much larger than the error in the redshift ratio.We take the PBH bubble with the same parameters setting in Fig. 1 as an example, see Fig. 2, which gives the calibrated redshifts z 1 ≃ 1.002 and z 2 ≃ 2.004.

C. Standard Timer via Redshift-Time Calibration
As we have discussed above, the redshift of the PBH bubble is calibrated, which can work as a redshift calibrator for other PBH bubbles.The redshift ratio η between a PBH bubble and the redshift calibrator can be determined by matching the redshifted mass functions.As a result, the redshift of a PBH bubble z PBH is given by where z 0 is the redshift of the redshift calibrator.Then physical time interval ∆t between z 0 and z PBH should equal to the physical evolution time between the redshift calibrator and the PBH bubble in Eq. ( 17), The accuracy of the calibrated physical time interval ∆t depends on the redshift ratio η.In order to calibrate the ∆t with high precision, we need to estimate the error in the calibrated time versus the error in redshift ratio, see Fig. 3, which produces the error in time calibration is around 0.1%, while the error in η is obtained around 0.03% and the redshift calibrator is set to z = 1.Consequently, the redshift-time calibration (z 0 , z PBH , ∆t) is obtained.After detecting a bunch of PBH bubbles in the Universe, a redshift-time calibration array can be constructed cover from the local Universe to the primordial Universe, which helps us constrain the cosmological models as the standard timer in the Universe.
-0.04 -0.02 0.00 0.02 0.04 The error in the time calibration in terms of the error in η measurement.We choose the lognormal mass function with M pk = 10 15 g and σ = 1 inside a PBH bubble and the redshift calibrator z = 1 as an example in calculation.

D. The Time-Dependent Mass Function
In the above discussions, we have presented the formal framework of a standard timer made of PBH bubbles.Here, we derive the approximated analytic formulas for the time evolution of PBHs' mass function due to their Hawking radiation, which enables us to extract the physical evolution time of PBH bubbles by using the methods studied in the above part.
We adopt the standard Hawking evaporation picture [65]: a black hole would directly radiate fundamental standard model particles whose de Broglie wavelengths are of the order of black hole size.Hence, PBHs would radiate heavier particles successively with time, as shown in the Table I.According to the Hawking radiation and energy conservation, we derive the mass-loss rate for a single PBH as follows, where ϕ(M ) measures the number of emitted particle degrees of freedom for a PBH with mass M and the factor "3" is included for convenience.The relativistic contributions to ϕ(M ) per degree of particle freedom are [65] f s=0 = 0.267, f s=1 = 0.060, f s=3/2 = 0.020, f s=2 = 0.007, where we introduced the dimensionless parameter f ≡ ϕ/(5.34 × 10 25 g 3 s −1 ).Hence, the PBH's mass at time t can be written as the following formula: where M f is the PBH mass at the formation epoch t f .According to this formula, we can also derive the lifetime of a single PBH with an arbitrary initial mass.
In order to make more precise analysis for the light PBHs which dominates EM signals from PBH bubbles, we extended the approximation for ϕ[M (t)] used in Ref. [70]: the small mass region is further separated, i.e., M ≤ M g , M g is the mass scale that the W, Z, Higgs bosons are directly emitted from such PBHs, while the large mass is also extended.We write: where κ ≃ 0.5, α ≃ 4 and ω ≃ 8 to sufficient precision, M * denote those PBHs' masses which are completing the evaporation at present, see the discussion below.When the PBH's mass goes down M q , the secondary emission would be triggered.Note that, in order to calculate the time-evolution of PBHs' masses, there are several characteristic mass scales of great interest.
• M a : First, let us see how long would a PBH with formation mass M f > 10M * fall into 10M * ≃ 5.1 × 10 15 g: where we also define M * ≡ (ϕ * t 0 ) 1/3 ≃ 5.07 × 10 14 f * 1.9 is the mass of a PBH currently evaporating if one neglects secondary emission once M (t) falls below M q , and we have used the approximation t f ≃ 0, and t 0 = 13.8Gyr is the lifetime of Universe.The value of f * is around 1.9 that can be derived by counting the quantum dofs of Hawking-radiated particles.Consider a formation mass M f = M a > 10M * , whose corresponding present mass is exact 10M * , we can derive • M c : Then, let us see how long would a PBH with formation mass M f > M q fall into M q ≃ 1.95 × 10 14 g: If the formation mass M f is larger than a critical mass scale M c , such a PBH would not reach M q at present, i.e., t q (M f ) > t 0 , which yields By using the range of M q , and then some intermediate value at the last step.So that only PBHs slightly larger than M * generate secondary emission at the present epoch.And the current mass is • M * : PBHs completing their evaporation today have M 0 = 0, the corresponding formation mass reads ≃ 5.14 × 10 14 g .
• M 1 : The formation mass whose current mass is exact M g :  4: The comparisons between approximations used in Ref. [70] and Eq. ( 25) for the mass evolutions of M f = 10 10 , 10 14 g, respectively.Our approximation ( 25) is more precise for the small formation mass M f .t=0 s t=10 5 s t=10 17   We can see that the mass scales M 1 and M * are very close to each other, and it is expected that the new threshold mass M g should not affect the large-mass region significantly, which is confirmed by the following results, see Fig. 4.
With above preparations, we can calculate the expressions for a PBH's mass at time t and its inverse M (M f , t) and M f (M, t), respectively, see Appendix C. The relation between M f and M is shown in Fig. 5, where we set the cosmic times t = {0, 10 5 , 10 17 } s after PBH formation.These curves are reasonable since those PBHs of M f ∼ 10 11 , 10 15 g are at the end of their lifetimes around t ∼ 10 5 , 10 17 s, respectively.
The time evolution of mass function of PBHs is the interplay of the Hawking radiation and the cosmic expansion.So, it is convenient to start with the comoving mass function which is defined as where dn c is the comoving number density over the mass interval (M, M + dM ).Note that the time-evolving mass function is implicitly dependent on time through the evaporating mass M (t).Using the chain rule for differentiation, we can relate the mass function at time t to the formation mass function n c (M f ) in the following form: and the explicit expressions are written as where the transfer function is calculated as Note that we do not need to take the differentiation to the Heaviside function Θ, as which is introduced to take into account the different evaporation stages, see Eq. (25).Then the physical mass function n(M (t)) is straightforward to calculate from the above, where is the formation physical mass function and λ(M, t) is shown in Eq. (37).
For the lognormal formation mass function, which is most common mass distribution of PBHs and the other types of mass functions (i.e. the power-law and the critical ones) can be rewritten in the form of lognormal distribution [73].The quantity β ≡ ρ PBH /ρ tot is the formation abundance of PBHs, and the total energy density of Universe at the formation epoch that can be calculated by the Friedmann equation, i.e., ρ tot (z here Ω r ≃ 10 −4 , Ω k ≃ 0, Ω m ≃ 0.315 and Ω Λ ≃ 0.6847 are normalized radiation, curvature, baryon and dark energy density parameters, respectively.z f is the formation redshift.For the purpose of numerical computation, we normalize the present scale factor a(t 0 ) = 1, and chose the Hubble parameter H 0 = 67.4km s −1 Mpc −1 .Using the redshift-time relation dt = − dz H(z)(1+z) , we can transfer the argument of any time-evolving function from the cosmic time t to the redshift z.By choosing parameters β = 10 −23 , σ = 1, M pk = 10 14 g, we plot the mass functions at different times t = 0, 10 5 , 10 17 s originated from the initial lognormal mass function, which is shown in Fig. 6.Fig. 7 plots the comparison between approximations used in [70] and Eq. ( 25) for the lognormal mass distribution, one can easily see that they are quite close to each other.M pk =10 14 g FIG.7: Two approximations used in Ref. [70] and Eq. ( 25) for the time-dependent mass function of M pk = 10 14 g, respectively.

IV. TEST OF COSMOLOGICAL MODELS VIA STANDARD TIMER
The cosmological redshift depends on the evolution of the scale factor a(t), and the cosmological time is the time elapsed from the primordial physical scale to the present physical scale which follows t ∼ H −1 0 .Hence, a cosmological model leaves the imprint on the redshift-time calibration which indicates that the cosmological parameters can be constrained by the standard timer.
From the definition of cosmological redshift 1 + z = a 0 /a and the Hubble parameter H(z) = ȧ/a, we obtain the redshift-time relation, Then, the redshift-time calibration (z 0 , z PBH , ∆t) between the PBH bubble and the redshift calibrator puts constraint on the Hubble parameter evolution as following, Considering the flat ΛCDM model as an example, H(z) can be expressed in following form, After extracting redshift-time calibration from a bunch of PBH bubbles covering from low redshifts to high redshifts, we can apply the Markov chain Monte Carlo (MCMC) simulation on flat ΛCDM model to constrains the cosmological parameters H 0 , Ω γ , Ω m and Ω Λ , In particular, we can evaluate the average Hubble parameter for two nearby redshifts z 1 and z 2 , where z 1 ∼ z 2 .In Eq. ( 40), we write Then we yield

V. CONCLUSION AND DISCUSSIONS
To conclude, we present a new approach to track the historical evolution of the Universe, which is named as the standard timer.We illustrate how the standard timer works by constructing the calibration between the redshift and physical time from the observed EM signals of PBH bubbles.The feasibility that PBH bubble can help track the evolution of the Universe is that their identical primordial origin produces the same initial mass function inside the bubble at the same epoch, then the redshift and the physical time calibration can be decoded from the observed EM signals and the internal mass function evolution.Also, the existence of PBH bubbles in the primordial Universe extends the narrow redshift windows in standard candle, standard ruler and standard siren, which makes the valid redshift range of standard timer cover from the primordial epoch to the present.As we have shown in Fig. 1 of [21], the EM signals from a PBH bubble with mass M bub = 10 38 g can be detected up z ∼ 10 3 .
In developing the method of standard timers, the key step is to calibrate the redshift and physical evolution time for two PBH bubbles.Following [63,64], we apply the inverse problem to Hawking radiation, by the least square method in Eq. ( 4), we can extract the redshifted mass function from observed EM signal spectrum.Then redshift ratio η can be obtained by matching redshifted mass functions for two PBH bubbles with error around 0.03%.Applying the obtained redshift ratio into evolution of mass function in Eq. ( 17), the calibration between redshift and physical time can be constructed with error around 0.2% in Fig. 2 and 0.1% in Fig. 3, respectively.This ensures the high precision of the standard timer.
Practically, various uncertainties should be taken into account in the observed EM signals from PBH bubbles, e.g., systematic uncertainties from the limited optical depth of the ultra-high gamma ray [74], measurement uncertainties in the gamma-ray spectrometry [75], etc.In data analysis, some regularization methods [69] in inverse problems should be applied to increase the accuracy of standard timer.In particular, under the low mass approximation, n(M ) ∼ M 2 in evolution of PBH mass function, which produces F (E) ∼ E −1 in the high-energy spectrum.This provides a specific template for filtering out the noise in gamma ray signals and evidence for potential PBH bubble candidate in gamma ray sources.
In general, there could exist other types of primordial phenomena that can track the evolution of the Universe, such as collision of cosmic strings [76] and domain walls [77].The connection between their primordial origin and the observed signals could build the one-to-one mapping between the internal physical evolution and the external evolution of the Universe, which are the potential tools in studying the history of the Universe.
Here, we define b(u) ≡ h ku 2 a( h ku ).Using the series expansion: we yield Here L(f ) is the Laplace transformation of f .The function f (u) is given by For a given PBH mass function n(M ), we can obtain photon number density emission rate as follows: Similar to the inverse black-body radiation problem, given an emission rate spectrum, finding the PBH mass function n(M ), is called the inverse Hawking radiation problem.The crucial difference between Hawking radiation and black-body radiation is the Γ s (E, M ) term, here we denote particle emission rate R(E) ≡ d 2 n/dtdE and consider the photon with s = 1, the (B3) can be written as where A and B are coefficients in the grey factor (B2).Therefore, we yield the relation: In order to obtain the mass function n(M ), we use the iteration method: Obviously , the iteration method works for the condition: So we need to check the convergence of the iteration relation Eq. (B6), we replace n m (M ) with n(M ) + ϵ m (M ) and yield Applying Eq. (B5) into Eq.(B8), we obtain Here, it is easy to find a positive real number ω < 1 to satisfy the inequality of Eq. (B9), we therefore prove that Hence, the mass function n(M ) can be solved by iteration method and Möbius inversion transformation we have discussed above.
Appendix C: The Analytic Formulas for M (t) and M f (M ) As discussed in Sec.III D, for the large formation mass M f , the time evolution remains unchanged, i.e., The next step is to express M f in terms of M and t.
for the condition − t × ωα −1 M 3 q + (1 − ωα −1 )M 3 g − ω M 3 function is written as [31,66] For simplicity, we set A = 1, ν = 0.35 and M peak = 10 14 g, respectively.Using the formulas ( 36) and ( 37), the time-evolving mass functions at various times are shown in Fig. 8.Note that the plots in Fig. 8 are more smooth than those of Ref. [70], since we have used the more precise approximation to the mass-loss rate, see Eq. ( 25).

Normalized mass function at z = 1 1 M
FIG.1:The normalized lognormal mass function with M pk = 10 15 g and σ = 1.The blue and red solid curves are the mass functions evolving from their formation to the redshifts z = 1 and z = 2, respectively.The blue and red dashed curves denote the normalized redshifted mass functions obtained from the inversion of Hawking radiation emitted at redshifts z = 1 and z = 2, respectively.

FIG. 2 :
FIG.2:The error in redshift calibration in terms of the error in η measurement.The redshift calibration is based on two PBH bubbles at different redshifts.We choose the lognormal mass function with M pk = 10 15 g and σ = 1 inside PBH bubbles and take z1 = 1 and z2 = 2 as an example in calculation.

FIG. 5 :
FIG. 5:The relations between the formation mass M f and the mass M at various cosmic times t = {0, 10 5 , 10 17 } s after their formations, denoted by the blue dashed, green dashed and red solid curves, respectively.The grey dotted line is obtained numerically by the public code BlackHawk[71,72].

FIG. 6 :
FIG.6:The mass functions at time t = 0, 10 5 , 10 17 s for the initial lognormal mass functions, denoted by the blue dashed, green dashed and red solid curves, respectively.We choose the parameters for the initial mass function: M pk = 10 14 g, β = 10 −23 .

3 ]
FIG.8:The mass functions at times t = 0, 10 5 , 10 17 s for the initial critical collapse mass function, denoted by the blue dashed, green dashed and red solid curves, respectively.We choose the parameters for the initial mass function: M pk = 10 14 g, β = 10 −23 .

TABLE I :
The approximate values of f (M ).