Revisiting the cosmic distance duality relation with machine learning reconstruction methods: the combination of HII galaxies and ultra-compact radio quasars

In order to measure luminosity distances in the Universe, we always turn to sources that have known (or standardised) intrinsic luminosity, such as type-Ia supernova (SN Ia) 9 , more distant quasars 13 ; 14 ; 15 ; 16 , and gamma-ray bursts (GRB) 17 , etc. In addition, the HII galaxies and extragalactic HII regions 18 ; 19 ; 20 constitute a large fraction of population that can be observed up to very high redshifts, beyond the feasible limits of supernova studies. It is well known that the luminosity $L($H$\beta)$ in H$\beta$ and the ionized gas velocity dispersion $\sigma$ of HII galaxies and extragalactic HII regions may have a quantitative relation (be known as “$L$–$\sigma$” relation). The physics behind this relation is based only on a simple idea, i.e., as the mass of the starburst component increases, the number of ionized photons and the turbulent velocity of the gas may both increase as well. Melnick et al. first found that the scatter of “$L$–$\sigma$” relation is very small and have the capability to determine the cosmological distance independent of redshift21 . More specifically, based on the measured flux density (or luminosity) and the turbulent velocity of the gas, one can infer the luminosity distance directly. Whereafter, the validity of the “$L$–$\sigma$” relation acting as the standard candle and its possible cosmological applications have been extensively discussed in the literatures 22 ; 23 ; 24 .

The “$L$–$\sigma$” relation between the luminosity $L($H$\beta)$ in H$\beta$ of a source and its ionized gas velocity dispersion can be expressed as 19

where $\alpha$ is the slope and $\kappa$ is the intercept. The $L(\textrm{H}\beta)$ is obtained from the reddening corrected flux density $F(\textrm{H}\beta)$ which only bases on a general equation $L(\textrm{H}\beta)=4\pi D^{2}_{L}F(\textrm{H}\beta)$. Thus, the equation above can be written as a relation of the observed flux density

Analogous to SN Ia applied in cosmology, the $\alpha$ and $\kappa$ parameters should also be optimized with the assumed cosmological model parameters. Fortunately, Wu et al. used the measurements of Hubble parameters from cosmic clocks (model-independent) to calibrate $\alpha$ and $\kappa$, and demonstrated that the calibrated values $\alpha=5.12\pm 0.08$ and $\kappa=33.08\pm 0.13$ are reliable for cosmological applications 24 . In this work, we will adopt these values with their corresponding uncertainties to get the luminosity distance.

The catalog of spectral and astrometric data from HII galaxies and extragalactic HII regions contain more than 100 sources by far, and its statistical properties can be preliminarily considered in cosmology. In this work, we will use the current observations of 156 HII objects compiled by Terlevich et al. 19 which contain 25 high redshift HII galaxies sources, 107 local HII galaxies sources, and 24 extragalactic HII regions sources covering redshift range $0<z<2.33$. This dataset is larger than the source samples used by Plionis et al. 22 and is more complete than the high redshift data used by Melnick et al. 21 . Full information (including name of the source, redshift, flux density, and turbulent velocity with corresponding observational uncertainties) about the sample of 156 HII regions can be found in Table 1 of the work 20 .

Quasars, among the most distant objects in the universe, have great potential as distance indicators. From an observational point of view, there are currently two types of quasar data that can be served as cosmological probes, i.e., the non-linear relation between the ultraviolet and X-ray fluxes of the quasar to construct the Hubble diagram 13 ; 14 ; 15 , and the angular size-distance relation of ultra-compact structures object in radio quasars (QSO) as the standard ruler of cosmology from the very-long-baseline interferometry (VLBI) observations 16 . The first type of quasar data provides the luminosity distance, but not the angular diameter distance, directly. Moreover, although the sample collected by 13 contained 1598 suitable quasars and redshift reaches to $z\sim 5.5$, the sample itself exhibits a large intrinsic dispersion. Take these factors into consideration, we will use the radio quasar sample to obtain the angular diameter distance information.

The angular size-distance relation in compact radio quasar for cosmological inference was first proposed by Kellermann et al. 25 , in which he tried to obtain the deceleration parameter with 79 compact radio sources from VLBI at 5 GHz. Whereafter, Gurvits 26 extended this method and attempted to investigate the dependence of characteristic size on luminosity and redshift based on 337 Active Galactic Nucleuses (AGNs) observed at 2.29 GHz 27 . In the subsequent analysis, the literature 26 adopted the modulus of visibility $\Gamma=S_{c}/S_{t}$ to redefine angular size of radio sources $\theta$, which can be expressed by $\theta(z)=2\sqrt{-\ln\Gamma\ln 2}/\pi B_{\theta}$, where $B_{\theta}$ is interferometer baseline measured in wavelengths, $S_{c}$ and $S_{t}$ are correlated flux density and total flux density, respectively. Based on a simple relation between the angle and distance with the intrinsic linear size $l_{m}$ of the compact structure in radio quasars, the angular size $\theta(z)$ can be written as

where $l_{m}=lL^{\beta}(1+z)^{n}$ describes the apparent distribution of radio brightness within the core, $l$ is the linear size scaling factor, $L$ is the intrinsic luminosity of the source, and $\beta$ and $n$ represent the possible dependence of the intrinsic linear size of the source on luminosity and redshift, respectively. With the gradually refined selection technique and observations, as well as the elimination of systematic errors caused by various aspects, Cao et al. 16 compiled milliarcsecond compact radio sample of 120 intermediate-luminosity $(10^{27}W/Hz<L<10^{28}W/Hz)$ quasars with reliable measurements of the angular size of the compact structure covering the redshift range $0.46<z<2.76$ from VLBI survey at 2.29 GHz. They showed that $l_{m}$ is independent of redshift and luminosity ($|\beta|\approx 10^{-4},|n|\approx 10^{-3}$), which suggests that it can be used to cosmological studies. However, the current problem is how to determine the value of $l_{m}$. In the subsequent analysis, the linear size, without pre-assuming a cosmological model, was determined to be $l_{m}=11.03\pm 0.25$ pc by Cao et al. based on the $D_{A}(z)$ reconstruction from $H(z)$ data obtained from cosmic chronometers16 . The calibrated intrinsic length and cosmological application of this sample had obtained stringent constraints on both the matter density parameter $\Omega_{m}$ and the Hubble constant $H_{0}$, which are consist with Planck 2018 observation 3 . The ultra-compact structures in radio quasars for exploring other cosmological models has been investigated in many literatures 28 ; 29 ; 30 ; Xu18 ; 31 ; 32 .Therefore, it is reasonable to ask whether the derived angular size depends on the intrinsic luminosity of the radio quasar, and consequently affecting testing the validity of CDDR. In fact, the derived angular size is obtained by a ratio of correlated and total flux densities, i.e., the modulus of visibility $\Gamma=S_{c}/S_{t}$. Therefore, from the perspective of observation, the intrinsic luminosity of the radio quasar does not affect the effectiveness of CDDR testing.

From a theoretical perspective, one can directly achieve CDDR testing by combining the $L$–$\sigma$ relation in HII regions with the angular size-distance relation of compact radio sources. From the observational point of view, however, there is currently a lack of data samples. Not only of the HII region, but also samples of quasars. Although their redshifts cover each other well, there are very few of them meeting the same redshift at the same time. In order to achieve the CDDR testing and obtain convincing results, we consider two non-parameterized technologies, Gaussian Process (GP) and Artificial Neural Network (ANN), to reconstruct the HII galaxy Hubble diagram and ultra-compact structure of radio quasar sources data, respectively. There is no reason to favor one technology over another, but mutually consistent results for different parameterized technique would strengthen the robustness of the conclusion.

Gaussian Process.— The Gaussian Process (GP) is a random process defined in the continuous domain, which can be regarded as a set of all random variables in the continuous domain, and any single or multiple random variables satisfy one-dimensional Gaussian distribution or multi-dimensional Gaussian distribution 33 . The GP can be determined by a mean function $m(\textbf{x})$ and a covariance function $k(\textbf{x},\textbf{x}_{\ast})$ (also called kernel function). The GP defines a priori function, one can assume that for a given x, y there follows a distribution $p(\textbf{y}|\textbf{x})=\mathcal{N}(\textbf{y}|m,\textbf{K})$, which $\textbf{K}=k(\textbf{x},\textbf{x}_{\ast})$. The purpose of GP is to learn a mapping function $f$ from x to y through x, y. Then, for the given new $\textbf{x}_{\ast}$, one can predict $\textbf{y}_{\ast}=f(\textbf{x}_{\ast})$. According to priori distribution of GP, the joint distribution $p(\textbf{y},\textbf{y}_{\ast}|\textbf{x},\textbf{x}_{\ast})$ of observed data y and forecast data $\textbf{y}_{\ast}$ are given by 34

where $\textbf{K}_{\textbf{y}}=k(\textbf{x},\textbf{x})$, $\textbf{K}_{\ast}=k(\textbf{x},\textbf{x}_{\ast})$ and $\textbf{K}_{\ast\ast}=k(\textbf{x}_{\ast},\textbf{x}_{\ast})$. According to the joint distribution of y and $y_{\ast}$, one can get the conditional distribution

where $\mu_{\ast}=\textbf{K}_{\ast}^{T}\textbf{K}_{\textbf{y}}^{-1}\textbf{y}$ and $\Sigma_{\ast}=\textbf{K}_{\ast\ast}-\textbf{K}_{\ast}^{T}\textbf{K}_{\textbf{y}}^{-1}\textbf{K}_{\ast}$ are the expectation vectors and covariance matrices of the posterior prediction distribution, respectively.

The kernel function has many choices, such as square exponential function. We take the Matérn ($\nu=9/2$) covariance function here, because it can provide more reliable results when using GP to reconstruct function 35

where $\ell$ denotes the characteristic length scale in $x$-direction and $\sigma_{f}$ is the signal variance in $y$-direction. One can see that, except for the hyper parameters in kernel function, there is no parameter estimation that was involved in the final prediction of the $\textbf{y}_{\ast}=f(\textbf{x}_{\ast})$. It should be emphasized here that whenever one performs the GP regression, the hyperparameters should be optimized by GP with the observed data set, along with other parameters of interest (cosmological parameters or not). Therefore, the proper way of performing such GP analysis is to treat the GP hyperparameters on the same footing as the cosmological parameters, i.e. varying all relevant parameters together and sampling the joint posterior. Such procedure, which guarantees that the reconstructed function is independent of the initial hyperparameter settings, has been extensively applied in different cosmological studies 36 ; 37 ; 38 ; 39 , especially high-fidelity constraints on the spatial curvature parameter and Hubble constant Colgain21 ; Dhawan21 .

In this analysis, we use Gaussian Processes in Python (GaPP) ¹¹1http://www.acgc.uct.ac.za/seikel/GAPP/index.html to realize the reconstruction of different functions. For the HII regions sample, the reconstructed logarithmic luminosity distance $\log D_{L,HII}(z)$, as a function of logarithmic redshift $\log z$, with the estimation of their $1\sigma$ confidence regions are shown in the top panel of Fig. 1. Meanwhile, we also reconstruct function $\theta(z)$ from compact radio sources observations, and final results with their corresponding $1\sigma$ uncertainties are shown in the bottom panel of Fig. 1. We reconstruct 1000 points for HII regions and compact radio sources, respectively. From this figure, one can see that the $1\sigma$ uncertainty from GP reconstruction is smaller than that of individual data points. Such issue has been extensively discussed in the recent works 33a ; 33b . More specifically, the final reconstructed confidence region depends on three factors, i.e., the observed errors of data, the optimization of hyper parameters of the GP method, and the product of the covariance matrixes $\mathbf{K}_{\ast}\mathbf{K_{y}}^{-1}\mathbf{K}_{\ast}^{T}$ between the predicted and current observed points. It should be noted that, if $\mathbf{K}_{\ast}\mathbf{K_{y}}^{-1}\mathbf{K}_{\ast}^{T}>\sigma_{f}$, the uncertainty of the predicted point will be less than uncertainties of observed points when there is a large correlation between the data. One can clearly see from Eq. (6) that the correlation between $z$ and $\tilde{z}$ will be large when $z-\tilde{z}$ less than $\ell$. Such condition, which is satisfied by most of the HII galaxy and quasar data points in our study, will result in smaller 1$\sigma$ confidence region from GP. We refer the reader to Ref. 33 for further details on this issue.

Artificial Neural Network.— For the second nonparametric approach, we turn to Artificial Neural Network (ANN) and show the reconstructed function from observational data. With the development of computer hardware in the recent ten years, machine learning technology has been gradually applied to many research fields in astronomy, and shown excellent potential for solving cosmological problems, such as analyzing gravitational waves 40 ; 41 and constraining cosmological parameters 42 ; 43 ; 44 ; 45 ; 46 ; 47 .

The main purpose of an ANN is to construct an approximate function or map that correlates the input data with the output data. The ANN has been shown to be ”universal approximator” that can represent a wide variety of functions 48 ; 49 . The ANN is made up of neurons, which are very simple elements that receive digital input. Generally speaking, the artificial neural network consists of an input layer, one or more hidden layers and an output layer. Each layer takes a vector from the previous layer as input, applies a linear transformation and a nonlinear activation function to the input, and propagates the current result to the next layer. Formally, in a vectorized way 50

where $\mathbf{x_{i}}$ is the input vector at the $i$th layer, $W_{i+1}$ and $\mathbf{b_{i+1}}$ are linear weights matrix and the offset vector which need to be optimized, $\mathbf{z_{i+1}}$ is the output vector after linear transformation, and $f$ is the activation function. Here, the Exponential Linear Unit (ELU) is acted as the activation function 51 , which is given by

where $\alpha$ denotes the hyper-parameter that controls the value to which an ELU saturates for negative net inputs. Compared to other activation functions (such as the rectified linear and the leaky rectified linear), when the network exceeds five layers, ELU can not only improve the learning speed, but also have better generalization performance 51 .

The ANN equals to a function $f_{W,\mathbf{b}}(\mathbf{x})$. The goal of ANN is to make its output to be as close as possible to the target value $\mathbf{y}$. Then, according to the difference between the predicted value $f_{W,\mathbf{b}}(\mathbf{x})$ of the current network and the target value $\mathbf{y}$, the weight matrix of each layer needs to be constantly updated for minimize the difference, which is defined by a loss function $\mathcal{L}$. The method used is gradient descent, that is, by constantly moving the loss value to the opposite direction of the current corresponding gradient to reduce the loss value. Formally, in a vectorized way 52

where the operator $\partial$ denotes partial derivatives, and $f^{\prime}$ is the derivative for the nonlinear function $f$.

According to the publicly released code by the work 50 , which explicitly describe the ANN method, we use the module called Reconstruct Functions with ANN (ReFANN) ²²2https://github.com/Guo-Jian-Wang/refann to perform the reconstruction of the HII regions and radio quasars data-sets. Similarly, we show the reconstructed function $\log D_{L,HII}(z)$ by using ANN method as a function of logarithmic redshift $\log z$ with the estimation of $1\sigma$ confidence region in the top panel of Fig. 2. For compact radio sources, the reconstructed function $\theta(z)$ with corresponding $1\sigma$ uncertainties by using ANN is given in the bottom panel of Fig. 2. Similarly, we also reconstruct 1000 points for HII regions and compact radio sources, respectively. Compared to GP technology, the ANN method does not assume random variables that satisfy the Gaussian distribution, which is a completely data driven approach. It is interestingly to note that the uncertainties of the data reconstructed by ANN are almost equal to that of the observations. Therefore, the $1\sigma$ confidence region reconstructed by ANN can be considered as the average level of observational error. We refer the reader to Ref. 46 for further details on this issue.

In order to test the validity of CDDR, we use the reconstructed data to directly perform a model-independent test of CDDR, i.e., we do not adopt any parameterized form to quantify the CDDR, which is given by the following form 53 ; 54

Note that any statistically significant deviation from $\eta(z)=1$ could indicate possible violation of the three basic CDDR assumptions. Furthermore, we turn to four parameterized forms of CDDR, which have been extensively discussed in the quoted papers Cao11b

Note that the first parameterized form is independent of redshift, therefore we can compare it to other parameterized forms to check its possible dependency on the redshift. In general, $\eta(z)$ can be treated as parameterized functions of the redshift. In this work, we also use other general parametric representations for a possible redshift dependence of CDDR including two one-parameter expressions and a two-parameter parametrization. The consistency of the results under different parameterized forms will enhance the robustness of the conclusion.

From the observational perspective, we obtain the luminosity distance $D_{L}(z)$ through the ”$L$–$\sigma$” of HII galaxies and extragalactic HII regions, and the diameter distance $D_{A}(z)$ can be derived from compact structure in radio quasars. For a given $D_{L,HII}$ data point, the angular diameter distance $D_{A,QSO}$ should be observed at the same redshift. To avoid introducing additional systematic errors, a cosmological model-independent selection criterion is considered. We take $|z_{HII}-z_{QSO}|<0.005$ in our analysis 6 ; 7 . If one only consider the actual observational sample, it is difficult to achieve a rigorous CDDR test and get convincing results. That is the reason why we used two non-parameterized techniques mentioned above to reconstruct the data. Using the reconstructed samples, we are able to have a one-to-one matching between HII regions and compact structure in radio quasars . After executing the redshift selection criterion by using reconstructed data samples, 892 data are remained. Subsequently, the observed $\eta_{obs}(z)$ can be represented by following form

Uncertainties have been assessed from the standard uncertainty propagation formula, based on (uncorrelated) uncertainties of observable quantities. The total uncertainty budget includes the ionized gas velocity dispersion $\sigma_{\log\sigma(\textrm{H}\beta)}$, flux density $\sigma_{\log F(\textrm{H}\beta)}$ and additional systematic errors introduced from the calibrations of $\alpha$ and $\kappa$ in HII regions data, the angular size $\sigma_{\theta}$, and additional systematic errors introduced in the calibrations of linear size $l_{m}$ in radio quasars. So the total uncertainty of $\eta_{obs}(z)$ can be expressed as

In order to determine the best fitting CDDR parameters and corresponding uncertainties, we use the Bayesian statistical methods to obtain the posterior probability density function of the CDDR parameters $\eta_{j}$ ($j=0,..,4$) corresponding to four parametrization forms of Eq. (12). The posterior probability density function is given by

where $\mathcal{L}$ is the likelihood function, and has a following form of

and $p(\eta)$ is the prior, and assumed the following uniform distribution: $p(\eta_{j})=U[-1,1]$. We use the Python module $emcee$ 55 to perform the Markove Chain Monte Carlo (MCMC) analysis.