Analysis of high-resolution FEROS spectroscopy for a sample of variable B-type stars assembled from TESS photometry

Spectroscopic data are necessary to break degeneracies in asteroseismic modelling of the interior structure in high- and intermediate-mass stars. We derive precise photospheric stellar parameters for a sample of 166 B-type stars with TESS light curves, suitable for detailed asteroseismic studies, through a homogeneous spectroscopic analysis. The variability types of these stars are classified based on all currently available TESS sectors. We obtained high-resolution spectra for all 166 targets with the FEROS spectrograph in the context of a large program. The Least-Squares Deconvolution method is employed to investigate spectral line profile variability and to detect binary systems. We identify 26 spectroscopic double-lined binaries; the remainder of the sample are 42 supergiants in the LMC galaxy and 98 Galactic stars. The spectra of the Galactic stars are analysed with the zeta-Payne, a machine learning-based spectrum analysis algorithm. We determine the five main surface parameters with average formal precisions of 70 K (Teff), 0.03 dex (logg), 0.07 dex ([M/H]), 8 km/s (vsini), and 0.7 km/s (vmicro). The average internal uncertainties we find for FEROS spectra with our spectrum analysis method are 430 K (Teff), 0.12 dex (logg), 0.13 dex ([M/H]), 12 km/s (vsini), and 2 km/s (vmicro). We find spectroscopic evidence that eight of the 98 Galactic targets are fast rotating g-mode pulsators occurring in between the slowly pulsating B (SPB) stars and delta Scuti instability strips. The g-mode frequencies of these pulsators are shifted to relatively high frequency values due to their rotation. Their apparently too low Teff relative to the SPB instability region can in most cases be explained by the gravity darkening effect. We also discover 13 new HgMn stars of which only one is found in a spectroscopic binary, resulting in a biased and therefore unreliable low binary rate of only 8%.


Introduction
While massive stars are important suppliers of chemical and mechanical feedback to the interstellar medium (e.g.Langer 2012), their structure and evolution is not yet well understood.Large theoretical uncertainties occur already during the early phases and have an impact on the further evolution (e.g.Maeder & Meynet 2000;Ekström et al. 2012;Aerts et al. 2019).To calibrate physical processes in these stars, observational constraints from their interiors are needed.With space asteroseismology, A&A 665, A36 (2022) profile of the convective boundary mixing region and the radiative envelope, the internal rotation rate, and the amount of chemical mixing both in slowly pulsating B (SPB) stars and β Cephei (β Cep) stars (e.g.Miglio et al. 2008;Szewczuk & Daszyńska-Daszkiewicz 2015;Moravveji et al. 2016;Pedersen et al. 2018;Michielsen et al. 2019;Moździerski et al. 2019) and in less massive γ Doradus (γ Dor) stars (e.g.Van Reeth et al. 2018;Mombarg et al. 2020;Angelou et al. 2020).These are three classes of pulsating main-sequence stars with predominant g-mode pulsations in γ Dor (1.3-1.9M ⊙ ) and SPB (3-8 M ⊙ ) stars and p-mode pulsations in β Cep stars (8-25 M ⊙ ).
There are not that many asteroseismic studies of massive stars as compared to low-mass stars, because long time-base, continuous and high-precision time series are needed to properly resolve the frequencies and identify the geometries of g-mode pulsations (Bowman 2020).This became possible with the start of science operations of space missions such as CoRoT (Auvergne et al. 2009) and Kepler (Borucki et al. 2010).Especially the Kepler mission revolutionised the field of g-mode asteroseismology with its four-year-duration light curves.However, Kepler was pointed to one fixed field of view for its entire nominal mission duration, where the field selection was driven by the primary science of the mission of detecting and characterising exoplanets around solar-like stars.Therefore, not many massive stars were observed by the mission.The ongoing Transiting Exoplanet Survey Satellite (TESS) mission (Ricker et al. 2015) is observing the whole sky to search for planets around bright stars, and it thereby also targets O-and B-type stars.The nominal mission divided the two ecliptic hemispheres into 13 sectors that are each observed for 27 days.The 13 sectors in the northern and southern ecliptic hemispheres overlap near the poles, where one-year-long observations are taken.In the southern hemisphere, the TESS continuous viewing zone (CVZ) includes the Large Magellanic Cloud (LMC) galaxy, which means that different metallicity regimes can be probed.Different studies have demonstrated the diverse variability among OBtype stars using Kepler/K2 and TESS light curves.Bowman et al. (2019b) and Burssens et al. (2019) found a number of coherent pulsators among more than 100 OB-type stars in the K2 mission and reported light curves dominated by stochastic low-frequency variability (SLF) for the most massive stars and dozens of blue supergiants observed with K2 and TESS.From the first two sectors of TESS, Pedersen et al. (2019) detected and classified the variability for a sample of 154 OB-type stars, including 40 LMC targets, showing that 90% have pulsational, binary, rotational, or SLF variability.Burssens et al. (2020) also found a large diversity of variability in a sample of OB-type stars observed within the first 13 sectors of TESS, and demonstrated the power of coupling TESS light curves with high-resolution spectroscopy.
Asteroseismic modelling of intermediate-mass stars entails a number of degeneracies due to strong correlations between the numerous parameters occurring in stellar models.The addition of spectroscopic, binarity, or astrometric data helps to lift these degeneracies, reduce uncertainties on the parameters, and supplies information that is required to establish a link between the interiors of stars and their atmospheric layers.With this in mind, many studies in the past decade have focused on obtaining spectroscopic parameters for samples of promising pulsating Kepler stars (e.g.Lehmann et al. 2011;Tkachenko et al. 2012Tkachenko et al. , 2013b,a;,a;Van Reeth et al. 2015;Niemczura et al. 2015Niemczura et al. , 2017)).During early phases of the Kepler mission science operations, ground-based spectroscopy was frequently used in tandem with Kepler photometry to aid the photometric classification of variable stars.Effective temperature (T eff ) and surface gravity (log g) estimates were used to place stars in the Kiel diagram (log g as a function of T eff ) and compare their positions with theoretical instability strips of different types of pulsators, thus confirming or rejecting the photometric classification (e.g.Uytterhoeven et al. 2011;Balona et al. 2011;Tkachenko et al. 2013a).The spectroscopic surface measurements were also used to check estimates of the parameters that were determined from photometry (e.g.Silva Aguirre et al. 2012;Thygesen et al. 2012;Pinsonneault et al. 2014).Conversely, asteroseismic values of log g are often more precise than spectroscopic ones and have thus been used in some studies to better constrain the spectroscopic parameters of stars (Bruntt et al. 2012;Thygesen et al. 2012).
In recent years, detailed asteroseismic modelling of intermediate-mass stars has made advances by including nonphotometric data.For example, Mombarg et al. (2019) added spectroscopic measurements of T eff and log g to the forward modelling of γ Dor stars to lift degeneracies and to estimate masses and ages of stars.Mombarg et al. (2020) studied atomic diffusion in AF-type stars not only by adding spectroscopic T eff and log g values to g-mode asteroseismic modelling, but also by comparing the predicted surface abundances to observational ones to evaluate models with and without atomic diffusion.Pedersen et al. (2018) studied different mixing profiles in more massive B-type stars and emphasised the importance of adding surface abundances to the asteroseismic modelling to constrain the shape of the mixing profiles in the radiative zones of these stars.Pedersen et al. (2021) performed forward asteroseismic modelling on a sample of 26 SPB stars.They applied spectroscopic parameters and an astrometric measurement of the luminosity to delimit the parameter space.In addition to spectroscopy, binary modelling can also add valuable constraints to asteroseismic modelling because it can deliver independent estimates of the radius and mass with very high precision.For example, Johnston et al. (2019) used binary information to reduce parameter uncertainties and correlations between parameters in the modelling of binaries with g-mode pulsators.Sekaran et al. (2020Sekaran et al. ( , 2021) ) combined Kepler photometry with high-resolution optical HERMES spectroscopy (Raskin et al. 2011) to perform a detailed asteroseismic study of a γ Dor-type g-mode pulsator in an eclipsing double-lined spectroscopic binary (SB2) system.The authors demonstrated that it is highly beneficial for asteroseismic modelling to add model-independent information on stellar mass and radius inferred from binary dynamics.
Various classes of chemically peculiar (CP) stars have been encountered among intermediate-mass B-type stars.One class are CP2 or ApBp stars, which show overabundances of Sr, Cr, and Eu.For a subset of ApBp stars, these overabundances, located in spots on the stellar surface, can be explained by the presence of a magnetic field (Smith 1996;Gray 2005).ApBp stars are typically slow rotators, and they are not likely to be in (close) binary systems (Abt & Snowden 1973;Mathys et al. 2020).The other group are CP3 or HgMn stars, which have enhanced abundances of Hg and Mn (and other elements e.g.Ga, Sr, and Y;Smith 1996).HgMn stars are slow rotators with v sin i ∼ 29 km s −1 (Abt et al. 1972) and have a high binarity rate with spectroscopic binarity percentages between 50-91% (e.g.Hubrig & Mathys 1995;Schöller et al. 2010).About half of the binaries are found in SB2 systems, mostly with short orbital periods.Although they are thought to be non-magnetic CP stars, they show rotational modulation consistent with spots of overabundant chemical elements on the surface (e.g.Korhonen et al. 2013;Kochukhov et al. 2021).Unlike for ApBp stars, the origin of these spots in HgMn stars is not yet fully understood.Possible explanations are atomic diffusion processes in the outer A36, page 2 of 17 S. Gebruers et al.: Analysis of FEROS spectroscopy for a sample of variable B-type stars with TESS photometry stellar layers (Alecian et al. 2011) or maybe weak magnetic fields (Hubrig et al. 2020).Another cause for the variability seen in HgMn stars in addition to rotational modulation are pulsations, but not many pulsating HgMn stars have been detected yet (e.g.Alecian et al. 2009;Kochukhov et al. 2021).
In this paper, we analyse the spectra of 166 B-type stars that exhibit photometric variability in their TESS light curves, with the purpose of future combined asteroseismic and spectroscopic modelling.For the sake of efficiency and consistency of the obtained results across the entire stellar sample, we employ the ZETA-PAYNE (Straumit et al. 2022) machine-learning-based spectrum analysis algorithm.The method is a generalisation of the originally proposed THE PAYNE algorithm (Ting et al. 2019) towards the inclusion of a model with an a priori unknown residual instrumental response function into the parameter vector.In this way, the parameters of the residual response function are optimised along with the atmospheric parameters of the star, thus removing the need to normalise stellar spectra manually in advance.
In Sect.2, we discuss the sample selection and observations of the targets.We describe some improvements made to the spectrum reduction pipeline in Sect.3. The spectroscopic classification is presented in Sect. 4 along with a summary of the ZETA-PAYNE framework and a few tests for the specific set-up in this paper.The spectroscopic and photometric results are given in Sect.5, while the internal uncertainties of the method and a subsample of peculiar stars are discussed in Sect.6.This section also presents the obtained results in the context of the spectroscopic Hertzsprung-Russell diagram (Langer & Kudritzki 2014).We conclude our study in Sect.7.

Sample selection and observations
We selected all the bright targets (stars with V-band magnitudes below 11.5 mag) in the samples of variable OB-type stars from Pedersen et al. (2019) and Bowman et al. (2019b).These samples were constructed following the asteroseismic potential of OB-type stars as determined from the first three TESS sectors that were available at that time.For the 166 brightest stars, we obtained spectroscopic observations with the Fibre-fed Extended Range Optical Spectrograph (FEROS; Kaufer et al. 1999), which is attached to the ESO/MPG 2.2 m telescope at La Silla, Chile.In an accompanying paper (Serebriakova et al., in prep.), the spectroscopic analysis of ESO UVES spectra for fainter stars (V > 10 mag) from the samples of Pedersen et al. (2019) and Bowman et al. (2019b) will be presented.
FEROS is a high-resolution spectrograph (R ∼ 48 000) that covers the wavelength range from 3600-9200 Å over 39 orders.The observations were taken during December 2019 and February 2020.In total, 166 objects were observed with a minimum of two epochs per target (separated by 2 months), except for seven targets that were not observable during the second run in February.The spectra were reduced with the publicly available Collection of Elemental Routines for Echelle Spectra (CERES; Brahm et al. 2017) 1 .This package consists of automated reduction pipelines for 13 different échelle spectrographs, among which FEROS.The pipelines are written in Python and complemented with C and Fortran routines to speed-optimise the time-consuming parts.It performs a bias subtraction, correction for bad columns, background subtraction, order extraction, wavelength calibration, barycentric correction, and a blaze correction.To achieve the goals of this paper, some changes to the CERES pipeline were necessary to obtain the best extracted 1D spectra.These changes are described in Sect.3.
To complement the spectroscopic data within this work and to classify the dominant variability type of each target, we reanalysed available two-minute-cadence photometry from TESS.We downloaded the available simple aperture photometry (SAP) and pre-data search conditioning (PDC-SAP) light curves output from the NASA SPOC pipeline (Jenkins et al. 2016) from the Mikulsi Archive for Space Telescopes (MAST2 ).Almost all bright (V ≤ 14 mag) O and early-B stars observed by TESS were prioritised for the two-minute cadence in at least one sector because of successful guest investigator proposals3 .Some stars still fall between the gaps in the CCDs and cameras, however, or are located close to the ecliptic.Cycle 4 of the TESS mission is currently in progress as of July 2021, which includes the second visit to the northern ecliptic hemisphere and parts of the ecliptic for the first time.
For (almost) all of our targets, at least one sector of TESS data is currently available, which is sufficient to classify its dominant variability.Because variability classification is sufficient for this paper as its focus is on the analysis of the newly obtained spectroscopic data for a large number of B-type stars, we focused our photometric analysis on the longer time spans of the two-minute TESS light curves that are available for each target.We checked both the available SAP and PDC-SAP light curves to compare and contrast the impact of instrumental systematics that may be present in either or both light curve formats.In most cases, the latter were sufficient for the variability classification.However, we note that correcting for instrumental effects (e.g.extracting custom light curves) to derive reliable lists of pulsation frequencies is the subject of future work.

Improvements of the CERES pipeline for the FEROS spectra
The CERES pipeline reduces 25 of the 39 FEROS orders, which corresponds to the wavelength range 3800-6800 Å, containing six Balmer lines and sufficient metal lines to determine stellar parameters (Brahm et al. 2017).We made four changes to the existing CERES pipeline in our quest to determine highprecision spectroscopic parameters of our targets.The first two were necessary to remove unphysical artefacts (hereafter referred to as wiggles) from the spectra and to eliminate bumps occurring at the edges of the orders.The other two changes concern an additional function that merges the individual orders and a method to remove cosmic hits.The original CERES pipeline only returns separate orders and not the fully merged spectrum because the main purpose of the developers was to determine precise radial velocities (RV), which is possible from individual orders.
When we obtained reduced spectra with the original CERES pipeline, we detected wiggles and a gap at certain positions in the 1D spectra (see the arrows in the top panel of Fig. 1).We traced these features back to bad columns in the 2D frames.The original pipeline corrects for bad columns by interpolating among the good pixels in each row at the positions of the bad pixels.However, because the orders are skewed, some parts of the bad columns are assigned too much or too little flux in this way, creating the observed wiggles in the 1D spectra.Instead of interpolating within each row, we found that almost all the bad A36, page 3 of 17 A&A 665, A36 (2022) Fig. 1.Comparison of the original CERES pipeline with the adapted pipeline.Top: spectrum of the B9V star HD 33244 reduced with the original CERES pipeline.Red and purple arrows indicate the positions of wiggles and a gap, and yellow circles show the positions of some bumps (see the text for an explanation of these features).Bottom: spectrum of HD 33244 reduced with the new version of the CERES pipeline that includes the four changes described in Sect.3. columns disappear when the bias frame is subtracted from the observed science frames.The bias frame contains the same bad columns as the science frames, therefore when the first is subtracted from the latter, the excess in flux is removed and the bad columns are no longer visible in the 2D frames.This is the case for all but one bad column.Column 1299 is the brightest column in the 2D frames and does not disappear when the bias is subtracted.For the pixels in this column we increased the variance by 10 000 such that they had an extremely small weight and contributed almost nothing to the extraction process.
The original CERES pipeline corrects for the blaze function by dividing each individual order by the normalised flat field in that order.The normalised flat field is assumed to represent the blaze function according to Brahm et al. (2017), and is defined as the flat field divided by its maximum in that order.However, we realised that this blaze correction induces large intensity jumps or bumps where successive orders overlap (see yellow circles in the top panel of Fig. 1).When instead the orders are divided by the non-normalised flat field, they connect much better and bumps are no longer visible (see bottom panel of Fig. 1).Therefore, we changed the blaze correction in CERES to a division of each spectral order by the non-normalised flat field.
The CERES pipeline only outputs the 1D spectrum for individual orders.While this is fine for RV determinations, for which the pipeline was designed, we needed the whole spectrum to perform a spectral analysis.Therefore, we implemented order merging into the code.First we created a unified wavelength grid in logarithmic scale that covers the full wavelength range of all the extracted orders and has a relative step size of 10 −5 .This step size is approximately equal to 1/(2R) with R the resolution of FEROS (R = 48 000).Each order was then shifted to the new wavelength grid by linear interpolation.We merged the orders by taking a weighted sum, with the square of the signal-to-noise ratio (w = (S /N) 2 ) as weight, where f i is the blaze-corrected flux in order i, and σ f,i is its uncertainty.If pure Poisson noise is assumed, (S/N) 2 i is equal to the number of photons from the star in a certain order so that the weighted mean in Eq. ( 1) is just a coaddition of photons over all orders.This means that in overlapping neighbouring orders, where part of the photons had ended up in one or the other order, all photons are added again by Eq. ( 1).
We added two steps to the pipeline to remove sharp lines originating from cosmic rays and chip defects.In the first step we only dealt with positive hits due to cosmic hits.At every position in the spectrum we computed the median flux and standard deviation within a window of 100 pixels (∼5 Å).When the difference between the flux at that point and the median flux was larger than four times the standard deviation, the flux value was replaced with the median value.This process was repeated five times.The second step also removed negative hits due to chip defects, but here the method was more conservative to prevent removing stellar absorption lines.This is especially difficult in slowly rotating stars that have many narrow absorption lines that must not be mistaken for cosmic hits.The method is the same as for the positive cosmic hits, only the window size was changed to 30 pixels (∼1.5 Å) and both positive and negative hits were removed.Although emission lines have been detected in B-type stars (Wahlgren & Hubrig 2000;Alexeeva et al. 2016), we are confident that our cosmic hit removal function does not remove these real emission lines because they are expected to have the same width as the absorption lines in the spectrum while only sharp and narrow features are being removed by our method.of Tkachenko et al. (2013c) 4 .The computation of LSD profiles requires normalised spectra and spectral line masks.Therefore, we ran each individual spectrum through the ZETA-PAYNE setup, which is described in Sect.4.2, and obtained normalised spectra from which the response function is removed.The line masks were computed with the GSSP code (Tkachenko 2015), which relies on atmosphere models obtained with the LLmodels code (Shulyak et al. 2004), the radiative transfer code SynthV (Tsymbal 1996), and the Vienna Atomic Line Database (VALD; Kupka et al. 1999) atomic data.For each star we computed a GSSP model for its preliminary T eff and log g value returned by ZETA-PAYNE, and we fixed [M/H] to 0 dex.
We identified SB2 systems as targets for which two components could be distinguished in the LSD profiles, while in some cases also revealing different relative positions of the spectral lines of the two components for the different epochs of that target.The LSD profiles were also used to detect line profile variability (LPV), which can have various physical origins, such as chemical inhomogeneities or pulsations.Examples of LSD profiles for an SB2 system and a star with LPV are shown in Fig. 2. In total, 26 targets were classified as SB2 systems or as targets for which it was not clear whether the variability in the LSD profiles is due to binarity or LPV.Some of these systems were already identified as (spectroscopic) binaries in SIMBAD.These 26 targets were excluded from the sample because they require the (simultaneous) analysis of both components, which is not yet possible with our current method.They are listed as SB2 systems in Table A.1.
The remaining targets were divided into two groups based upon the LPV in the different epochs of each target.The first group consists of 42 stars with LPV, either in Hα, in all Balmer lines or also in metal lines.Two of these stars are Be stars with not only Balmer lines in emission, but also metal emission lines.The other 40 stars are classified as (blue) supergiants in SIMBAD, and all of them are members of the LMC.For this type of star, non-local thermodynamic equilibrium (NLTE) analysis including wind physics is required to obtain atmospheric parameters.Our current set-up of ZETA-PAYNE is trained on LTE model spectra for log g > 3 dex, while supergiants typically have log g values between 1 and 3 dex.Therefore it could not be applied to the spectra of the supergiants.Instead, these 42 stars will be studied in the complementary paper by Serebriakova et al. (in prep.) that focuses on OBA-type supergiant LMC members.
The other group contains 92 stars without clear LPV or with variability in metal lines but stable Balmer lines, as well as six Be stars with Balmer lines in emission.These 98 (Galactic) stars are further studied in this paper.Among them, we identified singlelined spectroscopic binaries (SB1) using the RV determinations from ZETA-PAYNE.A target was classified as SB1 system when the RV difference between any two epochs of that target is larger than the uncertainties of the RV values.If the RV differences are smaller than the uncertainties or when only one spectrum is available, the star was assumed to be single.Twelve targets were found to be SB1 systems, and 86 are single stars.For each of the targets, we shifted all of its spectra to RV = 0 km s −1 and added them by taking a weighted sum, where the weight of each epoch is the square of the S/N value of that epoch computed following the algorithm from Stoehr et al. (2008).We subsequently used that same algorithm to compute the S/N of the averaged spectrum; these values are reported in Col. 3 of Table A.2.The number of epochs and binary information are given in Cols. 2 and 4 of the same table.
We computed the binary detection probability for all our targets.We distinguished four different observing strategies to obtain the observations required for our study: (1) two observations were obtained within 1 week, and one observation 6 weeks later, (2) two observations were obtained spread over 2 months, (3) four observations were obtained: two during 1 week, and two collected during another week 2 months later, and (4) two observations were obtained on consecutive nights and one observation 2 months later.Figure 3 displays the binary detection probabilities in these different campaigns.To compute these probabilities, we ran one million Monte Carlo simulations tuned for each observing strategy.We took a uniform period distribution in logarithmic space ranging from 1 to 1000 days and eccentricity values between 0 and 0.95.All orbits with eccentricities lower than 0.03 were considered as circular.The mass-ratio distribution was considered uniform from 0.01 to 1 and the initial mass for each star was taken between 2 and 5 M ⊙ , that is, in agreement with their spectral classification (Schmidt-Kaler 1982).We adopted a Salpeter initial mass function (IMF, with α = −2.35)for the primary masses.Our simulations assumed that the orbits are randomly oriented in 3D space, the time of periastron passage is uncorrelated with respect to the start of the RV campaign, and that the orbital parameters are uncorrelated.To detect a system as binary, it needs to fulfil the same criteria as we adopted for our analysis.
For orbital periods shorter than 10 days, the probability to detect binary systems is rather high with a probability of 60% or higher, regardless of the observing strategy.For systems with periods between 10 and 100 days, we still reach a detection rate between 30 and 60%.This rate significantly drops for longerperiod systems, and it is almost not possible to detect binary systems with orbital periods longer than 1000 days.

Spectrum analysis with ZETA-PAYNE
We used the ZETA-PAYNE (Straumit et al. 2022) spectrum analysis algorithm to analyse the averaged spectra of the 98 Galactic single stars or SB1 systems.The ZETA-PAYNE is a machinelearning framework that trains a neural network on a grid of synthetic spectra to obtain a spectrum interpolator.This interpolator predicts theoretical flux spectra for an arbitrary set of stellar labels and can be used to fit any observed spectrum and derive its surface parameters as long as they are within the parameter ranges of the training sample.
In this paper, the neural network was trained with 1D LTE synthetic models computed with GSSP.Our choice of the LTE formalism is well justified given that the T eff distribution of the sample stars is found to peak at some 11 500 K, with only one target being hotter than 20 000 K (see Sect. 5.1).In their recent study of OBA-type stars in the LAMOST survey, Xiang et al. (2022) demonstrated that the impact of NLTE effects on the spectroscopically inferred T eff values is negligible for stars cooler than 25 000 K. In particular, the authors showed an overall good consistency between the Kurucz LTE and Potsdam Wolf-Rayet (PoWR; Hainich et al. 2019) NLTE models at T eff ⪅ 25 000 K (their Figs.C.1 and C.2), and reported a good agreement between the parameters obtained with the LTE and NLTE models for stars below that same T eff cut-off value (their Fig. 11).Instead of using a regular grid in five dimensions, one for each of the parameters we wish to determine, namely T eff , log g, global metallicity ([M/H]), projected rotational velocity (v sin i), and microturbulent velocity (ξ), we created a grid that is quasirandomly spaced using Sobol numbers (Sobol 1967).It covers the whole parameter space for main-sequence BAF-type stars as listed in the second column of Table 1, but with fewer grid points than a regular grid would need.It is therefore easier to compute.Straumit et al. (2022) pointed out that for high-resolution spectra, the neural network performs worse in regions of the parameter space with low values of v sin i and temperatures in the late-A to F-type regime.This is due to the many narrow spectral lines

Parameter
Quasi-random grid Oversampled grid T eff (K) 6000-25 000 N(6000, 4000) log g (dex) 3-5 3-5 v sin i (km s −1 ) 0-400 −0.8 to +0.8 −0.8 to +0.8 that have to be resolved by the neural network in this regime.To improve the performance of the neural network for such complex spectra, the parameter space should be oversampled with model spectra computed for low v sin i and low T eff values.The final training sample therefore consisted of 5000 quasi-randomly spaced model spectra and 5000 additional random model spectra with T eff and v sin i drawn from the normal distributions N(6000, 4000) and N(0, 25) respectively.An overview of the parameter ranges for the whole training grid is given in Table 1.
All the model spectra were computed within the wavelength range from 3000 to 10 500 Å and for infinite resolution (no instrumental broadening) such that the grid can also be used for optical spectra from spectrographs other than FEROS.During the fitting procedure the spectra were convolved to the resolution of the respective spectrograph by either using a constant resolution value or a wavelength-dependent resolution (the line spread function, LSF).We used the constant resolution value of FEROS because the LSF computed from the ThAr frames that were taken during the observing run changes from night to night.The neural network was trained in the same way as described in Straumit et al. (2022); more information can be found in that paper.For every star in the sample, surface parameters and the RV were derived by fitting model spectra to the observed spectrum using the neural network interpolator, performing a Doppler shift, and minimising the χ 2 merit function.However, the neural network was trained on continuum normalised spectra while the observed spectra contain components from the instrument (the response function), the interstellar medium (ISM), and the Earth's atmosphere.Each observed spectrum had to be divided by its response function before it could be compared to the normalised model spectra.In the ideal case, the response function of a spectrograph is known and can be removed from any observed spectrum to obtain only the component from the star itself (and some additional features from the ISM and the Earth's atmosphere, e.g.telluric lines).This can be attempted by dividing an observed spectrum by a flat-field exposure, but some curvature always remains in the spectrum that is not related to the stellar atmosphere (Gray 2005).In practice, the response function is often removed by fitting a spline or a polynomial to manually selected points that are assumed to be part of the continuum.This is a rather subjective way of normalising a spectrum and can give different results when it is done by different researchers.
In order to minimise subjectivity, we followed Straumit et al. (2022) and normalised the spectra automatically during the fitting procedure in ZETA-PAYNE, such that the optimal response function was determined simultaneously with the best-fitting surface parameters.We assumed that the response function could be represented by a Chebyshev polynomial for which we optimised its coefficients (as described in Straumit et al. 2022).The number of coefficients is a free parameter that had to be fixed before starting the fitting procedure.Preferably, the number is as small as possible to ensure that no features are introduced into the spectrum, but it should also be high enough to capture the whole A36, page 6 of 17 shape of the response function.It also depends on the wavelength range, where generally more coefficients are needed when a longer wavelength range is used.For our analysis, we selected a wavelength range from 4000 to 5800 Å because it includes three Balmer lines (Hβ, Hγ, and Hδ) and many metal lines, but excludes the Hα region that complicates the optimal normalisation of FEROS spectra.Since all the response functions of the FEROS spectra have a similar shape, the same number of Chebyshev coefficients could be used for all the stars.
We tested the optimal number of Chebyshev coefficients for a typical FEROS response function by taking five random stars from the sample for which we fitted the spectra with 10, 15, 20, 25, 30 and 35 coefficients with ZETA-PAYNE.For each fit we computed the χ 2 value and plotted those values as a function of the number of coefficients (see Fig. 4).The χ 2 profile converges towards a certain value, and this plateau is reached around 25 coefficients.Therefore, we used 25 Chebyshev coefficients to represent the FEROS response functions.

Internal uncertainties
We determined the internal uncertainties of the stellar parameters that are inherent to the modelling set-up.We created 1000 artificial FEROS spectra in the wavelength range from 4000 to 5800 Å, quasi-randomly spaced within the parameter space, but not overlapping with the training sample.Each spectrum was computed with GSSP, given a random RV shift within ±50 km s −1 , and convolved to the resolution of FEROS.A random response function was introduced, represented by a Chebyshev polynomial with 25 coefficients.For every spectrum we obtained a noiseless version and a version with S /N = 150, a typical value of the real FEROS spectra.All the artificial spectra were analysed with the fitting routine of ZETA-PAYNE, and the resulting parameters were compared with their real values.Outliers were identified as spectra for which the T eff difference was larger than four times the 1σ standard deviation of the T eff difference distribution of all 1000 spectra.This is formulated as the reliability fraction of the neural network (Straumit et al. 2022) and is given in the last row of Table 2.
Using only the remaining spectra without the outliers, we computed the internal uncertainty for each parameter as the mean difference between the real values and those obtained with ZETA-PAYNE.These uncertainties for noiseless spectra and for spectra with S /N = 150 are given in Table 2, and should be taken into account when determining parameters for the observed FEROS spectra.
In addition to the total internal uncertainties for the full T eff range of the grid, we also computed the internal uncertainties for the noiseless spectra in bins of T eff and v sin i.An overview of this can be found in Table 3 and is plotted in Fig. 5.One of the interesting features in this figure is the tail of high T eff uncertainties in the ∼9000 K regime.A closer look at the properties of spectra forming the tail reveals low [M/H] and high v sin i values, which in turn implies an apparent dearth of spectral lines of metals.This significant loss of information, where one has to rely exclusively on the Balmer lines to estimate the T eff -log g pair of parameters, results in larger uncertainties in T eff and, to a lesser extent, log g of the star.There is also an increase in absolute uncertainty towards high T eff values, yet relative uncertainties remain at the level of 2-3%, similar to the lower T eff regime (see Table 3).We also note scatter in log g for spectra with T eff < 10 000 K, that is, for late A-and F-type stars.As discussed in Gebruers et al. (2021), the log g determination for these stars is degenerate with T eff , [M/H], and the continuum normalisation of the spectrum.Small offsets in the continuum result in large differences in log g, which explains the larger internal uncertainty in this T eff bin.The same trends in the internal uncertainties are reported by Straumit et al. (2022) for APOGEE (infrared) and BOSS (optical) spectra.

Effect of the selected wavelength region
The atmospheric parameters of the main-sequence stars in our sample were derived using the spectral wavelength range 4000-5800 Å.To test whether the analysis of a longer or shorter wavelength range has any effect on the obtained parameters, we selected six stars distributed over the parameter space and applied ZETA-PAYNE on three different wavelength ranges.The wavelength ranges are 3870-6750 Å (Hα, Hβ, Hγ, Hδ, and Hϵ), 3870-5800 Å (Hβ, Hγ, Hδ, and Hϵ), and 4200-5800 Å (Hβ and Hγ).For each of these wavelength ranges, we computed the difference in parameters with those obtained from the original wavelength range (4000-5800 Å).This is shown in Fig. 6.All the values lie within the internal uncertainties from Table 2.This demonstrates that the choice of wavelength region has no statistically significant effect on the surface parameters when the internal uncertainties are taken into account.
We realised that the uncertainties on the surface parameters increase with longer wavelength regions.There is also a trend of increasing log g when a longer wavelength region is A36, page 7 of 17 A&A 665, A36 (2022) analysed, although all values are still within the internal uncertainties.This is caused by the normalisation of the spectrum, which becomes more difficult when more parts of the spectrum are included.Especially when Hα is included, a wavelength range containing multiple telluric lines is covered and ZETA-PAYNE cannot model these.For the stars in Fig. 6, we found that for longer wavelength regions, the wings of the Balmer lines become more broadened than when shorter wavelength regions are analysed with ZETA-PAYNE.This is demonstrated in Fig. 7 and shows that our normalisation approach with a Chebyshev polynomial of 25 coefficients is not optimised to normalise a spectrum that is extended to Hα or to Hϵ, and that the deterioration in normalisation of the Balmer wings is compensated for by preferring a model with higher log g.Additional tests for the number of Chebyshev coefficients and methods to remove tellurics are needed to obtain a better normalisation in each specific wavelength range.Here we decided to limit our analysis to the 4000-5800 Å region because not that much information is added by extending to Hϵ and normalisation becomes much more uncertain when Hα is included.

Spectroscopy
The results for the 98 stars analysed with ZETA-PAYNE are given in Table A.2.The statistical uncertainties inferred from using ZETA-PAYNE (σ) and those that include the internal uncertainty (σ i ) are also provided.We advise to use the latter values since they are more representative for the typical precision with which the surface parameters of B-type stars can be derived.The distributions of the surface parameters are shown in Fig. 8.Most of the stars have parameters indicative of late B-type, with temperatures in the range 10 000-14 000 K. There are some cooler stars (early A-type) in the sample and a few hotter stars with T eff > 16 000 K. The log g distribution shows that the majority are main-sequence stars.However, there are six stars with log g < 3.5 dex, indicating that they are near the end or just beyond the main sequence.The v sin i values are fairly uniformly distributed within 0-350 km s −1 and [M/H] shows a normal distribution around the solar value with a mean of (−0.12 ± 0.27) dex.We note that the metallicity [M/H] in our case refers to all chemical elements with a mass higher than that of helium.All the fits with ZETA-PAYNE were visually inspected.For 20 stars we realised that the model predicted by ZETA-PAYNE as the best-fitting model does not fully correspond to the observed spectrum.In 13 of these cases, the observed spectrum contains more absorption lines than the predicted spectrum.These stars were identified as chemically peculiar, and specifically HgMn stars, which are further discussed in Sect.6.1.Three stars (HD 40435,HD 45527,and HD 63928) show strong LPV in their spectra, and four others (CPD-60 944B, HD 32493, HD 59426, and HD 6783) have very narrow spectral lines with strongly variable depths between different epochs.These stars are indicated as LPV and LPV* in Table A.2. From the TESS light curves discussed in the following section, we find that these seven stars are rotational variables.This means that they are likely chemically peculiar (i.e.ApBp stars), which explains the limited agreement with the best-fitting model spectra.In many of these stars Cr is overabundant, confirming the chemically peculiar nature.An indepth chemical abundance analysis is beyond the scope of this work.

Variability classification from TESS photometry
We analysed the SAP and PDC-SAP light curves independently for all stars to gauge the possible effect of systematics, and we verified that the allocated aperture mask and background pixels by the SPOC pipeline are reasonable.It has been shown that, for example, SAP light curves are preferable in the identification and analysis of bright eclipsing binaries (see e.g.Southworth et al. 2020Southworth et al. , 2021)).For each star, we calculated the amplitude spectrum of its SAP and PDC-SAP light curves using a discrete Fourier transform (DFT; Kurtz 1985) up to the Nyquist frequency, which for the two-minute TESS data is 360 day −1 .A high diversity and incidence of photometric variability occurs in massive stars, A36, page 8 of 17  with common causes including pulsations, rotational modulation, winds, and binarity (Pedersen et al. 2019;Bowman et al. 2019aBowman et al. ,b, 2020;;Burssens et al. 2020).Based on a group of expert classifiers, we visually inspected the light curves and amplitude spectra and determined the dominant variability type(s), which is included in Table A.1.
Stars pulsating in coherent p-and/or g-modes are particularly interesting for follow-up asteroseismic modelling (e.g.Moravveji et al. 2015Moravveji et al. , 2016;;Szewczuk & Daszyńska-Daszkiewicz 2018;Szewczuk et al. 2022;Pedersen et al. 2021), as are (pulsating) eclipsing binaries (e.g.Tkachenko et al. 2020;Sekaran et al. 2020Sekaran et al. , 2021) )   by gravity waves (e.g.Bowman et al. 2019bBowman et al. , 2020)).We list all variability types used in this work below, which are not mutually exclusive: -rot: rotational modulation; -SPB: slowly pulsating B star (i.e.predominantly g-mode pulsations); -const: constant star (i.e.no significant variability); -EB: eclipsing binary; -EV: ellipsoidal variable; -β Cep: early-B star (i.e.predominantly p-mode pulsations); -δ Sct: A-and early-F star with p-and g-mode pulsations; -SLF: stochastic low-frequency variability.Stars for which the variability type is also assigned a question mark have a tentative classification.In some cases, the inferred T eff and log g values from our spectroscopic analysis were also used to separate δ Sct and β Cep pulsators, for example, which can have similar amplitude spectra.
We show the TESS light curves and amplitude spectra of four example stars in Fig. 9 to demonstrate the power of including TESS light curves in our analysis.First is the SPB star HD 33599 (TIC 55295028), which has the typical signatures of multi-periodic g-modes appearing in distinct groups in its A36, page 10 of 17 amplitude spectrum (Kurtz et al. 2015;Van Beeck et al. 2021).Its light curve also has signatures of outbursts, and hence it is a candidate pulsating Be star.We also note that we have identified HD 33599 as an SB2 system with many emission lines in its spectra, which makes it interesting to follow this star up spectroscopically and understand the nature of its companion and evolutionary origin (see e.g.Shenar et al. 2020;Bodensteiner et al. 2020).Second, is the SPB/β Cep hybrid star HD 49193 (TIC 167523976), which has v sin i ≃ 40 km s −1 and T eff ≃ 21 000 K. It is the hottest star in our sample and is located close to the cool edge of the β Cep instability region.Third is the EB system AN Dor (TIC 220430912), which was recently revisited by Southworth & Bowman (2022) to derive masses and radii from archival RVs and new TESS light curves.Fourth is HD 40435 (TIC 31313111), which has the characteristic doublewave light curve and harmonic-series in its amplitude spectrum typical of rotation modulation (see e.g.Bowman et al. 2018).We identified LPV in the spectra of this star that is likely associated with an ApBp nature.
As a vast improvement over previous studies that used short TESS light curves, our reclassification has the main advantage of improved frequency resolution from long-term TESS light curves.Similar to previous work, we find a high variability fraction (>90%) in the TESS light curves of our sample.As mentioned previously, we did not extract pulsation, binary, or rotation frequencies for the follow-up modelling because this endeavour requires custom light curves and is the subject of future work.Our work has allowed us to rank the most promising stars with reliable spectroscopy and clear signatures of binarity and/or pulsations to follow them up with detailed seismic modelling.

Newly discovered HgMn stars
Thirteen of the main-sequence stars in our sample analysed with ZETA-PAYNE were found to be new HgMn stars.We optimised the Hg and Mn abundances of the 13 stars to confirm their HgMn classification.The wavelength range has been extended to 3930 Å in the blue part of the spectrum to include the Hg II 3984 spectral line.The differences between the observed and the model spectra around the Hg and Mn lines disappear when the abundances of these elements are optimised in the models (left versus right panels of Figs. 10 and 11).For some stars, such as HD 57808, some features are still present in the residuals after optimising for Hg and Mn.This is probably caused by overabundances of other elements such as Y, but a detailed abundance analysis is beyond the scope of this paper.
From the TESS light curves, we find that all 13 HgMn stars exhibit rotational modulation and none of them are pulsators.For a few targets, line profile variability is also found in the spectra (this is indicated in column 4 of Table A.2). Interestingly, only one of them is found in an SB1 binary: HD 44247.This would result in a binarity rate of only 8% (without correcting for observational biases) as opposed to what is found in other studies.Our observational biases include but are not limited to an intentional bias against binary stars at the very early stage of the sample selection, exclusion of all stars that were nevertheless identified as SB2 binary systems in our sample prior to its detailed spectroscopic analysis based on the obtained FEROS spectra, and our ability to detect (especially long-period) binaries is limited by the number of acquired epochs (two in most cases), the separation of some two months between them, and the S/N in the obtained individual epoch spectra.This last bias is represented by the binary detection probabilities quoted in Sect.4.1.Therefore, the above-quoted observed binary rate among HgMn stars in our sample is unreliable.None of these 13 stars have been classified as HgMn stars before, but µ Men was found to have an overabundance of Si (Houk & Cowley 1975)

(Spectroscopic) Hertzsprung-Russell diagram
The stars from Table A.2 are plotted in a spectroscopic Hertzsprung-Russell diagram (sHRD) in Fig. 12 in which the ordinate axis is defined as log (T 4 eff /g).In the top panel, the stars are colour-coded following their spectroscopic classification, and in the bottom panel, the symbols represent different photometric variability classes.Typical errors are shown, both for the uncertainties obtained from ZETA-PAYNE and when the internal uncertainties are taken into account.We used the non-rotating MESA evolutionary tracks for solar metallicity, an exponential core overshoot of f ov = 0.02, and the g-mode instability region from Burssens et al. (2020) to situate the stars in the sHRD.Most stars are intermediate-mass stars on the main sequence with masses around 2.5-4 M ⊙ .Two to six stars, however, depending on which errors are taken into account, appear to have evolved beyond the end of the main sequence.
We also include Fig. 13, which shows the same stars in the HRD, with luminosities calculated based on their Gaia G-band magnitudes, Gaia eDR3 photogeometric distances from Bailer-Jones et al. ( 2021), reddening values from the SFD dust map (Schlafly & Finkbeiner 2011), and bolometric corrections calculated with the prescriptions from Pedersen et al. (2020).In addition to the g-mode instability region from Burssens et al. (2020), we also plot the blue edge of the δ Sct instability strip from Dupret et al. (2005).Most of the stars classified as 'SPB' or 'SPB?' lie within the SPB instability strip.A large part of the sample is located outside of the SPB and δ Sct instability regions.This includes constant stars and stars with rotational modulation, but also 10 pulsators.We also added to Fig. 13 stars from the recent work by Sharma et al. (2022).These authors studied pulsating B-type stars in the Scorpio-Centaurus association using TESS light curves and also found pulsators below the SPB instability strip.Low-frequency pulsations in B-type stars cooler than the red edge of the SPB instability region have been reported before.Salmon et al. (2014) claimed that these are fast-rotating SPB stars that are seen equator-on, and due to gravity darkening, they appear cooler and fainter than if they were observed pole-on.A similar explanation is given for high-frequency pulsators found in the literature between the red edge of the β Cep instability strip and the blue edge of the δ Sct instability strip, although it is also possible that these are binaries or misclassified stars that are variable due to rotational modulation (Balona et al. 2015(Balona et al. , 2016)).
Our sample contains eight stars that are classified as SPB and are situated below the SPB instability strip.Their g-mode frequencies are higher than what is normally expected for these stars.These stars are HD 45835, θ Ret A, κ Men, HD 49531, µ Pic, HD 47478, 29 Dor, and HD 37027, and none of them are spectroscopic binaries.µ Pic, 29 Dor, and HD 37027 are known to be Be stars.We did not include the cores of the Balmer lines in the spectral analysis for Be stars because they are in emission.Therefore the T eff value that is derived from the Balmer lines might be subject to large uncertainties.θ Ret A is a relatively slow rotator with v sin i = 44 ± 13 km s −1 , but is situated relatively close to the red edge of the SPB instability strip.The other seven stars are rapid rotators and might thus be fast-rotating SPB stars seen equator-on, as explained by Salmon et al. (2014).Their high g-mode frequencies can be explained by fast rotation that has shifted their g-mode frequencies as observed in an inertial frame into a relatively high-frequency regime.This astrophysical interpretation has been proposed previously for bright SPB and Be stars observed from ground-based data (Aerts et al. 1999;De Cat & Aerts 2002;Aerts & Kolenberg 2005;Saesen et al. 2010Saesen et al. , 2013;;Mowlavi et al. 2013Mowlavi et al. , 2016;;Moździerski et al. 2014Moździerski et al. , 2019)), MOST space photometry (Walker et al. 2005;Aerts et al. 2006;Saio et al. 2007), CoRoT space photometry (Neiner et al. 2009;Diago et al. 2009), and Kepler/K2 space photometry (Pápics et al. 2017;White et al. 2017;Pedersen et al. 2021;Szewczuk et al. 2021).In fast-rotating stars, major frequency shifts are indeed expected for prograde and retrograde A36, page 12 of 17 sectoral modes (i.e.modes with an equal angular degree l and azimuthal wavenumber |m|) in the gravito-inertial regime.In this regime, g-mode oscillations are restored by both buoyancy and the Coriolis acceleration together (Townsend 2005;Salmon et al. 2014;Saio et al. 2017;Aerts et al. 2019).Our spectroscopic results and interpretation are fully in agreement with the findings by Gaia Collaboration (2022) based on Gaia DR3 time-series photometry.
Because of the rapid rotation of these SPB stars shown by the high v sin i values inferred from our FEROS spectroscopy, their apparently too low T eff relative to the SPB instability region is likely to be due to the gravity-darkening phenomenon.This phenomenon results from the variation in T eff with local gravity on the surface of a non-spherical star (von Zeipel 1924) and therefore implies a lower T eff at the equator than at the pole of a rapidly rotating star.As was recently demonstrated in Fabry et al. (2022, and private communication), the variation in T eff of a typical rapidly rotating B-type star is about 5% in the course of its main-sequence evolution.This T eff variation from stellar pole to equator can explain the apparent temperature offset from the red edge of the SPB instability strip in the HRD for all but three stars.The stars for which gravity darkening is not sufficient to place them within the SPB instability strip are two Be stars HD 37027 and 29 Dor and the star HD 47478.However, the 5% T eff difference derived in Fabry et al. (2022) is based on MESA models, which are limited to 1D approximations.A 2D study with the ESTER code (Bouchaud et al. 2020) has shown that for the A7-type star Altair, the temperature difference between equator and pole is approximately 2000 K. Thus 5% is only a lower limit, and gravity-darkening effects can be higher when the proper 2D structure of fast rotators is taken into account instead of 1D approximations.We also find four apparently p-mode pulsators (identified as δ Sct stars because their TESS light curves exhibit pulsation frequencies above ∼4 day −1 ) in between the β Cep and δ Sct instability regions (HD 44533, HD 55478, HD 67420, and η Phe).These four stars are all fast rotators, and HD 44533 and η Phe are Be stars.They are probably misclassified as δ Sct stars, just as the high-frequency mode near ∼4 day −1 of the rapidly rotating SPB pulsator HD 43317 was initially thought to be a p-mode (Pápics et al. 2012) while it actually concerns a rotationally shifted retrograde quadrupole mode (Buysschaert et al. 2018).In the case of the Be stars, their derived surface parameters are uncertain because of their fast rotation and circumstellar disk.

Multiple linear regression
We searched for correlations between all five spectroscopic parameters and the luminosity (L) via backward multiple linear regression.This was done by taking one of these six parameters as the dependent variable and performing a multiple linear regression with the other five parameters as independent variables.We computed the coefficient of determination (R 2 ) to determine how well the set of independent variables could predict the dependent variable, where R 2 = 1 would mean a perfect fit.For each of the independent variables we then verified if they were significant in this relation through a t-test.If the p-value from this t-test was above 0.05, the variable was assumed to be insignificant.The multiple linear regression was then repeated in the same way, but without the least significant variable, until only significant variables remained.The relations we found from this analysis are given in Table 4.
In the first three models, log L, T eff , log g, and [M/H] are related to each other.These are well-known relations expected from standard stellar structure and evolution theory.The fourth relation links ξ to log L, T eff , v sin i, and [M/H].It suggests a decrease in ξ with increasing luminosity within our sample.This contradicts the expectation that turbulent velocity fields increase as the log g (luminosity) value of star decreases (increases) (e.g.Gray 2005).This latter interpretation is supported by observational studies of single stars (e.g.Hunter et al. 2008) and binary systems (e.g.Tkachenko et al. 2014), as well as by theoretical studies of sub-surface convective zones in intermediateand high-mass stars (e.g.Cantiello et al. 2009).We currently lack a good physical interpretation of the (weak) opposite trend we found, but we speculate that it may have a methodological explanation in that ξ in B-type stars cannot be reliably inferred from the global fitting procedure of their entire optical spectra.Instead, a detailed analysis of carefully selected spectral lines that are sensitive to ξ might be required, such as the Si III triplet at 4552 Å, 4567 Å, and 4574 Å and/or He I lines in the optical part of the spectrum.This is different from cooler stars with many metal lines for which ξ can be determined from the entire spectrum because ξ values from all element lines, including possible incorrect values from lines that are less sensitive to ξ, are averaged out.
We also noted by eye a trend of decreasing v sin i with age in Fig. 13, especially for stars with a mass below 4 M ⊙ .However, by performing a linear regression analysis for v sin i and stellar age, which we obtained by interpolating the MESA evolutionary tracks at the positions of the sample stars, we concluded that this correlation is not statistically significant (R 2 = 0.02 and R 2 = 0.04 for M < 4 M ⊙ ).

Conclusions
We have analysed high-resolution FEROS spectra of a sample of 166 variable stars identified using TESS space photometry.Our follow-up spectroscopic data were reduced with the CERES pipeline, to which we added four new capabilities to obtain smooth reduced normalised spectra without artificial features and cosmic hits.
We have determined surface parameters with the ZETA-PAYNE framework for all Galactic single stars and SB1 systems in the sample, excluding SB2 systems and supergiants in the LMC galaxy.Except for the uncertainties deduced with ZETA-PAYNE, the internal uncertainties inherent to the neural network and fitting approach were also taken into account.It was shown that within the internal uncertainties, the parameters are not affected by the wavelength region that is analysed.If the neural network used in this paper is applied to longer wavelength regions, the number of coefficients for the Chebyshev polynomial representing the response function should be optimised for the respective wavelength range.Overall, in the considered wavelength interval from 4000 Å to 5800 Å and with 25 Chebyshev coefficients representing the residual response function of the instrument, the estimated internal uncertainties reach about 3% and 5% in T eff and v sin i, respectively, and ∼0.1 dex in [M/H] and log g.These values give an upper limit on the actual parameter uncertainties of the sample stars and reach the level required for high-precision asteroseismic modelling.For this, the parameter uncertainties are ideally as small as possible, but should definitely be below ∼3% in T eff , 5-10% in v sin i and ∼0.1 dex in log g and for individual abundances (Moravveji et al. 2016;Michielsen et al. 2021;Mombarg et al. 2022).To obtain better constraints on the spectroscopic parameters, information from different data sources should be combined and analysed together.For example, Gaia parallaxes and magnitudes can help constrain the otherwise degenerate parameters T eff and log g.This is beyond the scope of this paper, however, and we leave this endeavour for future work.
We have found that 13 of the stars in the sample are chemically peculiar HgMn stars.The best-fitting spectra returned by the ZETA-PAYNE framework do not fit the observed spectra of these stars well, specifically in spectral windows that contain lines of Hg, Mn, Y, and some other chemical elements.Our ad hoc abundance analysis of several chemical elements in the spectra of these HgMn stars with the GSSP software package revealed strong overabundances, confirming the peculiar nature of the stars.
The spectroscopic values and luminosities computed from Gaia eDR3 data were used to place the stars in a spectroscopic and classical HRD.In these diagrams, eight stars classified as SPB stars with relatively high g-mode frequencies are located below the (non-rotating) SPB instability region.In all cases, the location and variability of these stars can plausibly be explained by their fast rotation, which is expected to significantly affect the instability regions of early-type stars (see Szewczuk & Daszyńska-Daszkiewicz 2017) and shift their observed g-mode frequencies into the high-frequency regime.Therefore, we conclude that we found spectroscopic evidence of a group of fast-rotating g-mode B-type pulsators monitored by TESS, covering the area between the SPB instability strip and the blue edge of the δ Sct strip.For all but three stars, the lower temperatures relative to the SPB instability strip obtained for these stars could be related to the gravity-darkening effect.
The multiple linear regression analysis involving five spectroscopic parameters and stellar luminosity revealed well-known A36, page 15 of 17 A&A 665, A36 (2022) relations between log L, T eff , log g and [M/H] of the star as predicted by stellar evolution theory.A somewhat unexpected result is the weak negative relation found between log L and ξ, which we attribute to our method of inferring ξ from the entire optical spectrum instead of selecting specific lines sensitive to ξ.
Our spectrum fitting approach with the current neural network trained on LTE GSSP models is not suitable for the analysis of the spectra of supergiants.This will only be possible when NLTE models that can include physics of stellar winds are used to train the neural network.Similarly, it is not yet possible to analyse SB2 systems with this version of ZETA-PAYNE because this requires the analysis of either their composite or disentangled spectra.Analysis of the composite spectra of SB2 binaries is challenging and typically associated with large uncertainties in the inferred atmospheric parameters of the two binary components (see e.g.Tkachenko 2015).On the other hand, analysis of the disentangled spectra is as precise as spectrum analysis of spectra of single stars, provided a uniform orbital phase coverage is achieved for the SB2 system in consideration (see e.g.Hadrava 1995;Ilijic et al. 2004).Therefore, the detailed analysis of the spectra of evolved supergiants and of SB2 binary systems requires a dedicated effort and is beyond the scope of this paper.

Fig. 2 .
Fig. 2. Examples of LSD profiles.Top: LSD profiles for two epochs of the SB2 system HD 62153.Bottom: LSD profiles for two epochs of the single star HD 41297 with LPV.

Fig. 3 .
Fig.3.Binary detection probabilities as a function of the orbital period for four different observing strategies used to constitute our dataset: (1) two observations were obtained within 1 week and one observation 6 weeks later, (2) two observations were obtained spread over about 2 months, (3) four observations were obtained: two during 1 week, and two collected during another week 2 months later, and (4) two observations were obtained on consecutive nights and one observation 2 months later.

Fig. 4 .
Fig. 4. χ 2 values for best-fitting synthetic spectra to the observed spectra of five sample stars as a function of the number of Chebyshev coefficients.

Fig. 5 .
Fig. 5. Internal uncertainties for 1000 synthetic noiseless FEROS spectra.The left panels show the difference between parameters predicted with ZETA-PAYNE and the true values of the synthetic spectra as a function of T eff and the right panels show this difference as a function of v sin i.The orange dots are the average uncertainties with standard deviation in bins of T eff and v sin i.

Fig. 6 .
Fig. 6.Comparison of stellar parameters derived from different wavelength regions.From left to right the differences in T eff , log g, [M/H], v sin i and ξ are shown.Black, red, and yellow dots are the differences between the parameter value obtained from the 4000-5800 Å spectral range and the parameter value obtained from the wavelength range indicated with the respective colour in the legend of the figure, so p 4000−5800Å -p ∆λ .The pink region indicates the internal uncertainty interval from Table2.The systematic offset (insignificant within the internal uncertainty) in the log g parameter is discussed in detail in Sect.4.4.
Fig. 7. Comparison of the normalisation quality when different wavelength regions are analysed.Top panel: spectra computed in the 4200-5800 Å (black), 4000-5800 Å (purple), 3780-5800 Å (red), 3780-6750 Å (yellow) wavelength regions are shown.The bottom panel shows the differences between the spectrum computed in the 4000-5800 Å range and those computed in the other wavelength regions as indicated by the same colours as in the top panel.

Fig. 8 .
Fig. 8. Distributions of T eff , log g, [M/H], and v sin i for the sample of single and SB1 Galactic stars.

Fig. 9 .
Fig. 9. Subsection of the two-minute TESS light curve and amplitude spectrum of four example stars.From top to bottom: TIC 55295028 (HD 33599; SPB), TIC 167523976 (HD 49193; SPB / β Cep), TIC 220430912 (AN Dor; EB), and TIC 31313111 (HD 40435; rot).For the latter two stars, a phase-folded light curve is shown in the top right panel.
Fig. 10.Overabundance of Hg in HgMn stars.The left panel shows residuals around the Hg II 3984 line between the observed spectrum and best-fitting GSSP spectrum with the Hg abundance according to the metallicity of the star.The right panel shows the residuals when the Hg abundance is optimised for the star.The residuals in the right panel are colour-coded according to the Hg overabundance with respect to the abundance in the left panel.The value is also given for each star.The residuals of the different stars are offset with a constant value to plot all stars in the same figure.

Fig. 11 .
Fig. 11.Overabundance of Mn in HgMn stars.Similar as in Fig. 10, but for the Mn II 4252 and Mn II 4253 lines.

Fig. 12 .
Fig. 12. sHRD containing all stars from Table A.2. Top: stars are colour-coded according to their spectroscopic information.The hatched symbols are SB1 binaries.Bottom: same as in the top panel, but now the photometric variability is indicated by different symbols.They are colour-coded according to their v sin i values.Stars with dominant rotational variability are hatched with white horizontal lines, and symbols hatched with white vertical lines show stars that can either have variability according to the symbol or rotational variability.The right panels show a zoom in of the region below the SPB instability strip.Typical errors are plotted, both the statistical error returned by ZETA-PAYNE and the error including the internal uncertainties.Non-rotating MESA evolutionary tracks for different initial stellar masses (in grey) and the g-mode instability strip (red shaded area) from Burssens et al. (2020) are shown as well.
Fig. 13.HRD.Same figure as the bottom panel of Fig. 12, but with the luminosity on the y-axis.The small grey dots are B-type stars in the Scorpio-Centaurus association from Sharma et al. (2022).Non-rotating MESA evolutionary tracks from Burssens et al. (2020) are shown in grey, and the red shaded area is the g-mode instability strip from Burssens et al. (2020).The dashed purple line shows the blue edge of the δ Sct instability strip from Dupret et al. (2005).

Table 1 .
Parameter ranges of the training grid.

Table 2 .
Average internal uncertainties for the surface parameters and the RV computed from 1000 simulated FEROS spectra, without noise and for S /N = 150.

Table 3 .
Internal uncertainties for the surface parameters and RV in bins of T eff and v sin i computed from a sample of 1000 simulated noiseless FEROS spectra.

Table 4 .
A36, page 14 of 17 S. Gebruers et al.: Analysis of FEROS spectroscopy for a sample of variable B-type stars with TESS photometry Multiple linear regression results.