Constraining the Hubble constant and its lower limit from the proper motion of extragalactic radio jets

Approximately a century ago, the Universe was found to be expanding (Hubble, 1929). The Hubble parameter, $H(z)$, describes the expansion rate of the Universe. Its value in the current era ($z=0$), known as the Hubble constant, is denoted by $H_{0}$. Numerous measurements of the Hubble constant observed from different methods have been proposed to measure the values during the past few decades. For instance, Cosmic Microwave Background (CMB), remnant of the Big Bang in the early Universe, which is also a standard ruler as a cosmological probe. One of the latest measurements from Planck Collaboration et al. (2020) suggested $H_{0}=67.27\pm 0.60\,{\rm km\,s^{-1}\,Mpc^{-1}}$. Additionally, the method of the distance ladder, such as Cepheid variables in the late Universe, can be used to constrain the Hubble constant directly. The Supernova, H0, for the Equation of State (SH0ES; Riess et al., 2021), which adopted the method of distance ladder. They utilise Cepheids to calibrate 42 Type Ia supernovae in the same host galaxies, obtained a value of $H_{0}=73.04\pm 1.04\,{\rm km\,s^{-1}\,Mpc^{-1}}$. Other methods such as using baryon acoustic oscillations (e.g., Eisenstein et al., 2005; Cuceu et al., 2019), big bang nucleosynthesis (e.g., Cuceu et al., 2019; Seto & Toda, 2021), strong gravitational lensing of quasars (e.g., Wong et al., 2020), water masers (e.g., Herrnstein et al., 1999; Humphreys et al., 2013), gravitationally lensed supernova (e.g., Refsdal, 1964; Vega-Ferrero et al., 2018) were also utilised to measure the Hubble constant.

However, as we obtain more precise measurements, the results conducted from the early Universe further reveal an inconsistency with the Hubble constant from the late Universe. This discrepancy is known as "Hubble tension" and indicates a tension beyond $4.4\sigma$ (e.g., Verde et al., 2019). The reason for the tension is still under debate. One possible theory is that, there may be hints of new physics beyond the $\Lambda$CDM model. For instance, there are popular theories such as the early dark energy (e.g., Poulin et al., 2019), acoustic dark energy (e.g., Yin, 2020), or importing sterile neutrino (e.g., Adhikari et al., 2017) trying to tackle this issue. In terms of observation, aside from observing standard candles and standard rulers, several new methods have also been proposed to address this discrepancy for the next generation. For example, since gravitational waves (GWs) act as standard sirens, an independent way to infer luminosity distance from their amplitudes (Schutz, 1986), some studies suggested constraining cosmological parameters using GWs in the future (e.g., Abbott et al., 2017; Chen et al., 2018). Besides, fast radio bursts (FRBs), radio pulses with millisecond time scale, can also be used to constrain cosmology thanks to its unique observable quantities, dispersion measure (e.g. Li et al., 2018; Wu et al., 2021) and duration-luminosity relation (e.g., Hashimoto et al., 2019, Hsiao et al. in prep).

Relativistic jets are one of the energetic phenomena which produce piercing matter from regions near compact objects such as black holes and neutron stars. These robust plasma flows are expected to be along the axis of rotation of the host object. Although the relation between the accretion disc and the relativistic jets is still in debate, it is believed that the jets are driven by the tangled magnetic field (e.g., Blandford & Znajek, 1977; Hawley & Balbus, 2002; McKinney & Gammie, 2004). Previous studies revealed that it is possible to extract the information of the Hubble constant from the proper motion of jets (e.g., Lynden-Bell, 1977; Taylor & Vermeulen, 1997; Qin, 1999; Lu & Qin, 2021) and with GWs (e.g., Hotokezaka et al., 2019). Therefore, the observation of jets have potential to help us relieve the Hubble tension.

In this paper, we adopt the method proposed from Qin (1999) (see also Taylor & Vermeulen, 1997) to calculate the lower limit of the Hubble constant through the Friedman-Lemaître-Robertson-Walker (FLRW) cosmology and geometrical proper motion of extragalactic radio jets. With only the redshifts, and the receding and approaching angular velocities, the lower limit of the Hubble constant can be determined. The first measurement originated from Taylor & Vermeulen (1997), suggesting $H_{\rm 0,min}=37\,{\rm km\,s^{-1}\,Mpc^{-1}}$. Qin (1999) adopted a similar method and suggested $27.08\,{\rm km\,s^{-1}\,Mpc^{-1}}<H_{\rm 0,min}<53.15\,{\rm km\,s^{-1}\,Mpc^{-1}}$. Lu & Qin (2021) made use of an extragalactic radio source, NGC 1052, located at $z=0.005$. Their results suggested the lower limit of the Hubble constant is $H_{\rm 0,min}=51.5\pm 2.3\,{\rm km\,s^{-1}\,Mpc^{-1}}$, which successfully follows the method of Taylor & Vermeulen (1997) in constraining the Hubble constant. Apparently, the latest constraint from proper motion is not enough to relieve the Hubble tension. Therefore, in this paper, we investigate how well $H_{0}$ can be constrained from $H_{\rm 0,min}$ distributions by simulating samples of observed speeds in nearby bi-directional AGN jets. Furthermore, we perform Kolmogorov-Smirnov tests (K-S test; Massey, 1951) to compare the cumulative distribution function (CDF) of the mock data to theoretical distribution. K-S test aims to compare if two distributions are drawn from the same distributions, in order to obtain further information of $H_{0}$ and power-law index $k$ based on the cumulative distribution function of $H_{\rm 0,min}$.

The structure of this paper is as follows. We describe the theoretical framework and our method in Sec. 2. In Sec. 3, we present the results and, distributions of $H_{\rm 0,min}$ from our simulation according to two different cosmology, and further constrain $H_{0}$ and $k$ from K-S tests. The conclusion of this study is summarised in Sec. 4. We assume concordance cosmology with ($\Omega_{m}$, $\Omega_{\Lambda}$, $h$) = ($0.3$, $0.7$, $0.7$), unless otherwise mentioned.

According to FLRW cosmology, the Hubble law for a nearby extragalactic source with $z\ll 1$ can be written as follows:

where $D_{L}$ is the luminosity distance and $c$ is the speed of light.

In terms of the geometry, considering a bi-symmetric relativistic jet, the proper motion of the receding and approaching jets can be illustrated as (e.g., Rees, 1966; Behr et al., 1976; Blandford et al., 1977; Blandford & Königl, 1979; Mirabel & Rodríguez, 1994):

where $\mu_{r}$ and $\mu_{a}$ are the proper motions of receding and approaching jets, respectively. $\beta$ is the ratio between the jet velocity and the speed of light ($v/c$) while $\theta$ is the angle between the velocity vector of the approaching jet and the line of sight. Eqs. 1 and 2 yield (see also Qin, 1999; Lu & Qin, 2021):

Lastly, since velocity must be smaller than the speed of light ($\beta<1$), we obtain a relation among the lower limit of the Hubble constant, proper motion, and redshift:

We note that there are assumptions that the jets are straight and oppositely directed, and have identical bulk flow speed. If the velocity is higher (especially close to the speed of light), the $H_{\rm 0,min}$ will be close to the true $H_{0}$ at any fixed $\theta$. Contours of how $\theta$ and $\beta$ will change the fraction of $H_{\rm 0,min}$/$H_{0}$ is shown in Fig. 1. In principle, if there is a fraction of the speed of jets that is close to the speed of light, we expect to obtain a $H_{\rm 0,min}$ which is close to $H_{0}$ from those samples. Therefore, this method provides a model-free and distance-free way to constrain the Hubble constant only from geometry and FLWR cosmology, and may be able to alleviate $H_{0}$ tension.

We note that due to the relativistic beaming effect (or the Doppler boosting), the receding jets are not easily observed (e.g., Wilkinson et al., 1977; Bridle & Perley, 1984; Sparks et al., 1992; Laing & Bridle, 2002). For an ideal relativistic jet, the observed luminosity can be expressed as (e.g., Blandford & Königl, 1979; Cohen et al., 2007)

where $L_{\rm obs}$ is the observed luminosity and $L_{\rm int}$ is the intrinsic luminosity. $\delta$ is the Doppler factor ($\delta=\gamma^{-1}(1-\beta{\rm cos}(\theta))^{-1}$), and $\gamma$ is the Lorentz factor ($\gamma=1/\sqrt{1-\beta^{2}}$). $\alpha$ is the spectral index. $p$ is the Doppler boost exponent (for more detail, see Cohen et al., 2007). Contours of Dopper factor $\delta$ as a function of $\gamma$ and $\theta$ are shown in Fig. 2. For example, in the catalogue of Monitoring Of Jets in Active galactic nuclei with VLBA Experiments (MOJAVE) program (Lister & Homan, 2005; Homan et al., 2021), the two-sided jets account for $\sim 18\%$ of total sources. In this research, we simply assume the samples are those which have two-sided jets.

We propose a new method to test how well we can recover $H_{0}$ and the power-law index, $k$ (see the explanation in the next paragraph) from the distribution of $H_{\rm 0,min}$ gave a fixed sample size and a typical measurement error in the future. In this process, in addition to the high Lorentz factor jets, we also utilise lower Lorentz factor jets, in fact, we use the whole distribution to estimate $H_{0}$. In this way, the usable sample size becomes much larger. By conducting the K-S test between mock observational data (fixed $H_{0}=70\,{\rm km\,s^{-1}\,Mpc^{-1}}$ and $k=-1.5$) and the theoretical CDFs (with different prior $H_{0}$ and prior $k$), the p-value, a similarity of two profiles, can be represented as probabilities of $H_{0}$ and $k$. We simulate a large number of jets, i.e., $10^{5}$, with $H_{\rm 0,min}$, as theoretical distributions for each point between $65\,{\rm km\,s^{-1}\,Mpc^{-1}}<H_{0}<75\,{\rm km\,s^{-1}\,Mpc^{-1}}$ and $-1<k<-2$. Fig. 3 shows how the true $H_{0}$ and $k$ changes the features of CDFs. Figs. 4 and 5 show how our idea works. 1% error on the measurement of proper motions and 100 sources are assumed for the mock observational data here. Since the feature of a CDF of $H_{\rm 0,min}$ changes when $H_{0}$ changes, p-values from K-S tests will be small when $H_{0}$ of the observational data does not match that of the theoretical one, and vice versa. For example, in the upper panels of Fig. 4, we compare the similarity via the K-S test of the mock observational data (shown in blue curves) and the theoretical curves (the orange, the yellow, and the green curves are assumed $H_{0}=65$, $H_{0}=70$, and $H_{0}=75\,{\rm km\,s^{-1}\,Mpc^{-1}}$, respectively). The similarity between mock observational data and the yellow curve is the highest (p-value is the highest). Therefore, we infer the p-value as a probability as a function of $H_{0}$ to constrain $H_{0}$. After we conduct K-S tests between the mock observational data and the theoretical CDFs under different $H_{0}$ and $k$, we are able to constrain $H_{0}$ from the statistical p-value as a function of $H_{0}$ and $k$.

For both mock observational data and the theoretical data, we consider a distribution of the velocity following the power-law distribution (e.g., Lister & Marscher, 1997; Lister, 2003; Cara & Lister, 2008; Ajello et al., 2012; Lister et al., 2016; Yuan et al., 2018):

where $C$ is a normalised factor, $k$ is a power-law index. We assume $k=-1.5$ for the mock observational data, which is the best fit suggested by superluminal motions (Lister & Marscher, 1997). The interval of the velocity is set to be $1.01\leq\gamma\leq 100$ (e.g., Yuan et al., 2018). The samples are randomly distributed within the spatial volume (i.e., population has constant space density) between $0<z<0.02$ and $0<\theta<\pi/2$. Note that for 3-dimensionally oriented two-sided jets, random jet viewing angles in the population $P(\theta,\theta+d\theta)\propto\sin(\theta)d\theta$ (e.g., Law et al., 2009). The effects of Doppler orientation bias may be minimal for $z<0.02$ and the effects are ignored in the simulation for simplification. The $H_{0}$ of mock observational data is set to be $70\,{\rm km\,s^{-1}\,Mpc^{-1}}$. As a result, we calculate $\mu_{r}$ and $\mu_{a}$ for each jet according to Eq. 2. After the simulation, we construct $H_{\rm 0,min}$ based on Eq. 4. As for the uncertainty, it depends on not only the angular resolution and the distance to the object but also depends strongly on total time coverage and cadence, the brightness temperature of the moving feature, the relative stability of the core (reference point), confusion with nearby jet features, and possible accelerations/non-linear motions (e.g., Lister et al., 2021; Weaver et al., 2022). For simplicity, we simulate $\mu_{a}$ and $\mu_{r}$ with the observational uncertainty of 10%, 5%, and 1%, which is based on the uncertainties in the MOJAVE program (Homan et al., 2021). In addition, we assume three sets of jets: 100, 200, and 500 jets, which correspond to $\sim 20\%$, $\sim 40\%$ and $\sim 100\%$ of the number of the jets in the MOJAVE program. Furthermore, there are $\sim 120$ radio galaxies at $z<0.02$ with jets/lobes structure from 2 Micron All-Sky Survey (2MASS) Redshift Survey (van Velzen et al., 2012). We expect that observing the nearby sources with higher sensitivity radio telescopes will increase the sample size. To estimate errors of $H_{\rm 0,min}$, we conducted each simulation 1,000 times by adding randomized errors (i.e., 10%, 5%, and 1%) to $\mu_{a}$ and $\mu_{r}$.

Aside from the highest values of $H_{\rm 0,min}$, a comparison of the whole distribution between theoretical and mock data is also meaningful. Hence, we propose to perform the K-S test to confirm if the two distributions are statistically different. By assuming different cosmology and simulating a large enough number of jets (i.e., $10^{5}$), we perform K-S tests between the mock observational data with the theoretical CDF. Figs. 6 to 11 show examples of the p-values as a function of $H_{0}$ and $k$ from K-S tests between mock observational data and theoretical profile under different configurations. As expected, all of the constraints are centred on the value we assume for the mock observational data ($H_{0}=70\,{\rm km\,s^{-1}\,Mpc^{-1}}$ and $k=-1$). Figs. 6, 7, and 8 demonstrate the constraints of increasing sample sizes under a fixed accuracy (i.e., $1\%$, $5\%$, and $10\%$, respectively). Figs. 9, 10, and 11 demonstrate the constraints of improving accuracy under a fixed sample size (i.e., 500, 200, and 100 jets, respectively). Our result indicates that when $k$ and $H_{0}$ are close to the real physics in the universe, p-values are higher. However, there is a degeneracy when the accuracy of jets is high. Figs. 6 to 8 demonstrates that if we increase sample sizes, the constraint becomes tighter while the degeneracy remains. Interestingly, we also find that increasing sample sizes leads to tighter constraints on both power-law index and the Hubble constant at moderate accuracy (i.e., $10\%$ and $5\%$) while at $1\%$ accuracy, increasing sample sizes leads to tighter constraints on power-law index more. According to Figs. 9 to 11, improving accuracy results in better constraints in the Hubble constant compared with the power-law index in all cases but it also alleviates the degeneracy.

In Fig. 14, we show the expected error to be constrained of $H_{0}$ (the upper panel) and $k$ (the lower panel) under different sample sizes and the accuracy of the measurement. Compared with the uncertainties from the SH0ES (Riess et al., 2021) and Planck (Planck Collaboration et al., 2020) measurements, our result indicates that if we can meet the criteria of $1\%$ accuracy, the error of $H_{0}$ constrained from the K-S tests will be more accurate. It is also clear that at $10\%$ and $5\%$ uncertainties, increasing sample sizes will lead to tighter constraints. However, the constraint saturates at $1\%$ accuracy, which suggests increasing the sample size cannot improve the uncertainties as good as those at $10\%$ and $5\%$. As can be seen in Fig. 14, $H_{0}$ can be determined within $1\%$ with 500 jets with $1\%$ accuracy. Besides, increasing the sample size alone to 500 jets with $10\%$ accuracy meets the precision at the SH0ES level. As for the power-law index $k$, increasing sample sizes lead to tighter error compared with improving the accuracy of the proper motion.

We also note that due to the Doppler beaming effect (see the details in the Sec. 2), many receding jets are not easily observable. In our simulations, we simply assume that all jets among our samples have measurable receding jets. The Doppler de-boosting effect makes the receding jets dimmer. Improving the integration time and the sensitivity of the telescope will help us to address this issue in the future. Besides, our method requires a large enough sample size to have a meaningful statistical result from the K-S test. For instance, in the MOJAVE program, there are $\sim 90$ jets that are two-sided out of $\sim 500$ samples. Many of the samples are only detected one-sided due to the Doppler boosting effect. We suggest three ways to increase the sample size. The first way is to increase the total sample of radio sources. On-going and future radio sky surveys such as the Very Large Array Sky Survey (VLASS; Lacy et al., 2020) which is expected to detect 10 million radio sources including radio galaxies (with two-sided jets). The second way is to improve the sensitivity of the telescope. Once the limiting magnitude reaches the brightness of the receding jets, we will be able to measure the proper motions from both jets. For example, for an approaching jet with $\gamma=3$ and $\theta=\pi/6$, the Doppler factor is $\sim 1.82$. Therefore, the Doppler factor of the receding jet is expected to be $\sim 0.18$. Adopting $p-\alpha=3$ in the Eq. 5, the receding jet is $\sim 1000$ times fainter than the approaching jet due to the Doppler deboosting effect (and boosting effect for the approaching jet). If the SNR of the approaching jet is 1000 (e.g., Lister et al., 2019; Baczko et al., 2019), improving the sensitivity to $\sim 3$ to $\sim 5$ times better will be able to detect it. Hovatta et al. (2009) calculated the Doppler factors of quasars, BL Lacertae objects, and radio galaxies. The sources with low-Doppler factors are worth further observation on the receding jets when improving sensitivity in the near future. The third way is to monitor the two-sided jets which have not been monitored before. For instance, in the catalogue of FR I radio galaxies (with 219 FR I radio galaxies), FRICAT (Capetti et al., 2017), many jets are two-sided but only have image data, which suggests those jets have no measurement of proper motions.

In addition, it might be difficult to collect many proper motions of twin-jets with high Lorentz factor, especially the receding components. Therefore, we also test if the method still works under different maximum values of the velocity ($\gamma_{\rm upper}$). In our mock observational data, we adopt $\gamma_{\rm upper}=100$. We show how CDF changes under different $\gamma_{\rm upper}$ in Fig. 12. Regarding $H_{0}$ and $k$, $\gamma_{\rm upper}$ does not affect the CDF significantly. Here, we gradually reduce the maximum value of the Lorentz factor for the theoretical CDF and conduct K-S tests with the mock observational data. The assumption of the mock observational data is identical to that in Fig. 4, with 100 jets and $1\%$ accuracy. The result is shown in Fig. 13, which suggests that the critical value is at $\gamma_{\rm upper}\sim 30$. Values with $\gamma_{\rm upper}>30$ have p-values higher than 0.05. Once the maximum value of the Lorentz factor is smaller than the critical value, 30, our method fails with p-values smaller than 0.05.

Based on Fig. 2 and Tab. 1 of Pracy et al. (2016), there are $\sim 2200$ radio galaxies at $z<0.3$. Assuming the ratio of two-sided jets and total jet sample in the MOJAVE program ($\sim 90/\sim 500\sim 18\%$), we expect there are $2200\times 18\%\sim 400$ two-sided jets in the radio galaxy catalogue. Conducting VLBA to observe such samples will bring the two-sided jets to 500 including the 90 in the MOJAVE program. In the MOJAVE program, they monitored 447 AGN jets from 25-year multi-epoch observations (Homan et al., 2021). Thus, if aiming at these two-sided jets, $\sim 20$ years of observation at different epochs to measure the proper motion (e.g., $\sim 1\,$mas/yr for a jet of $\beta\sim 0.6$ at $z\sim 0.005$; $\sim 100\,\rm\mu as/yr$ for a jet of $\beta\sim 0.9$ at $z\sim 0.3$; Baczko et al., 2019) with VLBA ($\sim$mas resolution) may increase the sample size of two-sided jets to $\sim 500$ in the local Universe. The uncertainties of the proper motions reach $\rm\sim 5\mu as$ level with integration times of $30$ to $50$ minutes (e.g., Lister et al., 2019). Also, the radio galaxies included in Pracy et al. (2016) have a criteria with flux density at $1.4\,{\rm GHz}>2.8\,{\rm mJy}$. Therefore, MOJAVE VLBA with image rms levels of $\sim 0.1\,{\rm mJy\,beam^{-1}}$ is promising to detect and monitor these radio galaxies.

Last but not least, if there are too few samples in the nearby Universe, we can extend to higher redshift as long as the approximation of the FLRW metric still works. Moreover, even if the approximation fails at high redshift, if we assume cosmological parameters other than $H_{0}$, we can extend this method to higher redshifts. For instance, for a $\beta=0.9999$ ($\gamma\sim 70$), $\theta=\pi/2$, the $H_{0}$ is only $\sim 1\%$ more than the assumed $H_{0}$ at $z=0.05$. We can also revise the redshift distribution based on the distribution from observation for theoretical (simulated) CDFs.

Based on FLRW cosmology and geometrical relation of the proper motion, the lower limit ($H_{\rm 0,min}$) of the Hubble constant ($H_{0}$) can be determined only with three cosmology-free observables: redshift, approaching proper motion, and receding proper motion. We propose a new method, performing K-S tests between the observational (mock) data and the theoretical CDF of different $H_{0}$ and $k$, in order to constrain $H_{0}$ and $k$. We simulate $10^{5}$ jets as a numerical distribution of $H_{\rm 0,min}$ and $k$ between $65\,{\rm km\,s^{-1}\,Mpc^{-1}}<H_{0}<75\,{\rm km\,s^{-1}\,Mpc^{-1}}$ and $-1<k<-2$, as theoretical distributions and create mock observational data. The result shows that we can simultaneously constrain $H_{0}$ and $k$. If the values of $H_{0}$ and $k$ are close to the real universe, p-values are higher. We find that there is a degeneracy between $H_{0}$ and $k$ when the accuracy is high (e.g., $10\%$). Interestingly, we also find that increasing sample sizes leads to tighter constraints on both power-law index and the Hubble constant at moderate accuracy (i.e., $10\%$ and $5\%$) while at $1\%$ accuracy, increasing sample sizes leads to tighter constraints on power-law index more. Improving accuracy results in better constraints in the Hubble constant compared with the power-law index in all cases but it also alleviates the degeneracy. In the future, we will be able to alleviate the Hubble tension by the distribution of $H_{\rm 0,min}$ calculated from the proper motion and we could also further constrain $H_{0}$ and $k$.

We thank the anonymous referee for useful comments and constructive remarks on the manuscript. TG acknowledges the support of the National Science and Technology Council of Taiwan through grants 108-2628-M-007-004-MY3 and 111-2123-M-001-008-. TH acknowledges the support of the National Science and Technology Council of Taiwan through grants 110-2112-M-005-013-MY3, 110-2112-M-007-034-, and 111-2123-M-001-008-.

The data underlying this article will be shared on reasonable request to the corresponding author.