Model-independent Way to Determine the Hubble Constant and the Curvature from the Phase Shift of Gravitational Waves with DECIGO

In this Letter, we propose a model-independent method to determine the Hubble constant and curvature simultaneously by taking advantage of the possibilities of future spaceborne gravitational-wave detector DECIGO in combination with the radio quasars as standard rulers. Similarly to the redshift drift in the electromagnetic domain, accelerating expansion of the Universe causes a characteristic phase correction to the gravitational waveform detectable by DECIGO. Hence, one would be able to extract the Hubble parameter H(z). This could be used to recover a distance–redshift relation supported by the data not relying on any specific cosmological model. Assuming the FLRW metric and using intermediate-luminosity radio quasars as standard rulers, one achieves an interesting opportunity to directly assess the H 0 and Ω k parameters. To test this method, we simulated a set of acceleration parameters achievable by future DECIGO. Based on the existing sample of 120 intermediate-luminosity radio quasars calibrated as standard rulers, we simulated much bigger samples of such standard rulers possible to obtain with very long baseline interferometry (VLBI). In the case of (N = 100) of radio quasars, which is the size of the currently available sample, the precision of the cosmological parameters determined would be σH0=2.74 km s−1 Mpc−1 and σΩk=0.175 . In the optimistic scenario (N = 1000) achievable by VLBI, the precision of H 0 would be improved to 1%, which is comparable to the result of σH0=0.54 km s−1 Mpc−1 from Planck 2018 TT, TE, EE+lowE+lensing data, and the precision of Ω k would be 0.050. Our results demonstrate that such combined analysis, possible in the future, could be helpful in solving the current cosmological issues concerning the Hubble tension and cosmic curvature tension.


INTRODUCTION
Over the past few decades, the observations of type Ia supernova (SN Ia) revealed that the expansion of the universe is accelerating (Riess et al. 1998(Riess et al. , 2007)).Based on the cosmological principle (supporting homogeneous and isotropic universe) and Einstein's general theory of relativity (GR), it is generally accepted that this acceleration is caused by dark energy with negative pressure.The simplest solution is to use the cosmological constant term Λ as a model of dark energy driving the accelerated expansion of the universe.With addition of cold dark matter, it forms the ΛCDM model.This model has been supported by the vast majority of astronomical observations over the past few decades (Planck Collaboration et al. 2016;Caldwell et al. 1998;Anderson et al. 2014;Abbott et al. 2018).However, with the accumulation of precise astrophysical observations, some problems gradually emerged.For instance, there is a 4.4σ tension be-tween the Hubble constant (H 0 ) measurements inferred within the ΛCDM model from cosmic microwave background (CMB) anisotropy (both temperature and polarization) data (Planck Collaboration et al. 2020) and the local measurements of the Hubble constant by the Supernova H 0 for the Equation of State collaboration (SH0ES) (Riess et al. 2019).This inconsistency may be caused by unknown systematic errors in astrophysical observations or unrevealed new physics which is significantly different from ΛCDM (Freedman 2017;Di Valentino et al. 2021).Moreover, some works suggested that besides the H 0 tension, spatial curvature of the universe is still an open issue (Di Valentino et al. 2021, 2020;Handley 2021).More specifically, combining the Planck temperature and polarization power spectra data, Planck Collaboration et al. (2020) claimed that closed universe (Ω k = −0.044+0.018  −0.015 ) was supported by pure CMB data.However, the combination of Planck lensing data and low redshift baryon acoustic oscillations (BAO), revealed that flat universe with Ω k = 0.0007 ± 0.0019 was preferred.It should be emphasized that both methods depend, to some extent, on a particular cosmological model assumed, i.e. the non-flat ΛCDM.To better understand inconsistency between the Hubble constant and the curvature of universe inferred from local measurements and Planck observations, it is necessary to seek a completely model-independent method that can constrain both H 0 and Ω k simultaneously (Collett et al. 2019;Wei & Melia 2020;Qi et al. 2021;Liu et al. 2022Liu et al. , 2023)).
Traditionally, if the Hubble parameter H(z) and the luminosity distance D L (z) or angular diameter distance D A (z) are known at the same redshift z, one can directly measure H 0 and Ω k simultaneously.The difficulty, however, is that possibilities of measuring Hubble parameters in the electromagnetic window are limited.Current cosmic chronometer sample providing such measurements via differential ages of passively evolving galaxies consists of only of 31 H(z) data.Fortunately, gravitational wave (GW) observations provide such a chance to obtain H(z) directly.Recently, the discovery of the first GW event, GW150914 (Abbott et al. 2016), caused a great excitement in astronomy and physics communities.GW signals from inspiraling and merging compact binaries enable absolute measurements of the luminosity distance D L (z) (Schutz 1986), hence they are often referred to as standard sirens.Subsequent detection of the GW170817 event (Abbott et al. 2017) powered by the merger of binary neutron stars (NSs) was accompanied by electromagnetic counterpart detected in gamma-rays, X-rays, in the optical and radio bands.If the redshifts to standard sirens can be independently determined from their electromagnetic (EM) counterparts, such GW signals can be used for cosmology.More importantly, if the adiabatic inspiral phase could be monitored long enough, the accelerating expansion of the universe would cause an additional phase shift (more details in Sec.2.1) in the gravitational waveform, which would allow us to measure the cosmic acceleration directly (Yagi et al. 2012).This idea is similar to the redshift drift (Sandage 1962;Loeb 1998) but requires shorter cadence of observations.Current LIGO-Virgo-Kagra (LVK) network of detectors has negligible possibility to detect cosmic acceleration because chirp signals in their frequency range are present from seconds to a couple of minutes.The future generation of space-borne DECi-hertz Interferometer Gravitational-wave Observatory (DECIGO) will be a game changer, because the inspiral signals, ultimately to be detected by LVK ground based detectors, could be discovered and monitored a few years before the final coalescence.Thus, a cosmological model-independent way to provide H(z) with covering a wide redshift range in GW domain is possible (Seto et al. 2001).
Inspired by above mentioned ideas, we will determine the curvature and the Hubble constant simultaneously using the Hubble expansion rate from simulated future space-borne detector DECIGO in combination with the "angular size-distance" relation of ultra-compact structure in radio quasars (QSOs) from the very-long-baseline interferometry (VLBI) observations.The advantage of our work consists in assuming only the homogeneity and isotropy of the universe, thus, it is a geometric and cosmological model-independent way for the measurement of Hubble constant and spatial curvature simultaneously.The paper is organized as follows: In Sec. 2, we present the methodology and simulated data.The results and discussion are given in Sec. 3. Finally, we summarize our findings in Sec. 4.

DATA.
Assuming that the universe is homogeneous and isotropic on large scales, the space-time of our universe can be described by Friedmann-Lemaître-Robertson-Walker (FLRW) metric (Weinberg 2008): (1) where t is the cosmic time and r, θ, ϕ are the comoving spatial coordinates.The scale factor a(t) is the only gravitational degree of freedom and its evolution is determined by matter and energy content in the universe.The dimensionless curvature k = −1, 0, +1 corresponds to open, flat and closed universes, respectively.With such a metric, the angular diameter distance D A (z) can be expressed as (Weinberg 2008) where D H = c/H 0 is known as the Hubble distance, the dimensionless Hubble parameter E(z) is defined as H(z)/H 0 , where H(z) is the universe expansion rate at redshift z and H 0 = H(z = 0) is the Hubble constant.The curvature parameter Ω k is related to k as Ω k = −c 2 k/(a 0 H 0 ) 2 , where c is the speed of light.For convenience, we denoted S k (x) = sin(x), x, sinh(x) for Ω k < 0, = 0, > 0. If the angular diameter distance D A (z) and Hubble parameter H(z) are known, then we can constrain both H 0 and Ω k .Standard procedure is to assume the cosmological model like ΛCDM or wCDM (i.e. the cold dark matter model with dark energy not being exactly cosmological constant) and use H(z) formula suggested by the Friedman equation to calculate theoretical value of angular luminosity distance D th A (z; p) at the distance corresponding to the standard ruler, whose observed distance D obs A (z) can be measured.Then the chi-square based likelihood is maximized (or chi-square minimized) to obtain the parameters p, in such case p comprise also cosmological model parameters like Ω m , Ω Λ , etc..In this work, we use the opportunity to measure the acceleration parameter X(z) (more details see Sec.2.1, Eq. ( 5) and thereafter) from an additional phase shift in the gravitational waveform based on the expected performance of the future DE-CIGO space-based detector.This would enable to obtain directly the Hubble parameter H(z) at the redshift of NS-NS binary monitored with DECIGO.Having the large sample of acceleration parameters X(z i ), i = 1, ..., N we use Artificial Neural Network (ANN) (Wang et al. 2020b) to reconstruct the X(z) function supported only by the data, without assuming any particular cosmological model.This reconstructed relation (with the corresponding uncertainty band) can be used to calculate the distance Eq. (2) to every standard ruler used to measure the distance.In this case the set of cosmological parameters p comprise only H 0 and Ω k .In the role of standard rulers, we will use the ultra-compact radio quasars.Detailed description of GW data and radio quasars is given below.
] and underlies the phenomenon of the redshift drift.Its order of magnitude is the inverse of the age of the Universe, hence this correction is usually being neglected in most cosmological considerations.In the context of GW signals from inspiralling binaries, the relevant argument of the wave strain h(t) is the time to coalescence ∆t = t c −t with t c denoting coalescence time, measured in the observers frame.One can relate it to the time to coalescence ∆t e measured in the source frame by: ∆t = ∆T +X(z)∆T 2 , where ∆T = ∆t e (1+z) is the usual first order correction due to cosmological time dilation and the second order (Taylor expansion) factor has been denoted X(z) and is called the acceleration parameter (more details in the text following Eq.( 5)).Calculating Fourier transformed h(f ) and using stationary phase approximation (Yagi et al. 2012) it is straightforward to see that accelerating expansion introduces the phase shift Ψ acc = −2πf X(z)∆T 2 , which is additional with respect to the first order treatment neglecting the acceleration.1 See Fig. 1 of Kawamura et al. (2019) for an artistic illustration of this effect, i.e., GW signals coming from a neutron star binary at a given redshift with and without the acceleration of the expansion of the Universe.It was demonstrated (Seto et al. 2001;Nishizawa 2012) that in a gravitational wave signal emitted by a neutron star binary at the redshift of z = 1, accelerated expansion of the Universe would cause a phase delay of about l sec during the 10-year observing time.
Waveform templates used in matched-filtering technique do not account for acceleration.Hence, with ∼ 10 7 orbital cycles monitored during 10 years the accumulated phase shift due to acceleration would be detectable.As a second-generation space-based GW detector, DECIGO (Kawamura et al. 2011) is designed to improve the detection sensitivity of GW at lower frequencies, with its most sensitive frequencies between 0.1 and 10 Hz.It can be achieved by the proposed four clusters of spacecrafts, each of which having triangular 1000 km Fabry-Perot Michelson interferometer.This design improves the determination of GW sources' position in the sky, as well as the detection sensitivity of GW at lower frequencies.The expected sensitivity of DECIGO is ∼ 4 × 10 −24 Hz −1/2 for one cluster, which enables the early detection of inspiraling sources, and more importantly, the direct measurement of the acceleration of the expansion of the universe.The angular resolution is ∼ 1 arcsec 2 , which enables the unique identification of the host galaxy of the binary, and improves the precision of distance measurements (Seto et al. 2001).Therefore, DECIGO will provide us a new perspective on studying acceleration of cosmic expansion and on gravitational wave cosmology (?).
Besides, benefiting from the precise time prediction of coalescence time and high resolution of sky location of GW sources, redshift determination from EM follow-up observations would be more reliable and frequent.DE-CIGO is expected to observe 10 6 GW events coming from the neutron star binaries up to the redshift z ∼ 5 (Kawamura et al. 2019;Nishizawa 2012).For our purpose, we follow the work of Holz & Hughes (2005), where they assumed that the GW event randomly occurs in any one of the galaxies, and considering the large number of spectroscopic and photometric redshift measurements in future galaxy surveys (such as JDEM/WFIRST (Gehrels 2010;Levi et al. 2011), BigBOSS (Schlegel et al. 2010)/DESI (DESI Collaboration 2016, 2023) or 4MOST (de Jong et al. 2019)), they concluded that 10,000 GW events from DECIGO with redshift determinations could be obtained by EM follow-up observations.Considering the GW signal from NS-NS binary system whose components have masses m 1 , m 2 , one can define the redshifted chirp mass , where M = m 1 + m 2 is the total mass of the binary system, η = m 1 m 2 /M 2 is the symmetric mass ratio, and z is the source redshift.In the frequency domain, gravitational waveform including the effects of the cosmic acceleration can be written as where is the phase correction due to the cosmic acceleration, where x ≡ (πM z f ) 2/3 and Ψ N (f ) ≡ 3 128 (πM z f ) −5/3 .Terms with x −n represent the n -th Post Newtonian (PN) order, hence Ψ acc (f ) is of "4 -PN" order and becomes more important at lower frequencies.The acceleration parameter X(z) is defined as (Seto et al. 2001 Note such a parameter can be rewritten as X(z) = [ ȧ(0) − ȧ(z)]/2.Thus it directly reflects whether the universe expansion has accelerated (X(z) > 0) or decelerated (X(z) < 0).Angular diameter distance, given by Eq. ( 2) can be expressed in terms of the acceleration parameter in the following form: (6) The gravitational waveform without cosmic accelera-tion is given by h (7) where the amplitude includes the luminosity distance D L , which indicates that it can be measured directly.For the phase Ψ(f ), we use the restricted-2PN (Kidder et al. 1993) waveform including spin-orbit coupling at the 1.5PN order, which is given by where the direction of orbital angular momentum and the direction of the source from the Sun are also taken into account (Yagi & Tanaka 2010).The polarisation amplitude A pol,α (t) and the polarisation phase φ pol,α (t) represent the waveform measured by each detector (α = I, II).The Doppler phase φ D (t) denotes the difference between the phase of the wave front at the detector and the phase of the wave front at the solar system barycenter.All these quantities are explicitly given in Yagi & Tanaka (2010).
The waveform h(f ) acc in Eq. ( 3) depends on 11 parameters: where t c and ϕ c represent the time of and the phase at the final coalescence, β is the spin-orbit coupling parameter.( θS , φS ) indicates the direction of the source in the solar system barycenter frame, and ( θL , φL ) specifies the direction of the orbital angular momentum.For simplicity, we set m 1 = m 2 = 1.4M ⊙ and take t c = ϕ c = β = 0.For each fiducial redshift z, we randomly generate 10 4 sets of ( θS , φS , θL , φL ), and for each set, we calculate the uncertainties of the acceleration parameter X estimated, using the Fisher matrix formalism: where partial derivatives are with respect to the parameters and S h (f ) is the DECIGO noise power spectrum, which analytical expression is (Kawamura et al. 2006(Kawamura et al. , 2019;;Yagi & Seto 2011): + 4.94 The lower cut-off of frequency in Eq. ( 10) is related to the observation duration ∆T obs .The upper cut-off frequency of the detector f max is set equal to 100 Hz.In addition, since there will be eight uncorrelated interferometric signals from DECIGO, the Fisher matrix above should be multiplied by 8 (Kawamura et al. 2011;Yagi & Seto 2011).One can use Fisher analysis to estimate the measurement accuracy of the acceleration parameter X by marginalizing other parameters in the Fisher matrix where the σ X is the root-mean-square uncertainty for the j-th parameter (the acceleration parameter X), and the Γ −1 jj is the covariance matrix element for parameter X calculated from the Eq. ( 10).
Finally, the redshifts of GW sources are sampled according to the merger rate of double compact objects that reflects the star formation history (Dominik et al. 2013), which is called "rest frame rates" in the cosmological scenario2 .In our simulation, we assume flat ΛCDM as the proposed fiducial model with the cosmological parameters derived by Planck 2018 measurements (Planck Collaboration et al. 2020).From these simulated results ANN technique (Wang et al. 2020b) was used to reconstruct the X(z) function and the corresponding uncertainty band in the redshift range z ∈ [0, 5].Machine learning methods have become very popular in astrophysics and cosmology.Motivated by successful applications of ANN (Liu et al. 2021;Wang et al. 2020a;Qi et al. 2023), we chose this method, which is a general method that could reconstruct a function from any kind of data without assuming any particular parameterization of the function.Moreover, this non-parametric technique does not assume Gaussianity of measured random variables, hence it is a completely data driven approach.Based on the future generation of space-borne GW detector DECIGO, we simulated a dataset of 10 4 measurements of the acceleration parameter X(z).However, the acceleration parameters have large uncertainties, because there are too many unknown parameters in the waveform h(f ) acc , which leads to large uncertainties in the estimation from the Fisher matrix.Therefore, we divided this sample into 50 bins and adopted ANN to reconstruct the acceleration parameter as a function of redshift.We trained the ANN network on the simulated X(z) data, and then predicted the X(z) at other redshifts (where there were no GW data available).In general, the artificial neural network consists of an input layer, one or more hidden layers and an output layer.In this work, the input of the neural network is the redshift z, while the output is the corresponding cosmic acceleration parameter X(z) and its respective uncertainty σ X at that redshift.More details of data reconstruction using ANN can be found in (Wang et al. 2020a;Liu et al. 2021;Qi et al. 2023).Our strategy of non-parametric reconstruction of X(z) function implies that one does not need to match the QSOs and NS-NS sources by redshift, but rather for each QSO the corresponding value of acceleration parameter will be reconstructed from the total sample of GW signals.The values of the acceleration parameters X(z) and their corresponding uncertainties are shown in Fig. 1.Let us remark that individual acceleration parameters have considerable uncertainties, hence we binned the data in the redshift bins of width of width ∆z = 0.1 and in each redshift bin we calculated the weighted mean, . The ANN network was subsequently trained on these n = 50 bins.The final uncertainties level is similar as in the recent work (Sun et al. 2024), which used convolutional neural network (CNN) to derive the cosmic acceleration parameters.

Angular diameter distances from radio quasars
Quasars, as the brightest sources in the universe, have still a great potential to be used as useful cosmological probes, despite of high observed dispersion and extreme variability.Observations of these objects have long been sought to study cosmology beyond the limits of supernovae.With the advances in observational technology over recent decades, quasars can be observed at multiple wave-bands, and great efforts have been made to use them as standard candles or standard rulers in cosmology, using the Baldwin effect (Baldwin 1977), the Broad Line Region radius-luminosity relation (Watson et al. 2011), the properties of highly accreting quasars (Wang et al. 2013), etc. Recently, Risaliti & Lusso (2019) attempted to use quasars as standard candles through the non-linear relation between their intrinsic UV and X-ray emission.However, the large dispersion still remains the greatest problem in calibrating the QSOs as standard candles.In this work, we focus on the "angular size-distance" relation of ultra-compact structure in radio quasars both observed so far and those expected from the future VLBI observations.With the signal received at multiple radio telescopes across the Earth's surface, together with the registered correlated intensities considering the different arrival times at various facilities, the characteristic angular size of a distant radio quasar is defined as (Gurvits 1994) where B represents the interferometer baseline measured in wavelengths, and Γ = S c /S t is the visibility modulus which is defined by the ratio between the total and correlated flux densities.In a more detailed study involving the compact structure of radio quasars, Gurvits (1994) showed that the dispersion in linear size is greatly mitigated by retaining only the sources with a flat spectrum (−0.38 < α < 0.18).Subsequently, the application of compact radio sources to cosmology has made a breakthrough (Jackson & Dodgson 1997;Vishwakarma 2001;Lima & Alcaniz 2002).With gradually refined selection technique and the elimination of possible systematic errors, Cao et al. (2017) compiled a sample of 120 intermediate-luminosity quasars (10 27 W/Hz < L < 10 28 W/Hz) selected from 613 milliarcsecond ultra-compact radio sources covering the redshift range 0.0035 < z < 3.787 in 2.29 GHz VLBI survey.More importantly, they demonstrated that intermediateluminosity compact radio quasars can be used to obtain constraints on cosmological parameters.The details regarding the angular sizes θ(z) and redshifts z of 120 compact radio sources are listed in Table 1 of Cao et al. (2017), which extends the Hubble diagram to the redshift range 0.46 < z < 2.76.The angular size of the compact structure in radio QSOs can be written as the dimensionless expansion rate expansion rate E(z).
The parameter l m denotes the linear size scaling factor describing the apparent distribution of radio brightness within the core, which needs to be calibrated with external indicators such as SN Ia or cosmic chronometers.In this analysis, we adopt the calibrated value of this cosmological standard ruler at 2.29 GHz as l m = 11.03±0.25 pc, through a cosmological model-independent technique based on cosmic chronometer measurements (Cao et al. 2017).This sample of VLBI ultra-compact structure in radio quasars, which was extended to multi-frequency analysis (Cao et al. 2018), has been widely applied in various cosmological applications at high-redshifts (Li et al. 2017;Zheng et al. 2017;Qi et al. 2017Qi et al. , 2019a;;Xu et al. 2018;Liu et al. 2020bLiu et al. , 2021aLiu et al. , 2023)).
Considering that currently no other sample of calibrated radio quasars is available, we further simulate the possible future data on compact radio quasars from VLBI surveys based on better uv-coverage.Our simulation is based on the flat ΛCDM with H 0 = 67.4km s −1 Mpc −1 and Ω m = 0.308 adopted from the Planck results (Planck Collaboration et al. 2020), i.e. the same ΛCDM model as used for GW simulations.Taking the linear size of intermediate-luminosity quasars as l m = 11.03 ± 0.25 pc and following the redshift distribution of QSOs from (Palanque-Delabrouille et al. 2016), we simulate 1000 "angular size-redshift" data representative of compact radio sources in the redshift range 0.50 < z < 5.00.The fractional statistical uncertainty of the angular size θ was taken at a level of 3%, which is resonable for both current and future VLBI surveys based on better uvcoverage (Pushkarev & Kovalev 2015).Considering the intrinsic variance in the size of the AGN cores, derived from high angular resolution observations, an additional 10% systematical uncertainty in the observed angular sizes was also assumed in our analysis.In addition, we assume that the measured angular sizes follow a Gaus-sian distribution θ means = N (θ fid , σ θ ), where θ fid is obtained from Eq. ( 15) under assumed fiducial cosmological model.Simulated data are shown in Fig. 2. It should be remarked that sources with z < 0.5 are ignored, because the era of quasar formation has ended, and the nature of quasars has changed dramatically.This problem is further expanded and elaborated in (Gurvits 1994), which indicates that only the high redshift portion of the radio quasars can be used as a standard rulers.More details of the specific procedure of QSO simulation can be found in Cao et al. (2019); Qi et al. (2019a,b).

RESULTS AND DISCUSSION
In order to obtain the best-fit values of H 0 and Ω k with corresponding uncertainties, we used the Markov Chain Monte Carlo (MCMC) method (Foreman-Mackey et al. 2013) to minimize the χ 2 objective function: where the p represents the Hubble constant H 0 and cosmic curvature Ω k , the observed angular diameter distance is inferred from angular size-distance relation of ultra-compact structure in QSOs D obs A,i (z) = l m /θ i (z).Theoretical counterpart, D th A,i (z; p) is derived from Eq. ( 6) where the ANN reconstructed X(z) function was used to calculate respective integrals applying a simple trapezoidal rule.The uncertainties σ 2 Combining the simulated QSO data from future VLBI observations together with the acceleration parameters from DECIGO registered GW events, we obtained the best-fit values of H 0 = 67.19+0.77  −0.74 km s −1 Mpc −1 and Ω k = 0.022 +0.051 −0.046 at 68.3% confidence level.The posterior joint distribution function along with marginal distributions is shown in Fig. 3. Our results are consistent with the fiducial values of H 0 = 67.4km s −1 Mpc −1 and Ω k = 0 assumed, which proves the robustness of our simulation.Compared with the previous model-independent works based on the cosmic chronometer measurements H(z) and distances from popular cosmological probes, our results have a higher precision than the result of combined Pantheon SN Ia plus H(z) giving Ω k = 0.63 ± 0.34 (Wang et al. 2020a), and the result of UV+X-ray quasars plus H(z) yielding Ω k = −0.92± 0.43 (Wei et al. 2020).More importantly, such combination of radio and GW windows may achieve constraints with much higher precision on H 0 , namely with 1% uncertainty.Moreover, it creates an opportunity to estimate the Hubble constant and cosmic curvature using the objects at much higher redshifts (z ∼ 5) than currently available.Therefore the future measurements of this kind would shed light on the Hubble tension, which is recognized as an issue of local vs.CMB probes.
The above mentioned result was obtained assuming the optimistic scenario of having N = 1000 intermediate luminosity radio quasars calibrated as standard rulers.Therefore, we investigated the uncertainties of best fitted parameters as a function of QSO sample size N .Let us stress that the bottleneck is the QSO sample, while DECIGO is expected to provide enough data to extract acceleration parameters allowing for ANN reconstruction of X(z) relation.The results are displayed in Fig. 4-5, in which two cases are considered separately.The first is the case of H 0 and Ω k treated as free parameters.The second one is the precision of one parameter with a fixed prior on another one.Regarding the first case, one can see that even in the most conservative case (N = 50), the result of σ H0 = 4.36 km s −1 Mpc −1 i.e. 6% precision is encouraging for using the proposed cosmological model-independent method.When the number of standard rulers increases to N = 300, the precision of σ H0 = 1.36 km s −1 Mpc −1 is comparable to the result of σ H0 = 1.3 km s −1 Mpc −1 given by the SH0ES collaboration (Riess et al. 2019).If the sample size increases to N = 1000, the precision on H 0 will be improved to σ H0 = 0.77, which is comparable with the result of σ H0 = 0.54 km s −1 Mpc −1 from Planck 2018 TT, TE, EE+lowE+lensing data (Planck Collaboration et al. 2020).
The joint fit of the Ω k parameter, in the optimistic scenario (N = 1000) has the uncertainty of σ Ω k = 0.050.This is not competitive to the recent results from Planck lensing data, where a flat universe was precisely constrained with a 1σ uncertainty of σ Ω k = 0.002 (Planck Collaboration et al. 2020).However, some recent works (Di Valentino et al. 2021, 2020;Handley 2021) focused on the Hubble tension and other tensions in ΛCDM cosmology have raised the issue of the value of a spatial curvature of the Universe.More specifically, combining the Planck temperature and polarization power spectra data, these works showed that a closed universe (Ω k = −0.044+0.018 −0.015 ) was supported by CMB data, which provided an explanation of anomalous high lensing amplitude (Di Valentino et al. 2020) reported in Planck Legacy 2018.The robust conclusion of these works is to emphasize the importance of alternative and cosmological model-independent methods to constrain both H 0 and Ω k .Our method is one of this kind and the results are encouraging.The precision would be increased if instead of two-parameter fit, one would fit a single pa- rameter with the fixed prior on another one, which will be discussed below in more details.
Assuming the flat Universe, as shown by the red solid polyline in Fig. 4, the constraint on H 0 will be improved significantly.In this case, the constraint from only N = 100 QSOs would have uncertainty of σ H0 = 1.08 km s −1 Mpc −1 , which is slightly better than the result of σ H0 = 1.3 km s −1 Mpc −1 reported by the SH0ES collaboration.Only N = 500 QSOs are required to produce the results comparable with those obtained by Planck 2018 TT, TE, EE+lowE+lensing data.In the most optimistic scenario (N = 1000), we can get a result of σ H0 = 0.35 km s −1 Mpc −1 , which is at 0.5% level.All these results demonstrate that the combination of GW+QSO data proposed by us could become an important cosmological probe, precise enough to address the Hubble tension issue.With the fixed prior of H 0 = 67.4km s −1 Mpc −1 , the constraints on Ω k are also considerably improved.This can be understood in terms of degeneracy between Ω k and H 0 , which can also be seen from the joint constraint confidence regions of H 0 and Ω k shown in Fig. 3.In the most conservative case (N = 50) the uncertainty of Ω k is σ Ω k = 0.100, while in the most optimistic scenario (N = 1000), σ Ω k = 0.023.This is not as good as the result from the Planck 2018 data (Planck Collaboration et al. 2020), yet it refers to an alternative measurement using extragalactic objects in more local Universe than CMB.An interesting point is that as the number of QSOs increases, the constraint on Ω k does not improve significantly after reaching N = 500.From the statistical point of view, when the number of observations increases by N , the precision should increase by √ N .The saturation observed suggests the existence of a feature in the analysis that prevents from going lower.We have checked how the hypothetical assumption of zero uncertainty for the angular diameter distance influence the uncertainty of cosmological parameters inferred.It turned out that it has negligible effect σ H0 and σ Ω k were only slightly different from the above mentioned values (at the second significant digit).Better understanding of the origin of such an uncertainty floor would be de- sired and is left for future studies.Moreover, in future applications of our method with the real data, the influence of systematic effects cannot be ignored and should be studied carefully.
Our method is built on two pillars: GWs as a source of H(z) data and QSOs as source of D A (z) data.In order to gain some understanding of relative importance of uncertainties related to GW signals and QSO, we did the following test.Namely, based on N = 1000 QSO data paired with GW reconstructed X(z) values, we left the uncertainty of QSO data unchanged but doubled the uncertainty of GW inference.The results we obtained in this case were: σ H0 = 1.18 km s −1 Mpc −1 and σ Ω k = 0.070.Similarly, we kept the uncertainty of GW unchanged but doubled the uncertainty from QSOs.We obtained the precision of Hubble constant to be σ H0 = 1.04 km s −1 Mpc −1 and the precision of curvature parameter to be σ Ω k = 0.079.Fig. 6 displays the 1σ and 2σ confidence level contours for these two cases for estimated parameters.In conclusion, we found that the contribution of QSO and GW data uncertainties to the precision of parameter constraints were comparable.
Finally, we tested the influence of the GW data on the fits.For this purpose we simulated N = 1000 GW sources, reconstructed corresponding X(z) function and analysed jointly with N = 1000 QSOs.The result was H 0 = 65.69 ± 1.71 km s −1 Mpc −1 and Ω k = 0.03 ± 0.142 at 68.3% confidence level.One can see that with an order of magnitude smaller sample of GW signals the uncertainties raised 2 − 3 times relative to the main case studied.The precision would be at the present level attainable by local cosmological probes.Hence, the GW signal sample size and quality are crucial.

CONCLUSION
In this Letter, we proposed a cosmological modelindependent method to determine the Hubble constant and curvature parameter simultaneously based on the future generation of space-borne GW detector DECIGO in combination with the data regarding intermediateluminosity radio-quasars calibrated as standard rulers obtained from the VLBI observations.Since the accelerating expansion of the Universe would cause an additional phase shift in the gravitational waveform, we can extract the Hubble parameter H(z) useful for further cosmological studies.On the other hand, the future VLBI surveys will observe a large amount of compact structures in radio quasars, which would provide the angular diameter distances D A (z) based on a simple "angular size-distance" relation.Such combination of observations performed in the radio and GW windows creates a valuable opportunity to directly test H 0 and Ω k at much higher precision.
To fully demonstrate the potential of such combination of GW+QSO data-set, we simulated respective data representative of the expected yield of the DECIGO and VLBI.Simulation of the future VLBI data on QSOs was based on the properties of currently available sample of N = 120 intermediate-luminosity radio-quasars calibrated as standard rulers (Cao et al. 2017).The same fiducial cosmology was assumed in both GW and QSO simulations.In the most optimistic case of N = 1000 QSO calibrated as standard rulers, joint fits for H 0 and Ω k were consistent with the input simulation values and their 68.3%confidence level uncertainties were σ H0 = 0.75 km s −1 Mpc −1 and σ Ω k = 0.048, respectively.It means that achievable precision is at the level of 1%.Besides the joint fit on two parameters simultaneously, we considered fits of one parameter with fixed priors on the second one.In such cases precision improved as could be expected.In particular the precision of H 0 with a fixed prior on Ω k can reach 0.5% in the most optimistic case of N = 1000 QSOs.We also discussed the performance of the method on QSO samples of different size (from N = 50 to N = 1000).In particular regarding Ω k fit the uncertainty behavior did not reflect simple 1/ √ N behaviour, but saturated at N = 500.This phenomenon remains to be fully understood, but urges a careful consideration of systematic effects in future applications of the proposed method to the real data.
As a final remark, there are many potential ways to improve our results.In the future, multi-frequency VLBI observations will yield more high-quality quasar observations based on better uv coverage (Pushkarev & Kovalev 2015).Radio quasars having a compact structure, will be measured with higher angular resolution, and lower statistical and systematic uncertainty.On the other hand, we also look forward to a large amount of future data, not only from the radio quasars, but also from other astronomical probes, such as SNe Ia and BAO, allowing us to further improve the precision of H 0 and Ω k measurements.The results we obtained are very encouraging regarding the method proposed, which is alternative and complementary to CMB analysis.In light of current tensions in cosmology when local probes suggest different values of important parameters like the Hubble constant, from those obtained by CMB analysis, every cosmological probe alternative to already existing is valuable.

2. 1 .
The acceleration parameter from gravitational waveform detected by DECIGO In an expanding Universe, infinitesimal time intervals dt o measured in the observers frame are related to the analogous time intervals dt e in source (emitter) frame by the following (well known) relation: ) = 1 + z.This relation underlies the notion of the cosmological redshift.The next order derivative reads:

Figure 1 .
Figure 1.Simulated acceleration parameters X(z) based on the future DECIGO and the corresponding uncertainties σ X are shown in grey.Inset figure shows the effect of binning.Red bars represent weighted mean of X(z) in redshift bins of width ∆z = 0.1 with the corresponding uncertainties.Black solid line with corresponding blue uncertainty band display the ANN reconstructed X(z) relation.The advantage of binning is noticeable.

Figure 2 .
Figure 2. Simulated sample of 1000 angular size measurements of ultra-compact structure in radio quasars from the future VLBI observations.

Figure 3 .
Figure3.The 2-D plots and 1-D marginalized distributions with 1-σ and 2-σ contours of cosmological parameters (H 0 and Ω k ) from the simulated QSO and GW data.Dashed lines represent the fiducial parameters that we used in our simulation, i.e.H 0 = 67.4km s −1 Mpc −1 and Ω k = 0.

Figure 4 .
Figure 4.The precision on H 0 as a function of the number of simulated QSOs.The blue polyline represents the case of free Ω k parameter.The red polyline represents the case of fixed Ω k = 0 parameter.For the illustration of the uncertainty floor observed, standard statistical 1/ √ N curves have been overplotted.

Figure 5 .
Figure5.The precision on Ω k as a function of the the number of the simulated QSOs.The blue polyline represents the free H 0 parameter case.The red polyline represents the fixed H 0 = 67.4km s −1 Mpc −1 parameter case.For the illustration of the uncertainty floor observed, standard statistical 1/ √ N curves have been overplotted.

Figure 6 .
Figure6.1D and 2D marginalized probability distributions of H 0 and Ω k as in Fig.3but corresponding to doubling the uncertainties from GW or QSO data respectively, while keeping the other unchanged.Dashed lines represent the fiducial parameters that we used in our simulation, i.e.H 0 = 67.4km s −1 Mpc −1 and Ω k = 0.