⇦Chinese Physics Letters, 2024, Vol. 41, No. 5, Article code 059501 Cosmology-Independent Photon Mass Limits from Localized Fast Radio Bursts by Using Artificial Neural Networks Jing-Yu Ran (冉景遇)1,2†, Bao Wang (王宝)1,2†, and Jun-Jie Wei (魏俊杰)1,2* Affiliations 1Purple Mountain Observatory, Chinese Academy of Sciences, Nanjing 210023, China 2School of Astronomy and Space Sciences, University of Science and Technology of China, Hefei 230026, China Received 8 March 2024; accepted manuscript online 1 May 2024; published online 15 May 2024 †These authors contributed equally to this work.
*Corresponding author. Email: jjwei@pmo.ac.cn Citation Text: Ran J Y, Wang B, and Wei J J 2024 Chin. Phys. Lett. 41 059501    Abstract A hypothetical photon mass $m_{\gamma}$ can produce a frequency-dependent vacuum dispersion of light, which leads to an additional time delay between photons with different frequencies when they propagate through a fixed distance. The dispersion measure and redshift measurements of fast radio bursts (FRBs) have been widely used to constrain the rest mass of the photon. However, all current studies analyzed the effect of the frequency-dependent dispersion for massive photons in the standard $\Lambda$CDM cosmological context. In order to alleviate the circularity problem induced by the presumption of a specific cosmological model based on the fundamental postulate of the masslessness of photons, here we employ a new model-independent smoothing technique, artificial neural network (ANN), to reconstruct the Hubble parameter $H(z)$ function from 34 cosmic-chronometer measurements. By combining observations of 32 well-localized FRBs and the $H(z)$ function reconstructed by ANN, we obtain an upper limit of $m_{\gamma} \le 3.5 \times 10^{-51}$ kg, or equivalently $m_{\gamma}\le2.0 \times 10^{-15}$ eV/c$^2$ ($m_{\gamma} \le 6.5 \times 10^{-51}$ kg, or equivalently $m_{\gamma} \le 3.6 \times 10^{-15}$ eV/c$^2$) at the $1\sigma$ ($2\sigma$) confidence level. This is the first cosmology-independent photon mass limit derived from extragalactic sources.
cpl-41-5-059501-fig1.png
cpl-41-5-059501-fig2.png
cpl-41-5-059501-fig3.png
 DOI:10.1088/0256-307X/41/5/059501 © 2024 Chinese Physics Society Article Text 1. Introduction. As a fundamental postulate of Maxwell's electromagnetism and Einstein's special relativity, the principle of the constancy of the speed of light implies the masslessness of photons, as elucidated by the particle-wave duality. This principle is also embraced within the framework of general relativity (GR). However, even a very minute photon mass, if existent, would necessitate new physical theories, such as the renowned de Broglie–Proca theory,[1,2] the model of massive photons as an explanation for dark energy,[3] and other new ideas in the standard-model extension with massive photons.[4] Consequently, achieving precise constraints on the rest mass of photon, denoted as $m_{\gamma}$, remains an imperative. Numerous experimental and observational constraints on $m_{\gamma}$ have been derived from various effects resulting from the hypothetical nonzero photon mass. These include tests of Coulomb's inverse square law,[5] the Cavendish torsion balance,[6,7] gravitational deflection of electromagnetic waves,[8] mechanical stability of magnetized gas in galaxies,[9] magneto-hydrodynamic phenomena of the solar wind,[10-13] Jupiters magnetic field,[14] the spindown of a pulsar,[15] among others. However, these constraints are contingent upon specific theories of massive photons and are thus dynamic tests in nature. In contrast, a kinematic test based on the dispersion effect of the speed of light in vacuum is more inherently pure.[16,17] The dispersion relation for massive photons is governed by \begin{equation} E^2 = p^2c^2 + m_{\gamma}^2c^4, \tag {1} \end{equation} where $p$ represents the momentum. The group velocity of a photon with energy $E=h\nu$ can be expressed as \begin{equation} v = \frac{\partial E}{\partial p} = c\sqrt{1-\frac{m_{\gamma}^2c^4}{E^2}} \approx c\Big( 1-\frac{1}{2}\frac{m_{\gamma}^2c^4}{h^2\nu^2} \Big), \tag {2} \end{equation} where $h$ is the Planck constant, and the last term is valid when $m_{\gamma}\ll h\nu/c^2\simeq 7\times 10^{-42}(\nu/{\rm GHz})$ kg. Equation (2) suggests that if two photons with different frequencies are emitted simultaneously from the same source, they will be observed at different times due to their different velocities. Fast radio bursts (FRBs) provide an excellent celestial laboratory for detecting this dispersion effect, owing to their characteristics of (i) short time durations, (ii) long propagation distances, and (iii) low-frequency emissions. Based on the vacuum dispersion method, several studies have used FRBs to put stringent upper limits on the photon mass.[18-28] Since FRBs originate at cosmological distances, the cosmic expansion rate $H(z)$ has to be considered in constraining the photon mass $m_{\gamma}$. In all previous studies, the required $H(z)$ information is calculated within the standard $\Lambda$CDM cosmological model. It should, however, be emphasized that $\Lambda$CDM is rooted in the framework of GR, which also embraces the postulate of the constancy of light speed. Thus, there is a circularity problem in constraining the photon mass. To address this problem, one has to determine the Hubble parameter $H(z)$ in a cosmology-independent way. In this work, we propose a novel nonparametric approach, utilizing the artificial neural network (ANN) technology, to reconstruct a smooth $H(z)$ function that best approximates the discrete $H(z)$ measurements. By combining observations of well-localized FRBs and the $H(z)$ function reconstructed by ANN, we aim to establish a cosmology-independent constraint on the rest mass of the photon. In the following, we first introduce the theoretical framework and the observational data utilized for constraining the photon mass. Then, we present the cosmology-independent constraints on the photon mass and other relevant parameters. Finally, a brief summary and discussions are given. 2. Analysis Method and Data. 2.1 Theoretical framework. In our analysis, the observed time delay between different frequencies for FRBs is attributed to two main factors: (i) the propagation of photons through the plasma distributed between the source of the FRB and the observer, and (ii) the non-zero rest mass of the photon. (1) Due to the interaction between the plasma and electromagnetic waves, photons with higher frequency $\nu_{\rm h}$ travel faster than those with lower frequency $\nu_{\rm l}$, resulting in a dispersion effect. The time delay $\Delta t_{\scriptscriptstyle{\rm DM}}$ between photons of different frequencies caused by this dispersion effect can be described as[29,30] \begin{align} \Delta t_{\scriptscriptstyle{\rm DM}}&=\int \frac{\nu_{\rm p}^2}{2c}(\nu_{\rm l}^{-2}-\nu_{\rm h}^{-2})dl \notag\\ &=\frac{e^2}{8\pi^2m_{\rm e}\epsilon_0c}(\nu_{\rm l}^{-2}-\nu_{\rm h}^{-2}){\rm DM_{\rm astro}}, \tag {3} \end{align} where $\nu_{\rm p}\equiv \sqrt{\frac{n_{\rm e}e^2}{4\pi^2m_{\rm e}\epsilon_0}}$ is the plasma frequency with the number density of electrons ($n_{\rm e}$), the charge of an electron ($e$), the mass of an electron ($m_{\rm e}$), and the permittivity of vacuum ($\epsilon_0$). Here DM$_{\rm astro}$ is the dispersion measure (DM) contributed by the plasma, which is defined as the integral of the number density of electrons $n_{\rm e}$ along the line of propagation, given by DM$_{\rm astro} \equiv \int n_{\rm e} dl$. For an extragalactic source, DM$_{\rm astro}$ can be primarily divided into four components: the contributions from the Milky Way's interstellar medium (DM$_{\scriptscriptstyle{\rm ISM}}^{\scriptscriptstyle{\rm MW}}$), the Galactic halo (DM$_{\rm halo}^{\rm MW}$), the intergalactic medium (IGM; DM$_{\scriptscriptstyle{\rm IGM}}$), and the host galaxy (DM$_{\rm host}$). Therefore, DM$_{\rm astro}$ is given by \begin{align} {\rm DM_{\rm astro}}={}&{\rm DM_{\scriptscriptstyle{\rm ISM}}^{\scriptscriptstyle{\rm MW}}}+{\rm DM_{\rm halo}^{\scriptscriptstyle{\rm MW}}}\notag \\ &+{\rm DM_{\scriptscriptstyle{\rm IGM}}}+\frac{{{\rm DM}_{\rm host,0}}}{1+z}, \tag {4} \end{align} where the factor $(1+z)$ converts the DM component of the host galaxy in the rest frame, denoted as DM$_{\rm host,\,0}$, to the observed value DM$_{\rm host}$. (2) According to Eq. (2), the time delay $\Delta t_{m_{\gamma}}$ between massive photons of high and low frequencies can be expressed as \begin{equation} \Delta t_{m_{\gamma}}=\frac{1}{2}\Big(\frac{m_{\gamma}c^2}{h}\Big)^2(\nu_{\rm l}^{-2}-\nu_{\rm h}^{-2})H_{\gamma}(z), \tag {5} \end{equation} where $H_{\gamma}(z)$ is a newly defined function related to redshift, \begin{equation} H_{\gamma}(z)=\int_{0}^{z}\frac{(1+{z}')^{-2}}{H(z')}dz', \tag {6} \end{equation} where $H(z)$ is the Hubble parameter at redshift $z$. Note that here $H_{\gamma}(z)$ differs from the dimensionless redshift function defined in previous works, and it is in units of [km$\cdot$s$^{-1}\cdot$Mpc$^{-1}$]$^{-1}$. In our analysis, the total observed time delay is \begin{equation} \Delta t_{\rm obs}=\Delta t_{\scriptscriptstyle{\rm DM}}+\Delta t_{m_{\gamma}}. \tag {7} \end{equation} Observationally, the time delay of all FRBs exhibits a frequency dependence of $\nu^{-2}$, while both $\Delta t_{\scriptscriptstyle{\rm DM}}$ and $\Delta t_{m_{\gamma}}$ follow a $\nu^{-2}$ behavior. Hence, it is natural for us to analogously define the equivalent DM arising from the massive photons as[31] \begin{equation} {\rm DM}_{\gamma}(z)\equiv \frac{4\pi^2m_{\rm e}\epsilon_0c^5}{h^2e^2}H_{\gamma}(z)m_{\gamma}^2. \tag {8} \end{equation} Thus, the observed DM obtained from fitting the $\nu^{-2}$ behavior in the total frequency-dependent time delay can be written as \begin{equation} {\rm DM_{\rm obs}}={{\rm DM}_{\rm astro}}+{\rm DM}_{\gamma}. \tag {9} \end{equation} Once we are able to properly estimate the value of each DM term in Eq. (4), DM$_{\gamma}$ can then be effectively extracted from DM$_{\rm obs}$, thereby providing an upper limit on the photon mass $m_{\gamma}$. The DM$_{\scriptscriptstyle{\rm ISM}}^{\scriptscriptstyle{\rm MW}}$ term arising from the ionizing medium around our galaxy is well modeled by some galactic electron distribution models. Here we adopt the NE2001 model[32] for its wide application. The value of DM$_{\rm halo}^{\scriptscriptstyle{\rm MW}}$ is challenging to estimate accurately, but it has been expected to lie within the range of 50–80 pc$\cdot$cm$^{-3}$.[33,34] Here, we conservatively adopt a Gaussian prior of DM$_{\rm halo}^{\scriptscriptstyle{\rm MW}}=65\pm 15$ pc$\cdot$cm$^{-3}$.[27,35] The DM due to the IGM, i.e., DM$_{\scriptscriptstyle{\rm IGM}}$, depends on the cosmological model and is largely influenced by the number of halos intersected along the propagation path. Due to density perturbations on large-scale structures, it is challenging to calculate the precise value of DM$_{\scriptscriptstyle{\rm IGM}}$ for individual FRBs. Instead, we typically calculate the average value of DM$_{\scriptscriptstyle{\rm IGM}}$ using[36] \begin{equation} \langle {\rm DM_{\scriptscriptstyle{\rm IGM}}}(z) \rangle=\frac{21c\varOmega_{\rm b}H_0^2f_{\scriptscriptstyle{\rm IGM}}}{64\pi Gm_{\rm p}}H_{\rm e}(z), \tag {10} \end{equation} where $\varOmega_{\rm b}$ is the baryon density parameter at the present day, $H_0$ is the Hubble constant, $f_{\rm{IGM}}\simeq 0.83$ is the fraction of baryon in the IGM,[37] $G$ is the gravitational constant, and $m_{\rm p}$ is the mass of proton. The redshift-dependent function $H_{\rm e}(z)$ (in units of [km$\cdot$s$^{-1}\cdot$Mpc$^{-1}$]$^{-1}$) is defined by \begin{equation} H_{\rm e}(z)\equiv \int_{0}^{z}\frac{(1+z')}{H(z')}dz'. \tag {11} \end{equation} However, the actual value of DM$_{\scriptscriptstyle{\rm IGM}}$ may deviate from the average due to the inhomogeneity of IGM. To account for this variability, we employ a one-parameter model to simulate the probability distribution of DM$_{\scriptscriptstyle{\rm IGM}}$,[24,38,39] which is given by \begin{equation} P_{\rm{IGM}}({\rm DM_{\scriptscriptstyle{\rm IGM}}}|z)=A\varDelta ^{-\beta}\exp \Big[-\frac{(\varDelta^{-\alpha}-C_0)^2}{2\alpha^2s^2}\Big], \tag {12} \end{equation} where $\varDelta \equiv {{\rm DM}_{\scriptscriptstyle{\rm IGM}}}/\langle {{\rm DM}_{\scriptscriptstyle{\rm IGM}}}\rangle>0$, $A$ is a normalization constant to ensure the integral of $P_{\rm{IGM}}({{\rm DM}_{\scriptscriptstyle{\rm IGM}}}|z)$ to be unity, $C_0$ is chosen to satisfy $\langle \varDelta \rangle=1$, and $\alpha=\beta=3$.[39] The fractional standard deviation of ${\rm DM}_{\scriptscriptstyle{\rm IGM}}$ is approximated as $s=Fz^{-0.5}$, where $F$ is the free parameter that quantifies the strength of baryon feedback.[39] Both semi-analytic models and numerical simulations of the IGM and galaxy halos show that the probability distribution of DM$_{\scriptscriptstyle{\rm IGM}}$ can be well fitted by the quasi-Gaussian form [i.e., Eq. (12)].[39,40] Therefore, this analytic form is adopted in our analysis. The DM$_{\rm host}$ term is contributed by the source environment and the interstellar medium of the host galaxy, implying that the value of DM$_{\rm host}$ may vary significantly across different sources. To model this variability, we adopt a log-normal distribution for the probability density function of DM$_{\rm host}$ with the expression[39] \begin{align} &P_{\rm host}({\rm DM}_{\rm host}|\mu,\,\sigma_{\rm host})\notag\\ ={}&\frac{1}{\sqrt{2\pi}{\rm DM}_{\rm host}\sigma_{\rm host}} \times \exp \Big[-\frac{(\ln {{\rm DM}_{\rm host}}-\mu)^2}{2\sigma_{\rm{host}}^2}\Big], \tag {13} \end{align} where $\mu$ and $\sigma_{\rm{host}}$ are free parameters used to estimate the mean and standard deviation of $\ln {{\rm DM}_{\rm host}}$, respectively. With the above analyses, we can formulate the joint likelihood function for a sample of localized FRBs as \begin{equation} \mathcal{L}=\prod_{i=1}^{N}P_i({\rm DM}_{{\rm E},\,i}|z_i), \tag {14} \end{equation} where $N$ is the number of FRBs and $P_i({\rm DM}_{{\rm E},\,i}|z_i)$ denotes the likelihood for the $i$-th FRB with the corrected observable DM as \begin{align} {{\rm DM}_{\scriptscriptstyle{\rm E}}}&\equiv {\rm DM_{\rm obs}}-{{\rm DM}_{\rm halo}^{\scriptscriptstyle{\rm MW}}}-{{\rm DM}_{\scriptscriptstyle{\rm ISM}}^{\scriptscriptstyle{\rm MW}}}\notag\\ &={{\rm DM}_{\scriptscriptstyle{\rm IGM}}}+{{\rm DM}_{\rm host}}+{{\rm DM}_{\gamma}}. \tag {15} \end{align} For a certain FRB at the redshift of $z_i$, we can estimate the probability distribution of DM$_{\scriptscriptstyle{\rm E}}$ with the combination of Eqs. (12), (13), and (15), which gives the expression \begin{align} P(&{{\rm DM}_{\scriptscriptstyle{\rm E}}}|z_i) =\int_0^{{{\rm DM}_{\scriptscriptstyle{\rm E}}}-{{\rm DM}_{\gamma}}}P_{\rm host}({\rm {DM}_{\rm host}}|\mu,\,\sigma_{\rm{host}})\notag\\ &\times P_{\scriptscriptstyle{\rm IGM}}({{\rm DM}_{\scriptscriptstyle{\rm E}}}-{{\rm DM}_{\rm host}}-{{\rm DM}_{\gamma}}|z_i){{\rm dDM}_{\rm host}}. \tag {16} \end{align} Note that by Eq. (8), DM$_{\gamma}$ is a constant at a given $z_i$. Denoting $P_X(X)$, $P_Y(Y)$, and $P_Z(Z)$ as the probability density functions of random variables $X$, $Y$, and $Z$, respectively. The equality $P_Z(Z)=\int_{0}^{Z}P_X(X)P_Y(Z-X)dX$ holds for all non-negative independent random variables $X$, $Y$, and $Z$, such that $Z=X+Y$. As DM$_{\scriptscriptstyle{\rm IGM}}$ and DM$_{\rm{host}}$ are non-negative independent random variables, by setting $Z=\rm{DM}_{\rm{E}}-\rm{DM}_{\gamma}$, $X=\rm{DM}_{\rm{host}}$, and $Y=\rm{DM}_{\rm{IGM}}$, it is easy to derive Eq. (16). By maximizing the likelihood function $\mathcal{L}$, we can place constraints on the parameters in our model, including the rest mass of the photon $m_{\gamma}$. Nevertheless, there are two redshift-dependent functions, $H_{\gamma}(z)$ and $H_{\rm e}(z)$ [see Eqs. (6) and (11)], to be determined. In previous works,[18-27] the conventional approach is to adopt the flat $\Lambda$CDM cosmological model and to replace the denominator $H(z)$ with $H_0[\varOmega_m(1+z)^3+\varOmega_{\Lambda}]^{1/2}$. Subsequently, $H_{\gamma}(z)$ and $H_{\rm e}(z)$ are calculated by treating the cosmological parameters $\varOmega_m$ and $\varOmega_{\Lambda}$ as fixed values. However, it is crucial to note that the $\Lambda$CDM model is established within the framework of Einstein's theory of GR, which is derived from one of the fundamental postulates, i.e., the masslessness of photons. Consequently, adopting this model introduces a logical circularity problem in constraining the rest mass of the photon. To circumvent this issue, it is necessary to employ a cosmology-independent method to describe the relation between the Hubble parameter $H(z)$ and redshift $z$. 2.2 Artificial Neural Network. Gaussian process (GP) is a widely used approach for studying cosmology in a model-independent manner. This method assumes that the reconstructed value at a given point follows a Gaussian distribution, and the relationship between function values at different points is characterized by a selected covariance function.[41] In addition, the reconstructed function of $H(z)$ from GP tends to underestimate errors and is significantly influenced by the prior of the Hubble constant $H_0$.[42,43] In contrast, ANN, inspired by biological neural networks, is a mathematical model primarily used to discover intricate relationships between input and output data. The reconstructed function based on ANN makes no assumptions about the observational data and is entirely data-driven without parameterization of the function. Recently, researchers[44] developed a public code called ReFANN https://github.com/Guo-Jian-Wang/refann for reconstructing functions from data using ANN. Their reconstructed $H(z)$ function showed no sensitivity to the setting of $H_0$, indicating that reconstructing $H(z)$ from observational data using ANN may be more reliable than using GP. To make a cosmology-independent determination of $H(z)$, in this work we reconstruct a smooth curve of $H(z)$ from 34 model-independent measurements of $H(z)$ using ReFANN. The first step in using ReFANN to reconstruct a function from data is to find the best parameter configuration for the model used to train the data. Following the approach of Ref. [44], we determine the optimal network model by minimizing the risk:[45] \begin{align} {\rm risk}&=\sum_{i=1}^{N}{\rm Bias}_i^2+\sum_{i=1}^{N}{\rm Variance}_i\notag\\ &=\sum_{i=1}^{N}[H(z_i)-\bar{H}(z_i)]^2+\sum_{i=1}^{N}\sigma^2_{H(z_i)}, \tag {17} \end{align} where $N$ is the number of sample $H(z)$ and $\bar{H}(z_i)$ is the fiducial value of $H(z)$ at redshift $z_i$. For our purposes, we only consider the number of hidden layers and neurons in the hidden layer. Through experimentation, we determined that the optimal model configuration involves a single hidden layer with 4096 neurons. Additionally, to reduce sensitivity to initialization and stabilize the distribution among variables, we employ batch normalization[46] in our network model, with the batch size set to half of the number of $H(z)$ data points. We also use the Exponential Linear Unit[47] as the activation function, with the hyperparameter $\alpha$ set to be 1. Furthermore, we utilize the Adam optimizer[48] to update the network parameters in each iteration of the training process. Subsequently, we train this optimal network model with the observational $H(z)$ data. This training process provides us with a predicted value of $H(z_i)$ and its associated error $\sigma_{H(z_i)}$ at a given redshift $z_i$. This predicted function represents an approximate reconstruction of $H(z)$ based on the ANN model trained with observational data.
Table 1. The latest 34 $H(z)$ measurements obtained with the CC method.
$z$ $H(z)$
(km$\cdot$s$^{-1}\cdot$Mpc$^{-1}$)		 $\sigma$
(km$\cdot$s$^{-1}\cdot$Mpc$^{-1}$)		 Refs.
0.09 69 12  [49]
0.17 83 8  [50]
0.27 77 14
0.4 95 17
0.9 117 23
1.3 168 17
1.43 177 18
1.53 140 14
1.75 202 40
0.48 97 62  [51]
0.88 90 40
0.1791 75 4  [52]
0.1993 75 5
0.3519 83 14
0.5929 104 13
0.6797 92 8
0.7812 105 12
0.8754 125 17
1.037 154 20
0.07 69 19.6  [53]
0.12 68.2 26.2
0.2 72.9 29.6
0.28 88.8 36.6
1.363 160 33.6  [54]
1.965 186.5 50.4
0.3802 83 13.5  [55]
0.4004 77 10.2
0.4247 87.1 11.2
0.4497 92.8 12.9
0.4783 80.9 9
0.47 89 49.6  [56]
0.75 98.8 33.6  [57]
0.8 113.1 25.22  [58]
1.26 135 65  [59]
2.3 Hubble Parameter $H(z)$ Data. The Hubble parameter $H(z)$, which quantifies the expansion rate of the universe as a function of time (or equivalently redshift), has been extensively used for exploring the nature of dark energy and testing modified theories of gravity. $H(z)$ can be measured using two main methods. One approach is based on the detection of the radial baryon acoustic oscillation features.[60-62] However, this method relies on assumptions about the underlying cosmological model, typically the $\Lambda$CDM model. Therefore, in our analysis, we focus exclusively on $H(z)$ measurements obtained using another method, namely the cosmic chronometers (CC) method.[63] This method is model-independent as it relies on minimal assumptions, primarily the use of a Friedmann–Lemaître–Robertson–Walker metric. In this method, by combining the definition of $H(z)\equiv \dot{a}/a$ and the relation between the scale factor $a(t)$ and redshift, $1+z=1/a(t)$, $H(z)$ can be rewritten as \begin{equation} H(z)=-\frac{1}{1+z}\frac{\mathrm dz}{\mathrm dt}. \tag {18} \end{equation} This means that one can achieve a model-independent determination of $H(z)$ by calculating the differential age of the universe using passively evolving galaxies at various redshifts. We compile the latest 34 CC $H(z)$ measurements in Table 1, covering the redshift range of $0.07 < z < 1.965$.
cpl-41-5-059501-fig1.png
Fig. 1. Reconstruction of the Hubble parameter $H(z)$ from 34 CC $H(z)$ measurements using ANN. The shaded area corresponds to the 1$\sigma$ confidence region of the reconstruction. The blue dots with error bars depict the observational $H(z)$ data.
cpl-41-5-059501-fig2.png
Fig. 2. Redshift dependence of the reconstructed $H_{\gamma}(z)$ and $H_{\rm e}(z)$ functions (solid lines). The corresponding theoretical curves (dashed lines) derived from the flat $\Lambda$CDM model with $\varOmega_m=0.315$ and $H_0=67.4$ km$\cdot$s$^{-1}\cdot$Mpc$^{-1}$ are also shown.
Having obtained the dataset of $H(z)$, we adopt ANN to reconstruct the $H(z)$ function within the redshift range of $0 < z < 2$, and the results are shown in Fig. 1. The black line represents the best fit, and the blue shaded area indicates the 1$\sigma$ confidence region of the reconstructed function. The redshift-dependent $H_{\gamma}(z)$ and $H_{\rm e}(z)$ functions can then be derived by integrating the reconstructed $H(z)$ function with respect to redshift. As shown in Fig. 2, the reconstructed $H_{\gamma}(z)$ and $H_{\rm e}(z)$ functions exhibit significant differences in behavior as a function of redshift $z$. Such distinctions are crucial for breaking parameter degeneracy and enhancing sensitivity in testing the photon mass when a few redshift measurements of FRBs are available.[19-21,23,64] For comparison, we also plot $H_{\gamma}(z)$ and $H_{\rm e}(z)$ functions within the framework of the flat $\Lambda$CDM model using Planck 2018 parameters ($\varOmega_m=0.315$ and $H_0=67.4$ km$\cdot$s$^{-1}\cdot$Mpc$^{-1}$).[65] The corresponding theoretical curves are presented in Fig. 2 as dashed lines. It is evident that the reconstructions of $H_{\gamma}(z)$ and $H_{\rm e}(z)$ are roughly consistent with the predictions of the standard $\Lambda$CDM model, implying that the ANN method can offer a reliable reconstructed function from the observational data.
Table 2. Properties of 32 localized FRBs.
Name $z$ DM$_{\rm obs}$
(pc$\cdot$cm$^{-3})$		 DM$_{\scriptscriptstyle{\rm ISM}}^{\scriptscriptstyle{\rm MW}}$
(pc$\cdot$cm$^{-3})$		 Refs.
FRB 20121102 0.19273 557 188.0  [66]
FRB 20180301 0.3304 536 152.0  [67]
FRB 20180916 0.0337 348.76 200.0  [68]
FRB 20180924 0.3214 361.42 40.5  [69]
FRB 20181112 0.4755 589.27 102.0  [70]
FRB 20190102 0.291 363.6 57.3  [71]
FRB 20190523 0.66 760.8 37.0  [72]
FRB 20190608 0.1178 338.7 37.2  [73]
FRB 20190611 0.378 321.4 57.83  [74]
FRB 20190614 0.6 959.2 83.5  [75]
FRB 20190711 0.522 593.1 56.4  [74]
FRB 20190714 0.2365 504 38.0  [74]
FRB 20191001 0.234 506.92 44.7  [74]
FRB 20191228 0.2432 297.5 33.0  [67]
FRB 20200430 0.16 380.1 27.0  [74]
FRB 20200906 0.3688 577.8 36.0  [67]
FRB 20201124 0.098 413.52 123.2  [76]
FRB 20210117 0.2145 730 34.4  [77]
FRB 20210320 0.27970 384.8 42  [77]
FRB 20210807 0.12927 251.9 121.2  [77]
FRB 20211127 0.0469 234.83 42.5  [77]
FRB 20211212 0.0715 206 27.1  [77]
FRB 20220207C 0.043040 262.38 79.3  [78]
FRB 20220307B 0.28123 499.27 135.7  [78]
FRB 20220310F 0.477958 462.24 45.4  [78]
FRB 20220418A 0.622000 623.25 37.6  [78]
FRB 20220506D 0.30039 396.97 89.1  [78]
FRB 20220509G 0.089400 269.53 55.2  [78]
FRB 20220610A 1.016 1457.624 31  [79]
FRB 20220825A 0.241397 651.24 79.7  [78]
FRB 20220920A 0.158239 314.99 40.3  [78]
FRB 20221012A 0.284669 441.08 54.4  [78]
2.4 Fast Radio Burst Data. Expanding upon the sample of 23 localized FRBs utilized in Ref. [27], we incorporate 9 new localized FRBs recently detected by the 110-antenna Deep Synoptic Array.[78] Note that two new FRBs with redshift measurements, FRB 20220319D and FRB 20220914A, are not included in our sample. FRB 20220319D exhibits an observed DM$_{\rm obs}$ value of 110.98 pc$\cdot$cm$^{-3}$, while the DM$_{\scriptscriptstyle{\rm ISM}}^{\scriptscriptstyle{\rm MW}}$ value calculated using the NE2001 model stands at 133.3 pc$\cdot$cm$^{-3}$, leading to its exclusion from our analysis. Furthermore, FRB 20220914A is omitted due to its notable deviation from the expected DM$_{\scriptscriptstyle{\rm IGM}}$, i.e., $z$ relation outlined in Eq. (10). The total 32 localized FRBs, including their respective redshifts, DM$_{\rm obs}$, and DM$_{\scriptscriptstyle{\rm ISM}}^{\scriptscriptstyle{\rm MW}}$, are provided in Table 2. 3. Results. We explore the posterior probability distributions of the free parameters by maximizing the joint likelihood function $\mathcal{L}$ [Eq. (14)] using emcee, an affine-invariant Markov Chain Monte Carlo (MCMC) ensemble sampler implemented in Python.[80] The free parameters of our model now include the baryon density parameter $\varOmega_{\rm b} h_{0}^2$ [where $h_{0} \equiv H_0/$ (100 km$\cdot$s$^{-1}\cdot$Mpc$^{-1}$)], the DM contribution from the Milky Way's halo DM$_{\rm halo}^{\scriptscriptstyle{\rm MW}}$, the parameters related to the probability distributions for DM$_{\scriptscriptstyle{\rm IGM}}$ and DM$_{\rm host}$ (i.e., $F$, $\mu$, and $\sigma_{\rm{host}}$), and the photon mass $m_{\gamma}$. In our MCMC analysis, we set uniform priors on $\varOmega_{\rm b}h_{0}^2\in[0.01,\,\;1]$, $F\in[0.01,\,\;0.5]$, $e^{\mu}\in[20,\,\;200]$  pc$\cdot$cm$^{-3}$, $\sigma_{\rm{host}}\in[0.2,\,\;2]$, and $m_{\gamma}\in[0,\,\;10^{-49}]$ kg. For DM$_{\rm halo}^{\scriptscriptstyle{\rm MW}}$, we set a Gaussian prior, DM$_{\rm halo}^{\scriptscriptstyle{\rm MW}}=65\pm 15$ pc$\cdot$cm$^{-3}$, within a wide $3\sigma$ range of $[20,\,\;110]$ pc$\cdot$cm$^{-3}$. To incorporate the error of the ANN reconstruction into our analysis, at each MCMC step, we sample the Hubble rate $\tilde{H}(z)$ function according to the Gaussian distribution $\tilde{H}(z)=\mathcal{N}(H(z),\,\sigma_{H(z)})$ with mean $H(z)$ and standard deviation $\sigma_{H(z)}$. Here $H(z)$ is the best-fit function reconstructed by ANN and $\sigma_{H(z)}$ is the $1\sigma$ error of the reconstructed $H(z)$ function. The $1\sigma$ constraint results for these six parameters are summarized in Table 3. For the photon mass $m_{\gamma}$, both the $1\sigma$ and $2\sigma$ upper limits are presented. Posterior distributions and $1-2\sigma$ confidence regions for these parameters are displayed in Fig. 3. Notably, the constraints on $m_{\gamma}$ are determined as \begin{align} \tag {19} &m_{\gamma} \le 3.5 \times 10^{-51} \, \rm{kg} \simeq 2.0 \times 10^{-15} \, \rm{eV/c^2},\\ \tag {20} &m_{\gamma} \le 6.5 \times 10^{-51} \, \rm{kg} \simeq 3.6 \times 10^{-15} \, \rm{eV/c^2}, \end{align} at the $1\sigma$ and $2\sigma$ confidence levels, respectively. Meanwhile, we find that the baryon density parameter is optimized to be $\varOmega_{\rm b} h_{0}^2=0.030_{-0.007}^{+0.006}$, which is compatible with the value inferred from Planck 2018 ($\varOmega_{\rm b} h_0^2 = 0.0224 \pm 0.0001$) at the $1.2\sigma$ confidence level.[65] To investigate the impact of the prior assumption of DM$_{\rm halo}^{\scriptscriptstyle{\rm MW}}$ on our results, we also perform a parallel comparative analysis of the FRB data using a flat prior on DM$_{\rm halo}^{\scriptscriptstyle{\rm MW}}\in[20,\,\;110]$ pc$\cdot$cm$^{-3}$. The corresponding resulting constraints on all parameters are also reported in Table 3. Comparing these inferred parameters with those obtained from the Gaussian prior (see the first row in Table 3), it is clear that, except for the best-fit value of DM$_{\rm halo}^{\scriptscriptstyle{\rm MW}}$, the prior assumption of DM$_{\rm halo}^{\scriptscriptstyle{\rm MW}}$ only has a minimal influence on our results.
Table 3. The $1\sigma$ constraints for all parameters with different priors on DM$_{\rm halo}^{\scriptscriptstyle{\rm MW}}$. Both the $1\sigma$ and $2\sigma$ upper limits for the photon mass $m_{\gamma}$ are provided.
Priors of DM$_{\rm halo}^{\scriptscriptstyle{\rm MW}}$ ${\varOmega_{\rm b}}h_{0}^2$ $F$ DM$_{\rm halo}^{\scriptscriptstyle{\rm MW}}$
(pc$\cdot$cm$^{-3}$)		 $\mu$ $\sigma_{\rm{host}}$ $m_{\gamma}$
(10$^{-51}$ kg)
Gaussian prior $0.030_{-0.007}^{+0.006}$ $\geq 0.449$ $39_{-10}^{+10}$ $4.02_{-0.29}^{+0.46}$ $1.01_{-0.37}^{+0.22}$ $\leq 3.5\, (\leq 6.5)$
Flat prior $0.031_{-0.006}^{+0.006}$ $\geq 0.451$ $\leq 32$ $4.07_{-0.30}^{+0.48}$ $0.99_{-0.37}^{+0.21}$ $\leq 3.8\, (\leq 7.2)$
cpl-41-5-059501-fig3.png
Fig. 3. Plots of 1D marginalized probability distributions and 2D $1-2\sigma$ confidence contours for the parameters $\varOmega_{\rm b} h_{0}^2$, $F$, DM$_{\rm halo}^{\scriptscriptstyle{\rm MW}}$, $\mu$, $\sigma_{\rm{host}}$, and $m_{\gamma}$, constrained by 32 localized FRBs.
4. Discussion and Conclusions. As a matter of fact, it is impossible to prove experimentally that the rest mass of a photon is strictly zero. According to the uncertainty principle of quantum mechanics, the ultimately measurable order of magnitude for the photon mass is $m_{\gamma}\approx \hbar/(\Delta tc^2)\approx 10^{-69}$ kg, where $\hbar$ is the reduced Planck constant and $\Delta t\approx 10^{10}$ y is the age of our universe. However, due to the immense importance of Maxwell's theory and Einstein's theory of relativity, it is still intellectually stimulating and scientifically significant to approach this ultimate limit by various experiments. FRBs provide the current best celestial laboratory to test the photon mass $m_{\gamma}$ via the dispersion method. Constraining $m_{\gamma}$ with cosmological FRBs, however, one has to know the cosmic expansion rate $H(z)$. In all previous studies, the required $H(z)$ information is estimated within the standard $\Lambda$CDM cosmological model. Such $H(z)$ estimations would involve a circularity problem in constraining $m_{\gamma}$, since $\Lambda$CDM itself is built on the framework of GR and GR embraces the postulate of the constancy of light speed. In this work, aiming to overcome the circularity problem, we have employed an ANN technology to reconstruct a cosmology-independent $H(z)$ function from the discrete CC $H(z)$ data. By combining the DM–$z$ measurements of 32 localized FRBs with the reconstructed $H(z)$ function from 34 CC $H(z)$ data, we have placed the first cosmology-independent photon mass limit. Our results show that the $1\sigma$ and $2\sigma$ confidence-level upper limits on the photon mass are $m_{\gamma} \le3.5\times 10^{-51}$ kg (or equivalently $m_{\gamma} \le 2.0 \times 10^{-15}$ eV/c$^2$) and $m_{\gamma} \le 6.5 \times 10^{-51}$ kg (or equivalently $m_{\gamma} \le 3.6 \times 10^{-15}$ eV/c$^2$), respectively. Previously, under the assumption of fiducial $\Lambda$CDM cosmology, Ref. [24] obtained an upper limit of $m_{\gamma} \le 3.1 \times 10^{-51}$ kg at the $1\sigma$ confidence level by using a catalog of 129 FRBs (most of them without redshift measurement, and the observed DM$_{\rm obs}$ values were used to estimate the pseudo redshifts). The authors of Ref. [26] obtained $m_{\gamma} \le 4.8 \times 10^{-51}$ kg ($1\sigma$) by analyzing a sample of 17 localized FRBs in the flat $\Lambda$CDM model. Ref. [27] obtained $m_{\gamma}\le3.8\times 10^{-51}$ kg ($1\sigma$) for flat $\Lambda$CDM using a sample of 23 localized FRBs. Despite not assuming a specific cosmological model, the precision of our constraint from 32 localized FRBs is comparable to these previous results. Most importantly, this highlights the validity of our approach and suggests that as the number of CC $H(z)$ measurements increases, we can expect even more reliable model-independent tests of the photon mass. Acknowledgements. We are grateful to the anonymous referees for their helpful comments. This work was supported by the National SKA Program of China (2022SKA0130100), the National Natural Science Foundation of China (Grant Nos. 12373053, 12321003, and 12041306), the Key Research Program of Frontier Sciences of Chinese Academy of Sciences (Grant No. ZDBS-LY-7014), the International Partnership Program of Chinese Academy of Sciences for Grand Challenges (Grant No. 114332KYSB20210018), the CAS Project for Young Scientists in Basic Research (Grant No. YSBR-063), the CAS Organizational Scientific Research Platform for National Major Scientific and Technological Infrastructure: Cosmic Transients with FAST, and the Natural Science Foundation of Jiangsu Province (Grant No. BK20221562). References Rayonnement noir et quanta de lumièreSur la théorie ondulatoire des électrons positifs et négatifsMassive photon and dark energyCosmology and the massive photon frequency shift in the Standard-Model ExtensionNew Experimental Test of Coulomb's Law: A Laboratory Upper Limit on the Photon Rest MassExperimental Limits on the Photon Mass and Cosmic Magnetic Vector PotentialNew Experimental Limit on the Photon Rest Mass with a Rotating Torsion BalanceLimits on the Photon MassAstrophysical upper limits on the photon rest massThe role of finite photon mass in magnetohydrodynamics of space plasmasUsing plasma physics to weigh the photonRe-estimatation of the Upper Limit on the Photon Mass with the Solar Wind MethodSolar wind test of the de Broglie-Proca massive photon with Cluster multi-spacecraft dataLimit on the Photon Mass Deduced from Pioneer-10 Observations of Jupiter's Magnetic FieldTight Constraint on Photon Mass from Pulsar SpindownVariations of the Speed of Light with Frequency and Implied Photon MassThe mass of the photonCONSTRAINTS ON THE PHOTON MASS WITH FAST RADIO BURSTSPhoton mass limits from fast radio burstsFRB 121102 casts new light on the photon massBayesian framework to constrain the photon mass with a catalog of fast radio burstsLimits on the Weak Equivalence Principle and Photon Mass with FRB 121102 SubpulsesCombined limit on the photon mass with nine localized fast radio burstsBounding the photon mass with cosmological propagation of fast radio burstsBounding the photon mass with the dedispersed pulses of the Crab pulsar and FRB 180916BRevised constraints on the photon mass from well-localized fast radio burstsRevisiting constraints on the photon rest mass with cosmological fast radio burstsBounding the Photon Mass with Ultrawide Bandwidth Pulsar Timing Data and Dedispersed Pulses of Fast Radio BurstsThe physics of fast radio burstsTesting fundamental physics with astrophysical transientsNE2001.I. A New Model for the Galactic Distribution of Free Electrons and its FluctuationsProbing Galactic Halos with Fast Radio BurstsExploring the dispersion measure of the Milky Way haloAn 8 per cent determination of the Hubble constant from localized fast radio burstsCOSMOLOGICAL IMPLICATIONS OF FAST RADIO BURST/GAMMA-RAY BURST ASSOCIATIONSThe Cosmic Baryon BudgetReionization of the Inhomogeneous UniverseA census of baryons in the Universe from localized fast radio burstsLOCATING THE “MISSING” BARYONS WITH EXTRAGALACTIC DISPERSION MEASURE ESTIMATESReconstruction of dark energy and expansion dynamics using Gaussian processesAn Improved Method to Measure the Cosmic CurvatureModel-independent Constraints on Cosmic Curvature and OpacityReconstructing Functions and Estimating Parameters with Artificial Neural Networks: A Test with a Hubble Parameter and SNe IaNon-Parametric Inference in AstrophysicsBatch Normalization: Accelerating Deep Network Training by Reducing Internal Covariate ShiftFast and Accurate Deep Network Learning by Exponential Linear Units (ELUs)Adam: A Method for Stochastic OptimizationConstraints on the Equation of State of Dark Energy and the Hubble Constant from Stellar Ages and the Cosmic Microwave BackgroundConstraints on the redshift dependence of the dark energy potentialCosmic chronometers: constraining the equation of state of dark energy. I: H(z) measurementsImproved constraints on the expansion rate of the Universe up to z ∼ 1.1 from the spectroscopic evolution of cosmic chronometersFour new observational H ( z ) data from luminous red galaxies in the Sloan Digital Sky Survey data release sevenRaising the bar: new constraints on the Hubble parameter with cosmic chronometers at z ∼ 2A 6% measurement of the Hubble parameter at z∼0.45: direct evidence of the epoch of cosmic re-accelerationAge-dating luminous red galaxies observed with the Southern African Large TelescopeToward a Better Understanding of Cosmic Chronometers: A New Measurement of H(z) at z ∼ 0.7New Observational H(z) Data from Full-spectrum Fitting of Cosmic Chronometers in the LEGA-C SurveyA new measurement of the expansion history of the Universe at z = 1.26 with cosmic chronometers in VANDELSClustering of luminous red galaxies - IV. Baryon acoustic peak in the line-of-sight direction and a direct measurement of H ( z )The WiggleZ Dark Energy Survey: joint measurements of the expansion and growth history at z < 1The clustering of galaxies in the SDSS-III DR9 Baryon Oscillation Spectroscopic Survey: testing deviations from Λ and general relativity using anisotropic clustering of galaxiesConstraining Cosmological Parameters Based on Relative Galaxy AgesDispersion by pulsars, magnetars, fast radio bursts and massive electromagnetism at very low radio frequenciesPlanck 2018 resultsA direct localization of a fast radio burst and its hostCharacterizing the Fast Radio Burst Host Galaxy Population and its Connection to Transients in the Local and Extragalactic UniverseA repeating fast radio burst source localized to a nearby spiral galaxyA single fast radio burst localized to a massive galaxy at cosmological distanceThe low density and magnetization of a massive galaxy halo exposed by a fast radio burstThe Host Galaxies and Progenitors of Fast Radio Bursts Localized with the Australian Square Kilometre Array PathfinderA fast radio burst localized to a massive galaxyDissecting the Local Environment of FRB 190608 in the Spiral Arm of its Host GalaxyHost Galaxy Properties and Offset Distributions of Fast Radio Bursts: Implications for Their ProgenitorsA Distant Fast Radio Burst Associated with Its Host Galaxy by the Very Large ArrayThe host galaxy and persistent radio counterpart of FRB 20201124AA measurement of Hubble’s Constant using Fast Radio BurstsDeep Synoptic Array Science: First FRB and Host Galaxy CatalogA luminous fast radio burst that probes the Universe at redshift 1emcee : The MCMC Hammer
[1] de Broglie L 1922 J. Phys. Radium 3 422
[2] Proca A 1936 J. Phys. Radium 7 347
[3] Kouwn S, Oh P, and Park C G 2016 Phys. Rev. D 93 083012
[4] Spallicci A D A M, Helayël-Neto J A, López-Corredoira M, and Capozziello S 2021 Eur. Phys. J. C 81 4
[5] Williams E R, Faller J E, and Hill H A 1971 Phys. Rev. Lett. 26 721
[6] Lakes R 1998 Phys. Rev. Lett. 80 1826
[7] Luo J, Tu L C, Hu Z K, and Luan E J 2003 Phys. Rev. Lett. 90 081801
[8] Lowenthal D D 1973 Phys. Rev. D 8 2349
[9] Chibisov G V 1976 Sov. Phys. Usp. 19 624
[10] Ryutov D D 1997 Plasma Phys. Control. Fusion 39 A73
[11] Ryutov D D 2007 Plasma Phys. Control. Fusion 49 B429
[12] Liu L X and Shao C G 2012 Chin. Phys. Lett. 29 111401
[13] Retinò A, Spallicci A D A M, and Vaivads A 2016 Astropart. Phys. 82 49
[14] Jr Davis L, Goldhaber A S, and Nieto M M 1975 Phys. Rev. Lett. 35 1402
[15] Yang Y P and Zhang B 2017 Astrophys. J. 842 23
[16] Tu L C, Ye H L, and Luo J 2005 Chin. Phys. Lett. 22 3057
[17] Tu L C, Luo J, and Gillies G T 2005 Rep. Prog. Phys. 68 77
[18] Wu X F et al. 2016 Astrophys. J. Lett. 822 L15
[19] Bonetti L et al. 2016 Phys. Lett. B 757 548
[20] Bonetti L et al. 2017 Phys. Lett. B 768 326
[21] Shao L and Zhang B 2017 Phys. Rev. D 95 123010
[22] Xing N et al. 2019 Astrophys. J. Lett. 882 L13
[23] Wei J J and Wu X F 2020 Res. Astron. Astrophys. 20 206
[24] Wang H, Miao X, and Shao L 2021 Phys. Lett. B 820 136596
[25] Chang C M, Wei J J, Zhang S B, and Wu X F 2023 J. Cosmol. Astropart. Phys. 2023(01) 010
[26] Lin H N, Tang L, and Zou R 2023 Mon. Not. R. Astron. Soc. 520 1324
[27] Wang B, Wei J J, Wu X F, and López-Corredoira M 2023 J. Cosmol. Astropart. Phys. 2023(09) 025
[28] Wang Y B, Zhou X, Kurban A, and Wang F Y 2024 Astrophys. J. 965 38
[29]Lorimer D R and Kramer M 2012 Handbook of Pulsar Astronomy
[30] Zhang B 2023 Rev. Mod. Phys. 95 035005
[31] Wei J J and Wu X F 2021 Front. Phys. 16 44300
[32] Cordes J M and Lazio T J W 2002 arXiv:astro-ph/0207156
[33] Prochaska J X and Zheng Y 2019 Mon. Not. R. Astron. Soc. 485 648
[34] Keating L C and Pen U L 2020 Mon. Not. R. Astron. Soc. 496 L106
[35] Wu Q, Zhang G Q, and Wang F Y 2022 Mon. Not. R. Astron. Soc. 515 L1
[36] Deng W and Zhang B 2014 Astrophys. J. Lett. 783 L35
[37] Fukugita M, Hogan C J, and Peebles P J E 1998 Astrophys. J. 503 518
[38] Miralda-Escudé J, Haehnelt M, and Rees M J 2000 Astrophys. J. 530 1
[39] Macquart J P et al. 2020 Nature 581 391
[40] McQuinn M 2014 Astrophys. J. Lett. 780 L33
[41] Seikel M, Clarkson C, and Smith M 2012 J. Cosmol. Astropart. Phys. 2012(06) 036
[42] Wei J J and Wu X F 2017 Astrophys. J. 838 160
[43] Wang G J et al. 2017 Astrophys. J. 847 45
[44] Wang G J, Ma X J, Li S Y, and Xia J Q 2020 Astrophys. J. Suppl. Ser. 246 13
[45] Wasserman L et al. 2001 arXiv:astro-ph/0112050
[46] Ioffe S and Szegedy C 2015 arXiv:1502.03167 [cs.LG]
[47] Clevert D A, Unterthiner T, and Hochreiter S 2015 arXiv:1511.07289 [cs.LG]
[48] Kingma D P and Ba J 2014 arXiv:1412.6980 [cs.LG]
[49] Jimenez R, Verde L, Treu T, and Stern D 2003 Astrophys. J. 593 622
[50] Simon J, Verde L, and Jimenez R 2005 Phys. Rev. D 71 123001
[51] Stern D et al. 2010 J. Cosmol. Astropart. Phys. 2010(02) 008
[52] Moresco M et al. 2012 J. Cosmol. Astropart. Phys. 2012(08) 006
[53] Zhang C et al. 2014 Res. Astron. Astrophys. 14 1221
[54] Moresco M 2015 Mon. Not. R. Astron. Soc. 450 L16
[55] Moresco M et al. 2016 J. Cosmol. Astropart. Phys. 2016(05) 014
[56] Ratsimbazafy A L et al. 2017 Mon. Not. R. Astron. Soc. 467 3239
[57] Borghi N, Moresco M, and Cimatti A 2022 Astrophys. J. Lett. 928 L4
[58] Jiao K, Borghi N, Moresco M, and Zhang T J 2023 Astrophys. J. Suppl. Ser. 265 48
[59] Tomasetti E et al. 2023 Astron. Astrophys. 679 A96
[60] Gaztañaga E, Cabré A, and Hui L 2009 Mon. Not. R. Astron. Soc. 399 1663
[61] Blake C et al. 2012 Mon. Not. R. Astron. Soc. 425 405
[62] Samushia L et al. 2013 Mon. Not. R. Astron. Soc. 429 1514
[63] Jimenez R and Loeb A 2002 Astrophys. J. 573 37
[64] Bentum M J, Bonetti L, and Spallicci A D A M 2017 Adv. Space Res. 59 736
[65] Planck Collaboration, Aghanim N et al. 2020 Astron. Astrophys. 641 A6
[66] Chatterjee S et al. 2017 Nature 541 58
[67] Bhandari S et al. 2022 Astron. J. 163 69
[68] Marcote B et al. 2020 Nature 577 190
[69] Bannister K W et al. 2019 Science 365 565
[70] Prochaska J X et al. 2019 Science 366 231
[71] Bhandari S et al. 2020 Astrophys. J. Lett. 895 L37
[72] Ravi V et al. 2019 Nature 572 352
[73] Chittidi J S et al. 2021 Astrophys. J. 922 173
[74] Heintz K E et al. 2020 Astrophys. J. 903 152
[75] Law C J et al. 2020 Astrophys. J. 899 161
[76] Ravi V et al. 2022 Mon. Not. R. Astron. Soc. 513 982
[77] James C W et al. 2022 Mon. Not. R. Astron. Soc. 516 4862
[78] Law C J et al. 2023 arXiv:2307.03344 [astro-ph.HE]
[79] Ryder S D et al. 2023 Science 382 294
[80] Foreman-Mackey D, Hogg D W, Lang D, and Goodman J 2013 Publ. Astron. Soc. Pac. 125 306