Precision Ephemerides for Gravitational-wave Searches – IV: Corrected and refined ephemeris for Scorpius X-1

Low-mass X-ray binaries have long been theorised as potential sources of continuous gravitational-wave radiation, yet there is no observational evidence from recent LIGO/Virgo observing runs. Even for the theoretically ‘loudest’ source, Sco X-1, the upper limit on gravitational-wave strain has been pushed ever lower. Such searches require precise measurements of the source properties for sufficient sensitivity and computational feasibility. Collating over 20 years of high-quality spectroscopic observations of the system, we present a precise and comprehensive ephemeris for Sco X-1 through radial velocity measurements, performing a full homogeneous reanalysis of all relevant datasets and correcting previous analyses. Our Bayesian approach accounts for observational systematics and maximises not only precision, but also the fidelity of uncertainty estimates — crucial for informing principled continuous-wave searches. Our extensive dataset and analysis also enables us to construct the highest signal-to-noise, highest resolution phase-averaged spectrum of a low-mass X-ray binary to date. Doppler tomography reveals intriguing transient structures present in the accretion disk and flow driven by modulation of the accretion rate, necessitating further characterisation of the system at high temporal and spectral resolution. Our ephemeris corrects and supersedes previous ephemerides, and provides a factor three reduction in the number of templates in the search space, facilitating precision searches for continuous gravitational-wave emission from Sco X-1 throughout the upcoming LIGO/Virgo/KAGRA O4 observing run and beyond.


INTRODUCTION
With the completion of their third observing run, the catalogue of gravitational-wave (GW) signals detected by the Advanced Laser Interferometer Gravitational Wave Observatory (LIGO; LIGO Scientific Collaboration et al. 2015) and Advanced Virgo (Acernese et al. 2015) is growing rapidly (The LIGO Scientific Collaboration et al. 2021).These sources include the mergers of binary black holes (BHs) (Abbott et al. 2016), a binary neutron star (NS) (Abbott et al. 2017c,e), and NS-BH binaries (Abbott et al. 2021b).However, compact-object mergers are not the only sources expected to produced detectable GW emission.Unlike these transient signals, continuous GW (CW) sources present persistent quasi-monochromatic emission.This as-yet undetected type of gravitational radiation is emitted by rapidly rotating asymmetric NSs, whether found as isolated sources or within stellar binaries; see, e.g., Lasky (2015); Sieniawska & Bejger (2019); Piccinni (2022) for recent reviews.
Numerous mechanisms have been suggested to induce the time-★ E-mail: t.killestein@warwick.ac.uk varying quadrupole required for GW emission from spinning NSs, whose angular momentum loss would limit the spin rate and account for the observation that NS spins are measured at ≲ 700 Hz (Hartman et al. 2003;Hessels et al. 2006;Patruno et al. 2017), below estimated breakup frequencies (Chakrabarty et al. 2003).Low-mass X-ray binaries (LMXBs) with NS primaries have received particular attention as target sources (Abbott et al. 2007;Whelan et al. 2015;Aasi et al. 2015;Abbott et al. 2017aAbbott et al. ,a, 2019b;;Middleton et al. 2020;Zhang et al. 2021;Abbott et al. 2022b).In this scenario, a torque balance is achieved between the spin-up due to accretion from a stellar companion and spin-down due to GW emission (Papaloizou & Pringle 1978;Wagoner 1984;Bildsten 1998;Andersson et al. 1999).Such a rotational equilibrium leads to a characteristic strain that increases with increasing X-ray flux of the source (a proxy for the accretion rate; see Bildsten 1998).The most promising candidates are thus those that are most bright in X-rays.
The prototypical LMXB, Sco X-1 (Giacconi et al. 1962;Sandage et al. 1966;Shklovsky 1967), is composed of an accreting NS primary and donor star, and has been intensively studied since its discovery as among the closest X-ray binaries known (Gottlieb et al. 1975;Bradshaw et al. 1999;Fomalont et al. 2001).It shows strong emission across the electromagnetic (EM) spectrum, from radio (Fomalont et al. 1983) to gamma-rays (Brazier et al. 1990), powered by a near-Eddington accretion rate from the donor star.Intriguing null detections of very high energy (VHE, TeV) emission (Aleksić et al. 2011) suggest that the high-energy emission mechanism in Sco X-1 is markedly different to other systems.The donor star in the system remains enigmatic; the high accretion luminosity shrouds any stellar absorption features present in the near-infrared (Mata Sanchez et al. 2015), but these observations combined with dynamical constraints suggest a donor mass of 0.28 <  2 < 0.7  ⊙ (stellar type K4IV or later).As the strongest source of X-rays on the sky, Sco X-1 has been used as a 'lighthouse' to study intervening interstellar material (García et al. 2011), magnetic fields (Titarchuk et al. 2001), and even the Martian atmosphere (Rahmati et al. 2020).Very recently, polarised X-ray emission was detected by the PolarLight mission at keV energies (Long et al. 2022), providing a strong constraint on the system geometry and suggesting that the X-ray emission arises primarily from a compact, optically thin corona near the disc transition layer.
Despite being a cornerstone of our understanding of compact binary systems across the Universe, Sco X-1 still remains enigmatic, showing great complexity and variability.These complexities are simply not revealed in more distant systems where lower signal-tonoise ratio (SNR) limits their visibility.Putting this aside, Sco X-1 is the most luminous extra-solar X-ray source, which combined with its relative proximity (2.3 ± 0.1 kpc, Lindegren et al. 2021), implies it should be among the loudest CW sources detectable by current GW detectors (Watts et al. 2008) and has made it the target of a great number of search efforts.
No concrete detection of CW emission has yet been made of what we a priori expect is the strongest GW source, which remains puzzling.However, improvements in detector sensitivity have led to correspondingly stronger upper limits on the GW strain from Sco X-1, (Abbott et al. 2007;Aasi et al. 2015;Abbott et al. 2017a;Meadors et al. 2017;Abbott et al. 2017dAbbott et al. , 2019b;;Zhang et al. 2021;Abbott et al. 2022b), with the most recent analyses reaching ≲ 4 × 10 −26 (in the most sensitive frequency bands and assuming knowledge of the inclination angle; Abbott et al. 2022c).This level of strain begins to push below the torque-balance limit (Bildsten 1998), where the expected GW emission balances the spin-up torque from donor accretion.One complicating factor in the case of Sco X-1 is the absence of a measured spin period, unlike many other accreting NS systems (Abbott et al. 2022a).Sco X-1 has also been the target of directional, stochastic searches (Abbott et al. 2017b(Abbott et al. , 2019a(Abbott et al. , 2021a)), leading to further sensitive (yet also null) results.
As one of the primary uncertainties in any directed search, extensive X-ray timing observations have been performed in search of the spin period of the (assumed) NS primary, as revealed by X-ray pulsations or bursts (e.g.Galloway et al. 2010).The most recent placed an upper limit of 0.034% (90% confidence) on any putative X-ray variability using Rossi X-ray Timing Explorer (RXTE) data (Galaudage et al. 2022).It remains unclear whether Sco X-1 is intrinsically or intermittently variable in the X-ray, and whether any pulsed X-ray emission is being scattered away by surrounding material.
Another crucial component to facilitate CW searches are precise orbital constraints for the LMXB systems of interest (Watts et al. 2008).In order to combine long sections of data coherently, the motion of both the Earth relative to the Solar System and the orbital motion of the NS in the binary must be accounted for.It is computationally infeasible to marginalise over the vast, high-dimensional parameter space this presents, and therefore, it is critical to place con-straints on this orbital motion to reduce the search dimensionality and increase sensitivity (Dhurandhar & Vecchio 2001;Messenger et al. 2015;Leaci & Prix 2015).This presents a signficant observational challenge however; in a high accretion rate system like Sco X-1, the Balmer features are blurred significantly, limiting precision and making them unsuitable for precision radial velocity (RV) measurements.The disc structure is also dynamic and chaotic, leading to timeand phase-dependent changes in the emission line geometry.
Fortuitously, the irradiated donor star provides an alternative observational probe via narrow emission lines generated by Bowen fluorescence, which provide a remarkably precise probe of orbital motion, as first demonstrated by Steeghs & Casares (2002).The NS irradiates the face of the donor star with a strong X-ray flux, ionising He II which then de-excites, triggering a cascade of emission from C/N/O atomic orbitals (Kastner & Bhatia 1996).The (comparatively) compact emission geometry generates narrow emission lines (with little intrinsic broadening), which precisely trace the heated face of the donor in contrast to disc emission, which has contributions from either side of the disc and thus experiences significant Doppler broadening and complex geometric distortions from specific regions, such as the hot-spot/bulge or the gas stream.This tracer of orbital motion can be used to estimate the orbital parameters of the NS with RV measurements.Beyond Sco X-1, this technique has found broad applicability to LMXBs in general (e.g., Casares et al. 2003;Cornelisse et al. 2007;Brauer et al. 2018).To this end, Sco X-1 has benefited from an extensive spectroscopic monitoring campaign since 1999 as a part of the Precision Ephemerides for Gravitational-wave Searches (PEGS) project.This has lead to incremental improvements in the ephemeris accuracy (Galloway et al. 2014;Wang et al. 2018), with coverage continuing up to the present day through an extensive allweather campaign with the Very Large Telescope (VLT) Ultraviolet and Visual Echelle Spectrograph (UVES; Dekker et al. 2000) providing over 200 high-quality spectra.
In support of ongoing CW searches, in this paper, we present the most up-to-date ephemeris for Sco X-1 in the literature.We leverage 20 years of spectroscopic coverage with a homogeneous approach that accounts for potential systematic errors -with a view to maximising not only the precision of the ephemeris, but also the quality of uncertainty estimates -to correct and improve upon previous analyses.This is crucial when constraining the parameter space for CW searches; over-optimistic predictions lead to missing out valid areas of parameter space, whereas over-estimated errors greatly increase the computational burden of such searches.In Section 2, we describe the observational data used.In Section 3, we describe the method to infer the orbital properties of Sco X-1 via RV measurements and present our updated -and corrected -ephemeris.In Section 4, we further explore the uncertain variability in Sco X-1.We present our conclusions in Section 5

Spectroscopic observations
As part of the comprehensive analysis performed in this paper, we collate all spectroscopy presented in the previous PEGS project papers (Galloway et al. 2014;Wang et al. 2018) and combine this with more recent VLT/UVES data, extending the overall observational baseline to 22 years .We obtained (or retrieved) 264 spectra with the UVES (Dekker et al. 2000) instrument mounted on the 8.2m VLT Unit Telescope 2. All observations were obtained in service mode, across variable observing conditions but with a typical SNR in the range 50-100.We focus here on the spectra taken  with the blue arm using the CD#2 dispersive element, achieving a typical resolution Δ/ ∼ 40, 000 per spectral element with the 1" slit.Data were reduced using the v5.10.13ESO UVES pipeline.We also include all relevant legacy datasets presented in the previous PEGS releases -two runs of time-resolved spectroscopy taken with the ISIS double-beam spectrograph on the William Herschel Telescope (WHT) in 1999 and 2011, respectively.These observations were taken with the R400B grating and reduced using the molly software (Marsh 2019).All spectroscopic datasets are enumerated in Table 1.A pictorial summary of the VLT/UVES data is given in Figure 1, where we plot the phase-folded spectrogram of the Bowen line region implied by the binary constraints made in Section 3.2.For all datasets, we recompute both the barycentric time and velocity corrections at mid-exposure to reduce any scatter introduced by the different conversion routines used in the data processing pipelines.We use the astropy.timemodule and adopt the JPL DE405 ephemerides as our reference.To provide both consistency with previous ephemerides and an authoritative value, we give ephemeris times in both UTC Barycentric Julian Date (BJD) and GPS seconds.

Measuring Bowen line velocities
We obtain our Bowen radial velocities by fitting a constrained line model similar to the one used in Wang et al. (2018) to the reduced spectra.An example spectrum with the model is given in Figure 2.
The narrow Bowen components are modelled with Gaussians of fixed width (50 km/s) and variable amplitudes, offset from their rest frame with a common velocity.In the region of interest, we identify five narrow-line N III/C III/O II components that optimally constrain the Bowen RVs: 4634.13Å, 4640.64 Å, 4647.42Å, 4650.25 Å, and 4643.37 Å respectively.Some of these narrow lines are not visible at specific phases due to system geometry.To accommodate this, we constrain all line amplitudes to be ≥ 0, such that these emission lines cannot be forced by continuum modulation to unphysical nonnegative fluxes when they are not present.The broad-line component is modelled with a Gaussian of fixed width 1250 km/s.The amplitude and centroid of this component are fitted for to remove additional correlations with the centres of the narrow lines.The line widths above were measured by Gaussian fits to the individual components, and are kept fixed to limit the number of free parameters.Empirically, the narrow-line components do not change width significantly as a function of orbital phase, and any model-data mismatch is unlikely to cause large deviations in the fitted velocities as the emission-line peaks are typically clear and share a common velocity.
The final line model has eight free parameters that are well constrained by the dense spectral sampling (≈ 0.03 Å/pix) of the UVES data.Our specific line model mitigates correlations between parameters, such that uncertainties in the other fitted parameters do not skew the Bowen line velocities.Through experimentation, the UVES data could potentially benefit from the addition of more lines, but this leads to issues in the lower-quality WHT data with deviant line amplitudes so we opt for the above model for all datasets.
As a pre-processing step, we continuum-subtract the spectra by fitting a third-order Chebyshev polynomial to 10 Å regions at the red and blue ends of the Bowen line region (4605 Å-4675 Å).We fit the line model using nonlinear least squares, as implemented in the least_squares routine in the scipy package (Virtanen et al. 2020).The original ephemeris of Wang et al. (2018) is used to set the initial Bowen line velocity, and other parameters are initialised with sensible defaults from the highest SNR spectrum.The RVs are then computed from the mutual Doppler shift of the measured emission-line locations with respect to their known rest-frame values.Uncertainties are rescaled to enforce a unity reduced chi-squared test statistic.We add a constant value of 0.5 km/s in quadrature with the error values to ensure uncertainties are not underestimated due to poor conditioning of the least-squares fit, and to include uncertainties in the absolute velocity calibration of UVES.At this stage of reanalysis we are more concerned with correct relative error scaling between velocity measurements, as we rescale the errors between each dataset in latter analysis steps.

Corrections to previous Sco X-1 ephemerides
As part of our homogeneous reanalysis, we identified calibration errors that affect the timing parameters of both previously published ephemerides for the Sco X-1 system (Galloway et al. 2014;Wang et al. 2018).These were not readily apparent in previous work as the statistical uncertainties on the WHT-derived datasets masked their effect on the ephemeris.With the inclusion of four times as many VLT spectra, these discrepancies become significant and thus must be corrected.For full transparency, we elaborate in greater detail on the specific errors and their effect on the ephemeris: • The subset of VLT data used in PEGS I and PEGS III were inadvertently not corrected to the mid-exposure time, nor adjusted to the Solar System heliocentre.As a result, the values of  0 presented in these works bear systematic errors of ∼ 600s.
• Covariances between orbital parameters are underestimated in Wang et al. (2018) due to being taken from the initial least-squares fit to the data (using a two-point finite-difference Jacobian approximation), rather than being computed from the samples from which the ephemeris is derived.This changes the covariance by two orders of magnitude.
Discrepancies were identified by cross-checking the reduced data against the raw frames located at the European Southern Observatory (ESO) Archive, in particular the metadata stored in the file headers.Through a careful reanalysis of the datasets, we confirm that addressing these issues brings the data used in both G14 and W18 into strong agreement, and jointly agree with the VLT data presented in this paper (accounting for deviations in the period due to dependence on low-resolution WHT data).We plot marginal posteriors for each ephemeris in Figure 3 to illustrate how the ephemerides change with the corrections applied above.These issues underscore the importance of both robust archiving of astronomical data, such that raw frames are easily retrievable even decades later, and the importance of storing metadata alongside these frames and documenting reduction steps; without this, these timing issues could not have been easily diagnosed.With this additional scrutiny, we are now confident any calibration artifacts are accounted for in our reanalysis -our previously-published ephemerides are considered obsoleted by the ephemeris presented in this paper, and should not be used in CW searches going forward.

Bayesian modelling of the Keplerian orbit
Moving forward with the reanalysis, we apply here the velocity corrections derived during the previous steps to bring all spectra into  the Solar System barycentric frame, and convert observation times to BJD.We derive our ephemeris using a standard Keplerian orbit model (assuming zero eccentricity), using for the Bowen line RVs where  is the velocity semi-amplitude,  0 is the reference epoch,  is the period, and  is the systemic velocity.To ensure good initialisation prior to sampling, we first do a simple fit to the above model with nonlinear least-squares.To minimise covariance between  and  0 in our ephemeris, we use a simple bracketing line search to find the time of conjunction  ′ 0 =  0 +  that minimises cov(,  0 ) for an integer number of orbital cycles .This 'seed' ephemeris is then used to initialise parameters for the more complex modelling that follows, reducing burn-in time during sampling and ensuring our final ephemeris provides the minimal covariance.
To marginalise over systematics, we include two nuisance parameters per dataset: a constant offset term,   , to correct for differences in absolute wavelength calibration between the datasets, and an error scaling term,   , intended to correct for underestimated uncertainties.We assume a Gaussian likelihood on the RVs  and their uncertainties RV curve for Sco X-1.We show the 1 confidence intervals of the RVs measured from the 1999 WHT (blue), 2011 WHT (orange), and VLT (green) spectroscopic datasets, derived from the Bowen line model.We also plot the median posterior predictive RV curve implied by our ephemeris constraints (red) and the corresponding 1, 2, and 3 confidence regions (darker to lighter red shaded regions, respectively).Note the significantly reduced intrinsic scatter of the higher resolution VLT data with respect to the older WHT data.Bottom panel: Normalized residuals between the observational RV measurements with the median RV curve, scaled by the respective error factors    measured at time , given by log L (, , | 0 , , , ,   ,   ) ∝ ∑︁ where  represents the index of each individual dataset,   represents the per-dataset velocity offset,   represents the per-dataset error scaling.The VLT datasets have the most stable and accurate long-term calibration, and so the   value for this dataset is fixed to zero.It should be noted that this offset term effectively also absorbs longterm velocity variations in the system -e.g., due to perturbations from a distant companion.Over the observational baseline, these would be linear in order.Regardless, the focus of this work is to provide the most precise timing of the system possible, and we have no reason to expect significant secular motion.To avoid being biased by previous ephemerides calculated from these datasets, we assume uninformative uniform priors on all variables.To mitigate complications associated with multimodality arising from an unconstrained  0 , the prior range of  0 is centred on the seed value and bounded on either side within the seed period.We generate samples from the posterior distributions with Hamiltonian Monte Carlo (HMC) sampling (Duane et al. 1987;Betancourt 2017) and the No U-Turn Sampler (NUTS; Hoffman & Gelman 2011), using JAX (Bradbury et al. 2018) and NumPyro (Phan et al. 2019;Bingham et al. 2019).The number of burn-in and sampling steps are tuned to ensure convergence via the Gelman-Rubin split-r convergence diagnostic (Gelman & Rubin 1992).We sample four chains in parallel for 4000 steps in total, discarding the first 2000 samples of each chain as burn-in.This procedure takes just two minutes to complete on commodity hardware, running a single chain per core.
As a final step to further minimise cov(,  0 ), we repeat the process of shifting  0 by an integer number of periods, this time using the HMC samples to provide a more robust estimate of the covariances.We find a shift of -106 periods minimises this, and thus we adopt this value as  0 going forward.
Our updated ephemeris is presented in Table 2. Compared to (our corrected version of) the ephemeris from PEGS III, we achieve an factor 2.9 improvement in the uncertainty on the time of conjunction when propagated to the start of O4.In the top panel of Figure 3, we compare the two-dimensional marginal distribution of ( 0 , ) for our updated ephemeris with that for each of the previous (corrected and uncorrected) ephemerides, propagated to the epoch of the new ephemeris.Our corrections are clear in the now consistent values of the reference epoch,  0 .There are some deviations in the period , however this is expected due to the inclusion of more VLT data (with better resolution of the spectral lines) at later epochs.
In the bottom panel of Figure 3, we propagate the uncertainty on the PEGS IV  0 to past and future epochs.This again demonstrates the consistency of our corrected and new ephemerides with each other.Of course, the uncertainty in the reference epoch  0 grows in time, due to the posterior uncertainty in the ephemeris.However, our updated measurements reduce this rate of growth into the future observing runs, O4 and O5, of the LIGO/Virgo/KAGRA GW detectors, as is crucial for increasing the sensitivity of searches for CWs from Sco X-1.In Figure 4, we present the RV curve of Sco X-1 implied by our ephemeris inference.The high-resolution VLT spectra result in an excellent match with the fitted RV curve, with robust uncertainty quantification carried through our Bayesian measurements of the ephemeris.Finally, Figure 5 presents a corner plot of our samples, along with diagnostic statistics to illustrate the proper exploration of the parameter space.
For completeness, we present the ephemeris posterior and additional diagnostic quantities in Figure 5.To enable the principled calculation of search boundaries based on the full posterior distributions, we make high-quality samples available alongside this paper in the Supplementary Material section.
Examining the posteriors in Figure 5, our marginal distributions are largely Gaussian, and inter-parameter correlations are very weak.Our two-step line-search algorithm has successfully reduced the covariance between  0 and  significantly to an (absolute) value of 3.6 × 10 −15 d 2 -multiple orders of magnitude improvement over the value quoted in Galloway et al. (2014).It is further reassuring to see that both  1999  and  2011  , the systemic offset parameters for each dataset, are consistent with zero at the 1 confidence level, verifying the validity of the absolute wavelength calibration of both datasets -nevertheless their inclusion is important to account for a potentially non-trivial systematic effect on the ephemeris, and to provide correct uncertainty estimates on each parameter.
Our ephemeris provides an excellent (∼ 30s) precision on the time of conjunction T 0 , providing a significant (∼ factor 3) reduction in the required template search space (following Galloway et al. 2014).At these levels of precision, even subtle effects begin to have marked impacts on the ephemeris.It is therefore crucial to acknowledge that there are additional sources of statistical and systematic uncertainty on timings that may hamper efforts to push the ephemeris to even greater precision.Sco X-1 shows optical variability on the level of ∼ 0.5 mag, with typical variability timescales of ∼ 1 hour.As a result, we observe a flux-weighted average radial velocity over the length of our exposure, adding additional uncertainties to our velocity estimates.This is likely a small effect at the sub-km/s level owing to our comparatively short exposure times (≲ 5% flux variability over one ∼ 700s UVES integration), however may begin to be important when searching for deviations from Keplerian radial velocities -we discuss this further in Section 3.3.There are also potential uncertainties arising from the timestamps applied to each spectrum acquired -it is challenging to quantify the temporal accuracy as this is often poorly documented, and rarely validated experimentally.This is particularly true in the case of 'historic' datasets from many decades ago.In the ideal scenario, timestamps would be derived from GPS time (e.g.see Dhillon et al. 2007), or derived from system time that is updated regularly via Network Time Protocol (NTP), each of which can keep clock errors minimal.This also implies potential per-observatory time offsets, which may manifest in non-trivial ways.These effects are currently dwarfed by the statistical error present on our ephemeris, and are only likely to become problematic given a significantly increased volume of data -which will reduce the statistical errors on the ephemeris.We nevertheless caution that the 'true' uncertainties on our ephemeris may be larger than implied by our posterior samples, and encourage conservative allocation of CW parameter spaces to accommodate this.Future work will look to incorporate some of these effects into our Bayesian framework, as well as further extending the observational baseline to provide the volume of data required to model them.

Eccentricity constraints
We now consider potential deviations from the standard circular orbit model discussed above.As remarked upon in Wang et al. (2018), obtaining direct constraints on the orbital eccentricity of Sco X-1 is complicated by the Roche lobe geometry of the companion star.As it fills its Roche lobe, emission occurs from an extended surface of the donor, which continuously changes aspect ratio with orbital phase.This naturally leads to deviations from the pure sinusoidal RV curve expected from the base Keplerian orbit model, and adds additional periodic modulation on the ∼1 km/s level.This creates an apparent eccentricity dominating over -and potentially degenerate withany true orbital eccentricity, which makes it difficult to disentangle the effects of each.As with other LMXBs, there are strong physical reasons to expect the orbital eccentricity  of the system to be close to  = 0, due to tidal dissipation occurring during the main sequence lifetime of the NS progenitor (see, e.g., Tassoul & Tassoul 1992).As a further complication, the vanishing phase space (Lucy & Sweeney 1971) as  → 0 induces a bias towards non-zero orbital eccentricities, and therefore, any detection of  > 0 should be interpreted with caution mathematically also.
Bearing this in mind, we repeat the analysis of Section 3.2, but fit for a nonzero eccentricity, using a more general form of the Keplerian orbital model described previously.We adopt the common parameterisation ( √  cos , √  sin ) to sample the orbital eccentricty  and the argument of periapsis .This decorrelates the posterior at low eccentricities and make sampling easier, assuming uniform priors on both of these parameters in the range (-1, 1).We use the JAX-based kepler solver from the exoplanet package (Foreman-Mackey et al. 2021) for compatibility with the NUTS algorithm.
We empirically find an upper limit on the orbital eccentricity of Sco X-1 of 0.0132 (0.0161) at 90% (99%) confidence level.The commonly-employed significance test of Lucy & Sweeney (1971) is of limited utility here as it considers solely orbital eccentricity, and (falsely) suggests we should adopt an elliptical orbit as a result.We expect any non-zero eccentricity originates entirely from the Roche geometry, bearing in mind the tidal dissipation seen in similar LMXBs.This provides a conservative constraint of parameter space for potential CW searches; it is difficult to move beyond this without more intensive modelling efforts, involving full modelling of the binary system in both light curves and RVs.Some early progress is being made (see, e.g., Cherepashchuk et al. 2021), although it is important to note that light curve modelling carries potential systematics that must be carefully accounted for when combined with other independent constraints.Even with the high-quality UVES data here, the expected RV modulation from the distorted secondary (see Figure 6 of Wang et al. 2018, of order 100 m/s) is comparable to the statistical uncertainty, making the prospect of direct measurement unlikely -the greatest modulation occurs around phase zero, where the Bowen emission lines are weakest.Deviations from a pure sinusoid are directly informative on the inclination  (Masuda & Hirano 2021), making further reduction of RV uncertainties an important goal going forward.

BINARY PROPERTIES
As the prototypical LMXB, Sco X-1 has been the focus of intense study since its discovery in the early 1960s (Giacconi et al. 1962;Chodil et al. 1965).However, the high (super-Eddington) accretion rate and its effect on the Balmer lines means there is still uncertainty surrounding the structure of the accretion disc in Sco X-1  Corner plot of posterior samples for the  = 0 ephemeris, coloured according to each independent chain to confirm randomly-initialised chains converge on the same posterior mode.Inset: Histogram showing the distributions of marginal energies against the transition energies for each step taken by the NUTS sampler -the close match between these two distributions implies full and efficient (low autocorrelation) exploration of the parameter space.The estimated Bayesian fraction of missing information is close to 1 for all chains, providing a quantitative verification of this also.These diagnostics are unique to Hamiltonian Monte Carlo and provide an orthogonal check to the usual split-r diagnostics employed in MCMC that we use above.Refer to Betancourt (2016) for a full theoretical explanation, and further details.
-whether the observed form is stable over long timescales, and how this correlates with the various X-ray and optical states of the system (as identified in Scaringi et al. 2015).The long-term, highresolution dataset collected as part of the PEGS program provides the means to probe secular evolution of the disc, as well as search for period derivatives and other phenomena not detectable with the typical dense but short-baseline observational sampling presented in previous work, although one pays a penalty in the resolution of fine structure owing to the dynamic and changing disc state.We expand on this in the subsections below.

High-resolution, high-SNR line atlas for Sco X-1
With over 200 high-resolution spectra of Sco X-1, and a highprecision ephemeris, we construct a 'line atlas' by combining the spectra in the donor-star rest frame.By doing this, weak spectral lines not visible in the individual spectra can be detected, and the quality of all lines improved.The variability in the broad H/He line profiles leads to some smearing of these lines.However weaker, narrow lines are not affected as heavily.A significant number of VLT/UVES spectra taken as part of PEGS also have a simultaneous observations with the red arm, which, although of limited utility for constraining the binary motion due to the lack of strong and narrow spectral features, encodes information about the cooler material present in the Sco X-1 system.We include this in our line atlas primarily for completeness, although caution that telluric subtraction was not performed as part of data reduction, and so the red arm suffers from residual telluric contamination.This is somewhat suppressed by our shift-adding scheme.
We use the measured Bowen line velocities to shift all spectra into alignment, and normalise out the continuum flux by fitting Chebyshev polynomials (7th order for the UVES Blue arm, 5th for UVES Red arm CCD 1, 1st for UVES Red arm CCD 2), rejecting outlier points to avoid fitting spectral features.The heavy telluric contamination redwards of 9000 Å necessitates the use of a low-order polynomial to avoid erratic ringing effects.Spectra are then reinterpolated onto a common wavelength grid, and median-combined with a sigmaclipping outlier rejection scheme to remove any cosmic ray hits or defects.This process is performed separately for the blue arm and the red arm, and we take care to treat the two red CCDs separately to avoid discontinuities.Note that the continuum normalisation at the far-red end of the spectrum is hampered by the strong telluric bands, so we opt for a simple linear continuum only.
Our finalised line atlas is presented across Figure 7 (CD#2, blue) and Figure 8 (CD#4, red) -we achieve a median SNR of ∼ 1000 in the blue arm, marginally above the expected

√
scaling owing to our use of sigma-clipping, but likely affected by line variability.

Doppler tomography
We applied Doppler tomography (Marsh & Horne 1988;Marsh 2005) to our phase-resolved UVES spectra to probe the structure of the accretion disc surrounding Sco X-1.Our spectra are pre-processed to remove continuum flux using the same procedure as Section 4.1 and normalised in flux over the line region of interest to mitigate the impact of flux variability on the final Doppler map.We use the Python bindings of the DOPPLER1 code to produce the Doppler maps.The period and systemic velocity are fixed according to the median of our ephemeris posteriors (see Table 2), and we compute the Doppler map over a square grid 500 km/s in side length, with a per-pixel resolution of 1.5 km/s in the x-y disc plane velocities   and   .We neglect any velocities in the  plane for computational reasons, but note that our spectra are unlikely to be constraining along this axis (see Marsh 2022).Doppler maps are constructed for the H, H, He II and Bowen lines individually.These are shown in Figure 9, alongside the target  2 value used for the maximum-entropy optimisation procedure.Despite the marked resolution improvement afforded by VLT/UVES over previous datasets, the output Doppler tomograms are of similar quality to the previous maps presented based solely on WHT data in Wang et al. (2018).We speculate that this likely arises to long-term seasonal variations in the structure of the accretion disc around Sco X-1, and explore this possibility further in the following section.

Secular variability in the He II disc?
To explore the possibility of seasonal or secular variability in the disc structure of Sco X-1, we compute He II Doppler maps on subsets of the full UVES dataset presented in this paper.We manually group contiguous subsets of our data into five 'seasons' -selected roughly corresponding to the ESO observing semesters.We then compute individual Doppler maps for each season, to avoid averaging over any long-term variability.The target  2 is adjusted per-season to avoid over-resolving disc features.The resulting Doppler maps are shown in Figure 10.
We focus our attention on the He II Doppler maps, as these show a clearly-resolved disc structure.Some tenuous variability is present in the Bowen disc also, although this is hard to quantify or comment on phenomenologically, given that the donor star feature (at   = 0,   = 76.4 km/s) dominates the reconstruction.In the He II maps, we see a strong asymmetric emission component present in the accretion disc.Across multiple observing seasons, we see evidence for evolution in intensity of this 'rim'-like structure -with minimal intensity in Season 2 and maximal intensity in Season 5. We interpret this feature as an extended region of gaseous material suspended above the disc and corotating with it -potentially a signature of stream-disc overflow occurring in the system (see, e.g., Kunze et al. 2001).Given the high accretion rate of Sco X-1, this is not surprising, as the accretion process is stochastic and variable.As our Doppler maps are averaged over many orbital periods, this is likely to be a longer-term secular variation rather than an effect correlated with the disc state.This may potentially coincide with periods of enhanced accretion onto the NS, as implied by short-(∼ weeks-months) timescale fluctuations in the X-ray luminosity seen in MAXI light curves (Hynes et al. 2016).These disc states are wellknown (e.g., Scaringi et al. 2015), but we do not resolve these with our sparsely-sampled VLT/UVES dataset.We encourage additional observations and characterisation of these secular variations in the disc structure, preferentially via intensive spectroscopic observations over blocks of a few nights, every few months.This ensures each block has adequate phase coverage, and minimises the blurring of structure caused by a more coarse observing program such as our VLT/UVES data.Simultaneous Xray coverage is crucial to correlate structural changes with potential modulation of accretion rate, and to understand the disc state.As the prototypical LMXB, targeted studies of this phenomenon in Sco X-1 have the potential to provide insight into the longer-term dynamics of accretion across the domain of X-ray binaries and other high accretion rate compact binaries.

Behaviour of the Balmer lines
Complex absorption features are present in the Balmer lines -as reported by Scaringi et al. 2015 -reminiscent of P Cygni-style outflow profiles with red/blue absorption wings superimposed on a broad central emission feature.These occur antiphase with the overall optical brightness of the system (and with the Balmer line intensity itself), showing maximum absorption when the Balmer lines are strongest and no absorption when the Balmer lines are barely visible.Both red and blue components are visible, showing considerable variations in absorption depth.Given the strong time variability demonstrated in previous sections, it is challenging to probe this with our VLT/UVES dataset.Nevertheless, we present some broad overview properties below in the hope of encouraging further characterisation.Such profiles have also been noted previously in NIR spectra of the system on the Br line, e.g.see Bandyopadhyay et al. (1999).We focus on the H line as this shows the strongest absorption (see Figure 11)and thus provides the most robust velocity measurements, although this phenomenon is also visible in H  , .We apply a local continuum normalisation to bring all spectra onto a common baseline, then fit a Chebyshev polynomial to the line wing to provide a smooth approximation to the spectrum on this specific region.We find the minimum value of this polynomial and use bootstrap iterations to measure the uncertainty on our minimum wavelength value.To reject spectra where no absorption is present and thus ensure our velocity measurements are robust, we only consider spectra from this point on where there is ≥ 3% absorption in this red wing.
Under the assumption that this is an outflow, the 'apparent' velocity varies between 100 km/s-1000 km/s, with variable absorption depth as remarked upon above, and appears to show some correlation with orbital phase.However, such a phase-dependent morphology could equally be replicated with the addition of a broad absorption component superimposed on the emission component (see Figure 11).To illustrate this, we fit a toy model to one of our spectra showing high absorption, composed of a narrow emission line centred on H, and a broad (at least 2x wider than the emission component) absorption line with a free mean value.We note that the true structure of the Balmer lines is markedly more complex than this, and this simplified model is intended to be illustrative, rather than fully reproducing all features of the data.
Owing to degeneracies between the strength of the  line and this component it is not possible to constrain the amplitude of such absorption, but one solution suggests a width of around 8× that of .This could be created by optically-thick absorption of the inner disk surrounding the NS in Sco X-1, potentially probing some NS-driven outflow that may act as a crude tracer of the true NS radial velocity amplitude.This component appears to have some systemic velocity offset compared to the Balmer line in emission.It relies largely on serendipity to observe the disc in the right state to make these measurements, and would benefit strongly from simultaneous X-ray observations, to this end.Disentangling the origin is not possible with the sparse sampling of the VLT/UVES data, and we therefore advocate for additional intensive monitoring to investigate the nature of this absorption.

Constraints on other system parameters?
Despite the marked improvement in timing uncertainties offered by our reanalysis and extensive dataset, it is difficult to further constrain the system's orbital parameters, which are important in CW searches.The uncertainties on the other system parameters of interest are dominated by the poor constraints on the mass ratio  and the unknown neutron star radial velocity  1 , that are difficult to further improve on.The additional VLT/UVES data provide no stronger constraint on  1 owing to the variability we discuss in Section 4.3 blurring the Doppler tomograms.Similarly, we cannot tighten the lower bound on  as the lines are already fully resolved by our spectra.We recreated the Monte-Carlo analysis of Wang et al. (2018), casting the problem instead as a Bayesian forward-modelling approach (using the HMC methods discussed in Section 3.2), but this failed to provide any stronger constraints, owing to the above considerations.
Further observational work is required to better constrain these orbital parameters, yet the complexity of Sco X-1 makes this challenging -detection of donor star features in the near-infrared is hampered by the strong accretion flux (Mata Sanchez et al. 2015), and light-curve modelling is made more difficult by the distinct optical high and low states in the system.Only by jointly considering all available constraints, and pushing the capabilities of currentgeneration instruments, can we begin to more robustly constrain the orbital parameters of Sco-X 1 and further reduce the parameter space for CW searches going forwards.

CONCLUSIONS
As the prototypical LMXB, Sco X-1 continues to remain central to current searches for continuous GW sources.Tensions between the predicted GW emission and the GW strain limits obtained from previous observing runs are beginning to emerge, as specific subbands now reach below the torque-balance between accretion spinup and GW emission spin-down.However, this is largely contingent on the inclination constraint obtained from Fomalont et al. (2001), and further data is required to reach constraining upper limits whilst simultaneously marginalising over the unknown inclination; it may be the case, e.g., that the radio lobes do not trace the orbital inclination.
Although any GW emission thus far has eluded detection, the upcoming LIGO/Virgo/KAGRA O4 observing run promises to further improve instrument sensitivities.Furthermore, searches with increased sensitivity (e.g., Mukherjee et al. 2022) may yield more stringent constraints over previous results.With upcoming thirdgeneration GW detectors such as the Einstein Telescope (Maggiore et al. 2020) and Cosmic Explorer (Reitze et al. 2019) -bringing orders of magnitude increase in strain sensitivity and delivering high SNR GW detections that will unveil populations of compact objects currently out of reach for the current ground-based detectors (e.g., Cieślar et al. 2021) -likely including Sco X-1, alongside a wealth of other science outcomes (Kalogera et al. 2021).
From a compact object perspective, further studies of Sco X-1 focused on the transient structure of the disc and the nature of the companion star is crucial in underpinning our understanding of more distant LMXBs, and may yield even stronger constraints that may synergistically further constrain the parameter space for GW searches.It is clear many processes remain poorly understood.
Leveraging the improved search sensitivity afforded by the enhanced detectors and computational methods discussed above is contingent on a precise -and, critically, up-to-date -ephemeris for the donor star in Sco X-1.The ephemeris presented in this work will enable precision searches for CW emission from Sco X-1, which will further constrain any emission down to the torque-balance limit across the entire band of search frequencies -both in the upcoming LIGO-Virgo-Kagra O4 observing run and beyond.Sparsely-sampled, high-resolution observations of Sco X-1 over the coming years can efficiently keep this ephemeris current for the foreseeable future, underpinned and facilitated by the extensive and corrected VLT/UVES constraints that we present here.

Figure 1 .
Figure 1.Trailed spectrogram of the VLT/UVES data presented in this work, folded by the best-fit ephemeris derived in Section 3. Centred on 4640 Å is the Bowen line region of interest, and 4686 Å is the He II line.We restrict our analyses in this paper to this specific region of the spectrum, although other weaker Bowen lines are present in the spectrum.

Figure 2 .
Figure 2. Example plot of the Bowen region in one UVES spectrum after continuum subtraction, with the Bowen line model discussed in Section 2.2 overplotted along with the individual components.Note that one of the fitted line components is not visible in this spectrum so has zero amplitude.We take care to avoid including the bright HeII 4686 Å line.Individual spectra typically have SNR ∼ 50.

Figure 3 .
Figure 3. Marginal two-dimensional posteriors of  0 and  for all previous ephemerides, propagated forwards to the epoch of our ephemeris (see Table2).The contours show the 1-, 2-, and 3- confidence intervals respectively.

Figure 4 .
Figure 4. Top panel: RV curve for Sco X-1.We show the 1 confidence intervals of the RVs measured from the 1999 WHT (blue), 2011 WHT (orange), and VLT (green) spectroscopic datasets, derived from the Bowen line model.We also plot the median posterior predictive RV curve implied by our ephemeris constraints (red) and the corresponding 1, 2, and 3 confidence regions (darker to lighter red shaded regions, respectively).Note the significantly reduced intrinsic scatter of the higher resolution VLT data with respect to the older WHT data.Bottom panel: Normalized residuals between the observational RV measurements with the median RV curve, scaled by the respective error factors Figure 5. Corner plot of posterior samples for the  = 0 ephemeris, coloured according to each independent chain to confirm randomly-initialised chains converge on the same posterior mode.Inset: Histogram showing the distributions of marginal energies against the transition energies for each step taken by the NUTS sampler -the close match between these two distributions implies full and efficient (low autocorrelation) exploration of the parameter space.The estimated Bayesian fraction of missing information is close to 1 for all chains, providing a quantitative verification of this also.These diagnostics are unique to Hamiltonian Monte Carlo and provide an orthogonal check to the usual split-r diagnostics employed in MCMC that we use above.Refer toBetancourt (2016) for a full theoretical explanation, and further details.

Figure 6 .
Figure 6.Plot of uncertainty in T 0 as a function of epoch for all updated ephemerides presented in this work.The black bars overplotted denote the times of individual spectra contributing to the ephemeris.Our ephemeris provides a factor 2 improvement in uncertainty for the O3 observing run, growing out to later times.

Figure 7 .
Figure7.Co-addition of all VLT/UVES spectra in the binary rest-frame, to produce a line atlas that reveals the presence of many low-strength emission lines.Spectra are broken into chunks of 400 Å to aid visualisation.We annotate potential lines of interest to the LMXB community.

Figure 8 .
Figure8.Line atlas for Sco X-1 in the binary rest frame continued from Figure7-this figure shows the UVES red arm spectra, covering 5692 Å-9460 Å.Note that no telluric correction has been applied, and we instead use the binary motion to median the weaker features out.Some contamination persists in the red, hence we only attempt to identify the strong spectral lines of H and He, with the weaker C/N/O lines obscured by atmospheric absorption.A full treatment of telluric absorption is required to identify all faint lines present, but this is beyond the scope of this paper.UVES has a chip gap covering approximately 7540 Å-7660 Å in this configuration.Spectra are broken into chunks of 1200 Å to aid visualisation.

Figure 10 .
Figure10.Doppler tomograms for each observing 'season', using the method defined in Section 4.3.We compute season-averaged tomograms for the Bowen lines (top panel) and the He II lines(bottom panel), to search for long-term secular variations in the disc structure.The Bowen Doppler maps are largely dominated by the strong donor star signature, but some structure is visible.The red and white markers indicate the positions of the donor star, and the system centre-of-mass respectively.

Figure 11 .
Figure 11.Spectrum of Sco X-1 showing strong absorption wings on the H line.Overplotted is a two-Gaussian model that reproduces the line shape shown, with one positive component representing the broad H line, and one constrained-negative component reproducing the broad absorption wings present.

Table 1 .
Summary table containing all observations included in this ephemeris version, along with instrument and program IDs. 'Group' here refers to the fitting group introduced in Section 3.2 -datasets in the same group are fitted with common error scaling and velocity offset parameters.

Table 2 .
Tabulated posterior distribution summary statistics for our ephemeris.Values correspond to the posterior median, while uncertainties represent the marginal 1 confidence intervals.Explicitly, the value of  0 refers to the inferior conjunction of the companion star, and  asc,ns refers to the time of ascending node crossing for the NS, moving away from the observer.The top rows assume zero eccentricity.Our upper limit on the eccentricity is given in the bottom row, as the 90th percentile of the marginal eccentricity distribution.Accompanying this table is a corner plot of posterior samples in Figure5.