{'fraction_non_alphanumeric': 0.039264042573904014, 'fraction_numerical': 0.009020994905409353, 'mean_word_length': 4.720809355856692, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 1, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Cosmological Standard Timers from Primordial Black Hole Binaries

We propose that primordial black hole (PBH) binary systems can work as standard timers in tracking the evolution of the Universe. Through gravitational waves from monochromatic PBH binaries, the probability distribution on major axis and eccentricity from the same redshift is obtained. By studying the dynamical evolution of PBH binaries from the initial probability distribution to observed redshifted ones, the redshift-time calibration can be extracted, which can constrain cosmological models. A general formalism of the standard timer is further concluded based on the evolution of statistical distribution in dynamical systems.

For precise understanding in evolution of the Universe, various observational technologies have been proposed in studying the cosmological distance-redshift relation. Type Ia supernovae (SNe) produce consistent peak luminosity. The luminosity distance can be obtained by comparing the absolute and apparent magnitude of them, which works as the standard candle to provide luminosity distance-redshift relation [12,13]. Baryon acoustic oscillations (BAO) determine a fixed sound horizon, by observing sound horizons at different redshifts. The calibration between angular diameter and redshift is extracted, which serves as standard rulers [4,[14][15][16]. Gravitational waves (GWs) from binary systems and their electromagnetic counterparts provide the calibration between the luminosity distance and redshift as standard sirens [17]. However, there still exists puzzles in ΛCDM model, e.g., the Hubble tension [18][19][20] and the σ 8 tension [4]. Such tensions could be caused by unrecognized systematic uncertainty in measurements [21,22] or hidden new physics [23][24][25][26][27][28][29]. New observational methods are needed in cross checking the observed tensions.
Intuitively, the evolution of the Universe can be briefly characterized by the scale factor a(t) in the Friedmann-Robertson-Walker (FRW) metric, which is related to the cosmological redshift z by 1 + z(t) = a 0 /a(t), where a 0 is the present scale factor. Tracking z(t) provides another perspective in studying cosmological evolution (also used in the determination of the age of the Universe, see [30][31][32][33] for details), which could be achieved in the cosmological dynamic systems. Following their intrinsic dynamics, the physical evolution time from the initial state to the later state is attained. Meanwhile, redshift is decoded * qdingab@connect.ust.hk from their observable. Hence, redshift-time calibration could be constructed in cosmological dynamic systems. This approach is known as standard timers.
In the first study on standard timers [34], we have shown that through the Hawking radiation emitted from light primordial black hole (PBH) clusters, PBH stellar bubbles [35] can be used as standard timers. Due to the primordial origin of PBHs, the initial mass function of PBHs in clustering should be the same. With the emission of Hawking radiation, PBHs evaporate which deforms the mass function of PBHs. By studying the evolution of PBH mass function, its physical evolution time is attained. While observing the gamma ray spectrum from PBH stellar bubbles gives the redshifted PBH mass function where redshift is encoded. Hence, the calibration between redshift and physical evolution time of PBH mass function is constructed.
In this paper, we note that PBH binary systems can be used as standard timers. Under the assumption of random distribution of PBHs in space [36][37][38], the PBH binaries could decouple from the Hubble flow and form an identical initial probability distribution on the major axis a and the eccentricity e [39,40]. With the evolution of PBH binaries, the later probability distribution in PBH binary systems evolves accordingly. Through gravitational waves (GWs) emitted from PBH binaries, the probability distribution on major axis and eccentricity from the same redshift is obtained. By studying the evolution of PBH binaries from the initial probability distribution to later ones, cosmological redshift and physical evolution time in PBH binary systems are connected, which can constrain cosmological models. Considering the next generation of GW detectors, such as the Einstein Telescope [41], Laser Interferometer Space Antenna [42], can detect GWs from high redshift (z > 20), where PBH binaries may dominate the binary systems [43]. Detecting GWs from PBH binary systems at high redshift would own a pure GW background, which can construct the high-precision standard timer.
This paper is organized as follows. In Sec. II, we show how to construct the standard timer in PBH binary systems, including single parameter PBH binary systems in Sec. II A, multi-parameter PBH binary systems in Sec. II B and PBH binary systems without the initial probability distribution in Sec. II C. In Sec. III, the conclusion and discussions about the standard timer from PBH binary systems are given. In Appendix A & B, a formalism on constructing standard timers in general dynamic systems is shown.

II. THE STANDARD TIMER IN PBH BINARY SYSTEMS
In constructing the standard timer in PBH binary systems, two essential requirements are needed. One is the identical initial state of PBH binary systems. The other one is their later evolved state from the same redshift can be obtained. The former of requirements is achieved in the identical initial probability distribution on major axis and eccentricity [39,40] under the assumption of random distribution of PBHs in the space. The latter can be easily realized in local signal sources, such as the PBH cluster [35,44], where the redshift of signals from a local source should be the same. However, redshifts of PBH binaries are hardly classified, due to their being globally distributed in the Universe.
The potential monochromatic mass spectrum in PBH scenarios changes the story. It introduces PBH binary systems in the standard timer through GW channels. After collecting many GW signals from different redshifts, we can extract the redshifted chirp mass M z and mass ratio q from the GW waveform. Due to the unknown intrinsic chirp mass M, which follows M = M z /(1 + z) [45], the redshift of binary systems cannot be determined. However, if GW signals come from PBH binaries, under the assumption that the mass of PBHs is monochromatic, the mass ratio of PBH binaries follows q = 1 (also a number of BH binaries with q = 1 are found in GWTC-2, see [46] for details), and GW signals emitted from PBH binaries at the same redshift can give an identical redshifted chirp mass. Then, PBH binaries can be classified into different redshift shells based on their redshifted chirp mass. As the result, the standard timer can be constructed by comparing the initial probability distribution on major axis a and eccentricity e and later ones from the same redshift.

A. A toy model in single parameter PBH binary systems
By assuming the mass of PBHs is monochromatic, the state of PBH binary systems can be described by a probability distribution on major axis a and eccentricity e, which is dP/dade. For an intuitive understanding on how standard timers work in PBH binary systems, we start with circular binary systems. Then states of PBH binary systems only depend on the single parameter major axis a, which is dP/da. In studying the evolution of its probability distribution, we have where, a i and a t denote the major axis at the initial physical time t i and later physical time t, respectively. S(a; t) denotes the probability distribution on major axis a at physical time t and dP/da i is the initial probability distribution S(a; t i ). In connecting S(a; t) with S(a; t i ), we consider time evolution of the major axis, following [47] da dt = − 128 5 where M PBH is the monochromatic mass of PBHs, G and c are the Newton's constant and the speed of light, respectively. Integrating Eq. (2) from the initial physical time t i to later physical time t, we have a 4 i = a 4 t + δ 4 (∆t), where δ 4 is defined as δ 4 (∆t) ≡ 512G 3 M 3 PBH ∆t/5c 5 and ∆t is the physical evolution time ∆t = t − t i . Fixing the evolution time ∆t, the evolution of major axis da i /da t can be obtained, then the probability distribution of single parameter PBH binary systems from an identical redshift shell can be expressed as following With the expansion of the Universe, the observed major axis a z is redshifted by a z = (1 + z)a, which can be found in the Kepler's third law a z ∼ where M z and f z are the observed mass and GW frequency in binary systems, respectively. Therefore, the observed probability distribution follows Here, subscript o denotes the observational quantity. As above, we have the relation a 4 i (z) = a 4 z + δ 4 (∆t z ), where ∆t z depends on the redshift of PBHs binaries. Therefore, a correct redshift should be firstly obtained in calibrating the redshift-time relation from the observed probability distribution. We consider the condition a 4 z δ 4 (∆t z ), which infers a i (z) a z . In this large major axis limit, the observed probability distribution becomes Then the redshift can be numerically resolved from the where a is chosen in the large major axis limit. After obtaining the redshift, probability distributions can recover from redshifted ones by S(a; t) = (1 + z)S o (a(1 + z); t). Then, physical evolution time can be extracted in the condition a 4 t δ 4 (∆t), which indicates a i δ(∆t). In this small major axis limit, we have the log probability distribution from Eq. (3) as following log S(a; t) = log S(δ(∆t); t i ) Then δ(∆t) can be extracted from S(δ(∆t); t i )/δ 3 (∆t) in Eq. (6). δ(∆t) depends on the mass of PBHs M PBH and physical evolution time ∆t. In above calculation, redshift is obtained from the probability distribution in the large major axis limit, which helps determine the mass of PBHs from the redshifted mass M PBH = M z /(1+ z). Then, physical evolution time ∆t can be resolved from obtained δ(∆t). After obtaining the redshift and physical evolution time, the calibration between redshift and time can be constructed in single parameter PBH binary systems. In general, the standard timer requires two properties in the evolution of dynamical systems, we take the single parameter PBH binary system as an example in Fig The evolution of single parameter PBH binary systems dai/da. We assume the mass of PBHs is 30M and physical evolution time is the standard cosmological time from z = 3000 to the redshift denoted in the figure. The solid curve is the intrinsic evolution of dynamical systems and the dashed curve is the observed evolution of dynamical systems. the large major axis tail and the rapid evolution part in the small major axis tail. The flat constant part is the region where time evolution is negligible, the redshift can be extracted by comparing the redshifted and the initial probability distributions in this region. After obtaining the redshift, intrinsic evolution functions (solid lines) can recover from redshifted ones (dashed lines) in Fig. 1. The rapid evolution part describes the physical time evolution effectively changes the major axis of PBH binaries, then physical evolution time can be extracted in Eq. (6).

B. A practical model in PBH binary systems
A practical description of monochromatic PBH binary systems needs two parameters, major axis a and eccentricity e. In multi-parameter probability distributions, the evolution of probability distribution of PBH binary systems dP/dade from an identical redshift shell can be Here, dP/da i de i is the initial probability distribution of PBH binary systems, and J(a, e, ∆t) is the Jacobian of two-parameter PBH binary systems after the evolution of physical evolution time ∆t, which connects the initial and later probability distribution. In calculating J(a, e, ∆t), we consider the time evolution of parameters in PBH binaries, following [47] da dt = − 128 5 Due to the expansion of the Universe, the cosmological redshift is introduced in the observed probability distribution, which is where redshifted major axis becomes a z = (1+z)a. However, the redshift does not leave imprints on the eccentricity e, it can be found in time evolution of the eccentricity that the redshift effect in major axis, mass and time cancel with each other in Eq. (8), which results in no redshift effect appearing in observed eccentricity. In measuring eccentricity in binary systems, the precision is not very high due to the lack of suitable GW waveform templates [48]. Therefore, we only consider probability distributions on the major axis dP/da for a practical numerical solution, which can be obtained from Eq. (9) as following In the numerical solution of evolution in probability distributions, we consider two types of initial probability distributions in PBH binary systems. One is the Gaussian distribution localized at a specific major axis and eccentricity. After studying the Gaussian distribution, a general distribution can be decomposed into Gaussian distributions [49,50]. The other probability distribution in PBH binary systems is chosen from [39,51], where f PBH is the present energy density fraction of PBHs in the dark matter andx is the physical mean separation of PBHs at matter-radiation equality. Based on different types of initial probability distributions, a numerical study on the evolution of probability distribution of PBH binaries on the major axis a following Eqs. (8)(10) is shown in Fig .2.
As mentioned in Sec. II A, various regions of probability distributions behave differently. In the large major axis limit, the evolution of probability distribution is negligible, the redshift can be obtained by comparing the initial probability distribution (black line) and redshifted probability distributions (dashed lines) in Fig .2, which can be further numerically resolved from the equation (1 + z)dP/da z = dP/da i . However, the large major axis tail in observed probability distributions can hardly be obtained due to its extremely low gravitational wave frequency, as shown in the right panel of Fig. 2. In this case, we can use numerical solutions as templates of probability distributions to match observational results. After obtaining the redshift, observed probability distributions dP/da z can be restored to intrinsic probability distributions dP/da t following dP/da t = (1 + z)dP/da z . Then physical evolution time ∆t can be extracted from the numerical solution in the small major axis limit, where the evolution of probability distribution dominates. Accordingly, the redshift-time calibration is constructed in multi-parameter PBH binary systems.
Generally, the mass of PBHs is very essential in describing the evolution of PBH binary systems, when the mass function of PBHs is not monochromatic. The evo-lution of parameters follows [47] da dt = − 64 5 Here, m 1 and m 2 are the mass of PBHs in the binary. Therefore the probability distribution of general PBH binary systems should be described as dP/dadedM G , where M G ≡ (m 1 m 2 (m 1 + m 2 )) 1/3 . The evolution of probability distribution can also be numerically studied as we have discussed above. However, considering an extended mass function in PBH binary systems, it is hard to determine whether collected GWs come from the same redshift or not, because the redshifted mass depends on the mass of PBHs in binary and cosmological redshift, due to the mass uncertainty in extended mass spectrum of PBHs, the redshift is hardly determined. In inflationary scenario, the lognormal type mass function is a natural result of PBHs mass function [52][53][54], which is described by characteristic mass M c and width of mass spectrum σ. In the small σ limit, the mass function of PBHs can be approximated by a monochromatic mass function, the redshift can be determined as we have discussed in the monochromatic case. Otherwise, the determination of redshift needs further studies.
C. The standard timer without the initial probability distribution In constructing standard timers in PBH binary systems, the initial probability distribution of PBH binaries plays an important role in extracting the redshift and physical evolution time by comparing it with observed redshifted probability distributions. However, the initial probability distribution on major axis and eccentricity of PBH binaries is indeterminate, due to the unknown mass spectrum and space distributions of PBHs, etc. A number of previous works [40,51,55] have been done in calculating the initial probability distribution in different scenarios.
In practice, the information about PBHs is insufficient at present, such as the PBH mass, the initial probability distribution on major axis and eccentricity, etc. However, under the assumption of random distribution of PBHs in space, initial probability distributions of PBH binaries at different redshifts are same, which ensures one-to-one correspondence between later evolution probability distributions and physical evolution time. As the result, we cannot extract the redshift by comparing the initial and redshifted probability distribution in Fig. 2. Instead, the redshift ratio η can be obtained from two redshifted probability distributions, following dP/da i = (1 + z 1 )dP/da z1 = (1 + z 2 )dP/da z2 in large major axis limit, which further gives where η ≡ (1+z 2 )/(1+z 1 ). In order to obtain the redshift of observed GW signals, we can assume the cosmological evolution between two redshift shells follows the standard ΛCDM cosmology. Then, we assume the redshift of one observed probability distribution isz 1 and the redshift of the other one isz 2 = (1+z 1 )η −1. Under the assumption ofz 1 andz 2 , the PBH mass M PBH and physical evolution time ∆t can be obtained as we have discussed in Sec. II B. The cosmological time betweenz 1 andz 2 can be calculated as t z = z2 z1 dz/H(z)(1 + z), where H(z) is the Hubble parameter along the line of sight. Then proper redshiftz 1 is chosen such that physical evolution time is same as cosmological time between two redshifts as following After obtaining the redshift of PBH binaries from same redshift shell z 0 , the PBH mass can recover from the observed redshifted mass by M PBH = M z /(1 + z 0 ) and redshift of other probability distributions can be extracted by z PBH = (1 + z 0 )η − 1. Then we can numerically solve physical evolution time ∆t between z 0 and z PBH as discussed in Sec. II B. Consequently, the redshift-time calibration (z 0 , z PBH , ∆t) is obtained.
Furthermore, cosmological models can be tested in standard timers. Considering the cosmological redshifttime relation dz/dt = −(1 + z)H(z), we apply the obtained redshift-time calibration (z 0 , z PBH , ∆t) from standard timers, which gives Take the flat ΛCDM model as an example, H(z) = H 0 Ω γ (1 + z) 4 + Ω m (1 + z) 3 + Ω Λ . After constructing the redshift-time calibration in PBH binary systems from the primordial Universe to the present Universe, the Markov chain Monte Carlo (MCMC) simulation can be applied on the flat ΛCDM model in constraining the Hubble parameter H 0 , energy density fraction of radiation Ω γ , matter Ω m and cosmological constant Ω Λ .

III. CONCLUSION AND DISCUSSIONS
To summarize, we propose that PBH binary systems can work as standard timers to record the evolution of the cosmological redshift z(t). Under the assumption of random distribution of PBHs in space, PBH binary systems have an identical initial probability distribution on major axis and eccentricity. By studying the evolution of the probability distribution in binary systems, the physical evolution time between the initial and later probability distribution can be extracted. Then the redshift-time relation can be constructed by studying the probability distribution of PBH binary systems at different redshifts. In order to obtain the probability distribution on major axis and eccentricity from the identical redshift shell. We assume that PBH mass is monochromatic, through GWs produced from PBH binaries, its redshifted chirp mass can be obtained in GW waveforms, the PBH binaries from the same redshift have the same redshifted chirp mass and the mass ratio follows q = 1, then we extract the redshifted probability distribution on major axis and eccentricity from the same redshift.
For demonstrating how standard timers work in PBH binary systems, we perform an analytically studied toy model in single parameter PBH binary systems where eccentricity is set as e = 0 and a numerically studied practical model in non-circular PBH binary systems. We show that the redshift can be determined by comparing the initial and redshifted probability distribution at the large major axis limit and the physical evolution time can be obtained by comparing the initial and recovered redshifted probability distribution at the small major axis limit. Considering the initial probability distribution on major axis and eccentricity in PBH binary systems is indeterminate, the redshift of observed probability distributions cannot be directly obtained. We assume cosmological time between two redshifts follows the standard cosmology, then proper redshift of PBH binary systems should be chosen when the physical evolution time between two redshifted probability distributions equals its cosmological time, which further constructs standard timers in PBH binary systems without initial conditions.
In the above discussions, we mainly focus on PBH binary systems with a monochromatic mass spectrum, which helps classify the redshift of PBH binaries. In a general description of PBH binary systems, an extended mass spectrum should be taken into consideration, which is shown at the end of Sec. II B. The standard timer can be constructed by a numerical study in a probability distribution on major axis, eccentricity and mass of PBHs in binaries dP/dadedM G . However, an extended mass spectrum of PBHs could cause difficulties in redshift classification of PBH binaries, which needs further studies in redshift identification, e.g., [56][57][58] in standard sirens.
In general, the cosmological standard timer can be constructed based on dynamical systems in the Universe. Due to the same formation mechanism of dynamical systems, the statistical distribution of their initial states can be set as the standard reference, through the evolution mechanism of their statistical distribution, the elapsed time in the standard timer is evaluated. Meanwhile, cosmological redshift is encoded in the observable from dynamical systems. For signals from individual sources locally, the redshifted statistical distribution in dynamical systems from the same redshift can be obtained, and further gives their redshift by comparing with the initial state (see [34,35] for more details). For signals from global sources, as we have discussed in this article, GWs from PBH binaries globally, the redshifted statistical distribution from the same redshift can be extracted according to their redshifted parameters, and hence obtain their redshift. Consequently, the redshift-time calibration is constructed in a general dynamic system (see Appendix. A & B for the detailed formalism).
In constructing standard timers in dynamical systems, we need to set a particular condition of dynamical systems as a standard reference. For example, the standard reference in standard candles is consistent peak luminosity produced by Type Ia supernovae (SNe), and the standard reference in standard rulers is the fixed baryon acoustic oscillation (BAO) scale that the sound wave can travel before the recombination. Generally, initial states of dynamical systems are uncertainty under the Gaussian distribution of perturbations. However, the statistical distribution of initial states in dynamical systems can be unique due to the same physical mechanism behind them, which can be set as a standard reference. With the standard reference, the physical evolution time and redshift can be extracted by studying the evolution of observed dynamical systems, which can help calibrate z(t).
For simplicity, we start with a single parameter dynamical system, whose time evolution follows dM/dt = −f (M ). Here M is the observable physical parameter that characterizes the dynamical system and f (M ) is its time derivative function. The statistical distribution of the single parameter dynamical system S(M ; t) can be described as where M t ≡ M (t), N is the statistic of the distribution of dynamical systems. In order to trace the historical evolution of dynamical systems, a standard initial distribution is essential. Eq. (A1) can be written as Here, dN/dM i is the initial statistical distribution of dynamical systems S(M ; t i ). dM i /dM t describes the evolution of the dynamical system, which can be further expressed by its time evolution dM/dt = −f (M ), which gives Here, function g(M ) is an antiderivative of function 1/f (M ). Then, the evolution of dynamical systems can be written as .
Here, g (M ) ≡ dg(M )/dM and g −1 denotes the inverse function of g(M ). As the result, Eq. (A2) can be further expressed as .
The physical evolution time ∆t can be extracted by giving an initial statistical distribution dN/dM i in Eq. (A5).
In the observational aspect, the statistical distribution of the dynamical system is deformed, due to observed physical parameter M is redshifted by the cosmological expansion, which gives the observational distribution S o (M z ; t) (subscript o denotes the observational quantity) as .
Here, M z denotes the redshifted physical parameter, such as redshifted photon energy E z = E/(1 + z) and redshifted chirp mass in binary black hole systems M z = (1+z)M. dN/dM i (z) characterizes the cosmological redshift effect in the initial statistical distribution. Following Eq. (A3), we have g(M i (z)) = g(M z ) + ∆t z . In order to extract the redshift-time calibration, we consider two cases in Eq. (A6). For a fixed evolution time, the first case is g(M z ) ∆t z , which makes sure the time evolution is negligible and redshift can be extracted by comparing the redshifted physical parameter M z with the initial physical parameter M i in the initial statistical distribution. The second case is g(M z ) ∆t z , where ∆t z dominates in the redshifted physical parameter, which gives g(M i (z)) ∆t z . Then ∆t z can be extracted in following expression Above all, we have discussed the formalism of a standard timer in an observable dynamical system S(M ; t). However, this formalism does not apply to the case that M is not an observable of dynamical systems, meanwhile, the signals produced from them is an observable, e.g., electromagnetic waves and gravitational waves. In these cases, we consider the following integral equation, Here, K(E, M ) is the kernel function which transfers an unobservable distribution S(M ; t) to an observable distribution P (E; t). S(M ; t) can be extracted by an inverse integral equation where K −1 (E, M ) is the inverse kernel function of K(E, M ). Due to the cosmological expansion, the observed physical parameter E is redshifted to E z . Therefore, the observable becomes where Z 1 function describes the redshift effect in the observable E z . In order to construct the redshift-time relation, a redshift term need to appear in S(M ; t), which requires the connection between E and M in the kernel function, for instance, the primary Hawking radiation kernel follows H(E(1 + z), M ) = H(E, M (1 + z)) [34]. Therefore, we assume the kernel function follows Here, Z 1 and Z 2 function describe how the redshift term transfers from E z to M in kernel function. Then, Eq. (A10) can be written as As the result, S o (Z 2 (M ); t) is given by As we discuss in Eq. (A7), S o (Z 2 (M ); t) can be expressed in two conditions in Eq. (A14), which further gives the redshift-time calibration.
Appendix B: The standard timer from multi-parameter dynamical systems In general, we consider multi-parameter dynamical systems in the Universe, whose statistical distribution can be expressed as where M denotes the n-dimensional physical parameter vector which characterizes the dynamical system. By introducing the initial statistical distribution as the standard reference, Eq. (B1) can be written as Here, J is the Jacobian of the dynamical system which is defined as J ij ≡ ∂M i (t i )/∂M j (t). With the reference of the initial statistical distribution, the physical evolution time ∆t can be extracted from the determinant of the Jacobian det J(M, ∆t). However, due to strong coupling among different parameter components in its time evolution dM/dt = −f (M), the general analytical expression of the Jacobian element ∂M i (t i )/∂M j (t) can hardly be found, which indicates the numerical solution of det J(M, ∆t) is essential in extracting the physical evolution time ∆t.
In the observational perspective, the redshift term caused by the cosmological expansion also appears in the statistical distribution of multi-parameter dynamical systems as following where M z denotes the redshifted n-dimensional physical parameter vector and dN/d n M i (z) characterizes the redshifted initial statistical distribution. As we have discussed in Appendix. A, we consider two cases in extracting the redshift-time calibration. One case is that in the parameter space where the time evolution of parameters is negligible compared with their initial value, which gives det J(M z , ∆t z ) 1. Then the redshift z can be obtained by comparing the observed statistical distribution S o (M z ; t) dN/d n M i (z) with the initial one. The other case is that in the parameter space where the time evolution of parameters dominates their initial value, where physical evolution time ∆t can be extracted from the numerical solution.
In the scenario that S(M; t) is not observable, we consider the observable P (E; t) as following where, K(E, M) is the transfer kernel which transfers an unobservable S(M; t) to an observable P (E; t), V is the integral region of n-dimensional parameter M. With the expansion of the Universe, the redshift effect appears in the observable in the following form As we have shown in Eq. (A11), we introduce Z 1 and Z 2 function to transfer a redshift term from E z to M, which gives P o (E z ; t) = V K o (E z , Z 2 (M))S o (Z 2 (M); t)d n Z 2 (M) .

(B6)
As the result, the unobservable statistical distribution S o (Z 2 (M); t) can be obtained by an inverse kernel transformation as Eq. (A13), However, the analytical form of the inverse kernel in multi-parameter dynamical systems K o (E, M) −1 could hardly be found, which needs further numerical methods, e.g., the method for the least squares problem [59,60]. After obtaining S o (Z 2 (M); t), the redshift-time calibration can be extracted as Eq. (A14).