© The Authors 2025

1. Introduction

Directly measuring the expansion history of the Universe is among the most pressing tasks of observational cosmology, especially considering the evidence of its recent acceleration phase (Riess et al. 1998, 2021; Garnavich et al. 1998; Albrecht et al. 2006; Dawson et al. 2013, 2016; Martins 2024; Adame et al. 2025). Most studies aiming to do this include specific assumptions on an underlying model, or class of models, leading to cosmological constraints that are model-dependent. A conceptual alternative, not yet realised in practice but expected to be achieved by forthcoming astrophysical facilities, is the redshift drift measurement, also known as the Sandage test.

As was first explored by Sandage (1962), the evolution of the Hubble expansion causes the redshift of distant objects participating in the Hubble flow to change slowly with time. Just as the cosmological redshift, z, provides evidence of the expansion, the drift in this redshift provides evidence of the expansion's acceleration or deceleration between the epoch z and the present one. This implies that the expansion history of the Universe can be detected and mapped, at least in principle, by means of a straightforward spectroscopic monitoring campaign.

For simplicity, we take a generic metric theory of gravity with the further assumption of homogeneity and isotropy, in which case the evolution of the metric tensor is entirely specified by the global scale factor, a(t). A photon emitted by some object at comoving distance, χ, at the time, t_em, and observed by us at t_obs, will have a redshift of

$$1+z(t_{\mathrm{obs}},t_{\mathrm{em}})=\frac{a(t_{\mathrm{obs}})}{a(t_{\mathrm{em}})}.$$(1)

Now consider how the redshift of an object at a fixed comoving distance, χ, evolves with t_obs. Over a comparatively short timescale, Δt_obs, it will drift by

$$\frac{{\mathit{dz}}_{\chi }}{{\mathit{dt}}_{\mathrm{obs}}}\approx \frac{z_{\chi }(t_{\mathrm{obs}}+\mathrm{\Delta }{t}_{\mathrm{obs}})-z_{\chi }(t_{\mathrm{obs}})}{\mathrm{\Delta }{t}_{\mathrm{obs}}},$$(2)

which, differentiating Eq. (1) with respect to t_obs, leads to

$$\frac{{\mathit{dz}}_{\chi }}{{\mathit{dt}}_{\mathrm{obs}}}(t_{\mathrm{obs}})=\left[1+z_{\chi }(t_{\mathrm{obs}})\right]H(t_{\mathrm{obs}})-H(t_{\mathrm{em}}).$$(3)

Finally, we obtain, in a simplified form (McVittie 1962),

$$\dot{z}\equiv \frac{\mathit{dz}}{{\mathit{dt}}_{\mathrm{obs}}}=(1+z)\hspace{0.5em}{H}_0-H(z).$$(4)

Measuring ż for a number of objects at different z gives the function ȧ((z) and, given a(z), it is possible to reconstruct a(t). Therefore, a measurement of ż amounts to a purely dynamical and model-independent reconstruction of the expansion history of the Universe.

The effect is conceptually simple and provides a tool for key cosmological tests. Inter alia, one can easily verify that in decelerating universes, one expects a negative redshift drift at all redshifts. Conversely, a positive drift signal implies a violation of the strong energy condition, and hence the presence of some form of dark energy that accelerates the Universe (Liske et al. 2008; Uzan et al. 2008; Quercellini et al. 2012; Heinesen 2021). Mapping the redshift evolution of the drift can also lead to constraints on parameters in specific models, either alone or in combination with traditional cosmological probes (Alves et al. 2019; Rocha & Martins 2022; Marques et al. 2023, 2024).

The theoretical framework outlined above highlights the immense potential of redshift drift measurements to directly probe the expansion history of the Universe in a cosmological model-independent manner. However, achieving the sensitivity required to detect this effect is an observational challenge due to its exceedingly small expected magnitude. For instance, over a timescale of ∼10 years, the redshift drift, ż, in a Λ Cold Dark Matter cosmology (ΛCDM, Planck Collaboration VI 2020) is typically of the order of 10⁻¹⁰ to 10⁻⁹, or a few centimetres per second per decade in spectroscopic terms, for 0<z<6. Such precision necessitates both an adequate tracer population and a high-resolution spectroscopic capability that ensures the systematic uncertainties are smaller than the signal itself.

Among the possible tracers for this effect, the Lyman-α forest stands out as a particularly promising candidate (Loeb 1998). The Lyman-α forest consists of a series of absorption lines imprinted on the spectra of distant quasars due to intervening clouds of neutral hydrogen in the intergalactic medium (IGM). The suitability of the Lyman-α forest for redshift drift measurements arises from several key factors:

• Low cosmic overdensities: The Lyman-α forest arises from sparse intergalactic neutral hydrogen, which is physically unconnected with the background source, faithfully tracing the Hubble flow, apart from seldom strong absorbers due to dense intervening galactic environments (Rauch 1998).

• Abundance of tracers: The Lyman-α forest provides a dense sampling of absorption features along the line of sight (with number density dn/dz≈100−200 at z = 2−3, Kim et al. 2002), enabling the measurement of the redshift drift across numerous independent systems at varying redshifts.

• Spectroscopic accessibility: Even though the Lyman-α lines are not particularly narrow, with a typical linewidth of ∼30 km s⁻¹ (Kim et al. 2002), through the use of high-resolution spectroscopy their positions can be determined with sufficiently high precision to enable the measurement of the redshift drift. Moreover, the Lyman-α lines at z>1.5 are observable with large ground facilities.

• High-redshift sensitivity: At high redshifts, specifically deep in the matter era where the Universe is decelerating, the term (1+z) in Eq. (4) amplifies the magnitude of the drift, making the signal more detectable (see Fig. 2 in Liske et al. 2008).

Thus far, no detection of this signal is available, and existing upper limits based on other probes such as H I 21cm absorption lines (Darling 2012) and narrow metal absorption lines (Cooke 2020), have error bars three orders of magnitude larger than the expected signal. The first detections of the signal are expected to come from the ArmazoNes high Dispersion Echelle Spectrograph (ANDES) at the ESO's Extremely Large Telescope (ELT) (Marconi et al. 2024; Martins 2024) and by the full configuration of the Square Kilometre Array (SKA) (Klöckner et al. 2015; Kang et al. 2025), which are sensitive to the high redshift (matter era) and low redshift (acceleration era), respectively.

The Echelle SPectrograph for Rocky Exoplanets and Stable Spectroscopic Observations (ESPRESSO) spectrograph at the combined Coudé focus of ESO's VLT (Pepe et al. 2021) has a similar design to ANDES and enables us to bridge the gap between the present and the future. As part of the ESPRESSO Science Team's guaranteed time observations (GTO), we are carrying out the first dedicated redshift drift experiment, with a nominal experiment time of one year and an observation time of four nights (ca. 40 hours). This paper presents the first results, derived from 12h of observations of one quasar target.

Our two overarching goals are to improve the aforementioned upper limits and to start developing, testing, and validating a full end-to-end analysis pipeline, which will provide a baseline for the corresponding pipeline for the ANDES redshift drift measurement. This is the first of a series of reports on this experiment.

Specifically, in this work we set the first two epochs of observations with high-resolution, high-stability ESPRESSO spectroscopy and estimate the achievable precision. We check whether this aligns with previous estimates and if systematic effects of astrophysical or instrumental nature appear at the present level of signal-to-noise (S/N). We extrapolate our estimates to a realistic ANDES observational campaign to estimate the amount of time required to achieve a significant detection of the cosmological expansion.

In Sect. 2 we present the quasar targeted for our study. In Sect. 3 we describe the observation taken with ESPRESSO, the data reduction, the definition of the two observational epochs, and the recognition of metal lines and strong absorbers throughout the spectrum. The expected precision of a velocity drift measurement at the current level of S/N is presented in Sect. 4. In Sect. 5 we apply a direct pixel-by-pixel method, based on the work of Bouchy et al. (2001), to estimate the velocity shift, Δv, that has occurred between the two epochs. A different method of performing the measurement, based on modelling the Lyman-α forest, is presented in Sect. 6. In Sect. 7 we investigate the validity of our estimates and tackle the possible systematic effects arising in the measurement due to astrophysical or instrumental effects. In Sect. 8 we extrapolate our results to estimate the total observational time required to carry out a significant detection of the cosmological drift, either with ESPRESSO or with ANDES. Finally, we summarise and discuss our results in Sect. 9.

2. Target selection

The QUBRICS (QUasars as BRIght beacons for Cosmology in the Southern hemisphere) survey (Cristiani et al. 2023) provides a golden sample of seven bright quasars specifically selected to carry out the measurement of the redshift drift, observing their Lyman-α forest with an ESPRESSO- or ANDES-like spectrograph. For the present case study and the first epoch observation, we consider the brightest object from the golden sample, J052915.80-435152.0, hereafter called Super Bright 2 (SB2). In a parallel paper (Marques et al. in prep.) we shall discuss the observation of the second-most luminous object of the sample J212540.97-171951.4 (hereafter called SB1)¹. SB2 has the identifier QID 1128023 in the QUBRICS database (Calderone et al. 2019) with coordinates (J2000) RA 05:29:15.81 and Dec −43:51:52.1 and magnitudes r = 16.3005, i = 16.1184 (Onken et al. 2024) and G = 16.3452 (Gaia Collaboration 2021).

Besides being particularly suited for the Sandage test, SB2 is an outstanding object, the most luminous quasar known in the Universe. Medium resolution spectra taken with X-Shooter (Vernet et al. 2011) reveal that SB2 has a bolometric luminosity of log(L_bol/erg s⁻¹) = 48.27, supported by an accretion of 370 M_⊙ per year onto a supermassive black hole (SMBH) of mass log M/M_⊙∼10.24 or, in other terms, about one solar mass per day accreting onto a SMBH of ∼17 billion solar masses. From the same data, Wolf et al. (2024) report an emission redshift for SB2 of z_em = 3.962.

3. Data acquisition and treatment

In this work, we present high-resolution, high-S/N optical spectroscopy of this outstanding object obtained with ESPRESSO. ESPRESSO observations have been carried out using the single UT mode (Pepe et al. 2021) and achieve a resolving power of R=λ/Δλ∼135 000, where Δλ is the full width at half maximum of resolution element (∼2.2 km s⁻¹ in velocity space), with a 4x2 detector binning, chosen to optimise the S/N in the region covering the Lyman-α forest. A summary of the observations is given in Table 1.

The observations were reduced using the ESPRESSO Data Reduction Software² (DRS, Di Marcantonio et al. 2018), version 3.1.0. Laser Frequency Combs (LFC) frames were used for wavelength calibration of the spectra. We also produced spectra calibrated on the Fabry-Pérot (FP) etalon combined with a ThAr lamp frames to investigate the presence of possible instrumental systematics due to calibration stability and accuracy (see Sect. 7.5).

The final products of the pipeline were flat-fielded, blaze-corrected, sky-subtracted, and wavelength-calibrated. The spectra were also corrected to the Solar System barycentre using standard pipeline procedures, employing INPOP13c ephemeris (Fienga et al. 2014). Earth's precession is taken into account in the correction, consistently with the requirements for achieving sub-meter-per-second precision (Wright & Eastman 2014).

Order-by-order spectra were optimally extracted along detector rows (in the cross-dispersion direction), retaining for each pixel the original calibrated wavelength of the rows themselves, while the order-merged spectra were rebinned to a common wavelength grid.

The coaddition of the extracted spectra was performed using the ASTROCOOK software³ (Cupani et al. 2020), following the ‘drizzling’-like approach described in Cupani et al. (2016). A wavelength grid was defined with a fixed step of 1 km s⁻¹, and the pixels contributing to each bin in the grid were selected from the order-by-order spectra, to combine them all together in a single rebinning procedure. Contributing pixels were weighted by their associated error and by how much they overlapped with the wavelength range of the bins themselves. Order-by-order spectra were equalised prior to coaddition by multiplying their flux by scalars to get the same median flux level, to eliminate differences in flux levels due to varying observing conditions.

Coaddition of individual exposures resulted in a single spectrum, with a median S/N ∼ 86 per 1 km s⁻¹ pixel at the continuum level in the forest, and ∼105 redwards of the Lyman-α emission, as shown in Fig. 1. The Lyman-α forest highlighted in the figure is defined as the spectral interval between the Lyman-α and Lyman-β emission of the quasar, namely 509−603 μm, shown in more detail in Fig. 2. This section of the spectrum contains Lyman-α absorption lines with 3.18≤z≤3.96, and no Lyman-β absorptions.

Fig. 1. Combined spectrum of SB2. Top panel: Flux density in arbitrary units is shown in black, with the dashed purple line denoting the flux density error. Bottom panel: S/N per 1 km s−1 pixel. The shaded green area highlights the Lyman-α forest considered in the redshift drift measurement, bound by the Lyman-α and Lyman-β emissions of the quasar, namely between 509−603 m.

Fig. 2. Normalised transmitted flux of the Lyman-α forest in the combined spectrum of SB2 as a function of wavelength. The solid light blue line highlights the spectral regions masked due to metal absorption, the sub-DLA and bad pixels. The solid purple line reports the flux error.

Observations taken in the two ESO observing periods P110 and P112 (see Table 1), respectively, were also combined separately to provide two independent spectra that were used as the two epochs of the Sandage test. The first (P110) and second (P112) epoch spectra have a median S/N level at the continuum in the Lyman-α forest of ∼47 and ∼72, respectively. We defined the temporal baseline occurring between the two epochs’ spectra as the difference between the mean observational date of each period Δt = 0.875 yr.

3.1. Removing metal transitions from the Lyman-α forest

As has been shown by Cristiani et al. (2024), Lyman-α and metal lines have significantly different dynamical behaviours on the temporal and physical scales probed by the experiment. The authors show that absorption lines due to neutral hydrogen on adjacent sightlines (with sub-kiloparsec separations) appear identical within the noise, whereas metal lines show significant differences in their velocity structure. These intrinsic dynamical effects can be a source of strong systematic effects, inducing a velocity shift in the metal lines up to a few centimetres per second per year, of the same order of magnitude as the expected cosmic signal. We therefore removed the metal-polluted regions from our analysis to reduce possible systematics.

This process was done using the ASTROCOOK software (Cupani et al. 2020) to provide a secure identification of metal lines in the spectrum. First, the quasar continuum emission was estimated by applying an iterative sigma-clipping procedure to remove absorption features while masking sky emission lines, telluric absorption lines, or spurious, narrow spikes in flux (e.g. ‘cosmic rays’ that were not removed by the coadding process), and manually correcting the continuum level where necessary. Sky and telluric lines were masked using a reference SKYCORR model (Noll et al. 2014) and empirically defining proper cuts to conservatively remove even the line wings. We performed a search for metal absorption lines throughout the spectrum, starting from line doublets (e.g. C IV, Si IV and Mg II) redwards of the Lyman-α emission and following up with associated lines, i.e. other ionic transitions at the same redshift as the found doublets. Table 2 reports the absorption systems, i.e. the groups of ionic transitions with the same redshift, found along the spectrum of SB2.

Once identified, we fitted composite Voigt profiles to the metal absorption lines that fall redwards of the Lyman-α emission, estimating the absorber's redshift, column density, and Doppler broadening. Metal lines bluewards of the Lyman-α emission were not fitted, as this would require modelling multiple blended H I components, making the analysis significantly more complex. Since this is beyond the scope of the present study, we left it unaddressed and masked out such lines, removing the spectral regions in velocity space associated with higher wavelength transitions. A table reporting the fitted Voigt parameters for the metal transitions found redwards of the Lyman-α emission can be found in an electronic form at the CDS.

3.2. A sub-DLA at z_abs = 3.63

As is visible in Fig. 1, the spectrum of SB2 presents a strong H I absorption at λ∼563 μm or z∼3.63. We study this absorbing system in more detail to understand, based on its column density, whether it might belong to a dense environment prone to proper motions related to galactic dynamical processes, or it stems from the sparse cosmic web following the Hubble flow and thus can be considered in the redshift drift measurement. We found strong and complex associated metal absorption redwards of the Lyman-α emission, stemming from ions with both low (O I, Al II, C II, Fe II, Si II) and high (C IV, Si IV) ionisation states. We fitted the Lyman-α absorber within Astrocook using 3 components, tying the redshifts of H I and low ionisation species. The total H I column density of the system is log N_H I≈19.8, corresponding to a sub-Damped Lyman-α System (sub-DLA, Péroux et al. 2003).

Given this value, we decided to exclude the spectral range tied to this absorber from the analysis, as it is probably related to an intervening galactic system and could potentially induce systematics in the measurement due to its peculiar motions. We thus mask the region of the spectrum spanning ±2000 km s⁻¹ from the column density weighted average redshift position of the system. The final masked spectra used in our analysis have 41 737 pixels (1 km s⁻¹ wide) in the Lyman-α forest, compared to 50 805 pixels without mask. Meaning we consider only 82.15% of the total Lyman-α forest. Figure 2 shows the normalised flux of the Lyman-α forest in the spectrum of SB2, highlighting the masked-out regions.

4. Expected velocity precision

In a ΛCDM Universe (Planck Collaboration VI 2020), the expected cosmological redshift drift at redshift z = 3.573 (i.e. the middle of SB2's Lyman-α forest) is ż = −6.5 × 10⁻¹¹ yr⁻¹, equivalent to $\dot{v}=-0.43{\mathrm{cm}\mathrm{s}}^{-1}{\mathrm{yr}}^{-1}$ in velocity space. On the baseline of our experiment, the expected velocity shift between the two epoch spectra is therefore Δv=−0.38 cm s⁻¹.

Liske et al. (2008) provided a practical formula to assess the precision expected from an experiment performed on the ELT based on an ensemble of quasars at different redshifts. We adapted this formula to estimate the precision expected from our data. Given the SB2 ESPRESSO spectra at the two epochs, we estimated the predicted precision in the velocity shift measurement, adapting the scaling relation of Liske et al. (2008)

$${\sigma }_v=g\times 1.35{\left(\frac{\mathrm{S}/\mathrm{N}}{4075}\right)}^{-1}{\left(\frac{N_{\mathrm{QSO}}}{30}\right)}^{-0.5}{\left(\frac{1+z_{\mathrm{QSO}}}{5}\right)}^{-\gamma }{f}_{{\mathrm{Ly}}_{\alpha }}^{-0.5}{\text{cm s}}^{-1},$$(5)

where S/N = 86 is the median S/N at the continuum level per 1 km s⁻¹ pixel per object, accumulated over all observations (i.e. the S/N level of the combined spectrum). N_QSO is the number of quasars in the sample (N_QSO = 1 in our case), and z_QSO = 3.962 is the redshift of the quasar. The γ exponent is 1.7 for z_QSO≤4 and 0.9 above. The factor $f_{{\mathrm{Ly}}_{\alpha }}=0.8215$ is the fraction of the Lyman-α forest considered after masking. The factor g, called the ‘form factor’, depends on the observing strategy, being 1 if half of the exposures are taken at the beginning of the experiment and half at the end, and becoming larger if the measurements are spread out over time, reaching 1.7 for a uniform distribution (see Sect. 6 of Liske et al. 2008 for further details on the computation of the form factor). In our case, given the fact that 38% of the total integration time is collected at the first epoch, and 62% at the second, the form factor has a value of g = 1.03. From Eq. (5), we can estimate the expected maximum precision allowed by photon noise on our current data to be σ_v = 4.02 m s⁻¹.

Note that in this work we have analysed only one sightline, whereas the relation from which Eq. (5) is derived has been calibrated on a measurement based on an ensemble of spectra. Sightline-to-sightline deviations from the predicted measurement precision due to cosmic variance are expected due to the varying number of Lyman-α lines and strong absorbers found in the spectra of different objects.

5. Redshift drift measurement: Pixel-by-pixel method

We first measured Δv using a method developed for measuring radial velocity shifts between repeated observations of stars for exoplanet search (Bouchy et al. 2001), which we refer to as the pixel-by-pixel method. In summary, due to the drift of the absorption lines, the flux of the i-th pixel at the second epoch, F_2,i, can be expressed as a small perturbation of the flux in the same pixel at the first epoch, F_1,i:

$$F_{2,i}=F_{1,i}+\frac{{\mathit{dF}}_i}{d\lambda }\frac{\delta v_i}{c}{\lambda }_i,$$(6)

where λ_i is the observed wavelength of the i-th pixel, and dF_i/dλ is the slope of the spectrum at that pixel. The perturbation of the fluxes between two epochs defines, for each pixel, a small velocity shift, δv_i. Note that Eq. (6) holds only for δv_i much smaller than the typical line width (∼30 km s⁻¹). The single pixel's velocities can be averaged over the whole spectrum to estimate the velocity drift occurring between the two epochs.

$$\mathrm{\Delta }v=\frac{\sum \delta v_i{w}_i}{{\sum }_i{w}_i},$$(7)

where each pixel is weighted by $w_i={\sigma }_{v_i}^{-2}$, i.e. the inverse variance of δv_i. Accounting for the uncertainties in both epochs’ spectra and in their derivative, the variance becomes

$${\sigma }_{v_i}^2={\left[\frac{c}{{\lambda }_i({\mathit{dF}}_i/d\lambda )}\right]}^2\left[{\sigma }_{1,i}^2+{\sigma }_{2,i}^2+\frac{(F_{2,i}-F_{1,i}{)}^2}{({\mathit{dF}}_i/d\lambda {)}^2}{\sigma }_{F_i^{\prime }}^2\right],$$(8)

where σ_1,i and σ_2,i are the i-th pixel flux error in the two spectra, and ${\sigma }_{F_i^{\prime }}$ is the error on the flux slope at the i-th pixel. With this choice of weights, the uncertainty on the final Δv estimate is

$${\sigma }_v={\left[{\displaystyle \sum _i}{\sigma }_{v_i}^{-2}\right]}^{-1/2}.$$(9)

It is worth noting that, at first order, the flux slope term dF_i/dλ is epoch-invariant (under the assumption of very small shifts) and does not carry any epoch-dependent information. For this reason, we computed the flux derivative by means of finite differences of non-adjacent points, on the second epoch spectrum to exploit the higher S/N level and lower the noise in its evaluation.

We evaluated the velocity shift, δv, occurring between the two epochs, exploiting the normalised fluxes in the Lyman-α forest of the two spectra, where we have masked metal lines, the strong H I absorber at z∼3.63, and bad pixels. Our estimate yields Δv=−2.67±6.64 m s⁻¹. Taking into account the mean baseline between the two epochs of Δt = 0.875 years, we estimate a drift of $\dot{v}=-3.06\pm 7.59\mathrm{m}{\mathrm{s}}^{-1}{\mathrm{yr}}^{-1}$, equivalent to a redshift drift of ż = (−0.47 ± 1.15) × 10⁻⁷ yr⁻¹.

6. Redshift drift measurement: Model-based

The larger than expected uncertainty from the direct method is surprising and is not completely understood whether it stems from either instrumental systematics or problems in the measurement procedure, or whether it is intrinsic to our sightline. Loss of precision prolongs the time required for reliably detecting the redshift drift signal and raises the risk of critical failure. A secure estimate of the measurement uncertainty is therefore of crucial impact on the experiment timeline.

In this section, we develop a different approach to the measurement to investigate whether a higher precision can be achieved with the current data. Instead of directly comparing the two epochs’ spectra, we built an ensemble of analytical models of the forest of SB2 that are then correlated to the two epochs. The modelling procedure is calibrated on mock spectra extracted from state-of-the-art physical simulations of the IGM and is thus ‘informed’ on the spectral scales of the single Lyman-α lines and is able to de-prioritise pixels that might be affected by metal transitions or outlier pixels (e.g. non-removed cosmic rays). With this approach, we are capable of reducing the uncertainty in the measurement by a significant fraction and decreasing the total amount of integration time needed for the experiment.

Summarising, the main steps of the method are:

• Production of realistic mock sightlines that reproduce the properties of the ESPRESSO spectra of SB2 considered in the measurement.

• Calibration of the modelling procedure on the mock data.

• Construction of an ensemble of models of the Lyman-α forest of SB2.

• Evaluation of the velocity shift between the model ensemble and the first and second epoch data.

We stress that the mock data were only used to calibrate the modelling procedure and were not involved in the final measurement.

6.1. Mock spectra

We used realistic mock spectra that reproduce the properties of the three SB2 spectra (first epoch, second epoch, and combined) to calibrate the modelling pipeline and validate the measurement results. Starting from the hydrodynamical cosmological simulations of the IGM of the Sherwood Simulation Suite (Bolton et al. 2017) (with a box size of 10 Mpc/h and 2×1024³ dark matter and gas particles), we extracted the gas particle data and computed the H I optical depth along 5000 skewers piercing the box at every redshift step Δz = 0.1, within 3.1≤z≤4. The optical depth of the spectra was rescaled to match the mean optical depth observed at z∼3.57 (centre of the forest) from Becker & Bolton (2013). The sightlines (10 Mpc/h long) are randomly selected from the appropriate redshift interval, shifted so that they begin and end with regions of no absorption, and stitched together to form 100 spectra that match the length of the Lyman-α of SB2 (∼421 Mpc/h long, or 50 321 km s⁻¹).

The long spectra were then convolved with a Gaussian of FWHM ∼2.2 km s⁻¹ to mimic ESPRESSO's resolution and re-binned to a grid of pixels 1 km s⁻¹ wide (the average native pixel size extracted from the simulation was 0.593 km s⁻¹), matching the binning of our data. The noise properties of the original spectrum were added to the mocks following the prescription of Rorai et al. (2017) (see also Trost et al. 2025). Summarising, the original first epoch spectrum was divided into ten chunks, ∼10 nm wide. For each of these chunks, we:

• divided the pixels into 50 bins according to their normalised flux value, producing 50 corresponding noise distributions;

• assigned a noise value to each simulated pixel, where the value was randomly sampled from the distribution associated with the pixel's flux value, as defined in the point above.

This procedure was repeated for each chunk of the first epoch spectrum and its corresponding mock spectra, to preserve the relation between S/N at the continuum level and wavelength.

The second epoch spectra were generated in the same way from the same sightlines. However, before sightline extraction, we artificially added a velocity drift of −0.38 cm s⁻¹, i.e. the expected drift at z∼3.57 for a ΛCDM Universe during our baseline Δt = 0.875 yr, to the simulated particles’ velocities. These sightlines were then processed in the same way, and the noise of the second epoch was added to the flux, as explained above.

We combined the first and second epochs’ mock datasets via a weighted average, with weights proportional to the inverse of the pixels’ flux variance, to create the mock combined spectrum. Figure 3 shows the distribution of the S/N per pixel in the three mock datasets and the original spectra. To further validate the analysis (see Sect. 7), we also built second epoch mock datasets assuming either no drift between the two epochs or 10³ and 10⁴ years of temporal separation, adding a velocity drift to the simulation's particles of −3.8 m s⁻¹ and −37.6 m s⁻¹, respectively. On all spectra, we applied the same pixel masking of the original data, derived from metal, sub-DLA and bad pixels’ removal, to ensure the absence of systematics in the final measurement due to pixel masking and to conserve the total amount of pixels involved in the measurement.

Fig. 3. Distribution of the S/N per pixel in the Lyman-α forest of the three spectra of SB2 (dashed red) and of the corresponding mock spectra (solid black). Upper panel: First epoch. Middle panel: Second epoch. Bottom panel: Combined spectrum.

6.2. Model assessment

We relied on mock data to calibrate the modelling procedure before building the best possible model of SB2. Multiple sightlines were considered in the calibration to properly account for cosmic variance effects in the modelling procedure. Given the short timescales between the two epochs and that the total S/N is insufficient to detect the cosmic drift signal, we simplified the procedure by creating a fiducial model of the combined spectrum instead of modelling the two epoch spectra individually. A further simplification regards our decision to model the spectrum using cubic splines⁴ instead of the more traditional approach using Voigt profile decomposition. This was done because we were not interested in measuring physical properties of the gas but rather differences between the SB2's spectra in the two epochs, for which the spline method was sufficient.

In quasar spectra, Lyman-α lines have a typical width of ∼30 km s⁻¹ (Kim et al. 2002), while outlier features (e.g. cosmic rays) and metal transitions are typically narrower (down to a few kilometres per second). We calibrated the spline knot spacing on mock spectra containing only Lyman-α absorption, so that our procedure was ‘informed’ on the scales that we want to probe. Then, by generating an ensemble of models for SB2, we recovered a higher variance among models in the regions that might be affected by sharp non-Lyman-α absorption and could de-prioritise such features in the final measurement of Δv. The described procedure also prevents over-fitting of the noise.

In practice, let S^k(λ|n(A,ϕ)) be the spline model fitted on the k-th spectrum of the combined mock dataset described in the previous section, using the square inverse of the flux noise as weights in the fit, where n is a vector defining the positions of the spline knots along the spectrum. The knots are equispaced in velocity space, separated by a constant velocity, A, starting from an initial position (or phase) 0<ϕ<A with respect to the beginning of the Lyman-α forest.

Goodness of fit for the model was determined from its residuals to the data, $R^k(\lambda |A,\varphi )=\left(F-S^k(\lambda |\text{n}(A,\varphi ))\right)/{\sigma }_F$. We considered a model to be good if the residuals are centred at zero with a variance equal to unity, $R^k\sim \mathcal{N}(0,1)$.

One hundred mock spectra were fitted in this way, varying both the inter-nodal distance A and the initial phase ϕ, and their residuals were calculated as above. For each value of A and ϕ, we calculated the mean and the variance of the residuals for all 100 mock spectra. The optimal value of A was determined by examining which combination of parameters results in residuals most closely matching values drawn from $\mathcal{N}(0,1)$. Mathematically, we searched for a value of A that minimises:

$$\varepsilon (A)=\frac{{\left({\sigma }_{\mathrm{R}}(A)-1\right)}^2}{{\sigma }_{{\sigma }_{\mathrm{R}}}^2(A)},$$(10)

where σ_R(A) is the standard deviation of the fit residuals of splines with internodal distance A, averaged over all spectra and all considered initial node phases. ${\sigma }_{{\sigma }_{\mathrm{R}}}^2(A)$ is the variance of the quantity σ_R(A) due to sightline variance and variations in the initial node phase. Figure 4 shows σ_R−1 and its variance as a function of internodal distance A. Minimising Eq. (10) yields a best inter-nodal distance of $A=8.{37}_{-0.24}^{+0.26}\mathrm{km}{\mathrm{s}}^{-1}$. Note that this value of A is valid only for this precise S/N level, resolution, pixel size, and redshift. For any change in the spectral properties, A should be re-evaluated.

Fig. 4. Standard deviation of the fit residuals between the combined mock spectra's flux and the spline model with internodal distance A. The blue squares describe the σR of models on different mock spectra and randomly drawn initial node phase ϕ. The black error bars define the average and standard deviation of σR(A)−1. The vertical dashed purple line denotes the optimal value of A.

6.3. Building the forest model

Given the optimal knot spacing and its uncertainty, A and σ_A, we were able to construct a large number of models of the Lyman-α forest of SB2 that describe the data equally well in terms of a goodness-of-fit statistic (e.g. χ²). This allowed us to more thoroughly assess the impact of model non-uniqueness on the final result. This is conceptually similar to the approach taken by, for example, Lee et al. (2021), Milaković et al. (2024), Webb et al. (2025) in the context of fundamental constant and isotopic ratio measurements. Additionally, we identified parts of the spectrum where model non-uniqueness was minimal so that those regions could be prioritised in our analysis, whereas regions where model diversity (and hence uncertainty) was larger are considered ‘less trustworthy’. Moreover, we implemented a weighting scheme based on the model's gradient and its uncertainty, mimicking the pixel-by-pixel method (Sect. 5), to achieve a larger sensitivity to velocity shifts. Both aspects are major benefits of the model-based method, which results in a better handling of uncertainties.

We started by fitting the combined SB2 spectrum with N = 500 splines, S_j(λ | n_j). The node positions of the j-th model, n_j, were generated starting from an equispaced sequence of nodes in velocity space, n_j,0, with inter-nodal separation, A_j, drawn from a normal distribution, $\mathcal{N}(A,{\sigma }_A)$, and initial phase 0≤ϕ_j<A_j, drawn with uniform probability. The node positions were then perturbed by adding random noise in velocity space as n_j=n_j,0+m, where each element of the m array was randomly drawn from a normal distribution, $\mathcal{N}(0,{\sigma }_A)$.

From the ensemble of spline models, we evaluated the mean model, $\overline{S}$, and the variance between the models, ${\sigma }_S^2$.

Figure 5 shows a region of the first, second, and combined spectra, alongside the mean spline model. In the figure, flux noise and model variance are shown in the middle panel. The periodic pattern found in the flux noise results from the coaddition of several exposures and the subsequent rebinning to a fixed grid, which creates an aliasing effect in the error arrays, whose frequency depends on the position of the pixels within the spectral order. The bottom panel shows the pixel weights used in the measurement (see Eq. (12)), proportional to the local model derivative.

Fig. 5. Section of the SB2 spectrum. Top panel: Normalised flux of the first epoch (green), second epoch (orange), and combined (purple) spectra. The mean model, $\overline{S}$, is shown as a solid black line. The spectra are shifted vertically for clarity. Middle panel: Normalised flux errors, σF, of the three spectra shown in the panel above, with the same colour coding. The solid black line is the standard deviation of our model, ${\sigma }_{\overline{S}}$, i.e. it shows where our ensemble of models exhibits large variation, and hence uncertainty. Bottom panel: Normalised pixel weights, W, as defined in Eq. (12). Regions with higher values provide more constraints on our final ż measurement.

6.4. Velocity drift estimate

Strictly speaking, |ż| is larger at the red end of the Lyman-α forest than at its blue end, causing a redshift dependent compression of features in our spectrum. However, this effect is tiny and completely negligible at the S/N of our data. We therefore ignored the redshift dependence of ż and assumed that only a rigid translation Δv, in velocity space, occurred between the two epochs. Δv was measured indirectly, by first calculating the velocity shifts δv occurring between individual epochs and the model. The measured shift between two epochs is then: Δv=δv₂−δv₁, where subscripts identify the epoch.

The value of δv was obtained through a Monte Carlo Markov chain (MCMC) process with a modified Gaussian likelihood of the form:

$$ln\mathcal{L}(\delta v)=-0.5{\displaystyle \sum _i^{N_{\mathit{pix}}}}W({\lambda }_i|\delta v)\frac{{\left(F_i-\overline{S}({\lambda }_i|\delta v)\right)}^2}{{\sigma }_{F,i}^2+{\sigma }_{\overline{S}}^2({\lambda }_i|\delta v)},$$(11)

where the sum was performed over all non-masked pixels. F and σ_F are the normalised flux and its uncertainty, respectively, of either the first or second epoch spectra of SB2. N_pix is the number of non-masked pixels, $\overline{S}(\lambda |\delta v)$ is the mean spline model, and ${\sigma }_{\overline{S}}^2(\lambda |\delta v)$ is the variance within the ensemble of spline models, where the latter is included as an additional term at the denominator to effectively consider the pixels for which the ensemble's models are not consistent as less trustworthy. For both terms, the dependence on δv denotes the fact that the splines’ nodes have been solidly shifted by δv in velocity space, and the models are re-evaluated on the spectrum's pixel grid at each step of the MCMC process. The term W contains information on pixel weights.

Similarly to Bouchy et al. (2001), we wanted to assign a higher weight to pixels sensitive to velocity shifts, such as those with large flux gradients (i.e. the sides of the lines) and those for which variance in the flux gradient is minimal (i.e. for which models are consistent). The former aspect is similar to what is done in Sect. 5, but the latter aspect was only possible because we produced a large number of spline models. The two contributions were defined by evaluating the mean derivative (in absolute values) $\overline{S}\prime =\langle \left|{\mathit{dS}}_j/d\lambda \right|\rangle$ and the standard deviation within first derivatives’ magnitudes in the model ensemble σ_S′. The weights were therefore defined as

$$W({\lambda }_i|\delta v)=N_{\mathit{pix}}\frac{\overline{S}\prime ({\lambda }_i|\delta v)}{{\sigma }_{S\prime }({\lambda }_i|\delta v)}{\left[{\displaystyle \sum _j^{N_{\mathit{pix}}}}\frac{\overline{S}\prime ({\lambda }_j|\delta v)}{{\sigma }_{S\prime }({\lambda }_j|\delta v)}\right]}^{-1},$$(12)

where $\overline{S}\prime$ and σ_S′ are the mean absolute model derivative and its standard deviation among models. The normalisation of the W term is assumed such that $\sum _i^{N_{\mathit{pix}}}W({\lambda }_i)=N_{\mathit{pix}}$, as in the non-weighted case (W(λ_i) = 1 for all pixels). As $\overline{S}$ and ${\sigma }_{\overline{S}}$, W is also recalculated at each step of the MCMC process. In practice, the new scheme resulted in weighting line edges more strongly than the continuum and other flat regions, as seen in the bottom panel of Fig. 5.

We ran MCMC to estimate δv for each epoch independently. Figure 6 shows the posterior probability distributions probed by the MCMC analysis for δv₁ (filled green) and δv₂ (orange). The obtained results for the shifts are $\delta v_1=-2.{57}_{-3.73}^{+3.72}\mathrm{m}{\mathrm{s}}^{-1}$ and $\delta v_2=-3.{82}_{-2.44}^{+2.43}\mathrm{m}{\mathrm{s}}^{-1}$, where the uncertainties correspond to the 16th and 84th percentiles of the posterior distribution and should be considered equivalent to the 68% confidence interval. Taking the difference between the two values above and propagating their uncertainties, we obtained a velocity shift between the two epochs of $\mathrm{\Delta }v=-1.{25}_{-4.46}^{+4.44}\mathrm{m}{\mathrm{s}}^{-1}$, which translates into a drift of $\dot{v}=-1.{43}_{-5.10}^{+5.08}\mathrm{m}{\mathrm{s}}^{-1}{\mathrm{yr}}^{-1}$ or $\dot{z}=-2.{19}_{-7.78}^{+7.75}\times {10}^{-8}{\mathrm{yr}}^{-1}$.

Fig. 6. Posterior probability distributions obtained by the MCMC analysis comparing the mean model to the spectra of epoch 1 (filled green) and epoch 2 (orange) through the likelihood defined in Eq. (11) (see Sect. 6.4). Scatter points with error bars define the median values and the 68% confidence intervals.

Both estimates of δv_i have negative values. While this may first seem surprising, there is an easily understandable explanation related to the weighting scheme implemented in Eq. (11). To demonstrate this, we carried out a simple test of measuring δv on the combined spectrum, i.e. the same spectrum used to produce our 500 models. Naïvely, this should return a shift consistent with zero. Applying our measurement method to this data resulted in $\delta v_0=-2.{71}_{-2.02}^{+2.08}\mathrm{m}{\mathrm{s}}^{-1}$, a value that is negative and marginally inconsistent with zero. However, temporarily fixing W = 1 in Eq. (11) resulted in $\delta v_0=-0.{43}_{-3.51}^{+3.59}$, consistent with zero. Obviously, our choice of W has given more weight to pixels with negative contributions to δv while simultaneously decreasing measuring uncertainties. Another way of looking at this is through a comparison of δv₀ with δv₁ and δv₂. Comparing their values, one sees that δv₀ falls in between δv₁ and δv₂, as expected. In that sense, the obtained values of δv₁ and δv₂ are not surprising.

7. Validation and systematics

7.1. Measurement validation

We validated the measurement procedures by applying the pixel-by-pixel (see Sect. 5) and the model-based (see Sect. 6.4) methods to a synthetic sample of mock spectra for which we know a priori the velocity shift occurring between the two epochs. We thus produced mock datasets assuming the same S/N distribution, pixel size, spectral range, and resolution of the SB2 spectra, simulating a Δv of 0 cm s⁻¹, −0.38 cm s⁻¹, −4.3 m s⁻¹, and −42.9 m s⁻¹, corresponding to a temporal baseline of zero, 0.875, 10³ and 10⁴ years in a ΛCDM Universe.

We applied our methods to ten random sightlines from each mock dataset, with the results shown in Fig. 7. The top panel shows the measured velocity shifts, Δv, obtained with the pixel-by-pixel (blue) and model-based (red) methods in the four baseline cases. Both methods were applied to the same set of 10 sightlines. We found that the two methods are capable of recovering the imposed drift, with an uncertainty that depends solely on the S/N level of the two spectra (matched to the SB2 data, regardless of the simulated Δv) and on the number of Lyman-α lines contained in each sightline. The scatter among the ten measured velocity shifts is consistent with the estimated uncertainties. Notably, the pixel-by-pixel method shows a slight discrepancy between the expected signal and the mean measured Δv when applied to the longest baseline of 10 000 years. The cause of such bias is not clear, although it could stem from the perturbative assumption taken in Eq. (6) that might not hold any longer for the large shift imposed in this simulated case. Nonetheless, our aim is to measure very small shifts for which such an assumption is solid. The model-based method recovers the imposed shift well in all instances.

Fig. 7. Velocity shift, Δv, measured on ten random mock sightline pairs with both the pixel-by-pixel (blue square scatter points, see Sect. 5) and model-based methods (red diamond scatter points, see Sect. 6.4) where a different baseline is assumed between the two epochs. The horizontal blue dotted and dashed red lines report the average Δv measured with the pixel-by-pixel and the model-based methods, respectively. From left to right, the mock datasets are built assuming no shift (Δvexp = 0 m s−1), Δt = 0.875 yr (the same time elapsed between the two observed spectra, Δvexp=−0.38 cm s−1), Δt = 103 yr (Δvexp=−4.3 m s−1), and Δt = 104 yr (Δvexp=−42.9 m s−1). The dashed grey lines report the assumed shift.

7.2. Uncertainty across the spectrum

It is possible that a narrow segment of the Lyman-α forest contains unidentified contaminants that spoil the measurement on Δv. We investigated this possibility by dividing the Lyman-α forest into ten equally sized chunks (of width 9.4 nm) and measuring Δv over each individual one, using both the pixel-by-pixel and model-based methods, following the usual procedures described in Sects. 5 and 6.4.

Figure 8 shows the measured Δv for each chunk, as a function of the chunk's central wavelength λ_c, in the upper panel. The middle panel shows the velocity shift uncertainty, σ_v, computed on each chunk with the pixel-by-pixel (blue) and model-based (red) methods. In the absence of systematics, the measurement uncertainty is expected to depend solely on the median S/N at continuum and the number of pixels (or the fraction of the forest $f_{{\mathrm{Ly}}_{\alpha }}$) in each chunk, scaling as ${\sigma }_v\propto {\mathrm{S}/\mathrm{N}}^{-1}{f}_{{\mathrm{Ly}}_{\alpha }}^{-0.5}$. The expected measurement uncertainty is shown in the middle panel as a dash-dotted (dashed) grey line, scaled by the pixel-by-pixel (model-based) uncertainty obtained on the full spectrum. The estimated uncertainties closely follow the scaling expectations for both methods, highlighting the absence of wavelength-dependent systematic effects.

Fig. 8. Estimated velocity shift, Δv, computed on ten equispaced sections of the Lyman-α forest. Top panel: Δv estimated on each spectral chunk with the pixel-by-pixel (blue squares) and model-based (red diamonds) methods, as a function of the section's central wavelength λc. Middle panel: Measurement uncertainty, σv, for the two methods, shown with the same colour coding as the panel above. The dash-dotted and dashed grey lines report the expected uncertainty of the two methods, based on the median S/N and pixel number of each chunk. Bottom panel: Fraction of the Lyman-α forest $f_{{\mathrm{Ly}}_{\alpha }}$ considered in the measurement within each chunk (solid black line, left-hand axis), taking into account the masking of metals, strong absorbers and bad pixels. Note that the spectral chunk centred at 560.7 nm contains the sub-DLA absorber and has a significantly smaller fraction of usable pixels than the other regions. The right-hand axis of the same panel reports the median S/N at continuum level in each chunk in the first and second epoch spectra, shown with dashed green and orange lines, respectively.

7.3. Influence of the quasar

Given the extraordinary nature of SB2 (Wolf et al. 2024), with its high black hole accretion rate (Ṁ ∼ 370 M_⊙ yr⁻¹) and luminosity variability (∼15% over the last 6 years), one would expect to find a strong signature of IGM photoionisation in the Lyman-α forest close to the quasar emission redshift. However, our data do not show an important decrease in strength and number of Lyman-α lines close to z_em, nor a significant variability of such lines, which could be related to the quasar's strong activity. Therefore, in our previous analysis, we based the redshift drift measurement on the whole forest, without excluding the proximity region, about 5000 km s⁻¹ from the quasar emission redshift.

We checked for systematic effects due to the proximity to the quasar by performing the same measurement as before, but excluding such a region, and measured a velocity drift between the two epochs of Δv = 3.88±5.08 m s⁻¹. Having excluded the proximity region reduces the number of pixels used in the measurement to ∼37 400, or ∼74% of the total amount found in the forest. When ignoring the proximity region, we excluded the range close to the Lyman-α emission where the spectra have the highest S/N levels at the continuum, with a median of 89, 134 and 161 per pixel in the first epoch, second epoch, and combined spectra, respectively, and performed the measurement on spectra with effectively lower S/N. Therefore, the uncertainty grows not only due to a reduced pixel sample, but also due to a smaller median S/N level (∼82 per pixel at the continuum in the combined spectrum) outside of the proximity region. These two effects add up to a scaling factor ${\left(82/86\right)}^{-1}{\left(0.74/0.84\right)}^{-1/2}\sim 1.12$, explaining the new uncertainty. Still, the measurement is compatible with the expected value, and no systematic effect due to excluding or taking into account the proximity region is clearly visible.

7.4. Influence of local motions

As is shown by Inoue et al. (2020), the mass of the Milky Way (MW) and the other galaxies in the Local Group (LG), especially the Large Magellanic Cloud and M31, have a non-negligible impact on the proper motion of the Sun, providing a local source of acceleration along a certain sightline. Left uncorrected, this effect will impact the cosmological redshift drift signal measurement. The amplitude of this effect is sightline dependent and, in the case of SB2, induces an additional systematic drift in the Lyman-α lines of ${\dot{v}}_{\odot }=-0.16\pm 0.04\mathrm{cm}{\mathrm{s}}^{-1}{\mathrm{yr}}^{-1}$. This is of the same order of magnitude as the cosmological signal we want to measure. Thus, the expected velocity shift between two epochs of SB2 grows effectively to ${\dot{v}}_{\mathrm{SB}2}\sim -0.59\mathrm{cm}{\mathrm{s}}^{-1}{\mathrm{yr}}^{-1}$, implying it can be measured on a shorter temporal baseline, at fixed total integration time.

However, such a shift is the composite of the local and cosmological effects, where the former should be removed from the measured shift to infer the actual redshift drift. The required duration of the experiment is not affected by the local component, if not for the fact that its amplitude is inferred from measurements of the MW and LG's galaxies masses, and has an uncertainty that will be propagated to the redshift drift measurement's uncertainty. Such an uncertainty is one order of magnitude smaller than the expected redshift drift, of the order of 0.04 cm s⁻¹ yr⁻¹, where the cosmological signal is ∼−0.429 cm s⁻¹ yr⁻¹, and its propagation in the actual measurement is sub-dominant with respect to the measurement uncertainty related to the spectral analysis (see Sect. 8).

7.5. Instrumental systematics

ESPRESSO's wavelength calibration is a crucial aspect for robustly measuring the redshift drift signal. Previous studies identified differences as large as 20 m s⁻¹ between the two calibration methods available on ESPRESSO, the LFC and the Fabry-Pérot (FP) etalon combined with a ThAr lamp (Schmidt et al. 2021). Their impact is expected to decrease by considering a large number of lines over a broad spectral range, but the true quantification on a real measurement was so far lacking. We therefore repeated our analysis on SB2 observations calibrated using ThAr+FP frames instead of LFC, created specifically for this purpose following the same procedures used previously in Sects. 6.3 and 6.4. We measured $\mathrm{\Delta }v=-0.{81}_{-4.39}^{+4.47}\mathrm{m}{\mathrm{s}}^{-1}$ from these data, a change of +0.44 m s⁻¹ with respect to LFC-calibrated data. The two measurements agree within uncertainties.

Another source of uncertainty relates to assumptions on the shape of the line-spread function (LSF), used for determining line centres during wavelength calibration. The standard ESPRESSO DRS wavelength calibration procedures assume a Gaussian LSF shape and, because this cannot be modified without significant changes to the entire procedure, we were unable to directly quantify the impact of this assumption on our results. However, considering that using a Gaussian LSF should introduce wavelength calibration errors similar to the differences observed between using LFC and ThAr+FP (i.e. 20 m s⁻¹, Schmidt et al. 2021; Schmidt 2024), a reasonable expected additional systematic uncertainty is ∼0.5 m s⁻¹ through analogy with the previous paragraph.

8. Future perspectives

Based on our results, we estimated the time required for a statistically significant detection of the cosmological signal with ESPRESSO and ANDES. We envision an observational monitoring programme, carried out for N_e epochs, with a cadence of T hours of ESPRESSO integration of SB2 each year. The precision on the measurement of the cosmic acceleration ($\dot{v}\sim \mathrm{\Delta }v/\mathrm{\Delta }t$) scales in time, adapting Eq. (5) to account for the additional temporal dependency, as

$${\sigma }_{\dot{v}}\propto g(N_e)\times \mathrm{S}/{\mathrm{N}}_{\mathrm{tot}}(N_e,T{)}^{-1}\times {\left(N_e-1\right)}^{-1},$$(13)

where the experiment baseline is Δt=N_e−1 (e.g. 1 year of the programme has 2 epochs, 2 years of the programme have 3 epochs, and so on), and the total S/N obtained after N_e epochs scales as $\mathrm{S}/{\mathrm{N}}_{\mathrm{tot}}(N_e,T)\propto \sqrt{N_{e}T}$. The form factor, g(N_e), depends on the distribution across time of the observational epochs, which was 1.7 for the continuous monitoring programme that we considered.

Figure 9 shows the acceleration measurement uncertainty, ${\sigma }_{\dot{v}}$, achieved with ESPRESSO as a function of the temporal baseline of our experiment, starting from the data already collected, and following up with a monitoring campaign of 10, 100, or 1000 hours of integration time each year. The vertical solid lines denote at what epoch we can achieve a 3σ (or 5σ vertical dotted lines) detection of the cosmological signal in the three cases. The estimates of the required baseline for a detection are also reported in Table 3.

Fig. 9. Velocity shift uncertainty reached with ESPRESSO spectra of SB2, as a function of total temporal baseline of the experiment, assuming three different observational strategies with an integration time of 10 hours per year (blue), 100 hours per year (yellow), and 1000 hours per year (red). The horizontal dashed black line defines the magnitude of the cosmological signal expected for a ΛCDM Universe at the average redshift of SB2's Lyman-α forest. Vertical lines define the time needed to achieve a 3σ (solid) and 5σ (dotted) detection for the three observational strategies.

Clearly, the timescales for the experiment are long due to the limited collecting area of VLT. However, these results can be used to estimate, at first order, what the impact will be of moving the observational campaign to the ELT/ANDES spectrograph (Marconi et al. 2024). Assuming that the total efficiencies of the ESPRESSO and ANDES spectrographs are similar (about 10%), since the latter is fed by a telescope with a collecting area ∼20 times larger than the area of one VLT unit, the same accuracy in measuring Δv can be reached in 1/20-th of the total ESPRESSO integration time. By this logic, a 3σ (5σ) detection of the redshift drift can be carried out observing SB2 with ANDES in 68 years (95 years), with 50 hours of total integration time per year, instead of 1000 hours per year as in the case of ESPRESSO.

Similarly, assuming a fixed yearly integration time, the same level of precision in the measurement of the cosmic acceleration ${\sigma }_{\dot{v}}$ can be achieved with ANDES with a temporal baseline that is a factor 20^1/3∼2.714 shorter than in the case of an ESPRESSO experiment, due to the scaling of Eq. (13). In this case, a 3σ (5σ) detection of the redshift drift can be carried out by means of a monitoring campaign of SB2 with 100 hours of integration per year, in 54 years (75 years) with ANDES, instead of 145 years (203 years) with ESPRESSO. Table 3 reports the temporal baselines required for detection in an ANDES-based experiment, when allocating 10, 100, or 1000 total hours of exposure per year.

The ANDES spectrograph is yet to be constructed, and its efficiency could change in the future design phases. We can easily parametrise the required baseline of the measurement to reach a 3σ detection, letting the total efficiency of ANDES, ϵ, as a free parameter:

$$\mathrm{\Delta }{t}_{3\sigma }^{\mathrm{ANDES}}\approx 54{\left(\frac{T}{100\hspace{0.5em}\mathrm{h}{\mathrm{yr}}^{-1}}\right)}^{-1/3}{\left(\frac{ϵ}{0.1}\right)}^{-1/3}\hspace{0.5em}\mathrm{yr},$$(14)

where T is the total integration time of an ANDES-like spectrograph per year of the experiment. The total amount of integration hours needed for a 3σ detection is simply

$$\mathrm{Total}\hspace{0.5em}\mathrm{integration}\hspace{0.5em}\mathrm{time}=5.4\times {10}^3{\left(\frac{T}{100\hspace{0.5em}\mathrm{h}{\mathrm{yr}}^{-1}}\right)}^{2/3}{\left(\frac{ϵ}{0.1}\right)}^{-1/3}\hspace{0.5em}\mathrm{h},$$(15)

depending on the amount of telescope time that can be continuously allocated to this project.

9. Results and discussion

We set out to assess the feasibility of measuring the redshift drift signal using Lyman-α forest spectroscopy, based on the highest-quality data available prior to the commissioning of ELT/ANDES. We used new, high-resolution, and high-S/N spectra of the most luminous known quasar in the Universe, J052915.80-435152.0 (SB2; z_em = 3.962), obtained with ESPRESSO at the VLT. This exceptionally bright object is especially suited to perform the Sandage test of cosmic redshift drift, owing to its luminosity and emission redshift. We began this decades-long experiment with the ultra-stable spectrograph ESPRESSO over a one-year baseline. The primary goal is to investigate the measurement's intrinsic systematics and validate the predictions on the expected total integration time required to achieve a statistically significant detection. This work presents the first two epochs of the experiment, intended to serve as an anchor for comparison to future observations, initially carried out with ESPRESSO and subsequently with ELT/ANDES.

Our main results are:

• Applying the pixel-by-pixel method of Bouchy et al. (2001) (Sect. 5), we measured a velocity shift between two epochs of Δv=−2.67±6.64 m s⁻¹. These uncertainties are 65% above the prediction made by Liske et al. (2008) (σ_v,Liske = 4.02 m s⁻¹, see Eq. (5)).

• Using an ensemble of 500 models specifically tailored to the resolution, pixel size, and S/N of our data (Sect. 6), we identified spectral regions more sensitive to velocity shifts and gave them more weight in our analysis, obtaining $\mathrm{\Delta }v=-1.{25}_{-4.46}^{+4.44}\mathrm{m}{\mathrm{s}}^{-1}$. This is a 30% improvement in statistical uncertainty over the direct method, and only 10% greater than the value predicted by Liske et al. (2008).

• The redshift drift measured between the two epochs turns out to be $\dot{z}=-2.{19}_{-7.78}^{+7.75}\times {10}^{-8}{\mathrm{yr}}^{-1}$ at 〈z〉 = 3.57, which represents the main result of this paper. This value is fully consistent with the value expected within the ΛCDM model, −6.6×10⁻¹¹ yr⁻¹. The quoted uncertainties are only statistical and three orders of magnitude larger than the cosmological signal. Systematic uncertainties due to the motion of the Solar System, temporal variability of SB2's luminosity, and wavelength calibration are currently sub-dominant.

• Based on these findings, we estimate that a 3σ detection of the redshift drift could be achieved with ESPRESSO (ANDES) after 145 (54) years, assuming 100 hours of observation per year for SB2 (Table 3).

The best current upper limit on cosmic acceleration is $|\dot{v}|\le 2.2\mathrm{m}{\mathrm{s}}^{-1}{\mathrm{yr}}^{-1}$ at a 1σ confidence level (Darling 2012). Our measurement of $\mathrm{\Delta }v=-1.{25}_{-4.46}^{+4.44}\mathrm{m}{\mathrm{s}}^{-1}$ translates into $\dot{v}=-1.{43}_{-5.10}^{+5.08}\mathrm{m}{\mathrm{s}}^{-1}{\mathrm{yr}}^{-1}$ (Sect. 6), making our limit a factor of 2.3 larger than that of Darling (2012). However, their result is based on observations of H I 21 cm absorption line systems in ten objects over a baseline of 13 years, while we obtained a comparable result with only 12.4h of Lyman-α forest and a ∼1 year baseline. By expanding our dataset by another 12 hours over the next 1.5 years, we expect to match the uncertainty in the measurement of Darling (2012), showcasing the power of Lyman-α observations for ż measurements.

Interestingly, neither the pixel-by-pixel method (Sect. 5) nor our newly developed method based on splines (Sect. 6) precisely reproduce the uncertainties predicted by the scaling relation of Liske et al. (2008). This discrepancy may arise because the Liske et al. (2008) relation was calibrated assuming a measurement carried out on a sample of multiple objects, whereas we analysed a single sightline. Sightline-to-sightline fluctuations in the measurement uncertainties are expected due to cosmic variance, where the analysis of SB1 will provide further insights into this issue (Marques et al., in prep.).

We examined the validity of our methods and some potential sources of systematic uncertainties in Sect. 7. While both the direct and the spline methods recover simulated Δv (within uncertainties) for |Δv|<5 m s⁻¹, the pixel-by-pixel method shows a systematic underestimation of the shift when recovering the simulated |Δv|≈40 m s⁻¹ (Fig. 7). Such a result is not currently understood. Solar System accelerations add ∼1 cm s⁻¹ yr⁻¹ to the expected signal, being negligible at the current level of S/N. The presence of such an effect does not alter the experiment's expected temporal baseline (see Sect. 7.4). Excluding the quasar proximity zone from our analysis shifts Δv by 5.14 m s⁻¹, a value only marginally inconsistent (∼1.01σ) with the statistical uncertainty (Sect. 7.3). The shift is likely a consequence of a large (10%) decrease in the number of available pixels, and not a true estimate of uncertainty from quasar variability. More data and longer time baselines are needed to confirm whether temporal variability in quasar flux causes a systematic effect on Δv. Wavelength calibration uncertainties are estimated to be ≈1 m s⁻¹ (Sect. 7.5). Considering the latter is approximately one fourth of the statistical uncertainty, we expect that systematic effects will start playing a significant role in the error budget (in SB2 measurements) once the total S/N quadruples (see Eq. (5)).

A caveat of both methods used in this paper is that they are unable to disentangle a real velocity shift in the data from any changes to the instrument properties, such as its LSF. The LSF is known to change when important components are replaced (Lo Curto et al. 2015), and is likely to occur in experiments spanning several decades. There were no major interventions in ESPRESSO over the time period spanning observations presented here that would raise concerns about LSF or other instrument properties varying. However, removing instrumental imprints will be important in the future; for example, by employing advanced data reduction frameworks such as ‘spectro-perfectionism’ (Bolton & Schlegel 2010) or similar (Piskunov et al. 2021).

Based on considerations presented in Sect. 8, a 3σ detection of ż using SB2 will require 54 years of observations with ANDES when allocating 12.5 nights per year to this project (Table 3). Shorter timescales can be achieved with a more intensive campaign or by increasing the total efficiency of ANDES (see Eq. (14)). For obvious reasons, an intensive monitoring campaign of such magnitude focusing on only one object is not ideal. In a more realistic case, multiple bright objects will be observed simultaneously, spread over a large range of right ascensions. A golden sample of seven quasars compiled with this goal in mind has already been proposed by Cristiani et al. (2023). In such a case, SB2 being the brightest quasar of the sample, the experiment baseline will slightly increase with respect to Eq. (14) at fixed total observational cadence, T, with time shared among multiple targets. However, the other quasar's sightline might be better suited for the test, with more narrow lines, fewer strong absorbers, and a larger fraction of the forest available for the measurement than what we found on SB2. A measurement on such sightlines would have smaller uncertainties and counterbalance the lower S/N related to the fainter objects. At this level, it is not straightforward to formally predict the consequences of such effects on the baseline of an experiment involving all seven quasars of the golden sample. In a parallel paper (Marques et al, in prep.), we shall analyse the spectrum of SB1, the second brightest object of the sample, in order to address this point.

Data availability

Table A.1 is available at the CDS via anonymous ftp to cdsarc.cds.unistra.fr (130.79.128.5) or via https://cdsarc.cds.unistra.fr/viz-bin/cat/J/A+A/699/A159

Acknowledgments

AT is grateful to Vid Iršič for the prolific discussions. The INAF authors acknowledge financial support of the Italian Ministry of Education, University, and Research with PRIN 201278X4FL and the “Progetti Premiali” funding scheme. This work was financed by Portuguese funds through FCT (Fundação para a Ciência e a Tecnologia) in the framework of the project 2022.04048.PTDC (Phi in the Sky, DOI 10.54499/2022.04048.PTDC). CJM also acknowledges FCT and POCH/FSE (EC) support through Investigador FCT Contract 2021.01214.CEECIND/CP1658/CT0001 (DOI 10.54499/2021.01214.CEECIND/CP1658/CT0001). CMJM is supported by an FCT fellowship, grant number 2023.03984.BD. MTM acknowledges the support of the Australian Research Council through Future Fellowship grant FT180100194 and through the Australian Research Council Centre of Excellence in Optical Microcombs for Breakthrough Science (project number CE230100006) funded by the Australian Government. TMS acknowledges the support from the SNF synergia grant CRSII5-193689 (BLUVES). The work of KB is supported by NOIRLab, which is managed by the Association of Universities for Research in Astronomy (AURA) under a cooperative agreement with the U.S. National Science Foundation. EP acknowledges financial support from the Agencia Estatal de Investigación of the Ministerio de Ciencia e Innovación MCIN/AEI/10.13039/501100011033 and the ERDF “A way of making Europe” through project PID2021-125627OB-C32, and from the Centre of Excellence “Severo Ochoa” award to the Instituto de Astrofisica de Canarias. ASM and JIGH acknowledge financial support from the Spanish Ministry of Science, Innovation and Universities (MICIU) projects PID2020-117493GB-I00 and PID2023-149982NB-I00. NS and NN acknowledge financial support by FCT - Fundação para a Ciência e a Tecnologia through national funds by these grants: UIDB/04434/2020 DOI: 10.54499/UIDB/04434/2020, UIDP/04434/2020 DOI: 10.54499/UIDP/04434/2020. FP would like to acknowledge the Swiss National Science Foundation (SNSF) for supporting research with ESPRESSO through the SNSF grants nr. 140649, 152721, 166227, 184618 and 215190. The ESPRESSO Instrument Project was partially funded through SNSF's FLARE Programme for large infrastructures.

References

All Tables

All Figures

	Fig. 1. Combined spectrum of SB2. Top panel: Flux density in arbitrary units is shown in black, with the dashed purple line denoting the flux density error. Bottom panel: S/N per 1 km s−1 pixel. The shaded green area highlights the Lyman-α forest considered in the redshift drift measurement, bound by the Lyman-α and Lyman-β emissions of the quasar, namely between 509−603 m.
In the text

	Fig. 2. Normalised transmitted flux of the Lyman-α forest in the combined spectrum of SB2 as a function of wavelength. The solid light blue line highlights the spectral regions masked due to metal absorption, the sub-DLA and bad pixels. The solid purple line reports the flux error.
In the text

	Fig. 3. Distribution of the S/N per pixel in the Lyman-α forest of the three spectra of SB2 (dashed red) and of the corresponding mock spectra (solid black). Upper panel: First epoch. Middle panel: Second epoch. Bottom panel: Combined spectrum.
In the text

	Fig. 4. Standard deviation of the fit residuals between the combined mock spectra's flux and the spline model with internodal distance A. The blue squares describe the σR of models on different mock spectra and randomly drawn initial node phase ϕ. The black error bars define the average and standard deviation of σR(A)−1. The vertical dashed purple line denotes the optimal value of A.
In the text

Fig. 5. Section of the SB2 spectrum. Top panel: Normalised flux of the first epoch (green), second epoch (orange), and combined (purple) spectra. The mean model, $\overline{S}$, is shown as a solid black line. The spectra are shifted vertically for clarity. Middle panel: Normalised flux errors, σF, of the three spectra shown in the panel above, with the same colour coding. The solid black line is the standard deviation of our model, ${\sigma }_{\overline{S}}$, i.e. it shows where our ensemble of models exhibits large variation, and hence uncertainty. Bottom panel: Normalised pixel weights, W, as defined in Eq. (12). Regions with higher values provide more constraints on our final ż measurement.

In the text

	Fig. 6. Posterior probability distributions obtained by the MCMC analysis comparing the mean model to the spectra of epoch 1 (filled green) and epoch 2 (orange) through the likelihood defined in Eq. (11) (see Sect. 6.4). Scatter points with error bars define the median values and the 68% confidence intervals.
In the text

Fig. 7. Velocity shift, Δv, measured on ten random mock sightline pairs with both the pixel-by-pixel (blue square scatter points, see Sect. 5) and model-based methods (red diamond scatter points, see Sect. 6.4) where a different baseline is assumed between the two epochs. The horizontal blue dotted and dashed red lines report the average Δv measured with the pixel-by-pixel and the model-based methods, respectively. From left to right, the mock datasets are built assuming no shift (Δvexp = 0 m s−1), Δt = 0.875 yr (the same time elapsed between the two observed spectra, Δvexp=−0.38 cm s−1), Δt = 103 yr (Δvexp=−4.3 m s−1), and Δt = 104 yr (Δvexp=−42.9 m s−1). The dashed grey lines report the assumed shift.

In the text

Fig. 8. Estimated velocity shift, Δv, computed on ten equispaced sections of the Lyman-α forest. Top panel: Δv estimated on each spectral chunk with the pixel-by-pixel (blue squares) and model-based (red diamonds) methods, as a function of the section's central wavelength λc. Middle panel: Measurement uncertainty, σv, for the two methods, shown with the same colour coding as the panel above. The dash-dotted and dashed grey lines report the expected uncertainty of the two methods, based on the median S/N and pixel number of each chunk. Bottom panel: Fraction of the Lyman-α forest $f_{{\mathrm{Ly}}_{\alpha }}$ considered in the measurement within each chunk (solid black line, left-hand axis), taking into account the masking of metals, strong absorbers and bad pixels. Note that the spectral chunk centred at 560.7 nm contains the sub-DLA absorber and has a significantly smaller fraction of usable pixels than the other regions. The right-hand axis of the same panel reports the median S/N at continuum level in each chunk in the first and second epoch spectra, shown with dashed green and orange lines, respectively.

In the text

Fig. 9. Velocity shift uncertainty reached with ESPRESSO spectra of SB2, as a function of total temporal baseline of the experiment, assuming three different observational strategies with an integration time of 10 hours per year (blue), 100 hours per year (yellow), and 1000 hours per year (red). The horizontal dashed black line defines the magnitude of the cosmological signal expected for a ΛCDM Universe at the average redshift of SB2's Lyman-α forest. Vertical lines define the time needed to achieve a 3σ (solid) and 5σ (dotted) detection for the three observational strategies.

In the text