Measuring the speed of light with ultracompact radio quasars
Abstract
In this paper, based on a 2.29 GHz VLBI allsky survey of 613 milliarcsecond ultracompact radio sources with , we describe a method of identifying the subsample which can serve as individual standard rulers in cosmology. If the linear size of the compact structure is assumed to depend on source luminosity and redshift as , only intermediateluminosity quasars ( W/Hz W/Hz) show negligible dependence (, ), and thus represent a population of such rulers with fixed characteristic length pc. With a sample of 120 such sources covering the redshift range , we confirm the existence of dark energy in the Universe with high significance under the assumption of a flat universe, and obtain stringent constraints on both the matter density and the Hubble constant km sec Mpc. Finally, with the angular diameter distances measured for quasars extending to high redshifts (), we reconstruct the function using the technique of Gaussian processes. This allows us to identify the redshift corresponding to the maximum of the function: and the corresponding angular diameter distance Mpc. Similar reconstruction of the expansion rate function based on the data from cosmic chronometers and BAO gives us km sec Mpc. These measurements are used to estimate the speed of light: km/s. This is the first measurement of the speed of light in a cosmological setting referring to the distant past.
1 Introduction
Advances in cosmology over recent decades have been accompanied by intensive searches for reliable standard rulers in the Universe. Recently, attention of has been focused on large comoving length scales revealed in the baryon acoustic oscillations (BAO). The BAO peak location is universally recognized as a fixed comoving ruler of about 105 Mpc (where is the Hubble constant expressed in units of ). However, the socalled fitting problem (Ellis & Stoeger, 1987) still remains a challenge for BAO peak location as a standard ruler. In particular, the environmental dependence of the BAO location has recently been detected by Roukema et al. (2015, 2016). Moreover, Ding et al. (2015) and Zheng et al. (2016) pointed out a noticeable systematic difference between measurements based on BAO and those obtained with differential aging techniques. Recently efforts have also been made to explore the sizes of galaxy clusters at different redshifts, by using radio observations of the SunyaevZeldovich effect together with Xray emission (De Filippis et al., 2005; Bonamente et al., 2006). In similar spirit, powerful radio sources have served as a special population to test the redshift  angular size relation for extended FRIIb galaxies (Daly & Djorgovski, 2003), radio galaxies (Guerra, Daly & Wan, 2000), and radio loud quasars (Buchalter et al., 1998). More promising candidates in this context are ultracompact structures in radio sources (especially in quasars that can be observed up to very high redshifts), with milliarcsecond angular sizes measured by the verylongbaseline interferometry (VLBI) (Kellermann, 1993; Gurvits, 1994); for details regarding the respective definitions of angular size see Data. Possible cosmological applications of such compact radio sources as standard rulers have been extensively discussed in the literature (Jackson & Dodgson, 1997; Vishwakarma, 2001; Lima & Alcaniz, 2002; Zhu & Fujimoto, 2002; Chen & Ratra, 2003).
The idealized picture of standard rulers is that of a population of objects whose linear sizes are the same and either constant or change with redshift in a wellknown way. Of course, real objects like radio sources do not meet such stringent criteria. More realistically we hope for a population whose mean linear size satisfies these criteria. The extended radio sources whose linear extent is kpc are expected to depend strongly on the redshift evolution of the intergalactic medium. Compact radio sources depend more on the properties of the central engine (accreting central black hole) controlled by a limited number of parameters (like mass, accretion rate, magnetic field), and the ages of their jets are short compared with the age of the Universe. Consequently it is reasonable to suppose that such sources are much less affected by evolutionary processes. However, they are inherently a mixed population of AGNs (quasars, BL Lac objects etc.). There were suggestions (Gurvits, Kellerman & Frey, 1999; Vishwakarma, 2001) that the exclusion of sources with extreme spectral indices and low luminosities might alleviate the dependence of on the source luminosity and redshift. More recently, Cao et al. (2015) reexamined the same data in the framework of the CDM cosmological model assuming parametrized dependence of on luminosity and redshift (see Eq.(2) below). The general result was that whereas the evolution of with redshift is small, its dependence on luminosity can be substantial. We have divided the full sample into various subsamples according to optical counterpart and luminosity, with view to finding which subsample performs best in the role of standard rulers; quantitative details differ between the subsamples. In this paper we restrict our study to quasars and consider three luminosity ranges separately. Full details of this procedure will be the subject of a separate paper.
Here we focus on a very specific cosmological application of standard rulers. Namely, we reconstruct the angular diameter distance as a function of redshift up to and identify the redshift corresponding to the maximum of the function. According to Salzano et al. (2015), this allows us to calculate the speed of light , which is the first empirical assessment of the speed of light at an epoch much earlier than the present time.
2 Observational data
Cosmological parameter estimation from the angularsize/redshift relation for ultracompact radio sources was first discussed by Kellermann (1993), using a sample of 82 sources with 5 GHz VLBI contour maps; these images show a compact core with outlying weaker components, angular size being defined as the distance between the core and the most distant component with peak brightness greater than or equal to 2% of that of the core. Gurvits, Kellerman & Frey (1999) extended this work using a much larger sample, 330 VLBI images, taken largely from the CaltechJodrell Bank 5 GHz survey (Taylor et al., 1996). Their main conclusion was that their 5 GHz VLBI data are consistent with standard FLRW cosmologies with and . Their work has formed the basis of many subsequent investigations, early examples being Lima & Alcaniz (2002); Zhu & Fujimoto (2002); Chen & Ratra (2003).
The data used in the present paper were derived from an ancient 2.29 GHz VLBI survey undertaken by Preston et al. (1985) (hereafter called P85). By employing a worldwide array of dishes forming an interferometric system with an effective baseline of about wavelengths, this survey succeeded in detecting interference fringes from 917 radio sources out of a list of 1398 candidates selected mainly from the Parkes survey (Bolton, 1979). The results of this survey were utilized initially to provide a very accurate VLBI celestial reference frame, improving precision by at least an order of magnitude, compared with earlier stellar frames. An additional expectation was that the catalog would be “used in statistical studies of radiosource properties and cosmological models” (Preston et al., 1985). It should be noted that P85 does not list angular sizes explicitly; however, total flux density and correlated flux density (fringe amplitude) are listed; the ratio of these two quantities is the visibility modulus , which defines a characteristic angular size
(1) 
where is the interferometer baseline, measured in wavelengths (Thompson, Moran & Swenson, 1986, p. 13; Gurvits, 1994). It is argued in Jackson (2004) that this size represents that of the core, rather than the angular distance between the latter and a distant weak component. In this way Gurvits (1994) was able to construct an angularsize/redshift diagram for 337 sources with known redshifts; those with were ignored, and using just the highredshift data (258 objects) found marginal support for a lowdensity CDM cosmological model, but considered only those with . Using exactly the same data set, Jackson & Dodgson (1997) extended this work to the full plane, and concluded that the then canonical flat CDM model is ruled out at the 98.5% confidence level, and that the best values are if the Universe is spatially flat, later refined to (Jackson, 2004). This work was extended further by Jackson & Jannetta (2006), who updated the P85/Gurvits (1994) sample with respect to redshift, to include a total of 613 objects with redshifts . This extended data set is the one used in the present paper. Following Gurvits, Kellerman & Frey (1999), we exclude sources with inverted spectra and those with steep spectra: we apply an inclusion criterion on spectral index (), leaving 240 sources with relatively flat spectra. Appropriate information about all the 240 sources can be extracted from the source data set (Jackson & Jannetta, 2006). For each object, the corresponding optical counterpart and spectral index between 2700 and 5000 MHz can be found in P85. In order to enhance accuracy, in all cases we have derived the values for these quantities from the original references. We note that this final sample, which covers the redshift range , contains 181 sources identified as quasars (). By restricting the range of spectral indices and concentrating on milliarcsecond radio structures in quasars, we are better able to restrict the dispersion in intrinsic size in our analysis (Kellermann, 1993; Gurvits, Kellerman & Frey, 1999).
If we observe a given millarcsecond source at different reception frequencies , it is wellknown that its observed size falls as the frequency increases, being proportional to (Marscher & Shaffer, 1980; O’Sullivan & Gabuzda, 2009). At first sight an expectation would be that angular sizes in corresponding angularsize/redshift diagrams should fall by an additional factor , because the restframe emitted frequency has to be greater than the received frequency by a factor of : , thus masking other cosmological effects. However, observations show behaviour which is compatible with conventional cosmologies. The reason is almost certainly that there is a selection effect which is operating in our favour in this context; we are dealing with an ensemble of objects which may be intrinsically similar in their respective restframes, but appear to be very different in our frame. The underlying population consists of compact symmetric objects (Wilkinson et al., 1994), each comprising a central lowluminosity nucleus straddled by two oppositelydirected jets. Ultracompact objects are identified as cases in which the jets are moving relativistically and are close to the line of sight, when Doppler boosting allows just that component which is moving towards the observer to be observed, giving rise to the corejet structure observed in typical VLBI images, see for example Pushkarev & Kovalev (2012); the core is believed to be the base of the jet, rather than the nucleus (Blandford & Königl, 1979). Those jets which are closest to the line of sight appear to be the brightest. The radio sample is fluxlimited, so that as increases a larger Doppler boosting factor D is require (Dabrowski et al., 1995); it turns out that the ratio is approximately fixed, so that the restframe emitted frequency is also fixed. See Jackson (2004) for mathematical and astrophysical details. Note that this ratio is not necessarily unity, but the fact that it is fixed means that the interferometric angular sizes upon which this work is based are fit for purpose. This is another reason for ignoring sources with , which appear to be nonrelativistic. The above picture is based upon the unified model of active galactic nuclei and quasars (Antonucci, 1993, 2015).
3 Method
With respect to the milliarcsecond radiosources, our procedure follows that given in Sandage (1988). For a standard rod, the angular size at redshift can be written as
(2) 
where is the intrinsic metric linear size and represents the angular diameter distance. Moreover, we include two parameters and to respectively consider the “angular size  redshift” and “angular size  luminosity” relations. Our procedure follows the approach first proposed in Gurvits (1994) and later investigated in Gurvits, Kellerman & Frey (1999):
(3) 
where is the linear size scaling factor representing the apparent distribution of radio brightness within the core, is the luminosity normalized by WHz. In our analysis, it has been estimated from the measured flux density and the spectral index according to the formula: . In the latter we have used the socalled distance duality relation connecting luminosity distance with . This fundamental relation has triggered many studies in observational cosmology (Cao & Liang, 2011; Cao & Zhu, 2011; Cao et al., 2016). The parameter, which captures the dependence of the linear size on source luminosity, is highly dependent on the physics of compact radio emitting regions (Gurvits, Kellerman & Frey, 1999; Jackson, 2004). Besides cosmological evolution of the linear size with redshift, the parameter may also characterize the dependence of the linear size on the emitted frequency, as well as image blurring due to scattering in the propagation medium. Meanwhile, according to a recent analysis carried out by Koay & Macquart (2014), groundbased interferometers could hardly detect the scatter broadening of radio compact sources due to ionized intergalactic medium (IGM), which can only be detected by Space VLBI with baselines of 350,000 km at frequencies below 800 MHz. Considering that all VLBI images for our sample were observed at a frequency of 2.29 GHz, the effect of scattering in the propagation medium is not important in the present analysis (Gurvits, Kellerman & Frey, 1999).
It is obvious that, in combination with a reliable theoretical expression for , constraints on , and will help us to differentiate between subsamples and determine standard rods fulfilling the criteria of . To achieve this goal, we assume flat CDM model and use the cosmological parameters given by recent Planck Collaboration: and (Ade et al., 2014). Using routines available in the public CosmoMC package, we perform Monte Carlo simulations of the posterior likelihood for , , , in order to obtain their bestfit values, and determine the corresponding confidence regions for all three parameters characterizing compact radio structures. In computing we have assumed 10% uncertainties in the observed angular sizes, which in effect allows for both measurement errors and the intrinsic spread in linear sizes.
We should stress that, after identifying the best standard ruler population, the cosmologicalmodelindependent method will then be applied to determine the linear size of this standard rod. Next, using this value, the angular diameter distance function will be further investigated, searching for the redshift corresponding its maximum and the value of the speed of light at this redshift. It has been known for a long time (Weinberg, 1972) that theoretically, angular diameter distance will always ascend to a maximum value at a certain redshift and tend to decline at higher redshifts. However, it is very difficult to check this redshift evolution from astrophysical observations due to the socalled “redshift desert”. Moreover, the exact value for is highly dependent on the cosmological model adopted. In the framework of a flat Friedmann universe with no cosmological constant, the maximum redshift is fixed at . Numerical simulations by Salzano et al. (2015), who generated random cosmological models using the recent Planck+WMAP+highL+BAO best fitted cosmological parameters and their uncertainties, demonstrated that 95% of cases cover the range . In the present paper, considering the redshift range of our sample, especially that of the highredshift quasars we are able to trace the evolution of better and hence obtain more stringent constraint on . This is particularly interesting and important because the necessary condition for the extremum in the flat Universe leads to the formula for the speed of light:
(4) 
where is the Hubble parameter at redshift (Salzano et al., 2015). Therefore, having the real data on and one would be able to estimate the value of using extragalactic objects. This is the first such determination, since the previous papers on this subject (Salzano et al., 2015, 2016), although discussing the idea, referred only to simulated BAO data representative of future surveys like SKA or Euclid.
Quasar subsamples (Data number) 
[pc]  
Lowluminosity quasars (N=30) 

W/Hz  
Highluminosity quasars (N=31) 

W/Hz  
Intermediateluminosity quasars (N=120) 

W/Hz W/Hz  
Lowluminosity and highredshift quasars (N=13) 

W/Hz W/Hz;  
Highluminosity and highredshift quasars (N=31) 

W/Hz W/Hz;  
Combined High and Intermediateluminosity quasars (N=151) 

W/Hz  
Intermediateluminosity quasars (N=120) 

Cosmologyindependent method I  
Intermediateluminosity quasars (N=120) 

Cosmologyindependent method II  

4 Results and discussion
Let us start with the constraints on the parameters (,
, ) obtained for different samples corresponding to different
optical counterparts
Next, we consider only radio quasars, which choice is supported by
the following arguments. Firstly, quasars constitute the most
important part of our full radio sample. Secondly, the potential of
a quasar to act as a standard cosmological rod was already noticed
in Cao et al. (2015). However, the range of the parameter
suggests that the linear sizes of quasars are significantly affected
by their luminosity. Therefore, we further divide the full
quasar sample into three subsamples, according to their luminosity
. Quasars with luminosity smaller than W/Hz are
classified as relatively lowluminosity quasars, those with
luminosity larger than W/Hz are classified as relatively
highluminosity quasars, and intermediateluminosity quasars are
defined as having W/Hz W/Hz. The
intermediateluminosity range corresponds to the interquartile range
of luminosity distribution in our total sample of quasars. Since the
determination of luminosity involves cosmological model through
, we have crosschecked our sample selection by using two
different approaches: using a fiducial cosmology (flat CDM
model with parameters suggested by Planck) and
GPreconstructed from cosmic chronometers (passively
evolving galaxies). We present more details on cosmic chronometers
method later in this section. In both cases the sample contained the
same objects. So the sample selected by the luminosity seems not to
be affected much by the cosmological model assumed. Parameter
fitting results are summarized in Table 1 and displayed in
Fig. 1. The next issue is that our criterion seems uneven:
it cuts only one order of magnitude ( W/Hz
W/Hz) for intermediateluminosity quasars, allowing the other
classes – high and lowluminosity quasars – to vary over many
orders of magnitude. In order to check for this, we also made
similar fits to high and lowluminosity samples limited to one order
of magnitude. Strictly speaking, we included in our analysis two
following samples: 1) lowluminosity and highredshift quasars with
W/Hz W/Hz and ; 2) highluminosity
and highredshift quasars with W/Hz W/Hz
and . The results shown in Table 1, tend to support the
conclusion that only intermediateluminosity quasars exhibit no size
evolution and could be classified as standard rulers.
Fig. 2 shows the scatter plot of the three quasar
subsamples in the (, , ) space of observable
quantities. One can see that while lowluminosity and
intermediateluminosity samples are well separated, intermediate and
highluminosity ones are partly mixed. This motivated us to perform
fits for combined intermediate plus highluminosity samples with the
results displayed in respective rows of the Table 1. In conclusion:
one can see that quasars with intermediateluminosities meet the
requirement for a standard rod: , . We remark that the linear size is dependent on the
restframe frequency, thus the zero value of should be checked
with multifrequency VLBI data in the future. The bestfitted values
of and for this subpopulation are significantly
different from the corresponding quantities for other subsamples,
which supports the scheme of using a distinct strategy for treating
quasars with intermediate luminosities
The redshift of intermediateluminosity quasars ranges between and . In fact, the previous studies have given plausible reasons for ignoring lowredshift quasars with . The highredshift part () of the sample exhibits a smaller dispersion in luminosity , and the median angular size of quasars in this range is nearly independent of redshift (Gurvits, 1994; Jackson & Dodgson, 1997). In addition, there is no evidence of any selection bias concerning radio luminosity (Jackson, 2004) (weaker sources are distinctly smaller) when .
Another issue which might be raised regarding our selection of standard rods is the assumption of a particular cosmological model — CDM with the most recent parameters — in the course our estimation of , , and . Therefore, we have undertaken a similar analysis with different cosmological models and their best fitted parameters: WMAP9 observations in the case of CDM (Hinshaw et al., 2013), as well as Planck observations in the case of XCDM and CPL parameterizations (Ade et al., 2014). The results we obtained were very similar to those for the CDM model best fitted to Planck data. However, this approach still depends on the cosmological model assumed. Therefore, we also performed a model independent study. As is well known, assuming the FRW metric of a flat universe, the angulardiameter distance can be written as , where is the Hubble parameter. Following the recent works of Wei & Wu (2016) inspired by Gaussian processes (GPs) (Shafieloo et al., 2012; Seikel et al., 2012), we have reconstructed the function from the recent Hubble parameter measurements (Zheng et al., 2016) and then derived covering the redshift range of intermediateluminosity quasars. Such approach is independent of any specific cosmological model since the cosmic chronometers approach resulting in measurements is independent of any assumptions about cosmology. Moreover, in order to study the potential effect of the Hubble constant on the distance reconstruction with this method, two recent measurements of km (Bennett et al., 2014) and km (Riess et al., 2016) were adopted in the GP method. Distance reconstruction with the former prior on the Hubble constant is denoted as “Cosmologyindependent method I”, while “Cosmologyindependent method II” represents the distance reconstruction with the latter prior on the Hubble constant. Constraint results for the 120 intermediateluminosity quasars with the two modelindependent methods are shown in Fig. 3 and also in Table 1. They are consistent with those determined from the prior coming from the flat CDM model best fitted to Planck data, hence our results and discussions presented above are robust.
We invoked a specific cosmological model in order to identify the subsample which could then be used as standard rulers, and it turned out that intermediateluminosity quasars would serve this purpose the best. Therefore in order to calibrate further (i.e. determine more accurately) the linear size of the compact structure in the intermediateluminosity radio quasars we used the following method. Namely, the SN Ia are commonly accepted standard candles in the Universe and from their observed distance moduli we are able to recover luminosity distances covering the lower end of the quasar redshift range: . By virtue of the reciprocity theorem, they can be converted into angular diameter distances and used to calibrate the linear size of our milliarcsecond quasars. Applying the redshiftselection criterion, (according to (Cao & Liang, 2011)) to the Union2.1 compilation (Amanullah et al., 2010), we obtained 48 measurements of in the supernova/quasar overlap region, corresponding to the quasar redshifts. The intrinsic metric linear size distribution of the 48 intermediate luminosity quasars is also shown in the left panel of Fig. 4. Next we performed a similar fitting procedure, but now allowing the quantity to have only one free parameter, namely , so that the values of inferred from equation (3) match the supernova ones. As a result, we obtained the following:
(5) 
It should be noted that the 48 values of so derived are not used in any other context; specifically they are not the used directly in the preparation of the angulardiameterdistance/redshift diagram shown on Fig. 6, see discussion later. The likelihood function of is also shown in Fig. 4, from which one can see the perfect consistency between the result derived from the cosmologicalmodelindependent analysis and that determined from theoretical cosmological distances. In order to check the constraining power of quasars with the corrected linear size, using the “” relation for the full quasar sample, we were able to get stringent constraints on both the matter density parameter and the Hubble constant in the flat CDM cosmology, which are in reasonable agreement with those obtained from Planck observations. According to the unified classification of Active Galactic Nuclei (AGN), 10 pc is the typical radius at which AGN jets are apparently generated and there is almost no stellar contribution (Blandford & Rees, 1978). In the unified scheme this radius defines a compact opaque parsecscale core, which is located between the broadline region ( 1 pc) and narrowline region ( 100 pc) for QSOs. Recent AGN observations (Silverman et al., 2009) indicated that, within 10 pc around the central black hole, star formation rate is equal to the black hole’s mass accretion rate, which is consistent with the simulation results from Hopkins & Quataert (2010).
It is instructive to compare the estimates of angular sizes used here with those derived from more recent VLBI imaging observations based on better uvcoverage. Pushkarev & Kovalev (2015) presented the archival VLBI data of more than 3000 compact extragalactic radio sources (hereafter called P15) observed at different frequencies, GHz. Using these angularsize measurements of radio cores of active galactic nuclei (AGN) observed at 2 GHz, we can estimate the linear sizes of milliarcsecond sources covering redshift range , which provides a statistical value for the characteristic linear size, pc (=70 km sec Mpc). However, we should remark here that the P15 sample contains a wide class of extragalactic objects, including quasars (belonging to different luminosity categories), radio galaxies, and BL Lac objects, etc. Such measurement at 2 GHz gives a weighted mean size, taken over all the radio sources, rather the intermediateluminosity quasars. In fact, 58 intermediateluminosity quasars included in the modified P85 subsample used in the study, have also been observed by recent VLBI observations based on better uvcoverage and included in the P15 sample. Using these angularsize measurements of radio cores of intermediateluminosity quasars observed at 2 GHz, we estimate the statistical value for the characteristic linear size, pc, which is well consistent with our 2.3 GHz result derived from the P85 sample.
In the conical jet model proposed by Blandford & Königl (1979) [hereafter BK79], if we observe a given millarcsecond source at different frequencies , its observed size falls as the frequency increases, being proportional to (Marscher & Shaffer, 1980; O’Sullivan & Gabuzda, 2009). Therefore, VLBI results obtained at other frequencies on the same sources as used in the current study should be used for verification of the results. Based on the multifrequency angular size measurements from the P15 sample, we estimate characteristic linear sizes as pc at 5 GHz, pc at 8 GHz, pc at 15 GHz, pc at 24 GHz, and pc at 43 GHz. These figures are compatible with the expected behaviour based on the measurement from the P85 sample (See Fig. 5).
Next we attempted to construct an empirical relation extending to higher redshifts on the basis of individual quasar angular sizes, using Eq. (2); we obtained the measurement and corresponding uncertainty for each quasar. However, this procedure results in large uncertainties in – a problem which has been encountered previously (Gurvits, 1994; Gurvits, Kellerman & Frey, 1999), as can be seen from plots of the measured angular size against redshift therein. This problem remains even after 13 systems with very large () uncertainties are removed. In order to minimize its influence on our analysis, we have chosen to bin the remaining 107 data points and to examine the change in with redshift. The final sample was grouped into 20 redshift bins of width . Fig. 6 shows the median values of for each bin plotted against the central redshift of the bin. For comparison, the two curves plotted as solid lines represent theoretical expectations from the concordance flat CDM model and the Einstein  de Sitter model. One can see that the latter is disfavored at high confidence, more precisely with which corresponds to confidence level. More importantly, the angular diameter distance information obtained from quasars has helped us to bridge the “redshift desert” and extend our investigation of dark energy to much higher redshifts.
Finally, by combining the measurements of and the Hubble parameter at different redshifts, it should be possible to measure the speed of light . For this purpose one must measure the redshift at which reaches its maximum, then measure the angular diameter distance and the expansion rate at this particular redshift, and finally evaluate the speed of light through . Even though, from a purely theoretical point of view this idea is appealing (provided one is indeed able to measure and accurately and in the absolute sense) it turns out to be very difficult from the observational point of view. Therefore the papers (Salzano et al., 2015, 2016) by the proponents of this method refer only to simulated (fictitious) data. Our attempt is the first one performed on real data for this purpose.
Our quasar sample is sufficient to reconstruct the profile of up to the redshifts . However, large uncertainties of the measurements combined with the intrinsically large plateau about make it impossible to determine precisely. On the other hand, one can use a powerful reconstruction method (Shafieloo et al., 2012; Seikel et al., 2012) based on Gaussian Processes (GPs). We have employed these to reconstruct the function in order to find , and also to reconstruct function to find . The data on the function we used for the reconstruction, were acquired by means of two different techniques. The first was based on the so called cosmic chronometers (Jimenez & Loeb, 2002), i.e. the differential ages of massive, earlytype galaxies evolving passively on a timescale longer than their age difference. Recent analysis of the Baryon Oscillation Spectroscopic Survey (BOSS) Data Release by Moresco et al. (2016) enlarged the previous data set (Ding et al., 2015; Moresco et al., 2015) to a total number of 30 measurements of . The second technique relied on the analysis of BAO and provided us 6 further, precise measurements of at 6 different redshifts (see Table 1 of Zheng et al. (2016) for details and appropriate references to the original sources).
The function reconstructed from the binned quasar data is shown in Fig. 7 (upper left panel). From this reconstructed function we obtained the redshift at which reaches its maximum. The value is and the corresponding angular diameter distance Mpc. We have also enhanced the measurements with 38 galaxy clusters (Bonamente et al., 2006) and repeated the reconstruction of (see upper right panel of the Fig. 7). However, because of the redshift range of galaxy clusters (), their effect on the determination of turned out to be negligible. The lower left panel of the Fig. 7 shows the reconstructed function based on 36 Hubble parameter measurements. From this reconstructed relation we obtained Mpc/km/s at the maximum redshift. Therefore, our final assessment of the speed of light at is the following km/s. This is the first measurement of the speed of light in a cosmological setting referring to a distant past (at redshift the Universe was only Gyr old). However, the result is in perfect agreement with the value km/s measured “here and now”, which according to standard physics should be a universal constant of Nature. However, the agreement is only with our best fit central value, and the (inevitable) uncertainty in our result means that some window of opportunity is still open for exotic varyingspeedoflight scenarios. There are many ways in which such a constraint might be improved using our technique, which has paved one particular way in which one can derive the speed of light in the past. It should be stressed that our analysis and the assessment of were based on the radio source data from very old VLBI surveys. The prospect of future multifrequency VLBI surveys, comprising much more sources with higher sensitivity and angular resolution, especially at redshifts around where the redshift of maximum is very likely to be located, means that our constraints could be further improved by an order of magnitude. Moreover, one can imagine that a better reconstruction method might be devised, to process the raw radio quasar data at a preliminary stage. Hopefully the measurement of from quasars might be considerably improved by using future powerful probes such as BAO data from a survey like WFIRST (Spergel et al., 2013; Antonucci, 1993). Finally, we make a remark on the determination of the characteristic fixed length of our quasars. The distance moduli, used to obtain the 48 measurements of in the supernova/quasar overlap region, depend upon absolute magnitudes given in Amanullah et al. (2010) (Table 10), which in turn depend upon an assumed value of Hubble’s constant; the latter is given as km sec Mpc. Our values of scale as , as do the corresponding values of linear size. Hence our particular value should written as pc.
Acknowledgements
The authors would like to thank the referee for constructive comments, which allowed to improve the manuscript substantially. We are grateful to Zhengxiang Li, Yang Chen, and Meng Yao for helpful discussions. This work was supported by the Ministry of Science and Technology National Basic Science Program (Project 973) under Grants Nos. 2012CB821804 and 2014CB845806, the Strategic Priority Research Program “The Emergence of Cosmological Structure” of the Chinese Academy of Sciences (No. XDB09000000), the National Natural Science Foundation of China under Grants Nos. 11503001, 11373014 and 11073005, the Fundamental Research Funds for the Central Universities and Scientific Research Foundation of Beijing Normal University, China Postdoctoral Science Foundation under grant No. 2015T80052, and the Opening Project of Key Laboratory of Computational Astrophysics, National Astronomical Observatories, Chinese Academy of Sciences. Part of the research was conducted within the scope of the HECOLS International Associated Laboratory, supported in part by the Polish NCN grant DEC2013/08/M/ST9/00664  4. M.B. gratefully acknowledges this support. This research was also partly supported by the PolandChina Scientific & Technological Cooperation Committee Project No. 354. M.B. obtained approval of the foreign talent introducing project in China and gained special funding support from the foreign knowledge introducing project.
Footnotes
 In order to guarantee the appropriate size of each subsample, we decided to use, for the purpose of fitting, the full sample without the spectral index criterion.
 We are investigating the possibility that the distinction between mediumluminosity quasars and highluminosity ones has been overstated. There is a degree of degeneracy between the parameters and , which arises because the luminosity and redshift are correlated. Pending further investigations, we prefer to leave the 31 highluminosity quasars out of our considerations. Lowluminosity quasars are mostly at low redshifts, where the population is qualitatively different.
References
 Ade, P. A. R., et al. 2014, A&A, 571, A16
 Amanullah, R., et al. 2010, ApJ, 716, 712
 Antonucci, R. 1993, ARA&A, 31, 473
 Antonucci, R. 2015, arXiv:1501.02001
 Bennett, C. L., Larson, D., Weiland, J. L., & Hinshaw, G. 2014, ApJ, 794, 135
 Blandford, R. D. & Rees, M. J. 1978, “Some comments on radiation mechanisms in Lacertids, in BL Lac Objects” (ed. A.M.Wolfe), University of Pittsburgh, Pittsburgh, pp. 328¨C341
 Blandford, R. D., & Königl, A., 1979, ApJ, 232, 34
 Bolton, J. G., see Parkes Catalogue, 1990, Australia Telescope National Facility, Wright & Otrupcek, (Eds)
 Bonamente, M., et al. 2006, ApJ, 647, 25
 Buchalter, A., Helfand, D. J., Becker, R. H., & White, R. L., 1998, ApJ, 494, 503
 Cao, S., & Liang, N. 2011, RAA, 11, 1199
 Cao, S., & Zhu, Z.H. 2011, Science in China G: Physics and Astronomy, 54, 2260
 Cao, S., et al. 2015, ApJ, 806, 66
 Cao, S., Biesiada, M., Zheng, X., Zhu, Z.H., 2016, MNRAS, 457, 281
 Chen, G., & Ratra, B. 2003, ApJ, 582, 586
 Dabrowski, Y., Lasenby A., & Saunders R. 1995, MNRAS, 277, 753
 Daly, R. A., & Djorgovski, S. G. 2003, ApJ, 597, 9
 De Filippis, E., SerenoM, M., Bautz, W., Longo, G. 2005, ApJ, 625, 108
 Ding, X., Biesiada, M., Cao, S., Li, Z. X. & Zhu, Z.H., 2015, ApJL, 803, L22 [arXiv:1503.04923]
 Ellis, G. F. R. & Stoeger, W. 1987, ClassQuantGra, 4, 1697
 Guerra, E. J., Daly, R. A., & Wan, L. 2000, ApJ, 544, 659
 Gurvits, L. I., 1994, ApJ, 425, 442
 Gurvits, L. I., Kellerman, K. I., & Frey, S. 1999, A&A, 342, 378
 Hinshaw, G., et al. 2013, ApJS, 208, 19
 Hopkins, P. F., & Quataert, E. 2010, MNRAS, 407, 1529
 Jackson J. C., & Dodgson M. 1997, MNRAS, 285, 806
 Jackson, J. C. 2004, JCAP, 11, 7
 Jackson J. C., & Jannetta A. L. 2006, JCAP, 11, 002; see also http://nrl.northumbria.ac.uk/13109/
 Jimenez, R. & Loeb, A. 2002, ApJ, 573, 37
 Kellermann, K. I. 1993, Nature, 361, 134
 Koay, J. Y. & Macquart, J.P. arXiv:1410.7612v1
 Lima, J. A. S., & Alcaniz, J. S. 2002, ApJ, 566, 15
 Marscher A. P., & Shaffer D. B. 1980, AJ, 85, 668
 Moresco, M., 2015, MNRAS, 450, L16
 Moresco, M., Pozzetti, L., Cimatti, A., et al. 2016, arXiV:1601.01701
 O’Sullivan, S. P., & Gabuzda, D. C. 2009, MNRAS, 393, 429
 Preston R. A., Morabito D. D., Williams J. G., Faulkner J., Jauncey D. L., & Nicolson G. D., 1985, AJ, 90, 1599
 Pushkarev, A. B., & Kovalev, Y. Y. 2012, A&A, 544, 34
 Pushkarev, A. B., & Kovalev, Y. Y. 2015, MNRAS, 452, 4274
 Riess, A. G., Macri, L. M., Hoffmann, S. L., et al. 2016, ApJ, 826, 56
 Roukema, B., Buchert, T., Ostrowski, J. J. & France, M. J. 2015, MNRAS, 448, 1660
 Roukema, B., Buchert, T., Fuji, H. & Ostrowski J. J. 2016, MNRAS, 456, L45
 Salzano, V., Dabrowski, M. & Lazkoz, R. 2015, PRL, 114, 101304
 Salzano, V., Dabrowski, M. & Lazkoz, R. 2016, PRD, 93, 063521
 Sandage, A. 1988, ARA&A, 26, 561
 Seikel, M., Clarkson, C., & Smith, M. 2012, JCAP, 06, 036
 Shafieloo, A., Kim, A. G., & Linder, E.V. 2012, PRD, 85, 123530
 Silverman, J. D., et al. 2009, ApJ, 695, 171
 Spergel, D., et al. 206, WFIRSTAFTA Final Report, arXiv:1305.5422
 Taylor, G. B., Vermeulen, R. C., Readhead, A.C.S., et al. 1996, ApJS, 107, 37
 Thompson, A.R., Moran, J.M., & Swenson G.W. 1986, Interferometry and Synthesis in Radio Astronomy (John Wiley & Sons, New York)
 Vishwakarma, R. G. 2001, Classical Quantum Gravity, 18, 1159
 Wei, J. J. & Wu, X.F. 2016, arXiv:1611.00904v1
 Weinberg, S. Gravitation and Cosmology: Principles and Applications of the General Theory of Relativity (John Wiley & Sons, New York, 1972), p.423
 Wilkinson, P. N., et al. 1994, ApJ, 432, L87
 Zheng, X., Ding, X., Biesiada, M., Cao, S. & Zhu Z.H. 2016, ApJ, 825, 17
 Zhu, Z.H., & Fujimoto, M. K. 2002, ApJ, 581, 1