Identification and Spectroscopic Characterization of 128 New Herbig Stars

We present optical spectroscopy observations of 145 high-mass pre-main-sequence candidates from the catalog of Vioque et al. 2020 From these, we provide evidence for the Herbig nature of 128 sources. This increases the number of known objects of the class by ∼50%. We determine the stellar parameters of these sources using the spectra and Gaia EDR3 data. The new sources are well distributed in mass and age, with 23 sources between 4 and 8 M ⊙ and 32 sources above 8 M ⊙. Accretion rates are inferred from Hα and Hβ luminosities for 104 of the new Herbigs. These accretion rates, combined with previous similar estimates, allow us to analyze the accretion properties of Herbig stars using the largest sample ever considered. We provide further support to the existence of a break in accretion properties at ∼3–4 M ⊙, which was already reported for the previously known Herbig stars. We re-estimate the potential break in accretion properties to be at 3.87−0.96+0.38 M ⊙. As observed for the previously known Herbig stars, the sample of new Herbig stars independently suggests intense inner-disk photoevaporation for sources with masses above ∼7 M ⊙. These observations provide robust observational support to the accuracy of the Vioque et al. 2020 catalog of Herbig candidates.


INTRODUCTION
The group of intermediate-to high-mass pre-main sequence (PMS) sources (M>1.5 M ) play a particularly important part in understanding the differences in the formation of, and accretion of matter onto, low-and high-mass stars.Whereas the formation of lowmass stars is widely accepted to be due to magnetically controlled accretion (Bouvier et al. 2007), higher mass stars are generally non-magnetic, and it would be expected that this scenario does not apply to them (see Mendigutía 2020).Recent studies (e.g.Wichittanakom et al. 2020;Grant et al. 2022) have provided evidence for the change in accretion mechanism happening at around 4 M .
In addition, the protoplanetary disks of intermediateto high-mass PMS stars show significant differences with respect to the disks around low-mass stars.The more luminous stars promptly photoevaporate their inner disks with FUV photons (Kunitomo et al. 2021), causing large inner cavities.The impact this has on the thermal and chemical evolution of the disks may significantly impact the formation and evolution of planets (Panić & Min 2017;Miley et al. 2021).Mm-wavelength observations have also shown that the disks around intermediatemass stars are more massive than those around low-mass stars (Andrews et al. 2013;Pascucci et al. 2016;Stapper et al. 2022) and have a higher fraction of detected structures in the dust disk (van der Marel & Mulders 2021; Stapper et al. 2022).Furthermore, there is evidence pointing toward particular disk structures being favored in more massive objects (e.g., spirals have been mostly found in early spectral type stars, Garufi et al. 2018).Hence, it has been theorized that they have a higher incidence of giant planets (e.g.Reffert et al. 2015;van der Marel et al. 2021).Therefore, the disks around intermediate-to high-mass PMS stars play a key role for understanding the structures and evolution processes of protoplanetary disks, and their link to planet formation.
The intermediate-to high-mass PMS regime comprises different types of young stellar objects (YSOs).The Herbig Ae/Be group contains stars at the latest stages of pre-main sequence evolution.In total, around 255 Herbig Ae/Be sources have been historically considered and studied (e.g.Vioque et al. 2018), although a smaller fraction has been properly characterized and is free of contaminants (see the seminal papers of Thé et al. 1994, Vieira et al. 2003, and Hernández et al. 2004).Recent compilations and studies about the general properties of Herbig Ae/Be stars can be found in Vioque et al. (2018) and Guzmán-Díaz et al. (2021).The cooler predecessors of the Herbig Ae/Be stars are occasionally referred to as Intermediate-Mass T Tauris (IMTTs).The defining difference between both groups is arbitrary in the literature and different thresholds in spectral type and mass have been used (e.g.Calvet et al. 2004;Povich et al. 2016;Villebrun et al. 2019;Nuñez et al. 2021;Valegård et al. 2021).Valegård et al. (2021) compiled most IMTTs within 500 pc, describing the general properties of a sample of 49 sources.Higher-mass PMS objects (M>8-10 M ) are generally called Massive Young Stellar Objects (MYSOs).However, the difference between the former two categories and the MYSOs is also ambiguous, often depending on the optical visibility of the sources.The most comprehensive catalogue of MYSOs to date can be found in Lumsden et al. (2013) with several hundred sources, although only a very small fraction of those have been characterized in great detail (e.g.Frost et al. 2021;Koumpia et al. 2021).
A large caveat of all the aforementioned results is that there is no homogeneous or complete survey of intermediate-to high-mass PMS stars.The existing samples are small heterogeneous collections of often randomly discovered objects, with a non-negligible amount of contaminants.This is a direct consequence of the Initial Mass Function (IMF): there are few nearby high-mass stars in any given molecular cloud (e.g., there are three stars with M > 1.5 M in the ALMA survey of Lupus, Ansdell et al. 2016).Any global conclusions that have been drawn are therefore subject to an unknown bias.Furthermore, as more massive stars evolve faster, the younger, high-mass PMS stars have been barely considered by the vast majority of studies (e.g. in the samples of van der Marel et al. 2021 andStapper et al. 2022 all sources are older than ∼ 4 Myr).This has important consequences on the conclusions obtained so far for these objects.For example, at those late stages, planet formation is mostly over in the low-mass regime (Cieza et al. 2021), so it is probable that in general we have not been tracing actively planet-forming stars.
To properly investigate the evolution of intermediateand high-mass PMS stars, what is needed is a wellselected sample covering a large range in age and mass.Vioque et al. (2020) produced a large homogeneous catalogue of a few thousand intermediate-to high-mass PMS ('Herbig star'1 ) candidates by combining large scale photometric surveys (Gaia, 2MASS, WISE, IPHAS, and VPHAS+; covering from the optical to the mid-IR and including Hα photometry) with machine learning techniques.
In this paper, we present spectroscopic observations for 145 Herbig candidates from the catalogue of Vioque et al. (2020).We discuss the results of the observations and present a comprehensive list of 128 newly confirmed Herbig detections.We start the paper with a description of the observations in Sect. 2. For each observed object, we characterize the extinction and the stellar parameters in Sect.3.After discarding some contaminants in Sect.4, in Sect. 5 we present accretion rates for the fraction of the sources with Hα or Hβ emission line measurements.In Sect.6 the Herbig nature of the observed sources is assessed.In Sect.7 the derived accretion rates are analyzed, and we compare them with the accretion rates of previously known Herbig stars.We discuss our results in Sect.8 and conclude in Sect.9.

OBSERVATIONS
A total of 145 Herbig candidates from the Vioque et al. (2020) catalogue were observed in low-to mediumresolution optical spectroscopy during three different ob- ).These 145 sources were selected because their absolute magnitudes suggest that their stellar masses cover the Herbig mass regime in a representative manner.In order to obtain precise stellar parameters, we targeted sources with accurate parallaxes ( , all but five have σ( )/ ≤ 0.1).This is because the parallax error dominates the uncertainty of the stellar parameters when these are obtained from the locations of the objects in the HR diagram (see e.g.Vioque et al. 2018).
The signal-to-noise ratio of these spectra are typically on the order of 100.Observing dates and instrumental setups are detailed in Table 1.The three observing runs can be summarised as follows: • 56 Herbig candidates were observed with the Intermediate Dispersion Spectrograph (IDS) instrument which is at the Cassegrain focus of the 2.54metre Isaac Newton Telescope (INT).The INT is located at the Roque de los Muchachos Observatory in the island of La Palma, Spain.Two different configurations were used.One block used the R900V diffraction grating, which covers the ∼ 3600 − 5000 Å spectral range.The other block used the R1200R grating (∼ 5700 − 6700 Å).
• 50 Herbig candidates were observed with the ESO Faint Object Spectrograph and Camera (v.2) or EFOSC2 in two settings.EFOSC2 is installed at the Nasmyth B focus of the 3.58-metre New Technology Telescope (NTT) at La Silla Observatory, Chile.The first block used the G7 diffraction grism (∼ 3300−5300 Å) and the second block used the G20 grism (∼ 6000 − 7200 Å).
A log of the observations is presented in Table B1.Bias, flat and arc frames were taken each night for the reduction of the observations.Standard procedures were used in order to process the data, which was reduced using the Image Reduction and Analysis Facility (IRAF).We started with bias subtraction.Next, flat field division was used to correct the pixel-to-pixel variation of the CCD signal.Then, one-dimensional spectra were extracted and sky subtracted from science frames.Last, arc frames of Cu-Ar and Cu-Ne comparison lamps were used to obtain the wavelength calibrated normalized spectra.

Comments on the blue spectral range
The blue spectral region considered in the three observing runs (∼ 3300 − 5400 Å) covers the main wavelength range to determine spectral types.An example of the normalized spectra obtained in this region is shown in Fig. 1.This region is especially useful for the earlier spectral type stars such as A and B stars, as their spectrum in the wavelength range beyond 5000 Å is fairly line free.This region also covers the Hβ line, which is a common tracer of circumstellar activity in YSOs.
The chosen grisms allow for efficiently obtaining spectral types and effective temperatures (T eff ).These are tabulated in Table B2.For some objects the temperature was estimated directly from model fitting the spectra (see Wichittanakom et al. 2020).Otherwise spectral types were obtained by comparison with model spectra (BOSZ models, Bohlin et al. 2017) and published spectral standards (Digital Atlas by Gray2 ), and the T eff in Table B2 are the values that correspond to those spectral types according to Pecaut & Mamajek (2013).In this latter case, the T eff uncertainties are of one subspectral type.The procedure followed for each object is detailed in Table B2.T eff could not be estimated for VOS 2164 and VOS 603 because of strong emission line spectra (see Sect. 4).'VOS' names refer to sources from the catalogue of 'Vioque, Oudmaijer, Schreiner, et al. 2020'.

Comments on the red spectral range
At the INT and NTT telescopes additional observations at a higher resolution were performed for each source, covering a redder spectral range (∼ 5800 − 7000 Å).The exceptions to this are VOS 50 and VOS 4463 (i.e.104 sources in total, see Table 1).This red range covers the important diagnostic Hα line, which enables the determination of accretion rates.An example of the normalized spectra obtained in this region around the Hα line is shown in Fig. 1.
The measured Hα equivalent widths (EW obs , observed above the continuum) are tabulated in Table B3.The Hα line profile was classified into single-peaked, doublepeaked, or P-Cygni profile (regular or inverse), following the classification scheme of Vioque et al. (2018).In addition, in Table B3 we state whether the Hβ line covered in the blue range (Sect.2.1) appears in emission.

STELLAR PARAMETERS
In this section we use the Gaia Early Data Release 3 (EDR3, Gaia Collaboration et al. 2016, 2021) to derive the stellar luminosity of the observed sources and place them in the Hertzsprung-Russell (HR) diagram.The HR diagram in combination with theoretical tracks provide us with estimations of the sources' stellar mass and age.

Data acquisition and calibration
We obtained EDR3 source identifications by using the DR2 source identifications provided in Vioque et al. (2020) and the gaiaedr3.dr2neighbourhood table of the Gaia Archive (see Torra et al. 2021).
The EDR3 parallax (Lindegren et al. 2021a) to distance conversion was done as follows.For the 140 sources with σ( )/ ≤ 0.1 we obtained the distance by simple inversion of the parallax.For the five sources with 0.1 < σ( )/ ≤ 1 we used the geometric prior of Bailer-Jones et al. (2021).In all cases the correction to the zero point parallax bias of Lindegren et al. (2021b) was applied.To trace problematic parallaxes that may lead to spurious or heavily inaccurate distances, we used the Fidelity parameter of Rybizki et al. (2022).Only three observed sources have fidelities below 90% (VOS 1385, VOS 1440, and VOS 2158).
The Gaia photometry of the sources (described in Riello et al. 2021) is presented in Table B1.To obtain the Gaia EDR3 intrinsic G BP -G RP colors of the observed sources, we converted the intrinsic Johnson-Cousins colors for dwarf stars of Pecaut & Mamajek (2013) to Gaia G BP -G RP colors with the polynomial equation presented in Table C.2 of Riello et al. (2021).In order to evaluate the validity of this approach, we compared the intrinsic colors derived this way with those obtained from GALAH+ spectra (Casagrande et al. 2021, see their Figure 3), from RAVE DR6 data (Steinmetz et al. 2020) and from the Tycho-2 Spectral Type Catalogue (Wright et al. 2003).The correspondence in all cases is within the 0.03 mag error.
Extinctions were obtained using the effective temperatures (derived in Sect.2.1), the intrinsic colors, and the color dependent extinction coefficients of Casagrande et al. (2021, see their Figure 1).In all cases R V = 3.1 was assumed.The median error for the derived A V values is 0.17 mag.
Observed sources By fitting an atmosphere model from Castelli & Kurucz (2003) of the corresponding T eff to the dereddened Gaia G RP photometry we derive the total stellar flux for each source.Combining this flux with the distance we obtain the total luminosity (L, in a procedure similar to that of Vioque et al. 2018).We have assumed solar metallicity and log(g) = 4.0.The effect of these parameters in the derived luminosities is negligible.
The luminosities of all sources are presented in Table B2, together with their distances and A V extinctions.

SEDs and IR excess
Fitting an atmosphere model to the dereddened photometry (Sect.3.1) allows us to generate Spectral Energy Distributions (SED) to estimate the amount of infrared (IR) excess.For doing this we use the 2MASS and WISE passbands (from AllWISE, Cutri et al. 2013) which are available for all sources (see Vioque et al. 2020).These passbands range from J band (1.24 µm) to W 4 band (22 µm).The derived IR excesses (L IR /L * ) appear tabulated in Table B3.
We note that WISE bands, especially W3 (12 µm) and W4, have a Point Spread Function that might lead to contaminated photometry in crowded regions or with bright backgrounds (see e.g.Koenig & Leisawitz 2014, Ribas et al. 2014).In general, the IR excesses derived this way should be considered as indicative, as a fraction of them might be affected by this caveat.A warning flag was included in Vioque et al. (2020) catalogue indicating potential problematic W3 and W4 passbands.In Table B3 this warning flag has been refined by examining the images at W1, W2, W3, and W4 as a set.In total, we found 29 sources where contamination is suspected.
There are six sources (VOS 76, VOS 78, VOS 491, VOS 1240, VOS 1385, and VOS 2161) with an IR excess over 70% of the stellar bolometric luminosity (L IR /L * > 0.7), which is the maximum excess typically observed in Herbig stars (e.g.Pascual et al. 2016;Banzatti et al. 2018).This implies that the IR excess luminosity of those sources likely includes a large amount of contamination from extended background emission.

Hertzsprung-Russell diagram
We present the 143 observed sources with T eff determinations in an HR diagram in Fig. 2. From this HR diagram we derive masses and ages by using the PAR-SEC 1.2S pre-main sequence tracks (Bressan et al. 2012;Marigo et al. 2017).These masses and ages are tabulated in Table B2.We note that the ages are very model dependent, are based on an arbitrary decision of the age 'zero', and are very susceptible to the HR diagram location uncertainties.
There are nine sources which are inconsistent with masses M>1.5 M .Therefore, they rather belong to the T Tauri regime.These are highlighted in Table B2.The remaining 134 sources are compatible with the Herbig regime (M>1.5 M ).

CONTAMINANTS
The spectral analysis of Sect. 2 and the derivations of Sect. 3 allow us to identify some evolved contaminants among the observed PMS candidates.VOS 1634 and VOS 1806 have spectral types corresponding to an M and a K star, respectively.Their Gaia EDR3 parallaxes assign them luminosities over 400 L .Hence, these are probably post-MS giants rather than PMS sources, as YSOs this massive would not be optically visible at such early stages of evolution.VOS 2164 spectra clearly correspond to a planetary nebula (and it appears as so in previous literature; e.g.Kaler et al. 1976).In addition, previous literature allowed us to identify VOS 603 as a cataclysmic variable star (dwarf nova; e.g.Otulakowska-Hypka et al. 2016), VOS 1240 as a carbon star (e.g.Groenewegen et al. 2002), VOS 1380 as a Type II Cepheid (e.g.Schmidt et al. 2004), and VOS 1385 as a RV Tau variable (Braga et al. 2018).VOS 458 was identified as an AGB candidate in Robitaille et al. (2008).However, it is a bit hot and under-luminous for an AGB star.We included VOS 458 in the list of contaminants because of its uncertain nature, although we note that the stellar parameters derived for VOS 458 are not incompatible with a YSO nature.
In total, we have identified eight evolved stars within the 145 observed PMS candidates.These appear marked in the HR diagram of Fig. 2 and are removed from the rest of the analysis.Hence, the number of observed sources we consider in what follows is 137 (145 − 8).
We should point out that separating Herbigs from classical Be stars (very similar emission-line non-PMS sources, Rivinius et al. 2013) was the main task of the machine learning algorithm used in Vioque et al. (2020).Hence, the sample observed in this work is already filtered of classical Be stars.

MASS ACCRETION RATES
For the sources with Hα and Hβ lines in emission we corrected the measured EWs for the underlying line absorption.To do this we used the typical EW absorption values of each spectral sub-type (Joner & Hintz 2015).These corrected equivalent widths (EW cor ) are tabulated in Table B3, together with the observed ones (EW obs ).Following, for example, Fairlamb et al. (2017): where F λ is the continuum flux density corresponding to the central wavelength of the Hα or Hβ line.We obtained this flux by using atmosphere models from Castelli & Kurucz (2003), which were scaled to the dereddened G RP flux of each star.Then, L Hα,β = 4πd 2 • F Hα,β , where d is the distance to the sources.Following Fairlamb et al. (2017, see also Mendigutía et al. 2011, Wichittanakom et al. 2020, and references therein) we can derive the accretion luminosity (L acc ) from the lines as: where A and B are constants.For Herbig stars, Fairlamb et al. ( 2017) determined these constants to be A = 2.09 ± 0.06 and B = 1.00 ± 0.05 for the Hα line and A = 2.60 ± 0.09 and B = 1.24 ± 0.07 for the Hβ line.
Finally, the mass accretion rate ( Ṁ acc ) can be derived as: using the stellar parameters derived in Sect.3. In case information about both the Hα and Hβ line is available the derived accretion rates come from the Hα line EW cor (this was done for consistency, as stars with Hα emission often have Hβ in absorption).We derive accretion rates for 92 sources using the Hα line and for 12 other sources using the Hβ line.The mass accretion rates derived this way are tabulated in Table B3.We note that Eqs. 1 and 2 assume a magnetospheric accretion mechanism (Hartmann et al. 2016 and references therein).

PMS NATURE OF THE OBSERVED SOURCES
In this section, we assess the PMS nature of the observed sources by means of their location in the HR diagram, IR excesses, line profiles, and accretion rates.We note that it is beyond the capabilities of the data presented in this paper to assert with absolute certainty whether all the observed sources are indeed new Herbig discoveries.In fact, it even proved difficult and controversial for much more intensely studied objects (e.g.HD 45677, Oudmaijer & Miroshnichenko 2017).Nevertheless, in this section we provide ample and independent evidence to conclude that the vast majority, if not all, of the 137 considered sources (145 − 8, Sect.4) are of a Herbig nature (see Appendix A for a description on the most dubious sources).

HR diagram, stellar masses, and IR excesses
In the HR diagram of Fig. 3 it can be seen that most of the 137 observed sources are massive hot objects.In addition, none of the observed sources are located outside PMS locations in the HR diagram, although the proximity of some of them to the Zero Age Main Sequence (ZAMS) hampers a clear PMS identification.In the HR diagram of Fig. 3 we also show previously known Herbig Ae/Be stars and IMTTs (from the compilations of Vioque et al. 2018, Guzmán-Díaz et al. 2021, and Valegård et al. 2021) with good astrometric solutions [σ( )/ ≤ 1, RUWE < 2, and Fidelity > 0.95].The stellar parameters of these previously known sources were rederived following Sect.3 procedures to compare HR diagrams which are affected by the exact same systematics and uncertainties.We conclude that the observed sources are similarly distributed in the HR diagram to the previously known Herbigs.
In Fig. 4 we present the mass distribution of the 137 sources.In this figure we also show the mass distribution of the previously known Herbig Ae/Be stars and IMTTs.We note that the observed sources and the previously known Herbigs cover a similar mass range and have an analogous mass distribution.This is likely due to the fact that both groups are tracing the massive end of the IMF (see Vioque et al. 2018;Guzmán-Díaz et al. 2021).It is noteworthy that only three sources with M>15 M were observed, whereas there are several previously known Herbig stars over that mass.This is because PMS sources with those masses, in addition New and previously known Herbig stars to being rare and often optically obscured, are typically at large distances and thus tend to have poor parallaxes.This caused those sources to be systematically excluded from the target selection of Sect.2, which sampled the catalogue of Herbig candidates of Vioque et al. (2020).
In Fig. 5 we plot the IR excesses derived in Sect.3.2 for the observed sources as a function of mass.In this figure we also show the values obtained for the previously known Herbigs in Vioque et al. (2018).The observed sources have similar IR excesses to the previously known Herbigs, and the break in inner disk dispersion efficacy at 7 M discussed in Vioque et al. (2018) is also present for the observed objects.

Emission lines and accretion rates
Regarding the presence of emission lines, 92 of the 100 sources with Hα line observations show the line in emission.In addition, as discussed in Appendix A, of those eight sources without Hα emission, seven show other emission lines.Of the 37 sources for which no Hα line information is available 12 have emission in Hβ.We note that Hβ emission may not be present even if Hα emission is.This high percentage of sources with hydrogen emission supports the PMS nature of the group.Regarding the Hα line profiles, of the 92 stars with Hα emission, 37 show single-peaked emission, 44 doublepeaked emission, and 11 P-Cygni emission (of which five have inverse P-Cygni profiles).Therefore, 40% are single-peaked, 48% are double-peaked and 12% are P-Cygni.These percentages are similar to those observed for known Herbigs (31%, 52%, and 17%, respectively; see Vioque et al. 2018).We suspect that the small difference between both groups is caused by the lower resolution of our observations, which moved many P-Cygni and double-peak profiles to the 'single-peaked' group.We point out that similar percentages were found for the Brγ line by Grant et al. (2022) in known Herbig stars.
The high fraction of double-peaked line profiles is also suggestive of the PMS nature of the observed sources; stellar activity would not necessarily result in doublepeaked emission profiles and stellar winds would have resulted in a larger fraction of P-Cygni like shapes.The fraction of single-peaked Hα profiles is consistent with the fraction of pole-on accretion disks that would be expected in a random distribution.
The accretion rates of this group of 92 sources with Hα emission and 12 sources with Hβ emission are presented in Fig. 6 as a function of mass.In Fig. 6 we compare these accretion rates with those obtained via the same procedure and assumptions in Wichittanakom et al. (2020) for a set of 163 previously known Herbig Ae/Be stars.The overlap between both sets in this parameter space is consistent with the observed sources being of a Herbig nature.This further supports that the Hα and Hβ emission used to derived the accretion rates is originated in a PMS accretion disk.There are, however, a few outliers in the trend of Fig. 6.The accretion rates of the observed sources are analyzed in more detail in Sect.7.
In this section we have shown evidence for the Herbig nature of most of the observed sources.This can be summarised by their PMS location in the HR diagram (Fig. 2), the amount of IR excess and its correlation with stellar mass (Fig. 5), the presence of emission lines in their spectra, the shapes of the Hα line, and the derived mass accretion rates (Fig. 6).Therefore, of the 145 observed objects, we propose 128 sources to be new Herbig identifications, nine sources to be new massive T Tauri identifications (Sect.3.3), and eight sources to be evolved stars of non-PMS nature (see Sect. 4).A closer look to the 128 proposed Herbig stars allowed us to determine 20 less secure identifications.These 20 stars are discussed in Appendix A.

ACCRETION PROPERTIES
The mass accretion rates derived in Sect. 5 from Hα and Hβ luminosities, together with the similar results of Wichittanakom et al. (2020), allow us to construct the largest sample of Herbig stars with mass accretion rate determinations to date.In total, we compile 258 Herbig sources with derived accretion rates between both works.We note that 48 of these sources have accre- tion rates that were derived in Fairlamb et al. (2015) by measuring the UV excess over the Balmer jump.This is a more direct measurement of the accretion rate, and it is free of the assumptions made when correlating accretion luminosity and emission line luminosities (see Mendigutía et al. 2015).
With this enhanced sample, we revisit the accretion rate properties of Herbig stars.In Fig. 6 we show the accretion rate of this combined sample as a function of stellar mass.We note that for both new and previously known sources the mass accretion rate increases with stellar mass.As Wichittanakom et al. (2020) discussed, we also find that the accretion rate decreases with time during the PMS phase.However, it is not trivial to disentangle this effect from the dependence of the age on the mass.
In Wichittanakom et al. (2020) it was concluded that lower mass Herbigs have a dependence of the mass on the accretion rate, characterised by a gradient that matches the gradient observed in the T Tauri regime (e.g.Calvet et al. 2004;Natta et al. 2006).However, higher mass Herbigs show a smaller gradient in the mass vs. accretion rate relation.The break between both groups was set at 3.98 +1.37 −0.94 M .This accretion break and similar accretion gradients to those found in Wichittanakom et al. (2020) have also been identified by Grant et al. (2022), using a similar sample and the Brγ line as the accretion tracer.
In a similar way to what was done in Section 5 of Wichittanakom et al. ( 2020), we looked for the mass value where the difference between the gradients of the lowand high-mass regimes maximizes.We note that this approach does not take into account the uncertainties in the accretion rate, nor the caveats of estimating the accretion rate from emission lines (see e.g.Fairlamb et al. 2017;Mendigutía 2020).To study this difference between the observed gradients we use the S parameter.This parameter represents the significance of the difference of the slopes.It is defined as: where b is the slope in each regime of the mass vs. accretion rate linear fit in log space, and σ 2 is the variance in that slope.The results of this study on the whole sample of 258 sources are illustrated in the top panels of Fig. 7.In the top left panel we show the S parameter as a function of the mass used to separate the low-from the high-mass regime.We calculated the S parameter for the sample of new Herbigs, the sample of previously known Herbigs, and the combination of those two samples.In all cases the S parameter increases with mass, up to a maximum at around 3-4 M , and then monotonically decreases.Using a simple t-test approach we conclude that the difference between gradients is significant to within 95% confidence if S 2. Hence, it is evident from Fig. 7 that, within the Herbig regime, there is a change in the gradient of the accretion rate as a function of mass (see, however, Sect.8.2 for a description of the caveats to consider when using line luminosities as accretion tracers).
The maximum S value found for all sources corresponds to a stellar mass for the break in gradient of 3.26 M (top right panel of Fig. 7).This number is consistent with the value found in Wichittanakom et al. (2020) of 3.98 +1.37 −0.94 M .However, in the top right panel of Fig. 7 it is noticeable that several sources are far from the general correlation.If we remove those 'outliers' from the analysis (reducing the sample to 198 sources, bottom panels of Fig. 7) we obtain a maximum S value at 3.87 M , much closer to the central value of 3.98 derived in Wichittanakom et al. (2020).
Although there is a range of masses in which the gradient difference as traced by the S parameter is significant; there is a clear peak for S values in all samples at around 3-4 M .This peak stands after removing the sources that deviate the most from the correlation.Thus, we confirm that the break in accretion rate detected in Wichittanakom et al. (2020) in the 3-4 M range holds with the sample of new Herbigs.By using the maximum S values obtained for the different Herbig subsamples, we further constrain the potential break in accretion rate to a mass of 3.87 +0.38  −0.96 M (corresponding to the mass of a B7-B8 MS star).

DISCUSSION
In the previous sections we present and discuss 128 new Herbig stars homogeneously selected and observed, for which we provide accurate stellar parameters.In this section we put these sources in context with the historically considered Herbigs, and explain why these new Herbigs provide interesting insights to the intermediateto high-mass star formation scenario.

General remarks
Among the main caveats of previous studies dedicated to intermediate-to high-mass PMS stars is the low number of sources in any given mass or age range (e.g.Vioque et al. 2018;Guzmán-Díaz et al. 2021;van der Marel & Mulders 2021).The sample of 128 new Herbig stars contains both low-mass objects at the boundary with the T Tauri regime and very-massive PMS objects (see Fig. 3 and 4).In the mass range of 1.5 to 4 M we present 73 new sources increasing by 42% the number of known objects.In the mass range of 4 to 8 M we present 23 new sources, increasing by 70% the number of known objects (55 sources are above the 4 M threshold of Wichittanakom et al. 2020 for the break in accretion properties).In the mass range of PMS stars above 7-8 M (typically considered the MYSO regime) we present 32 new sources, increasing by 80% the number of known objects.This is the threshold of Vioque et al. (2018) from which very effective inner disk-dispersal mechanisms are acting.These statistics, summarised in Table 2, 2021) with mass determinations.However, other less famous objects of the class do exist in the literature, like the recently proposed 58 Herbig Ae/Be sources from the LAMOST survey (Shridharan et al. 2021;Zhang et al. 2022), the 13 proposed Herbig Ae/Be stars in the Small Magellanic Cloud (Keller et al. 2019), or the 77 IMTTs found in the Carina Nebula (Nuñez et al. 2021).
The most massive of the newly discovered Herbig stars appear to overlap with the class of Massive Young Stellar Objects (Frost et al. 2021;Koumpia et al. 2021).More than 300 such objects with luminosities larger than 5000 L , corresponding to masses larger than around 8 − 10 M , are listed as MYSO in the RMS catalogue (Lumsden et al. 2013).These massive young stars are infrared bright due to the large amounts of dusty material obscuring them from sight and, as a result, they were historically assumed to be optically invisible.Intriguingly however, some were already reported to be visible in the optical.For example, the objects PDS 27 and PDS 37 appear in the RMS catalogue (see also Koumpia et al. 2019), but with V band magnitudes of ∼13 mag they are optically bright enough to have been recognised as Herbig Be stars (Ababakr et al. 2015;Vioque et al. 2018).Given that Gaia, with its faint magnitude limit of ∼20 mag, pushes the definition of optically visible to much fainter magnitudes, the new Gaia-discovered Herbig stars may well bridge the gap between the (optically bright) Herbig Be stars and the optically faint MYSOs.The fact that this project already adds 32 new sources to this mass regime evidences this hypothesis, while pre- Regarding the age of the sources, the sample of 128 new Herbig stars contains both sources close to the ZAMS and at earlier stages of evolution.Although some stars could be considered IMTTs (see e.g.Valegård et al. 2021), most sources clearly belong to the Herbig Ae/Be regime (see Sect. 1 and footnote 1).With the exception of the four sources discussed in Appendix A, the age range covered at every mass bin is equivalent to the age range covered by previously known and analyzed objects.Given that the Vioque et al. ( 2020) catalogue is HR-diagram independent, this is probably not caused by a selection bias but by the fact that younger objects are too embedded to appear in the Gaia survey.
There are some biases in the catalogue of Vioque et al. (2020) that affect the sample of 128 new Herbig stars.These are mainly that sources with strong IR excesses and Hα emission were favored in the selection.This biases the sample toward the more 'active' PMS objects.We refer the reader to Vioque et al. (2020) for more details.
In Fig. 8 we plot the distances and galactic coordinates of the 128 new Herbig detections.To these we add all previously known Herbigs (Table 2) with Gaia EDR3 parallaxes and σ( )/ < 1.The new Herbigs are confined to the Galactic plane because of selection effects in Vioque et al. (2020).They are generally fainter than the previously known Herbigs (see Table B1), and hence they are typically further away.There are only 5 new Herbigs within 500 pc, 80% of the new sources being beyond 1 kpc.

Interpretation of the change in accretion gradient
In this section we interpret the break in accretion gradient at 3 − 4 M identified in Wichittanakom et al. (2020) and Grant et al. (2022), that we extend to the new Herbig stars in Sect.7 and constrain to 3.87 +0.38   −0.96   M .The most accepted interpretation for this change in accretion gradient is that low-mass objects are subjected to magnetospheric accretion (see Bouvier et al. 2007), whereas a fraction of the more massive objects are accreting through a different mechanism; possibly the boundary layer accretion mechanism (Mendigutía 2020).
We should caution, however, that the accretion rates derived in Sect. 5 assume a magnetospheric accretion scenario.If a different accretion mechanism applies for some objects their derived accretion rates might be highly inaccurate (see Section 6 of Wichittanakom et al. 2020).For example, if a boundary layer accretion mechanism (Mendigutía 2020) is applying for sources more massive than 3 − 4 M , the accretion luminosity to accretion rate relation of Eq. 2 would need to be corrected by the relative difference in rotational velocities between the star and the gas contact phase (see Sect. 1.3 of Wisniewski et al. 2021 for more details).In addition, the relation between the line luminosities and the accretion luminosity (Eq. 1) would have to be revisited for the boundary layer scenario.Therefore, we advise the reader to treat the accretion rates derived in this work for massive stars with caution.
In addition, there are alternative explanations to the change of trend in hydrogen line luminosity properties, apart from a change in accretion mode (see Mendigutía et al. 2015;Marcos-Arenal et al. 2021).For example, outflowing material could be dominating the line emission (caused by e.g.disk photoevaporation, Guzmán-Díaz et al. 2021).Indeed, it is not straightforward to characterize the origin of the hydrogen emission lines in high-mass Herbig stars (Mendigutía et al. 2017;Mendigutía 2020, and references therein).This contrasts to what we observe for the lower mass objects, where the hydrogen emission lines seem to originate in the magnetosphere (Bouvier et al. 2020;Gravity Collaboration et al. 2020).
Furthermore, there is an observational bias to consider.High-mass PMS stars are optically visible for a much shorter time than low-mass PMS stars, and they are often heavily obscured, especially at the younger ages.Hence, in our optical analysis we might be biased against the strongest accretors in the high-mass regime.
Because of the reasons stated in this section, it is not straightforward to deduce a break in accretion properties only from the break in the behaviour of the hydrogen lines.However, this break at 3 − 4 M has been independently detected using different techniques (e.g.near-IR interferometry: Monnier et al. 2005; optical-  et al. 2015).It is the combination of those results together with the accretion rates derived from emission lines that lead us to conclude that the break in accretion properties of Fig. 7 is likely due to a change between accretion mechanisms; from magnetospheric accretion applying to the lower-mass objects to a boundary layerlike accretion mechanism acting in some, or most, of the more massive objects.

Evaluation of Vioque et al. (2020) results
The target list of the observations presented in this work was extracted from the catalogue of new intermediate-to high-mass PMS candidates of Vioque et al. (2020).In this section we reevaluate the accuracy and quality of that catalogue.
In the Vioque et al. ( 2020) catalogue 2226 new Herbig candidates were presented (with σ( )/ ≤ 0.2, this number remains similar when the astrometry is updated to EDR3).We note that, although the catalogue is astrometry-independent, arbitrary cuts to the astrometric quality are necessary to select massive objects with a certain degree of confidence.In this work we have observed 145 objects, which is roughly 6% of the whole catalogue.As mentioned in Sect.2, the target selection was based only on the absolute magnitude and on the parallax quality of the sources.Hence, the target list is representative of the whole catalogue of new intermediate-to high-mass PMS stars of Vioque et al. (2020).
The number of contaminants we have found in this work (8/145 or 5.5%) is consistent with the estimated lower-limit precision of the Vioque et al. (2020) catalogue (P 81%).This affirmation holds true even when we consider as contaminants the more dubious sources of Appendix A (28/145 or 19%).We note that the Vioque et al. (2020) catalogue is HR diagram independent and thus the high proportion of massive stars targeted (134/145 sources, 92%, are above 1.5 M , see Sect.3.3) is a compelling positive assessment of that catalogue.Therefore, we conclude that the observations presented in this work give independent support to the quality and robustness of the Vioque et al. ( 2020) catalogue.
In addition, 14 classical Be candidates from the Vioque et al. (2020) catalogue were observed (these are non-PMS stars which are typical contaminants in Herbig samples, see Vioque 2020).None of these classical Be candidates could be identified as a misclassification (e.g. by showing a PMS nature).The discussion on the observations of these sources will be presented in an independent paper.
Because of sensitivity limitations, mostly candidates at the bright end of the Vioque et al. (2020) catalogue were observed (90% of the observed sources are in the 12<G<14 mag range).One could therefore argue that the observed sample is biased.However, the Vioque et al. (2020) catalogue is distance independent.Hence, the conclusions that arise from these observations can be extrapolated to the fainter objects of the catalogue, given that they were all selected homogeneously by the machine learning algorithm.

CONCLUSIONS
In this work we discuss the results of the spectroscopic observations of a sample of 145 Herbig candidates from the catalogue of Vioque et al. (2020).The main results and conclusions of these observations are the following: • We propose 128 sources as new 'Herbig' identifications (i.e.PMS stars with M > 1.5 M ).We provide ample evidence supporting this classification.This evidence is based on their PMS location in the HR diagram, the amount of IR excess and its correlation with stellar mass, the presence of emission lines in their spectra, the shapes of the Hα line, and the derived mass accretion rates.Only 5 sources lie within 500 pc, whereas 75% of the stars are between 1 and 4 kpc.Twenty sources were flagged as less secure identifications.
• We derive extinctions and accurate stellar parameters for all sources, placing them in the HR diagram by means of Gaia EDR3 parallaxes.
• This sample of 128 new Herbig stars increases the number of known objects of the class by ∼ 50%.The sources are distributed over a representative range in mass and age when compared to previously known Herbig stars.According to classical definitions, most of the observed sources fall within the Herbig Ae/Be or the Massive Young Stellar Objects regime, but some stars can be considered Intermediate-Mass T Tauris.In particular, 23 of the new sources have masses between 4-8 M and 32 sources have masses > 8 M .
• Four sources were identified as new 'unclassified B[e]' (FS CMa) discoveries.Nine other sources were also identified as having a PMS nature, but their masses assign them to the T Tauri regime.
• We derive accretion rates for 104 of the new Herbig stars by using hydrogen emission lines (Hα and Hβ luminosities).This is a 60% increment to the number of Herbig stars with derived accretion rates.
The change in accretion gradient as a function of mass in the 3-4 M range described in Wichittanakom et al. ( 2020) is also present for the new Herbig stars.This provides further support to a change in accretion mechanism happening withing the Herbig regime.We constrain the mass for this possible change to 3.87 +0.38 −0.96 M (the mass of a B7-B8 main sequence star).
• There are four sources (VOS 63, VOS 67, VOS 821, and VOS 1635) of 5-6 M which are younger than previously known PMS stars of that mass.
• The sudden decrease in the amount of near-and mid-IR excess at ∼ 7 M described in Vioque et al. (2018) for the historically considered Herbig stars is also present for the new sources.For this group of Herbig stars M > 7 M corresponds to sources with T eff 15000 K.This further supports the idea of very effective inner-disk dispersion mechanisms acting on massive stars, like the disk photoevaporation mechanism proposed in Guzmán-Díaz et al. (2021).
• The observations described in this work provide independent support to the accuracy and highquality of the catalogue of new intermediate-to high-mass PMS stars presented in Vioque et al. (2020), of which these observations constitute a mere 6%.
These observations yield a well-defined set of new intermediate-to high-mass PMS stars.Contrary to previous samples, these new Herbig stars were homogeneously identified and observed.Therefore, this set of objects will be the basis for future surveys and followup observations dedicated to the Herbig group and their protoplanetary disks.The sample of new Herbig stars presented in this work will be complemented by an X-Shooter Very Large Telescope survey focusing on newly identified intermediate-mass T Tauri stars (Iglesias et al. in prep.).All together, the two surveys will increase by a factor of two the number of known intermediateto high-mass PMS objects, covering representatively all the stages of the optical evolution of massive forming stars.SIMBAD database developed and operated at CDS, Strasbourg, France.This article is based on observations made in the Observatorios de Canarias del IAC with the Isaac Newton Telescope (INT) operated on the island of La Palma by the Isaac Newton Group of Telescopes in the Observatorio del Roque de los Muchachos.In addition, this article is based on observations collected at the Centro Astronómico Hispano-Alemán (CAHA) at Calar Alto, operated jointly by Junta de Andalucía and Consejo Superior de Investigaciones Científicas (IAA-CSIC).Finally, this article also used observations collected at the European Organisation for Astronomical Research in the Southern Hemisphere under ESO programmes 0104.C-0937(A) and 0104.C-0937(B).

APPENDIX
A. LESS SECURE IDENTIFICATIONS Among the 128 new Herbig stars proposed there are 20 sources which, for different reasons, have a less secure PMS nature than the rest.These sources are discussed in detail in this appendix.
Eight sources have Hα and Hβ fully in absorption (VOS 209, VOS 448, VOS 491, VOS 495, VOS 668, VOS 854, VOS 1922, andVOS 2060).However, three of these sources (VOS 448, VOS 495, and VOS 1922) show an asymmetric Hα absorption profile, which might hint some Hα emission.These eight sources are shown in the HR diagram of Fig. 2.Although emission in hydrogen lines is historically one of the defining properties of Herbig stars, some intermediate-mass PMS stars lack hydrogen emission.Hence, this fact alone is inconclusive for removing these sources from the PMS category.In addition, of the eight sources for which we do not detect clear hydrogen emission, seven show other emission lines (the exception is VOS 854).However; VOS 448, VOS 495, VOS 1922, andVOS 2060 do not show any significant level of IR excess, and the IR excess of VOS 491 is clearly spurious (see Sect. 3.2).Because hydrogen emission lines and IR excess are the main indicators of YSO nature, we label these latter five sources plus VOS 854 as less secure identifications.We note that 37 observed sources lack Hα information, and hence the aforementioned analysis could not be applied to them.Of these, seven sources have Hβ in absorption and display little IR excess (L IR /L * <0.03).These are VOS 821, VOS 879, VOS 1225, VOS 1276, VOS 1771, VOS 1913, and VOS 2051, which are also labeled as less secure identifications.
We now consider the nature of VOS 209 and VOS 668, which have emission lines and IR excess but lack hydrogen emission.van den Ancker et al. (2021) found that the intermediate-mass YSO HD 152384 has all hydrogen lines strictly in absorption, but has refractory lines in emission.This led van den Ancker et al. ( 2021) to suggest that HD 152384 is at the late stages of the PMS phase and is surrounded by a tenuous circumstellar disk caused by the collision of rocky planets (see also the extreme debris disks described in Moór et al. 2021).We pose VOS 209 and VOS 668 as PMS sources of a similar nature to HD 152384.Indeed, the age estimates derived in Sect.3.3 imply that both VOS 209 and VOS 668 are compatible with being close to the main sequence (7.27 +0.10 −0.19 and 6.61 +0.16 −0.09Myr, respectively).These two sources are marked in Fig. 2.
In addition, we found four sources (VOS 1405, VOS 1440, VOS 2158, and VOS 2161) that have so many permitted and forbidden emission lines that the underlying photospheric absorption spectrum is hardly visible.They are reminiscent of the 'unclassified B[e]' objects (Lamers et al. 1998), which constitute a class that includes both evolved stars and young Herbig stars.However, it is often hard to decide on the evolutionary nature of the objects.For example, the archetypal unclassified B[e] star HD 45677 may or may not be a young star (see Oudmaijer & Miroshnichenko 2017;Hofmann et al. 2022).Such sources are also referred to as FS CMa objects (Miroshnichenko 2007).We therefore propose VOS 1405, VOS 1440, VOS 2158, and VOS 2161 to be new 'unclassified B[e]' (FS CMa) discoveries.The absence of clear photospheric hydrogen absorption lines combined with multiple emission lines led to this FS CMa classification.The nature of these sources is unclear, but they could still be of a YSO nature.These new 'unclassified B[e]' discoveries are highlighted in the HR diagram of Fig. 2.
Finally, there are four sources (VOS 63, VOS 67, VOS 821, and VOS 1635), which are younger than all previously known Herbig stars of a similar mass (5 -6 M , see Fig. 3).PMS objects of this high-mass and young age are expected to be quite embedded.Thus, these sources require a closer look at their nature.VOS 821 was already mentioned in this appendix regarding its lack of observed hydrogen emission and IR excess.However, we do not have any reason to suspect of the PMS nature of the other three objects, although a post-MS nature can neither be entirely discarded.It has been proposed that the FU Ori outbursting phenomena, which can cause T Tauri stars to get bluer and more luminous (e.g.Vorobyov et al. 2017;Kuffmeier et al. 2018), might explain the position of PMS sources in this region of the HR diagram.However, should that be the case we would expect to measure large Hα EWs signposting high accretion rates, and that is not the case for these sources.
The 20 sources discussed in this appendix are flagged in Tables B2 and B3.

B. TABLES
This appendix contains Tables B1, B2, and B3. or spectral types derived from spectra (Sect.2.1) appear in boldface.T eff values or spectral types not in boldface were approximated from the spectral derivations using the conversions of Pecaut & Mamajek (2013).See Sect. 3 for a description on the derivation of the stellar parameters, distances and extinctions.Stellar radii were derived following Eq. 2.

Figure 1 .
Figure 1.Example normalized spectra of VOS 140 (spectral type B9.5), observed at INT. Top: Blue spectral range covered by our observations.Bottom: Portion of the red spectral range covered by our observations centered around the Hα line.

Figure 4 .
Figure 4. Histogram of the number of stars observed in this work per 1 M bin.In blue the 110 observed stars with stellar mass determinations larger than three times the uncertainty [M > 3σ(M )].In orange the 27 stars with M < 3σ(M ).Contours trace the previously known Herbig Ae/Be stars and IMTTs with good astrometric quality.

Figure 5 .
Figure 5. Mass vs. IR excess (LIR/L * ).In blue the stars observed in this work with M > 3σ(M ).Orange crosses show the IR excesses derived in Vioque et al. (2018) for the previously known Herbigs with a good astrometric solution and M > 3σ(M ).The gray line traces the 7 M break in inner disk dispersion efficacy.

Figure 6 .
Figure 6.Mass vs. mass accretion rate.Blue filled circles indicate the 92 observed stars with accretion rates determined from Hα measurements.Orange crosses indicate the 12 observed stars with accretion rates determined from Hβ measurements.Black plus symbols trace the 163 previously known Herbig stars with mass accretion rates derived in Wichittanakom et al. (2020).

Figure 7 .
Figure 7. Left plots: Stellar mass used to separate the low-mass regime from the high-mass regime as a function of the S parameter.Red line traces the maximum S value when all sources are considered, with its estimated uncertainty (gray shaded region).Right plots: Mass vs. mass accretion rate for the low-and high-mass regimes, which were defined by the maximum S value obtained in the plots on the left.The fits that gave the gradients with the maximum S value are shown, with a 95% confidence interval.

Figure 8 .
Figure 8.The 128 newly confirmed Herbig stars are presented as blue dots.Orange crosses trace all 232 previously known Herbigs in the literature with Gaia EDR3 parallaxes ( ) and σ( )/ < 1. Panels on the right are a close up look at the panels on the left.Top left: Galactic longitude vs. distance, each circular grid line is 1 kpc, up to 5.5 kpc.Top right: Galactic longitude vs. distance, each circular grid line is 200 pc, up to 1000 pc.Bottom left: Galactic longitude vs. Galactic latitude.Bottom right: Galactic longitude vs. Galactic latitude limited to −5.5 < b < 5.5 deg, where the newly identified stars are confined.andnear-UV spectropolarimetry:Ababakr et al. 2017;  spectro-photometry: Mendigutía et al. 2011, Fairlamb  et al. 2015).It is the combination of those results together with the accretion rates derived from emission lines that lead us to conclude that the break in accretion properties of Fig.7is likely due to a change between accretion mechanisms; from magnetospheric accretion applying to the lower-mass objects to a boundary layerlike accretion mechanism acting in some, or most, of the more massive objects.

Table 1 .
Observing dates and instrumental setups for the three observing runs; including telescope, spectrograph, CCD detector, grating or grism used, spectral range in Å, reciprocal dispersion in Å/pixel, and spectral resolution in Å.Note-The signal-to-noise ratio of these spectra are typically in the order of 100.A total of 145 different sources were observed.

Table 2 .
Number of known Herbig stars per stellar mass range.

Table B1 .
Log of observations.Column 1 contains the name of the sources as assigned in this work.Column 2 contains the Gaia EDR3 source id of each source.Column 3 shows other literature names.Columns 4 and 5 are Gaia EDR3 right ascensions and declinations.Column 6 and 7 list the observation dates (for the blue and red settings, see Sect.2).Columns 8 and 9 present the exposure times.Column 10 states the spectrograph used in each case.Columns 11 to 13 show the Gaia magnitudes.See Table1for details on the instrumental setups used at each telescope.

Table B2 .
Stellar parameters, distances and extinctions (AV) for the 145 observed sources, ordered by name.Effective temperatures (T eff )