© The Authors 2024

1 Introduction

Low-metallicity massive stars remain poorly understood. The lack of detailed understanding of the global, photospheric, and wind properties of these massive stars, and of the physical mechanisms operating in their interiors and atmospheres directly propagate to topics as diverse and important as the Universe’s chemical enrichment in heavy elements (Kobayashi et al. 2006; Hansen et al. 2014), the formation of high-ionization emitting gas (Kashino et al. 2023), the feedback of massive stars on the local environment (Doran et al. 2013), host galaxies and reion-ization of the early Universe (Codoreanu et al. 2018), and on the astrophysical interpretation of gravitational wave detections (Belczynski et al. 2010; Abbott et al. 2020). One of the reasons of the current limitations stems from the lack of a sufficiently high-quality observational database of a broad set of massive stars in various low-metallicity environments.

With typical metallicities of 0.5 Z_⊙ (Hunter et al. 2009) and 0.2 Z_⊙, respectively (Bouret et al. 2003; Rickard et al. 2022), the Large and Small Magellanic Clouds (LMC and SMC) are our closest windows to the subsolar metallicity environment. Furthermore, both galaxies are close enough for individual massive stars to be studied in great detail. In this context, the ULLYSES Director’s Discretionary program (Roman-Duval et al. 2020, 2023) dedicated 500 orbits of the Hubble Space Telescope (HST) to obtain high-quality far-ultraviolet (FUV) spectroscopy of a representative set of OB and WR stars in the Magellanic Clouds. The initial ULLYSES sample was further complemented by archival HST observations to reach a grand total of 130 LMC and 124 SMC targets.

The X-Shooting ULLYSES project (XShootU) is an accompanying ESO Large Program whose aim is to provide ground-based spectroscopy of all ULLYSES targets (Vink et al. 2023, hereafter Paper I). Using the ESO X-shooter spectrograph (Vernet et al. 2011), XShootU complements the FUV HST observations with homogeneous intermediate-spectral resolution spec-troscopy, covering the optical and near-infrared (NIR) domains. While a subset of the ULLYSES targets had existing optical spectroscopy prior to XShootU, the spectral coverage and data quality was far from homogeneous.

Optical spectra are needed to enable the full potential of the ULLYSES campaign. Stellar parameters such as surface gravity (log g) or projected rotation rates (υ sin i) are difficult to reliably constrain from FUV spectra alone, which translates into significant uncertainties, for example on stellar wind parameters and surface abundances. These, however, are important quantities to constrain theories of stellar structure, evolution, and outflows which constitute some of the most fundamental building blocks of our current understanding of the multi-messenger Universe.

Similarly, high-quality libraries of massive star spectra and the physical parameters of the associated stars, are crucial ingredients for population synthesis tools that are used to analyze distant populations of massive stars observed in integrated light. Specifically, they are needed in the context of the current and upcoming generation of new telescopes that probe the distant, low-metallicity universe, such as the James Webb Space Telescope (JWST) and upcoming 30 m class telescopes.

The goals of the XShootU project have been described in Paper I, including a top-level description of the target sample and the observing campaign. In this paper we present the uniform reduction and calibration of the optical X-shooter data of the Magellanic Cloud stars. We provide advanced data products that can be used for numerous science applications, from detailed atmospheric analyses to the creation of stellar template libraries. The NIR data reduction requires specific attention and will be presented in a subsequent paper in the XShootU series.

We would like to provide a note on the naming convention used in this paper. While Table 1 in Paper I provides all target names, these are often not suitable for electronic filenames (e.g. they contain spaces or special characters). For this reason and for consistency with the FUV ULLYSES dataset, we opted to use the MAST target names for all fits files in the XShootU data release. An overview of Paper I names and corresponding MAST target names is given in Table A.1.

This paper is organized as follows. Section 2 provides a detailed description of the observing campaign and of the collected RAW data. Section 3 describes the data reduction process, including nonstandard procedures for the flat-fielding and instrument response correction. Slit-loss correction, absolute flux calibration, normalization, telluric correction, resampling and co-addition of the final science products are described in Sect. 4. Section 5 provides the radial velocities (RVs) while Sect. 6 describes the advanced data products that are made available to the community.

2 Observation campaign

2.1 Instrumental setup

The ESO X-shooter spectrograph uses two dichroics located after the focal plane. These dichroics split the light into three arms (UVB, VIS, and NIR), with a separation between the arms at λ ≈ 560 nm and 1024 nm (Vernet et al. 2011). In this paper, we focus exclusively on the spectra obtained with the UVB and VIS arms. Both arms are equipped with 4k × 2k detectors, covering 12 and 15 spectral orders, respectively.

The X-shooter spectral resolving power (R = λ/Δλ) is set by the width of the entrance slits. The XShootU campaign used slits of 0.″8 for the UVB arm and 0.″7 for the VIS arm, corresponding to R ≈ 6700 and 11 400. In this setup, the point spread function (PSF) along the dispersion direction is well sampled, with 5 and 4.4 RAW pixels across the full width at half maximum of the UVB and VIS arms PSFs, respectively. Similarly, the UVB and VIS detectors provide a RAW pixel scale of ≈ 1.″6 pix⁻¹ in the spatial direction. The slow, high-gain readout mode was chosen for both detectors with no binning (readout mode: 100,lxl,hg; gain ≈ 0.63 e⁻ ADU⁻¹), delivering a low readout noise of 2.4 e⁻. The dark current level is extremely low (<2 e⁻ pix⁻¹ h⁻¹) and is neglected. Both UVB and VIS detectors have excellent cosmetics.

Each XShootU target was observed in STARE mode using sequences of observation (OBs) composed of one exposure or, for the fainter objects, two back-to-back exposures with identical exposure times. The individual exposure times vary between 100 s and 2955 s according to the brightness of the stars, with the aim of obtaining a signal-to-noise ratio per resolution element (S/N) of at least 100 at 490 nm and 600 nm for all LMC and SMC targets. No dedicated night-time calibrations were requested by the XShootU program, such that we rely on the standard X-shooter calibration plan. The latter typically offers one flux standard star and one telluric standard star observed during the same night as the scientific observations.

An overview of the observing campaign is provided in Table 1 while the full journal of observations is in Table 2. The observations were spread over 102 different nights from October 2020 to January 2022 and include 129 LMC and 103 SMC stars. In total, 358 UVB and 357 VIS individual exposures were obtained. Seventy-three objects have multi-epoch observations (i.e., two or more exposures), of which 34 have these observations spread over different nights. The main reasons for these repeated observations over two nights or more, are due to operators declaring that the weather conditions were below requirements or because part of the spectrum was saturated. These repeated observations provide a limited insight into the variability of these objects.

2.2 Quality control of the RAW data

We inspected each RAW data frame and checked for (near) saturation (>60 000 ADU). We identified one target with very strong saturation (VFTS 482) and 21 other targets where the Hβ emission line or nebular lines were saturated. For those targets, observations were repeated with reduced exposure times to record the full spectrum without saturation. However, all data were reduced in the same way and saturated parts of the spectra were masked in the co-addition (see Sect. 4.5).

During this quality control process, we identified the presence of an undocumented ghost that affects about 15% of our observations. The ghost appears as a bright, circular spot of typically 7 pix in diameter (Fig. 1, left panel). In most cases, the spot appears on an unused part of the detector (outer part, inter-order regions). However, in 3% and 1% of the frames, it landed, respectively, on a region used for the sky or overlapped with the science spectrum. In the latter case, the spot produces a clear artifact in the reduced 1D spectrum (Fig. 1, right panel). In most cases, the ghost is saturated, so that it can easily be identified in the cumulative distribution of counts in each pixel. However, we have also visually identified frames where the ghost was not saturated, which made it more difficult to automatically flag. The origin of the ghost has not been identified. Some back-to-back exposures separated by 10–15 min reveal that the ghost moves on a short timescale or even disappears. In a handful of cases, multiple ghosts could be identified in the RAW image.

Fig. 1 Example of the ghost artifact identified in various 2D RAW images (left panel). Whenever the ghost position overlaps with the extraction region of the science target, it produces an artifact in the extracted ID spectrum (right panel, at λ ≈ 910 nm in this example). The insets are in the same units as the main panel.

Fig. 2 Comparison of a reduced spectrum of SK-68D137 without (left panels) and with (right panels) response-curve correction.

3 Data reduction

3.1 General overview

The data reduction was performed using the ESO X-shooter pipeline v3.5.0 (Goldoni 2011). The following steps were applied using dedicated routines of the pipeline: bias subtraction, order guess position estimate, order tracing, flat-fielding, wavelength calibration, spectral extraction, sky-subtraction and response-curve correction. No dark correction was applied given the very low dark current of the UVB and VIS detectors. We decided early on to work with response-curve-corrected spectra as these provide better correction for the blaze functions and an overall smoother spectrum, which eases the rectification to the continuum (Fig. 2). For the vast majority of the observations, RAW calibration frames obtained in the mornings after each observation night were adopted. When these were not available, suitable calibration frames nearest in time were used.

While the general data reduction strategy closely follows the standard recommendation of ESO, the following sections detail issues that were identified and the mitigation strategies implemented. This section focuses on the data reduction with the ESO CPL pipeline recipes. Additional calibration steps, such as slit-loss correction, absolute flux calibration, telluric correction, and normalization, are addressed in the next section.

Fig. 3 Comparison of reduced spectra of SK-68D137 obtained with the different flat-fielding strategies when the SCIENCE target and the FLUX standard star have been observed on different nights. Red indicates when different FF frames were used to reduce the observations of the SCIENCE and FLUX stars. Black indicates when the same FF frames were used.

3.2 Flat-fielding

Flat-field (FF) frames are used in the X-shooter calibration chain to correct for the pixel-to-pixel response variation and to determine the blaze functions of the échelle orders. Under the standard X-shooter calibration plan, the validity of these FF calibration frames is three days and they are thus not taken daily. The default reduction procedure further uses FF frames with matching slit widths to reduce the SCIENCE and FLUX observations (i.e., 0.8″/70.7″ for the UVB/VIS SCIENCE frames in our case, and 5″ for FLUX frames).

While this is a standard approach, the net result is that different FF frames are used to reduce the FLUX and the SCIENCE observations, so that any difference in the pixel-to-pixel or blaze-correction between the FLUX and the SCIENCE data reduction leaves a residual imprint in the flux-calibrated spectrum. On average, the effect is negligible but might become noticeable for some orders close to the edge of the detectors, especially when the FF-frames used to reduce the FLUX and SCIENCE observations are taken on different days (Fig. 3).

To investigate the most suitable flat-fielding strategy, we compared different approaches. In addition to the standard approach described above, we used the same broad-slit FF frames to reduce both the FLUX and the SCIENCE frames. We also used the narrow-slit FF frames corresponding to the SCIENCE observations to reduce the FLUX data. The size of the slit has a negligible effect on the outcome of the procedure. However, using the same FF frames to reduce the FLUX and SCIENCE observations improves the determination of the blaze function (Fig. 3). Our final approach used the broad-slit 5″ FF frames obtained for the FLUX observations to reduce both the FLUX and the SCIENCE frames. From our tests, using FF-corrections obtained with different slit widths has no significant impact as long as the FF slitwidths is at least as large as the smallest one used for either FLUX or SCIENCE. This is logical as the pixel-to-pixel sensitivity correction is performed at the pixel level on the 2D images. However, forcing the FLUX and SCIENCE frames to be corrected by the same FF calibration guarantees that the imprints of the FF correction residuals of the FLUX and SCIENCE frames cancels out when correcting for the response curve of the instrument (see Sect. 4.3).

Fig. 4 Radial velocity difference measured between the UVB and VIS arms after correction for known wavelength calibration offsets (Sect. 3.3). The panels show the values for different data quality (increasing from top to bottom). The UVB+VIS quality of a star is given by the sum of all points listed in Table 5. The median and 68.3% quan-tile are listed at the top of each plot.

3.3 Wavelength calibration

The wavelength calibration performed by the standard ESO pipeline uses a physical model (Moehler et al. 2010). The transformation from pixel to wavelength space is optimized through the analysis of multi(9)-pinhole mask ThAr (UVB, VIS) lamp frames, instead of ThAr observations using the science slits. The predicted positions of the lines obtained through the multi(9)-pinhole mask are fitted using a two-dimensional Gaussian to recover the actual positions on the frame. The resulting wavelength calibration is slightly offset in wavelength compared to observations through the science slits as different optical components are in the light path. These offsets can be recovered using quality control information of ThAr slit observations, present in the ESO QC1 database¹. We found that we have to adjust² the UVB spectra by adding 4.477 km s⁻¹, and the VIS spectra by adding 1.001 km s⁻¹.

After correction for these offsets, RVs measured from the UVB and VIS arms come into much better agreement (see also Sect. 5). Figure 4 is based on a preliminary determination of the radial velocities that processed the UVB and VIS arms separately. It indeed shows that the systematic residual offset is well below 1 km s⁻¹, independent of the quality of the spectral lines of the object and of the S/N of the data. This amounts to about 10% of the statistical uncertainties on these measurements. The larger radial velocity differences for the lower-quality data are mainly due to emission lines, which are either stellar (in WR stars), or nebular.

3.4 Updated flux standard models and instrument response

Following the flat-fielding, the extracted spectra are corrected for the instrumental response curve using a set of seven x-shooter spectro-photometric white dwarf standard stars (GD 71, GD 153, EG 274, FEIGE 110, EG 21, LTT 3218, and LTT 7987) observed during the same or adjacent nights. The standard reduction procedure in the public ESO pipeline (Moehler et al. 2014) resulted in inconsistent spectral energy distribution slopes between the standard stars and small “wiggles” around the Balmer lines (on the order of a few percent), which can affect the measurement of vital stellar parameters such as temperature and surface gravity. We updated the reference models and the knot anchors points (Fig. 5) used in the corresponding pipeline recipe to reduce these wiggles and boost the consistency in the relative flux calibration (Fig. 6).

Of the seven stars used as X-shooter standards, two (GD 71 and GD 153) are also used as standards for the HST. For consistency with the flux-calibrated HST spectra, we used the HST reference models for GD 71 and GD 153 as inputs for the X-shooter pipeline (Bohlin et al. 2020). To correct the models of the remaining five X-shooter standards, we collected all their observed data between October 12, 2020 and April 21, 2021, and reduced them as if they were science frames, using the nearest GD 71 observations as the flux calibrator and the HST reference models for the reduction. In total, 110 such observations were available (4 for EG 274, 23 for FEIGE 110, 33 for EG 21, 38 for LTT 3218, and 12 for LTT 7987). We co-added all reduced spectra for each of these standards and performed an interpolated grid-based chi-square minimization to obtain the best-fitting synthetic white dwarf model for each.

Given the parameter range covered by these five standard stars (T_eff ≈ [9, 45] kK, log(g) ≈ [5, 8.5]), the master grid used for the fitting consisted of models taken primarily from the WDGrid (Bohlin et al. 2020), supplemented by models from TMAPGrid2 (Werner & Dreizler 1999; Werner et al. 2003; Rauch & Deetjen 2003), and from a grid computed in Koester (2010). After the best fitting model was determined, we followed the procedure outlined in Moehler et al. (2014) to remove the large-scale bumps in the residuals that result from missing physics in the synthetic models using the SciPy UnivariateSpline method (Virtanen et al. 2020). Finally, we adjusted the knot points used in the flux calibration pipeline so that they are more symmetric around the Balmer lines in order to reduce the wiggles initially seen in the reduction. Figures 6 and 7 illustrate the improved consistency obtained with the new flux standard models. The updated static calibration files including the new flux models and anchor points is provided on Zenodo (files xsh_response_fit_waves_cat_UVB, and xsh_star_catalog_UVB and _VIS).

Fig. 5 Comparison of knot points between the original knots (in gray) and the new knots (in red) used for interpolation in determining the response curve for GD 71 in the UVB arm. The top panel shows the observed spectra as a function of wavelength for one observation, while the bottom panel shows the reference model divided by the observation in black and the resulting response curve after interpolation in orange.

Fig. 6 Comparison of new and old reference models for the EG 274 and EG 21 FLUX standard stars used by the XShootU program. The top panel of each star shows the observed co-added spectrum (in black) and the new and old reference models (in red and blue, respectively). The sub-panels show the residuals between the observation and each of the reference models using the same color scheme.

3.5 Sky-subtraction

The advanced products that we provide are sky-subtracted: the average flux level of an ‘empty’ region of the detector has been subtracted from the flux of the science object. According to this definition, the ‘sky signal’ contains the illumination of the night sky (moon, dark sky brightness) as well as any nebulosity or extended emission that will fall on the slit. We effectively work under the assumption that these latter components are not variable. While all the advanced products provided are sky-subtracted, we provide the sky spectrum used so that it is trivial for any user to reconstruct the original, non-sky-subtracted signal. We briefly describe below the sky-subtraction algorithm as well as the procedure applied when multiple objects were visible in the 2D non-sky-subtracted spectra.

The sky level is estimated from regions outside the object mask and away from the detector edges by more than 0.″5 (approx. 3 pix). The object mask extends over ±2″, thus ≈25 pixels in width around the centre of the slit. This leaves two regions of approximately 3″ above and below the science spectrum used estimating the level of the sky. The selected sky pixels are then tabulated as a function of wavelength and the 2D spectrum is then smoothed along the wavelength axis using a running median over 51 pixels.

For 10 targets, one or more additional objects are present in the slit beside the main target. These contaminants were identified by visual inspection of the RAW images and their cross-dispersion profile. For these objects, the sky regions were selected manually upon visual inspection of the reduced 2D spectra. The adopted sky regions are given in Table 3. An example is shown in Figs. 8 and 9. Additional comments are provided in Appendix B.

Fig. 7 Comparison of the Hγ (top) and Hβ (bottom) line profile of the target AV-238 for three different observations (Dec. 2020, May 2021, and June 2021) and calibrated with three different FLUX standard stars (black: EG 21, blue: LTT 9789, red: GD 153). The residual plots show the ratio of data calibrated with LTT 7987 (blue) and GD 153 (red) and those calibrated with EG 21. The 1 and 2% inconsistencies in the blue wing of Hβ and Hγ, respectively, disappeared when using the new flux models. The remaining differences located in the core of the lines are likely of astrophysical origin (e.g., line profile or RV variations, or nebular correction residuals).

Fig. 8 2D reduced spectrum of AV-251 (left) and AV-6 (right), revealing a faint companion 1″ below the spectrum of AV6. The overlaid grid indicates the wavelength axis (horizontal, in nm) and the spatial axis (vertical, in arcsec from the center of the slit). The default science and sky extraction windows and those adopted for AV-6 in this work are overlaid on the 2D spectrum. The color bar indicates the ADU level in the images (in log scale). The cross sections along the spatial direction are in units of ADU (in linear scale).

Fig. 9 Comparison of the extracted spectra of AV-6 with the default (black) and manually selected extraction and sky regions for the components centered at +1.05″ (blue) and −0.8″ (red). The left and right panels correspond to the UVB and VIS arms, respectively. The insets focus on various spectral lines of interests and are in the same units as the main panels.

3.6 Cosmic ray hits

The ESO CPL pipeline allows cosmic ray hits to be flagged by using a Laplacian edge detection algorithm (van Dokkum 2001). The default detection threshold (–valueremovecrhsingle-sigmalim=20) resulted in a significant number of cosmic ray hits retained in the extracted 1D science spectra. Experimenting with different values, we opted for a detection threshold at 10σ (–valueremovecrhsingle–sigmalim=l0). Figure 10 shows a comparison of the extracted 1D spectra with the default and with our adopted value.

Fig. 10 Comparison of reduced spectra of AV-251 obtained with the different cosmic rejection thresholds. Black: –valueremovecrhsingle– sigmalim=20 (default). Red: –valueremovecrhsingle-sigmalim= 10 (adopted). The spectra are slightly shifted for clarity. The insets focus on spectral lines of interest and are in the same units as the main panel.

3.7 Order re-connection

Order re-connection in échelle spectrographs are a frequent issue and we paid particular attention to check the quality of the extracted and merged 1D science spectra. Poor order re-connection not only complicates the normalization process, but may also induce significant biases in the science analysis if the impact on the shape or strength of spectral lines is unidentified. In the X-shooter pipeline, order reconnection is implemented as a weighted mean in the overlapping parts of the orders. The weight are defined as the inverse of the variance, hence 1/σ² where σ is the signal-to-noise ratio. This approach is statistically correct in the case of a normally distributed noise and no systematic biases. In this context, we visually inspected each re-connection region (see Fig. 11 for an example) and concluded that the order re-connection in the response-corrected spectra was excellent and void of any systematics.

4 Advanced data products

4.1 Slitloss correction

The XShootU spectra were obtained with narrow slits (0.″8 and 0.″7 for the UVB and VIS arms, respectively). The width of the entrance slits is thus of the same order of magnitude as the median seeing value at Cerro Paranal (Sarazin et al. 2008) and thus a fraction of the star PSF does not enter the spectrograph slits. The FLUX standard star observations used to determine the response curve correction (Sect. 3.4) were observed with 5" wide slits such that the entire star PSF is transmitted through the slit. To obtain absolute flux-calibrated spectra, we corrected the response-curve-corrected spectra for slit losses due to seeing and image quality across the detector. The correction for slit-losses is wavelength dependent, because the seeing value (σ_λ) changes with wavelength.

The seeing in Paranal is measured in real time by a differential image motion monitor (DIMM). Using the DIMM values recorded in the ESO archive, we estimated the average seeing during each XShootU observation. As the DIMM seeing (σ_DIMM) is defined at 500 nm at zenith, this value needs to be corrected for the wavelength dependence and the airmass of each observation (Fried 1966; Goodman 2000): $${\sigma }_{\lambda }={\sigma }_{\text{DIMM}}{\left(\frac{\lambda }{500\text{\hspace{0.17em}nm}}\right)}^{-0.2}\times {\text{airmass}}^{0.6}.$$(1)

The slit transmission fraction (T) at each wavelength is given by: $$T=\text{erf}(\sqrt{\text{ln\hspace{0.17em}}2}\frac{w_{\text{slit}}}{{\sigma }_{\lambda }}),$$(2)

where w_slit is the slit width and erf(x) is the error function: $$\text{erf}(x)=\frac{2}{\sqrt{\pi }}{\displaystyle {\int }_0^x{\text{e}}^{-t^2}}\text{d}t.$$(3)

We then corrected the (relative) response-corrected spectrum for the slit losses (1 – T) at each wavelength bin. An example is given in Fig. 12.

Fig. 11 Examples from the reduced UVB spectrum of 2DFS-3694 illustrating the absence of order re-connection artifacts. The blue and orange lines correspond to the wavelength regions of the first and second orders (as presented at each panel). The vertical dashed lines indicate the locations of the order re-connections (highlighted as red lines, which is the overlapping of the blue and orange regions). The order numbers and the common wavelength range of the orders are indicated at the top of each panel.

4.2 Telluric correction

The slit-loss-corrected spectra of the VIS arm were corrected for telluric absorption using the MOLECFIT tool v3.0.3 (Smette et al. 2015; Kausch et al. 2015). The telluric model includes O₂ and H₂O absorption lines, appropriate for the VIS wavelength region. The standard wavelength regions used by MOLECFIT were adjusted to provide the best overall results for the full sample of stars. In particular, all wavelength regions with strong telluric features are included, while regions with stellar spectral lines are excluded. The adopted wavelength ranges are given in Table 4.

Generally, the spectrum can be well fitted by the telluric model, with typical fit residuals on the order of 1–2% (Fig. 13). This results in a good telluric correction in most of the wavelength regions with weak telluric lines, with residuals well below 1% (see Fig. 14). Notable exceptions are the strongest molecular bands of O₂ at 760–770 nm and H₂O at 930–980 nm, both of which are close to or fully saturated (depending on weather conditions). Residuals in these regions typically reach up to 10% (see Fig. 14, insets), but can be much higher in rare cases with particularly poor weather conditions.

Fig. 12 Comparison of the spectrum of AV-43 before slit-loss correction (lower spectrum), after slit-loss and telluric absorption correction (middle spectrum), and after the scaling to the B magnitude and lining up the VIS to the UVB (upper spectrum). The corresponding V and R magnitudes are indicated with star symbols. The transmission curves are from Bessell (1990).

4.3 Absolute flux calibration

After the response-curve and slitloss correction steps (Sects. 3.4 and 4.1), one might expect the relative scaling between the UVB and VIS arms to be correct, and the integration of absolute fluxes (e.g., over the B or R-band) to match the corresponding observed photometric magnitudes. However, despite considerable improvements, the agreement is generally not perfect and differences of about 5%, on average, between the two arms and with respect to published photometry remain. These likely arise from residuals from arm-dependent flux losses, sky-subtraction, and variations in the sky transparency, seeing, and airmass during the integration time. It might also result from differences between the DIMM seeing and the true seeing. We thus implemented two additional steps to further improve on the absolute flux calibration of the XShootU spectra.

In a first step, we scaled the slit-loss-corrected spectra to published photometric measurements of our targets. To do this, we collected B and R magnitudes for all our targets. All used B magnitudes and their references are listed in Table 1 of Paper I, with the exception of the magnitudes for AV-232, HV5622, LH58-496, NGC346-ELS-25, NGC346-ELS-51, VFTS-169, and VFTS-586, which come from Massey (2002). All R-band photometry comes from Massey (2002), with the exception of AV264, which comes from Zacharias et al. (2013). The scaling factor was calculated using the B filter for the UVB arm and the R filter for the VIS arm as these filter’s bandpasses are almost entirely contained within the wavelength range of the two arms (see Fig. 12). Transmission curves were adopted from Bessell (1990). The U- and V-band values were used to check the consistency of the applied scaling. Figure 15, top panel, shows the distribution of scaling factors for the UVB arm, per spectral type.

After this step, small shifts between the UVB and VIS arms remained for many observations (3% on average), for example due to uncertainties in the published photometry or target variability. To create a smooth SED across the two arms, we applied a second scaling. We calculated the average flux in each arm from 545 nm to 555 nm (the overlapping wavelength region), and applied a global scaling factor to the VIS arm to align it with the flux of the UVB arm. For most observations, this correction was less than 10%. For targets with no published photometry, we only applied this second step. For some targets, no R-band magnitudes are available. In this case, we scaled only the UVB spectra to published photometry and the second scaling step was used to align the UVB and VIS arms. Figure 15, bottom panel, shows the distribution of the additional scaling factors for the VIS arm, per spectral type. This figure only includes the stars for which R-band photometry was available. For both panels, no trend with spectral type is visible, and variations must be connected to seeing variations during observations, variable cloud coverage or the uncertainties in the photometry values. A comparison of spectra before and after slit-loss correction, and after the final scaling to the B magnitude, is shown in Fig. 12.

4.4 Rectification to continuum

The response-curve and slit-loss-corrected individual exposures were individually normalised by adjusting the continuum with a functional form. The rectification procedure follows the general algorithm described in Sana et al. (2013) but uses different functional forms. In short, the algorithm works in a pseudo wavelength range (centered on 0 and ranging from −1 to +1 for numerical stability) and proceeds in three steps: (i) an initial detrending of the flux-calibrated spectrum using a linear fit; (ii) the identification and removal of absorption and emission line regions using the ratios of a minimum or a maximum filter with to a median filter, and (iii) an iterative κ − σ clipping to adjust a pre-defined function to the remaining continuum points. Three iterations were usually adopted, but this could be modified upon a visual inspection of the results.

To make the automatic rectification process more robust, we selected regions where no continuum points were allowed (e.g., regions around He II λ4686 or Hα) but also anchor continuum regions where the σ-clipping could not remove continuum points (see Fig. 16). These anchor regions help the rectification process to more robustly converge in case the first step in the iterative cleaning of spectral lines is too aggressive.

The full spectral range was divided in four regions, two in each arm, that were separately rectified (see Fig. 16). For the UVB arm, the region blueward of the Balmer jump (310365 nm) was rectified using a quadratic function while the region redward of the Balmer jump (390–555 nm) uses a generalized Planck-function: $$f_{\text{c}}(\lambda )=a{\lambda }^b{\left({\text{e}}^{c/\lambda }-1\right)}^{-1},$$(4)

with a, b and c being variables of the fit. Similarly, the VIS arm wavelength range was split into two sections. The range blueward of the Paschen jump (555–820 nm) was adjusted with Eq. (4) while the range redward of the Paschen jump (800–1020 nm) was adjusted with a quadratic function.

The choice of the functions was motivated by the respective sizes of the wavelength ranges to fit and by our intention to avoid introducing high-order frequencies in the rectified signal. Such high-frequency fluctuations could impact shallow absorption or emission profiles that are of astrophysical significance (e.g., extended Balmer wing absorption, diffuse interstellar bands or extended wind line emission) that would be overcorrected using a local approach or a more aggressive normalization function (e.g., high-order polynomials, splines). In doing so, we provided a global high-quality normalization (within 0.5% of the continuum in most of the spectral range). Further local normalization may need to be applied, depending on the science case.

The regions around the Balmer and Paschen jumps were excluded from the rectification procedure because they cannot be adjusted without an atmospheric model. As described above, the overall result is quite satisfying given the large spectral ranges covered (up to 265 nm) and the limited number of free parameters (only three) in each region. Some limitations are, however, apparent:

• The bluest region (< 350 nm) displays only a pseudo continuum due to the high density of spectral lines. This prevented us from retaining continuum points so that the adopted continuum is an interpolation of the quadratic function used in the 324.64–364.0 nm region.

• There is a tendency to underestimate the continuum below 410 nm, while overestimating the continuum just redward of Hδ (by ~0.5%); see Fig. 16, top-right panel.

• Small (~0.5%) dips in continuum blueward of the Balmer lines are often present.

• Broad bumps (~ 1%) in flux around 570 nm are caused by the O₂ model of the telluric absorption correction (see Sect. 4.2).

• There is a broad bump in flux from 900–960 nm. This is a result of the telluric absorption correction as the telluric lines are (close to) saturated in this region.

Fig. 13 Example of a telluric model fit to the strong O2 (left, top) and H2O (right, top) bands for one of the exposures of AV-43. The residuals of the fit (observation – model in 10−16 erg s−1 cm−2 Å−1) are shown in the sub-panels.

Fig. 14 Comparison of the spectrum of AV-43 before (SLC) and after (TAC) using the model from Fig. 13 to correct for telluric absorption. The insets show the two regions with the strongest absorption features; they are in the same units as the main panel.

Fig. 15 Scaling factors for absolute flux correction, as a function of spectral type. Top panel: scaling factors applied to the UVB arm to make the spectra match the B magnitude photometry. Bottom panel: extra scaling factors applied to the VIS arm, after matching them to the R photometry, to align the UVB and VIS arms. The dotted lines in both panels correspond to ±3σ ranges.

4.5 Co-addition and final products

The 1D flux-calibrated spectra corresponding to individual exposures of each arm were corrected for barycentric motion. They were also re-sampled to a common wavelength grid starting at 300 nm and 545 nm for the UVB and VIS arm, respectively, adopting a constant wavelength step of 0.02 nm (see Sect. 3). The re-binned spectra of all exposures were co-added using a weighted mean where the weight is given by the inverse of the square error spectrum provided by the pipeline. Regions of 1D spectra that were identified as saturated (Sect. 2.2) were masked in the co-addition process.

We did not co-add the normalised spectra of each epoch but re-did the normalization after co-addition of the flux-calibrated UVB and VIS spectra using the same procedures as described in Sect. 4.4. This prevents that any uncertainty on the rectification of the individual epochs from being propagated to the co-added data.

Finally, the co-added flux-calibrated spectra of the UVB and VIS arms were merged by stitching them together at λ = 550 nm. We chose to adopt a strict transition from one arm to another rather than, for example, a progressive transition over a small wavelength range. This choice was driven by the fact that each arm has a different spectral resolving power, such that it is clear for the user at which wavelength the change of resolution occurs. Data above 550 nm for the UVB arm and below 550 nm for the VIS arm are still made available in the data products of the individual epochs (see Sect. 6.3), but are ignored in the co-added data products (Sect. 6.4).

4.6 Overall quality: S/N distribution

We estimated the S/N obtained for each observation using continuum windows in the reduced and normalised spectra. For the UVB arm, the continuum region was taken from 475 to 480 nm while it ranged from 675 to 685 nm in the VIS arm. An exception occurs for broad emission line stars for which no automatic normalization was applied. In these cases, the VIS continuum region was restricted to the range of 683 to 685 nm and a local normalization was applied before computing the S/N. In each case, the S/N was computed per pixel as the ratio of the median to the 1σ-dispersion. The resulting distributions are displayed in Fig. 17 and show a median (per pixel) S/N above 200 in the UVB arm and around 150 in the VIS arm of the co-added spectra. These correspond to a median S/N per resolution elements of above 350 and 250 in the UVB and VIS co-added spectra, respectively.

5 RV measurements

The radial velocities were determined by cross-correlating the observed spectra from the DR1 release with theoretical template spectra calculated with the PoWR model atmosphere code (Hamann & Gräfener 2003; Sander et al. 2015). We used template spectra from the publicly available model grids³ (Hamann & Gräfener 2004; Todt et al. 2015; Hainich et al. 2019). In total, 3610 OB and WN model spectra were used. We did not employ the WC models, as their strong emission lines make radial velocity determinations difficult. For the SMC and LMC targets, we used the individual 1D rebinned spectra from the DR1. For each spectrum we take the column WAVELENGTH_AIR for the wavelength and the column SCI_FLUX_NORM for the normalized flux. The wavelength region over which we cross-correlate is limited to 400–890 nm. The blue limit avoids the Balmer jump, while the red limit avoids the increasing level of noise.

We exhaustively compared the observed spectra to all theoretical spectra, to determine the radial velocity for each observing epoch. A grid of ten projected rotational velocities was also tried (υ sin i = 0, 20, 50, 100, 150, 200, 250, 300, 400, 500 km s⁻¹) and all PoWR spectra were convolved with the rotational profile (Gray 2005). They were also convolved with a Gaussian instrumental profile (R = 6700 for the UVB spectrum and R = 11 400 for the VIS spectrum).

For each observed target epoch, a loop is run over all 3610 × 10 rotationally broadened PoWR spectra. Each PoWR spectrum is cross-correlated with the observed spectrum. We masked out a number of wavelength regions: the region around He II λ4686, the interstellar Na I lines, the regions corrected for telluric lines, the gap around 636 nm (which has been interpolated over in the normalized spectra), and the Paschen P16 line (which has an incorrect wavelength in the PoWR spectra for grids computed before 2021; see bug warning on the PoWR website). We further masked the pixels that have a nonzero value in the SCI_FLUX_QUAL column. When there are repeated observations of one target, all spectra of that star are fit with the same theoretical spectrum and the same rotational broadening, but allowing for different radial velocities for each epoch. Our cross-correlation technique uses a velocity grid from -200 to +600 km s⁻¹ with a step-size of 10 km s⁻¹. The final value for the radial velocity is then found by a parabolic interpolation around the maximum of the cross-correlation function.

The cross-correlation function is computed as (e.g., David et al. 2014) $$C({\upsilon }_{\text{s}})=\frac{{\displaystyle {\sum }_i}(F_i-\overline{F})(T_i-\overline{T})}{\sqrt{{\displaystyle {\sum }_i}{\left(F_i-\overline{F}\right)}^2}\sqrt{{\displaystyle {\sum }_i}{\left(T_i({\upsilon }_{\text{s}})-\overline{T({\upsilon }_{\text{s}})}\right)}^2}},$$(5)

where F_i are the observed fluxes, with the average observed flux $\overline{F}$, and T_i(υ_s) the theoretical fluxes, shifted by a velocity v_s, with the average theoretical flux $\overline{T({\upsilon }_s)}$.

A parabolic fit is applied to find the maximum of the cross-correlation function, and thus the corresponding radial velocity and radial velocity uncertainty. For the uncertainty σ_RV, we use: $${\sigma }_{\text{RV}}^2=-{\left[N\frac{C^{″}C}{1-C^2}\right]}^{-1},$$(6)

where C is the cross-correlation function at maximum, C″ its second derivative, and N the number of points in the spectrum (Zucker 2003).

Among all the PoWR spectra we chose the one that had the highest maximum in the cross-correlation function and used that one to determine the radial velocity. When there were repeated observations, we chose the PoWR spectrum that had the highest sum of the cross-correlation maxima. The resulting radial velocities are given in Table C.1. Table C.2 lists the stars for which we could not determine the radial velocity, including the reason.

Although we list the chosen PoWR spectrum for each star in Table C.1, we caution that its parameters are not necessarily the ones that best represent the properties of the star. The radial-velocity determination procedure focuses on other features of the spectrum than a stellar parameter determination procedure does. As one example, we recall that the He II λ4686 line was excluded from the procedure. In addition, most of the chosen PoWR spectra we used turned out to have SMC metallicity, even for the LMC stars.

The formal uncertainty derived from Eq. (6) is small for most of the stars. As an alternative estimate of the uncertainty, we measured the short-term radial velocity behavior. Figure 18 shows the difference in radial velocity between observations of the same star obtained up to 1 h apart. The differences are ~2 km s⁻¹, which is a more representative value of the radial velocity uncertainty.

As a further indication of the reliability of our radial velocity determination, we also determined a quality index. We introduced different quality criteria that give a weight to how well the theoretical spectrum fits certain observed spectral lines. For the VIS spectra, this includes judging the depth of the Paschen lines, which play an important role in the radial velocity determination. We judged by eye whether the observed lines are broad or narrow, and how well the theoretical spectrum fits the observed one. For each of these criteria we assigned points (positive or negative) and determined the final quality indicator by adding up the points. The details are given in Table 5.

For approximately half of our stars, Simbad also lists a radial velocity. Figure 19 compares our values to those of Simbad. The results are color-coded according to the quality of the RV determination. The comparison shows how the agreement between Simbad and our results improves as we go to higher quality measurements as defined in Table 5.

Fig. 16 Example of automatic normalization of 2DFS-163 for each range in the UVB and VIS arms. Regions shaded in red and green are regions that are rejected or forcefully included as continuum regions.

Fig. 17 Distribution of S/N for the UVB (upper panel) and VIS (lower panel) part of the spectra.

Fig. 18 Radial velocity difference between observations of the same star (listed on the x-axis) taken 1 h or less apart. The title lists the median and the dispersion. The dashed lines indicate the ±2 km s−1 range.

Fig. 19 Comparison between our radial velocities and those from Simbad. The color scale indicates the quality of our radial velocity determination. The formal errors on our results are, in most cases, smaller than the symbol.

6 Data release

The XShootU pipeline reduced data, and the 2D and 1D data products are released to the community using Zenodo⁴. A part of these products are also ingested into the ESO phase 3 archive⁵ and on the ULLYSES STScI website⁶. The UVB/VIS data release associated with this paper is dubbed DR1. It replaces nonpublic early releases (eDR) that have been distributed in the consortium for test and development purposes. We do not envision significant changes of the UVB/VlS data reduction beyond DR1, but the numbering will allow us to keep track of changes should they be required. DR2 should include the NIR but is not anticipated before 2025 because the ESO pipeline is currently undergoing significant improvement of the wavelength calibration using telluric lines. In addition, the HST ULLYSES target list has recently been expanded using archival data. The XShootU collaboration is currently seeking to obtain the corresponding X-shooter observations. Should they be acquired, they will be included in a subsequent data release.

The data format and the various data products that we provide to the community are extensively described in dedicated documentation available on Zenodo. We provide here a brief overview of the available content.

6.1 Pipeline reduction products

Pipeline reduction products (PRPs) are produced by running the esorex pipeline cascade (Sect. 3.1). We collected all intermediate PRPs and master calibrations in addition to the reduced 1D and 2D spectra. For example, we obtained 20 PRP files from the reduction of each scientific frame. These PRPs are described extensively in the XSHOOTER pipeline manual that can be retrieved (e.g., from the instrument webpage⁷) and is not be repeated here.

6.2 Consolidated 2D data products

We provide a subset of the most important 2D PRPs as consolidated 2D reduction products (CRP_2D). These are formed by one set of three multi-layer fits files for each star:

• XSHOOTU_<target>_FLUX_MERGE2D: flux-calibrated 2D science spectra;

• XSHOOTU_<target>_FADU_MERGE2D: non flux-calibrated 2D science spectra;

• XSHOOTU_<target>_SKYADU_MERGE2D: non flux-calibrated 2D sky spectra.

These files provide, for each target, the reduced (flat-fielded, wavelength-calibrated, straightened, resampled and merged orders) 2D spectra corresponding to each science exposure. These are produced, respectively, from the SCI_SLIT_FLUX_ MERGE2D_<ARM>, SCI_SLIT _MERGE2D_<ARM> and SKY_SLIT _MERGE2D_<ARM> PRPs. Aside from the main extension (empty, though fits descriptors allow indentification of the target and DR version), each file thus contains N_epoch × N_arm × 3 layers, where N_epoch is the number of individual observing epochs and N_arm is the number of arms in the DR. The factor 3 comes from the fact that we provide FLUX, ERRS and QUAL layers for each exposure. Each extension is tagged by the ARM, observing time and nature (FLUX, ERRS, QUAL) of the data. ESO headers are preserved in the headers of the individual extensions. These CRP_2D files should be the entry points to investigate the spatial surroundings of the XShootU targets (e.g., nebulosity, nearby objects in the slit), or to perform a custom-made sky-subtraction or extraction of the science spectrum.

6.3 Consolidated 1D data products

All extracted 1D science and sky spectra (epochs, arms), and the advanced data products of each exposure are consolidated into a set of two fits files for each target:

• XSH-U_<target>_INDIV1D: non-barycentric-corrected 1D science and sky spectra;

• XSH-U_<target>_INDIV1D_REBIN: barycentric-corrected and resampled 1D science and sky spectra;

The INDIV1D and INDIV1D_REBIN files have the same structure and contain similar information, but the wavelength calibration scale (WLC) of the INDIV1D files is imported directly from the PRPs while the WLC of the INDIV1D_REBIN files has been corrected for the barycentric motion and the 1D data has been resampled onto a common WLC grid using a quadratic interpolation. This common grid starts at 300 and 545 nm for the UVB and VIS arm, respectively, and has a step of 0.02 nm, thus providing about 2 to 4 pixels per resolution elements. The WLC is provided in air (λ_air) and in vacuum (λ_vac) so that the data can be easily combined with the HST ULLYSES UV data that are in vacuum. Vacuum wavelengths were computed using the PyAstronomy⁸ (Czesla et al. 2019) routine airtovac2, which iteratively solves the equation: $${\lambda }_{\text{vac}}=n{\lambda }_{\text{air}},$$(7)

where the refractive index n is taken from Ciddor (1996): $$n=\frac{0.05792105}{238.0185-{\lambda }_{\text{vac}}^{-2}}+\frac{0.00167917}{57.362-{\lambda }_{\text{vac}}^{-2}}+1.$$(8)

Here λ_vac is the vacuum wavelength in µm, and assuming standard air⁹. The flux-calibrated spectra were not resampled when moving from air to vacuum, leading to a small systematic error ≲10⁻³ relative to the flux level. This error is much smaller that the precision of our flux calibration.

As for the CRP_2D, separate extensions provides access to the spectra obtained at different epochs and for the different arms. Each extension is tagged by the ARM and observing time. ESO headers of the SCI_SLIT_FLUX_MERGE1D_<ARM> PRPs are preserved in the headers of the individual extensions. Each extension provides a fits table containing, among others, the response-curve-corrected science spectra, the science and sky spectra in ADU (non-response-corrected), as well as spectra corrected for slit-losses, telluric absorption and absolute flux-calibration. Correction functions (e.g., response-curve, wavelength-dependant slit-correction and normalization functions) are not included to limit the size of the data but can be reconstructed from the different columns. These files containing the individual epochs are thus the entry point to study variability, and to access the extracted but unprocessed 1D spectra of the science target or the sky.

6.4 Co-added one-dimensional (1D) data products

Our final data products result from the co-addition of all collected 1D spectra, and the connection of the UVB and VIS arms into a single spectrum. The result is provided in a single fits file per target:

• XSH-U_<target>_COADD1D: barycentric-corrected, resam-pled and co-added 1D science and sky spectra.

As for the consolidated 1D data products, the COADD1D file provide the co-added and response-curve-corrected science spectra, the co-added science and sky spectra in ADU (non-response-corrected), as well as the co-added spectra corrected for slit-losses, telluric absorptions and absolute flux-calibration. In case of multi-epoch observations, the co-addition ignores variability. The co-added 1D products should be suitable for atmosphere analysis or spectral library applications depending on whether the target is suitable for the science goals of the analysis (e.g., in terms of variability vs. average spectrum). In particular, we do not apply any correction for RV variations of the source, such as those expected from spectroscopic binaries.

6.5 Quality control figures and atlas

Our data release is accompanied by a number of figures that allow the community to quickly assess the quality of the spectra, of the telluric correction and of the normalization. Examples are given in Figs. 12, 14 and 16 and are available on Zenodo. An atlas of the flux-calibrated and normalized co-added spectra is also provided.

7 Summary

We reduced the UVB and VIS X-shooter data obtained by the XShootU large program in a consistent and systematic way, improving upon the standard reduction in a number of ways that include the flat-fielding strategy and the use of improved FLUX standard models. We provide advanced reduction products, including the correction for slit losses and telluric lines, and an absolute flux calibration. We also provide normalized spectra using smooth functions that should not have introduced high-order frequencies in the resulting spectra as the latter might modify the quantitative science content of the data. Various data products are provided to the community, including 2D and 1D consolidated reduction products as well as barycentric-motion-corrected and co-added spectra for which all individual epochs and the data from the UVB and VIS are combined to optimize the S/N. We also provide quick-look figures that enable a first assessment of the quality of various steps in the data processing and of the final data products. The DR1 data should be suitable for immediate scientific exploitation.

Acknowledgments

We thank the referee for their constructive comments. This project has received funding from the KU Leuven Research Council (grant C16/17/007: MAESTRO). M.A. acknowledges support from “la Caixa” Foundation (ID 100010434) under the fellowship code LCF/BQ/PI23/11970035. L.M. thanks FAPESP (grant 2022/03703-1) and CNPQ (grant 307115/2021-6) for partial funding of this research. AACS and VR are supported by the Deutsche Forschungsgemeinschaft (DFG - German Research Foundation) in the form of an Emmy Noether Research Group - Project-ID 445674056 (SA4064/1-1, PI Sander). AACS and VR are further supported by funding from the Federal Ministry of Education and Research (BMBF) and the Baden-Württemberg Ministry of Science as part of the Excellence Strategy of the German Federal and State Governments. JSV gratefully acknowledges support from STFC via grant ST/V000233/1. CJKL gratefully acknowledges support from the International Max Planck Research School for Astronomy and Cosmic Physics at the University of Heidelberg in the form of an IMPRS PhD fellowship. GM acknowledges funding support from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (Grant agreement No. 772086). F.N. acknowledges funding by grants PID2019-105552RB-C41 and MDM- 2017-0737-19-3 Unidad de Excelencia “María de Maeztu”. C.K. acknowledges financial support from the Spanish Ministerio de Economía y Competitividad under grants AYA2016-79724-C4-4-P and PID2019-107408GB-C44, from Junta de Andalucía Excellence Project P18-FR-2664, and from the grant CEX2021-001131-S funded by MCIN/AEI/ 10.13039/501100011033. PAC and JB are supported by the Science and Technology Facilities Council research grant ST/V000853/1 (PI. V. Dhillon). R.K. acknowledges financial support via the Heisenberg Research Grant funded by the German Research Foundation (DFG) under grant no. KU 2849/9.

Appendix A Target names

Appendix B Individual notes

AV-324

Observations in poor seeing conditions (up to 2″). We adopted a broader extraction window for the central object and smaller windows for sky regions.

BAT99 105

The faint companion on the lower part of the slit has no impact on the sky estimate, probably because of the median approach used in the pipeline to construct the sky spectrum.

HD269927C (= BAT99 120)

We note a very faint companion in the lower part of the slit, contributing less that 1% of the flux of the central object. We also note variable nebulosity across the slit. The standard extraction windows, where the companion is included in the sky regions, result in a better nebular subtraction than when only using the upper sky region (‘Default Sky 2’ in Fig. 8). We opted for the former approach even though the companion probably contributes to ~ 0.5% of the extracted flux. The shapes of the Balmer lines (strong wind lines) in the extracted spectrum of HD269927C do not seem to be different in the two approaches.

SK-67D106

Observations in poor seeing conditions (up to 2″). We adopted a broader extraction window for the central object and smaller windows for the sky regions.

SK-191

Observations in poor seeing conditions (up to 2″). We adopted a broader extraction window for the central object and smaller windows for the sky regions.

VFTS-482

The slit contains at least three other objects, at approximate positions of −2.″4, +2.″2 and +5.″0. Only a ≈1.″2 window toward the bottom of the slit is available for the estimate of the sky background. The brightest of the three companions is too close to the edge of the slit for a meaningful extraction while the other two are too faint.

VFTS-506

Both edges of the slit seem contaminated by nearby sources, but this is not affecting the 2D sky-subtracted spectra thanks to the limited contamination and the median approach used. The sky level seems higherthan forotherstars in the 30 Dorregion, possibly due to nearby faint companions that are not clearly resolved. The level of contamination, if present, is likely faint or even negligible so that we used the standard extraction window forVFTS 506.

VFTS-542

A nearby visual companion is identified at +4.5″ in the upper part of the slit. Another much fainter companion is also present at + 2.2″ so that we do not use the upper part of the slit to estimate the background and sky levels. Finally, there is a possible very faint companion at about −3″ (i.e., within the adopted sky region). Its trace is only visible in the reddest part of the VIS arm on the raw image and its contribution is probably completely negligible, and so we have ignored it. The companion at +4.5 presents clear signatures of Balmer lines, and of He I λλ4443, 4471 and He II λ4686 absorption.

Appendix C Radial velocity results

References

All Tables

All Figures

	Fig. 1 Example of the ghost artifact identified in various 2D RAW images (left panel). Whenever the ghost position overlaps with the extraction region of the science target, it produces an artifact in the extracted ID spectrum (right panel, at λ ≈ 910 nm in this example). The insets are in the same units as the main panel.
In the text

	Fig. 2 Comparison of a reduced spectrum of SK-68D137 without (left panels) and with (right panels) response-curve correction.
In the text

	Fig. 3 Comparison of reduced spectra of SK-68D137 obtained with the different flat-fielding strategies when the SCIENCE target and the FLUX standard star have been observed on different nights. Red indicates when different FF frames were used to reduce the observations of the SCIENCE and FLUX stars. Black indicates when the same FF frames were used.
In the text

	Fig. 4 Radial velocity difference measured between the UVB and VIS arms after correction for known wavelength calibration offsets (Sect. 3.3). The panels show the values for different data quality (increasing from top to bottom). The UVB+VIS quality of a star is given by the sum of all points listed in Table 5. The median and 68.3% quan-tile are listed at the top of each plot.
In the text

Fig. 5 Comparison of knot points between the original knots (in gray) and the new knots (in red) used for interpolation in determining the response curve for GD 71 in the UVB arm. The top panel shows the observed spectra as a function of wavelength for one observation, while the bottom panel shows the reference model divided by the observation in black and the resulting response curve after interpolation in orange.

In the text

	Fig. 6 Comparison of new and old reference models for the EG 274 and EG 21 FLUX standard stars used by the XShootU program. The top panel of each star shows the observed co-added spectrum (in black) and the new and old reference models (in red and blue, respectively). The sub-panels show the residuals between the observation and each of the reference models using the same color scheme.
In the text

Fig. 7 Comparison of the Hγ (top) and Hβ (bottom) line profile of the target AV-238 for three different observations (Dec. 2020, May 2021, and June 2021) and calibrated with three different FLUX standard stars (black: EG 21, blue: LTT 9789, red: GD 153). The residual plots show the ratio of data calibrated with LTT 7987 (blue) and GD 153 (red) and those calibrated with EG 21. The 1 and 2% inconsistencies in the blue wing of Hβ and Hγ, respectively, disappeared when using the new flux models. The remaining differences located in the core of the lines are likely of astrophysical origin (e.g., line profile or RV variations, or nebular correction residuals).

In the text

Fig. 8 2D reduced spectrum of AV-251 (left) and AV-6 (right), revealing a faint companion 1″ below the spectrum of AV6. The overlaid grid indicates the wavelength axis (horizontal, in nm) and the spatial axis (vertical, in arcsec from the center of the slit). The default science and sky extraction windows and those adopted for AV-6 in this work are overlaid on the 2D spectrum. The color bar indicates the ADU level in the images (in log scale). The cross sections along the spatial direction are in units of ADU (in linear scale).

In the text

	Fig. 9 Comparison of the extracted spectra of AV-6 with the default (black) and manually selected extraction and sky regions for the components centered at +1.05″ (blue) and −0.8″ (red). The left and right panels correspond to the UVB and VIS arms, respectively. The insets focus on various spectral lines of interests and are in the same units as the main panels.
In the text

	Fig. 10 Comparison of reduced spectra of AV-251 obtained with the different cosmic rejection thresholds. Black: –valueremovecrhsingle– sigmalim=20 (default). Red: –valueremovecrhsingle-sigmalim= 10 (adopted). The spectra are slightly shifted for clarity. The insets focus on spectral lines of interest and are in the same units as the main panel.
In the text

Fig. 11 Examples from the reduced UVB spectrum of 2DFS-3694 illustrating the absence of order re-connection artifacts. The blue and orange lines correspond to the wavelength regions of the first and second orders (as presented at each panel). The vertical dashed lines indicate the locations of the order re-connections (highlighted as red lines, which is the overlapping of the blue and orange regions). The order numbers and the common wavelength range of the orders are indicated at the top of each panel.

In the text

	Fig. 12 Comparison of the spectrum of AV-43 before slit-loss correction (lower spectrum), after slit-loss and telluric absorption correction (middle spectrum), and after the scaling to the B magnitude and lining up the VIS to the UVB (upper spectrum). The corresponding V and R magnitudes are indicated with star symbols. The transmission curves are from Bessell (1990).
In the text

	Fig. 13 Example of a telluric model fit to the strong O2 (left, top) and H2O (right, top) bands for one of the exposures of AV-43. The residuals of the fit (observation – model in 10−16 erg s−1 cm−2 Å−1) are shown in the sub-panels.
In the text

	Fig. 14 Comparison of the spectrum of AV-43 before (SLC) and after (TAC) using the model from Fig. 13 to correct for telluric absorption. The insets show the two regions with the strongest absorption features; they are in the same units as the main panel.
In the text

	Fig. 15 Scaling factors for absolute flux correction, as a function of spectral type. Top panel: scaling factors applied to the UVB arm to make the spectra match the B magnitude photometry. Bottom panel: extra scaling factors applied to the VIS arm, after matching them to the R photometry, to align the UVB and VIS arms. The dotted lines in both panels correspond to ±3σ ranges.
In the text

	Fig. 16 Example of automatic normalization of 2DFS-163 for each range in the UVB and VIS arms. Regions shaded in red and green are regions that are rejected or forcefully included as continuum regions.
In the text

	Fig. 17 Distribution of S/N for the UVB (upper panel) and VIS (lower panel) part of the spectra.
In the text

	Fig. 18 Radial velocity difference between observations of the same star (listed on the x-axis) taken 1 h or less apart. The title lists the median and the dispersion. The dashed lines indicate the ±2 km s−1 range.
In the text

	Fig. 19 Comparison between our radial velocities and those from Simbad. The color scale indicates the quality of our radial velocity determination. The formal errors on our results are, in most cases, smaller than the symbol.
In the text