Discovering new B[e] supergiants and candidate Luminous Blue Variables in nearby galaxies

Mass loss is one of the key parameters that determine stellar evolution. Despite the progress we have achieved over the last decades we still cannot match the observational derived values with theoretical predictions. Even worse, there are certain phases, such as the B[e] supergiants (B[e]SGs) and the Luminous Blue Variables (LBVs), where significant mass is lost through episodic or outburst activity. This leads to various structures around them that permit dust formation, making these objects bright IR sources. The ASSESS project aims to determine the role of episodic mass in the evolution of massive stars, by examining large numbers of cool and hot objects (such as B[e]SGs/LBVs). For this, we initiated a large observing campaign to obtain spectroscopic data for $\sim$1000 IR selected sources in 27 nearby galaxies. Within this project we successfully identified 7 B[e] supergiants (one candidate) and 4 Luminous Blue Variables of which 6 and 2, respectively, are new discoveries. We used spectroscopic, photometric, and light curve information to better constrain the nature of the reported objects. We particularly note the presence of B[e]SGs at metallicity environments as low as 0.14 Z$_{\odot}$.


Introduction
How exactly single massive stars, born as O/B-type main-sequence stars, progress to more evolved phases and eventually die remains an open question. Binarity, which has an important implication in the evolution, even further complicates the quest for an answer. Observational data has revealed a number of transitional phases in which massive stars can be found, also known as the massive star "zoo". Whether they pass through certain phases or not depends on their initial mass (≥ 8 M ), metallicity (Z), rotational velocity (v rot ), mass loss properties and binarity [1][2][3][4]. Although some of them are quite distinct (e.g. Wolf-Rayet stars as opposed to Red Supergiants -RSGs), there are phases which display common observables, such as B[e] supergiants (B[e]SGs) and Luminous Blue Variables (LBVs).
The B[e] phenomenon is characterized by numerous emission lines in the optical spectra [5]. In particular, there is strong Balmer emission, low excitation permitted (e.g., Fe II), and forbidden lines (of [Fe II], and [O I]), as well as strong near-or mid-IR excess due to hot circumstellar dust. However, this can be observed in sources at different evolutionary stages (such as in Herbig AeBe stars, symbiotic systems, and compact planetary nebulae, see [5] for detailed classification criteria). The B[e]SGs form a distinct subgroup based on a number of secondary criteria. They are luminous stars (log L/L 4.0), showing broad Balmer emission lines with P Cygni or double-peaked profiles. They may also display evidence of chemically processed material (e.g., 13 CO enrichment, TiO) which points to an evolved nature, although it is not yet certain if they are in pre-or post-RSG phases [6,7]. The presence of the hot circumstellar dust is due to a complex circumstellar environment (CSE) formed by two components, a stellar wind radiating from the poles and a denser equatorial ring-like structure [8][9][10][11][12]. However, the formation mechanism of this structure remains elusive. A variety of mechanisms have been proposed, such as the following: fast rotation [13], the bi-stability mechanism [14], slow-wind solutions [15], magneto-rotational instability [16] , mass transfer in binaries [17], mergers [18], non-radial pulsations or the presence of objects that clear their paths [19]. Although poorly constrained, their initial masses range from roughly 10 M to less than 40 M (Mehner 2023, IAU S361, subm.).
The LBVs are another rare subgroup of massive evolved stars, considered to represent a transitional phase from massive O-type main-sequence to Wolf-Rayet stars (e.g., [3,[20][21][22]). They experience instabilities that lead to photometric variability, typically referred to as S Dor cycles [22], as well as outbursts and episodic mass loss, similar to the giant eruption of η Carina that resulted in large amounts of mass lost through ejecta (e.g., [23]). It is not yet fully understood whether these two types of variability are related (e.g., [24]). Apart from the evident photometric variability, their spectral appearance changes significantly during their outburst activities (S Dor cycle). It is typical to experience loops from hot (spectra of O/B type) to cool states (A/F spectral types while in outbursts). Depending on the luminosity, the brightest LBVs (log L/L > 5.8) seem to directly originate from main-sequence stars (with mass > 50 M ), while the less luminous ones are possibly post-RSG objects that have lost almost half of their initial masses (within the range of ∼25-40 M ) during the RSG phase (Mehner 2023, IAU S361, subm.). Currently, various mechanisms have been suggested, such as radiation and pressure instabilities, stellar rotation, and binarity (see the reviews on the theory and observational evidence in [22,24], Mehner 2023, IAU S361, subm., and the references therein) and, as such, no comprehensive theory exists to explain them.
Therefore, if and how these two phases are linked remains an open question. B[e]SGs tend to have initial masses with a wide range below the most luminous LBVs, and in accordance with the less luminous ones. The presence of similar lines in their spectra points to similarities in their CSEs, with shells and bipolar nebulae observed in both cases [22,25,26].
Due to their photometric variability, LBVs are more commonly detected in other galaxies compared to B[e]SGs, which generally display less variability. 1 . Therefore, B[e]SGs need to be searched for to be discovered. This has only been successful for 56 (candidate) sources in the Galaxy and for the Magellanic Clouds (MCs), M31 and M33, and M81 [7], and only recently in NGC 247 [27]. On the other hand, LBVs have been found in more galaxies (additional to the aforementioned), such as IC 10, IC 1613, NGC 2366, NGC 6822, NGC 1156, DDO 68, and PHL 293B [22,[28][29][30][31], summing up to about 150 sources (including candidates).
This paper presents the discovery of new B[e]SGs and LBV candidates found with a systematic survey to identify massive, evolved, dusty sources in nearby galaxies (≤5 Mpc), as part of the ASSESS project 2 (Bonanos 2023, IAU S361, subm.). In Section 2 we provide a short summary of the observations and of our approach, in Section 3 we present the new sources, and in Sections 4 and 5 we discuss and conclude our work.

Galaxy Sample
For the ASSESS project, a list of 27 nearby galaxies (≤5 Mpc) was compiled (see Bonanos 2023, IAU S361, subm.). In this paper, we present our results from a sub-sample of these galaxies (Table 1) for which the spectral classification is final, while for another set we have scheduled observations in queue and have submitted proposals. For some galaxies (e.g., MCs) data have been collected through other catalogs/surveys and are presented separately (e.g., [32][33][34][35]). Table 1. Properties of galaxies examined in this work: galaxy ID (column 1), sky coordinates (columns 2 and 3), galaxy type (column 4), distance (column 5), metallicity (column 6), and radial velocity (RV, column 7).  [57] 227 [46] 1 The numbers presented here reflect the mean value per galaxy. 2 The RV errors correspond to the statistical error and not the systemic one, which is (typically) larger.

ID
The aim of the ASSESS project is to determine the role of episodic mass loss by detecting and analyzing dusty evolved stars that are primary candidates to exhibit episodic mass loss events (Bonanos 2023, IAU S361, subm.). This mass loss results in the formation of complex structures, such as shells and bipolar nebulae in Wolf-Rayet stars and LBVs (e.g., [25,58]), detached shells in AGBs and RSGs (e.g., [59]), disks and rings around B[e]SGs (e.g., [7,11], or even the dust-enshrouded shells within which the progenitors of Super-Luminous Supernovae lay (e.g., [60][61][62]). The presence of these dusty CSEs makes these sources bright in mid-IR imaging. Therefore, we based our catalog construction on published point-source Spitzer catalogs [63]. Since IR data alone cannot distinguish between these sources, the base catalogs were supplemented with other optical and near-IR surveys (Pan-STARRS1; [64], VISTA Hemisphere Survey-VHS; [65], Gaia DR2; [66,67]). Gaia information was also used to remove foreground sources when possible (see [68], and Tramper et al., in prep., for more details).
Given this data collection, we performed a selection process to minimize contamination by AGB stars and background IR galaxies/quasars. An absolute magnitude cut of M [3.6] ≤ −9.0 [34] and an apparent magnitude cut at m [4.5] ≤ 15.5 [69] were applied to avoid AGB stars and background galaxies, respectively. In order to select the dusty targets we considered all sources with an IR excess, defined by the color term m [3.6] − m [4.5] > 0.1 mag (to exclude the majority of foreground stars, for which this is approximately 0, and to select the most dusty IR sources). The three aforementioned criteria served as a minimum to consider a source as a priority target. Consequently, the reddest and brightest point-sources in the Spitzer catalogs were given the highest priority. An extensive priority list/system was constructed by imposing certain limits for the color term, M [3.6] , and the presence of an optical counterpart (for more details, see Tramper et al., in prep.). Depending on the galaxy size we ended up with a few tens to hundreds of targets per galaxy.
To obtain spectroscopic data for such a large number of targets we required instruments with multi-object spectroscopic modes. With these we could allocate up to a few tens of objects per pointing. Multiple pointings (with dithering and/or overlap) were applied to cover more extended galaxies and when the density of the target was high. Therefore, when we were creating the necessary multi-object masks we were forced to select sources based on the spatial limitations (e.g., located out of the field-of-view or at the sensor's gap) and spectral overlaps. Consequently, some priority targets were dropped and, additionally, non-priority targets ("fillers", i.e., sources dropped through the target selection approach described previously) were added to fill the space.

Observations
To verify the nature of our selected targets we needed spectroscopic information. Since this is not available for the majority of the ASSESS galaxies, we initiated an observation campaign to obtain low resolution spectra. Given the large number of targets, along with the sizes of the galaxies, we used the multi-object spectroscopic modes of the Optical System for Imaging and low-Intermediate-Resolution Integrated Spectroscopy (OSIRIS; [70]), on the 10.4 m GTC ( [71], for the galaxies visible from the Northern hemisphere, i.e., IC 10 and NGC 6822). We used the FOcal Reducer/low dispersion Spectrograph 2 (FORS; [72]), at 8.2 m ESO/VLT (for the Southern galaxies, i.e., the rest of Table 1). The resolving power and wavelength coverage was similar for both instruments, at ∼500-700 over the range R ∼5300-9800 Å for GTC/OSIRIS and R ∼1000 over the range ∼5200-8700 Å for VLT/FORS2. Details for the observations and data reduction can be found at Munoz-Sanchez et al., in prep., for the GTC/OSIRIS campaign and Tramper et al., in prep., for the VLT/FORS2 campaign. Here we provide only a short overview of the data reduction followed.
For the OSIRIS data we used the GTCMOS package 3 which is an IRAF 4 . This pipeline for spectroscopic data combines (for each raw exposure) the two CCD images from the detector (correcting for geometric distortions) and performs bias subtraction. Although it can perform the wavelength calibration and can correct the curvature across the spatial direction in 2D images, we noticed that it was not perfect. For this reason we opted to perform a manual approach and extracted a small cut in the image around each slit. We performed the wavelength calibration individually for each of these images (slits) and tilt was corrected when necessary. The science and sky spectra were extracted (in 1D), and followed by flux calibration. We used IRAF to extract the long-slit spectra for standard stars, and then the routine standard and sensfunc to obtain the sensitivity curve. This was applied through the calibrated routine to the science spectrum.
For the FORS2 data, we used the FORS2 pipeline v5.5.7 under the EsoReflex environment [74]. This resulted in flux-calibrated, sky-subtracted 1D spectra for each slit on the mask. However, for some slits the pipeline did not produce suitable spectra, due to multiple objects in the slit, strongly variable nebular emission, slit overlap, and/or strong vignetting at the top of the CCD. For this reason, we also performed the reduction without sky subtraction and manually selected the object and sky extraction regions from the 2D spectrum. For each slit, the automatically and manually extracted spectra were visually inspected, and the best reduction was chosen.

Spectral Classification
The resolution and wavelength range (as described in the previous section) provide access to a number of spectral features, such as Hα (a mass loss tracer for highṀ stars), the TiO bands (present in cool stars), He I and He II lines (indicative of hot stars), various metal lines (notably Fe lines), and the Ca triplet (luminosity indicator). Therefore, we were able to effectively classify the vast majority of our targets.
Both B[e]SGs and LBVs are characterized by strong emission lines, indicative of their complex CSEs. Hα is usually found in very strong emissions and is significantly broadened in the presence of strong stellar winds and/or the presence of a (detached) disk (e.g., [7,12]). There were a number of He I lines (at λλ5876.6, 6678.2, 7065.2, 7281.4) within our observed range, which manifest in the hottest sources. In the quiescence state of LBVs, the presence of He lines indicates hotter sources (of B/A spectral type, which can be observed even with P-Cygni profiles when stellar winds are strong, such as, for example, in [21]). However, when an outburst is triggered and evolves outwards, the temperature temporarily decreases until the ejecta become optically thin. As a result of this temperature shift, the spectral lines typical for the quiescent LBV weaken and metal emission lines strengthen (e.g., [75]). During this phase, and depending on the temperature and density conditions of the circumstellar material, they may also display some forbidden Fe lines. B[e]SGs display additional forbidden emission lines, due to their more complex CSEs, with typical examples being [O I] λλ5577, 6300, 6364 and [Ca II] λ7291, 7324. The latter is more evident in the more luminous sources (e.g., [10,21]).
Therefore, among all sources identified with strong Hα emissions, we classified as being B[e]SGs those with evident [O I] λ6300 [5,21], and as being LBVs those without. Both classes may display forbidden emission lines from Fe and Ca (e.g., [21,76,77]), while all of them display Fe emission lines. We note here that these LBVs are candidate sources, since there is no absolute way to characterize an LBV from a single-epoch spectrum (in contrast to B[e]SGs). It has to be supplemented with more spectroscopic or photometric observations that reveal variability (and possibly the return to a hotter state). We also note that our sample contained more interesting sources that displayed Hα in emission (i.e., main sequence O-stars and blue supergiants), but these were left for future papers (e.g., Munoz-Sanchez et al. 2023, IAU S361, submission).

Statistics
From our large observational campaign, we were able to robustly classify (after careful visual inspection) 465 objects in the 12 targeted galaxies (see Table 1). Only 11 out of all of these (∼3%) contained features in their optical spectra that indicated a B[e] SG/LBV nature (which was the subject of the current work, with the rest being left for future papers). Other stellar sources related to massive stars included mainly RSGs (∼37%), other Blue Supergiants (∼7%), and Yellow Supergiants (∼ 5%). There was a small number of emission objects (∼2%), carbon stars (∼6%), and AGN/QSO and other background galaxies (∼4%), while another bulk of sources were classified as H II regions (∼22%) and foreground sources (∼14%). In Table 2 we present the identified objects. We note that, although the same approach was followed for all 12 galaxies, we obtained null results for five of them: IC 10, NGC 1313, Sextans A, M83, NGC 6822. In addition, there were only four objects (∼36%) with previous spectral information, for which we confirmed or updated classification. It is also interesting to note that ∼64% of these sources were considered priority targets in our survey ( Table 2, col. 4), while the rest failed to pass our selection criteria (see Section ??). We further discuss these facts in Section 4. Table 2. Properties of the sources identified in this work: source ID in this work (column 1), sky coordinates (columns 2 and 3), priority target (column 4), source ID in Spitzer (base) catalog (column 5), SNR (column 6), spectral type from this work and literature (columns 7 and 8), and radial velocity from this work (RV, column 9).

Spectra
All spectra showed a strong, broadened Hα component, accompanied by several other characteristic emission lines. We present their spectra in Figures 1 and 2, where the strength of the Hα emission for all objects is highlighted in the right panel. The order of the spectra (from top to bottom) was one of decreasing Hα strength.
We identified a series of Fe II emission lines in the left wing of Hα (∼6200-6500 Å), and, when the spectrum extended far enough to bluer wavelengths, we identified another series ranging from roughly ∼5100-5400 Å. Figure 3 showcases these lines in a zoom-in on the ∼6200-6500 Å region. We used the Fe II emission lines in this region to correct for the radial velocity (RV) shift. The obtained RV values are shown in column 9 of Table 2. Therefore, we verified that the RVs were in agreement with the motion of their host galaxies, confirming that these stars were, indeed, of extragalactic origin.
According to the classification criteria presented in Section 2.3, we robustly identified 6 sources as being B[e]SGs: WLM-1, NGC55-1, NGC247-1, NGC253-1, NGC300-1, and NGC300-2. Figure 1 presents the full spectra for the B[e]SGs, while Figure 3 shows the characteristic [O I] λ6300 line. It is particularly interesting to note the very strong He I lines of NGC300-1. These emission lines require a hotter formation region, such as a spherical or a bipolar shell formed by a strong stellar wind, in addition to the structures that give rise to the forbidden emission features. We also note the absence of [Fe II] lines for the WLM-1, NGC253-1, and NGC300-2 sources. Half of the sources (NGC55-1, NGC247-1, NGC300-1) displayed strong [Ca II] emission lines, while for one source (NGC253-1) they were very faint (limited by the noise), and were totally absent for two of the sources (WLM-1, NGC300-2; see Figure 4). These lines were stronger in luminous sources (e.g., [76,77]). The very low SNR for the NGC253-1 and NGC300-2 (see Table 2, column 6) justified the lack of Fe and Ca lines. In the case of WLM-1, the SNR was sufficiently good that the lack of forbidden Fe lines should be considered a real non-detection (similar to source WLM 23 from [78]. We further discuss this in Section 4.2). Unfortunately, due to overlapping slits in the mask design, some of these spectra suffered from artifacts from the reduction processing (in particular, NGC300-2  [80,81]), O VI Raman-scattered lines (λλ6830, 7088) of symbiotic systems (e.g., [81,82]), or even the absorption lines of Li I 6708 present in young stellar objects (e.g., [83]). Moreover, during the visual screening of all these spectra, objects with such characteristic lines would be classified differently, as all possible objects were considered. Additionally, at the distances we were looking at, we were mainly probing the upper part of the Hertzsprung-Russell diagram, while their RVs were relatively compatible (within their error margins) with those of their host galaxies. Our Gaia cleaning approach removed the majority of the foreground sources (naturally, a small fraction remained hidden in our target lists). Therefore, we consider these objects to be strong supergiant candidates.
We characterized as LBVc the following 4 sources: NGC55-2, NGC55-3, NGC247-2, and NGC3109-1 (see Figure 2). NGC55-2 was the hotter of all these sources as it was the only LBVc with all He I lines in emission. NGC3109-1 displayed He I lines in absorption, while the rest did not show any of these lines. During the outbursts the He I lines decrease and vanish, as the temperature and the density (due to the expanding pseudo-photosphere) drop significantly to allow for other lines to form. It is during these cooler states that Fe lines become evident in LBVs. Depending on the conditions, forbidden emission lines may form. This was the case with NGC55-3, which displayed the [Ca II] lines in emission, along with a few [Fe II] lines. The other sources did not show any forbidden lines. Similar to the B[e]SG spectra, there were unavoidable residuals and artifacts, due to the slit overlap and reduction issues.
Of these cases, NGC7793-1 was the most extreme example 5 . The region at [O I] λ6300 was highly contaminated with a sky residual line from another source in the slit. Therefore, we could not conclude whether this line existed or not. We noticed the presence of some [Fe II] and the [Ca II] lines, but a B[e]SG or LBV classification solely from this spectrum was not possible. However, additional information could be retrieved from photometry (see Section 4.3), so that we could propose a B[e]SG candidate (B[e]SG c) classification for NGC7793-1.
The final classification for each star is provided in column 7 of Table 2.

Light Curves and Variability
We collected variability information for all targets from both Pan-STARRS DR2 6 and the VizieR 7 services. We found four sources (WLM-1, NGC247-1, NGC274-2, and NGC3109-1) with data in the Pan-STARRS DR2 release 8 (with an approximate coverage between 2010 and 2014). We considered only values with psfQfPerfect>0.9 to select the best data. For three sources (NGC55-2, NGC55-3, and NGC247-1) we found additional data in the catalog of largeamplitude variables from Gaia DR2 (covering 2014 to 2016; [84]), and NGC3109-1 had already been reported as a variable [85]. In Table 3 we summarize the collected information for all sources and their corresponding magnitude differences (peak-to-peak) for all (5) Pan-STARRS filters, the two Gaia filters (for which we doubled the quoted values in the catalog to match the Pan-STARRS definition of magnitude difference, and some additional variability studies).
In total, we found light curves for two B[e]SGs (WLM-1 and NGC247-1) and four LBVc (NGC55-2, NGC55-3, NGC247-1, and NGC3109-1). We show the Pan-STARRS light curves in Figure 5 and 6, where we plot the magnitude difference at each epoch with the mean for the particular filter (indicated on the y-axis label). For the B[e]SGs we noticed a (mean) variability of 0.25-0.3 mag, while for the LBVs it was slightly larger, at 0.3-0.44 mag. There were no obvious trends in the B[e]SG light curves, while, in the case of NGC3109-1, a dimming across all filters was observed. Menzies et al. [85] also detected such a trend, although smaller, for this target, due to the different filters used. Limited by the photometric data, they argued that a background galaxy or AGN was not excluded, but, given our spectrum and its consistent RV value with its host galaxy, we could actually verify its stellar nature. For NGC247-2, the light curves were generally flatter. There was a noticeable peak present in the y light curve (at MJD∼ 56300 days), which was not evident in the other filters (although we note that there were no observations around the same epoch). The quality flags corresponding to these particular points did not show any issue. However, we should be cautious with this, as further mining of the data is needed to reveal if this is a real event or an artifact.
NGC247-1 was the only source for which we had multiple sources of variability information. Very good agreement between the Pan-STARRS and Gaia data is evident, and consistent with the value quoted by Solovyeva et al. [27] (∆V=0.29 ± 0.09 mag). Although Davidge [86] quoted a smaller value (∆g ∼ 0.1 mag), their time coverage was limited to about 6 months, a time frame that definitely does not cover the whole variability cycles for these sources.
Traditionally, LBVs are considered variable at many scales (e.g., [22,24,87,88]). The (optical) S Dor variability is of the order of 0.1 mag to about 2.5 mag with cycles ranging from years to decades. The giant eruptions, although much more energetic (∼5 mag) are less frequent events (a time frame in the order of centuries), and, therefore, a smaller subgroup of LBVs have been observed to display such events. On the other hand, the B[e]SGs are considered more stable, with variability that does not exceed ∼0.2 mag (optical; [5]). However, this is changing and significant variability is observed, due to binary interactions and possible pulsations (e.g., [12,19,89]). Therefore, it is not surprising to observe similar magnitude differences between the two classes.

Demographics
As mentioned already in Section 3.1, we did not detect B[e]SGs or LBVs in the following five (out of 12) galaxies: IC 10, NGC 1313 Sextans A, M83, NGC 6822.
M83 and NGC 1313 are the most distant galaxies (at 4.9 and 4.6 Mpc, respectively) and confusion becomes an important issue (unsurprisingly, M83 is the galaxy for which we detected the most H II regions; see Tramper et al., in prep.). Due to the spatial resolution of Spitzer and the increasing distance of some of our target galaxies, H II regions or other point-like objects (e.g., clusters) were included in the point-source catalogs and, therefore, considered to be viable targets in our priority system. The farthest galaxies, for which the majority of observed targets were, indeed, resolved point sources and at least one was either an LBVc or a B[e]SG, were NGC 7793 and NGC 253 (at ∼3.4 Mpc). Therefore, the null detections for IC 10 and NGC 6822 (less than 1 Mpc) and for Sextans A (at 1.34 Mpc) were not due to distance and confusion.
Massey et al. [90] detected one LBV in NGC 6822 (J194503.77-145619.1) and three in IC 10 (J002012.13+591848.0, J002016.48+591906.9, J002020.35+591837.6). Our inability to recover these targets was due to two reasons. Firstly, we imposed strict criteria to prioritize our target selection (see Section ??) based on relative strong IR luminosity and color. Almost all of these targets (except for IC 10 J002020.35+591837.6) had m [4.5] > 15.5 mags, which directly excluded them from further consideration. This was further supported by the fact that four out of or our 11 discoveries initially did not pass as a priority target (see Table 2), but were observed as "filler" stars (see Section ??). This was particularly important for galaxies with smaller sizes, where only one (IC 10, Sextans A) or two pointings (NGC 6822) were performed. Therefore, the second reason was the limitations that arose from the particular pointing(s) to the galaxy, as targets might have been located out of the field-of-view or at a sensor's gap (which was the case for IC 10 J002020.35 + 591837.6), and therefore not be observable. Other reasons (not corresponding to the aforementioned targets) that could impact the selection of a target or render its spectrum useless include overlapping slits, a poor wavelength calibration and/or SNR, or other reduction issues.
In total, our approach was successful in detecting these populations, and it was mainly limited by technical issues.

Comparison with Previous Classifications
Four of our sources had previous classifications (see Table 2). WLM-1 had been identified as an Hα source previously [90], through a photometric survey, and identified as an Fe line star through spectroscopic observations (WLM 23 in [78]). Even though the presence of the [O I] λ6300 line was noted, the source was not classified as a B[e]SG, due to the lack of forbidden Fe lines (see e.g., [20,21,92] on Fe stars). Therefore, we updated its classification to a B[e]SG from an Fe star. We also noted that our spectrum (obtained on November 2020) was very similar to theirs (obtained on December 2012), which might indicate that the star was rather stable over this eight-year period (however, this should be treated with caution due to the lack of systematic observations).
NGC55-2 and NGC55-3 had been identified as candidate LBV/WN11 (ids B_34 and B_13, respectively), with both Balmer and He I lines in emission and with P-Cygni profiles [79]. Their spectra were within the 3800-5000 Å range and outside ours. However, given that the diagnostic [O I] line was not present, we classified both of these sources as LBVc, consistent with the previous results 9 .
For NGC247-1 we provided a classification of B[e]SG, similar to what was suggested by Solovyeva et al. [27]. We note here that their spectral coverage was ∼4400-7400 Å which overlapped with our observed range. Hence, we can also comment that no significant differences existed between the two observations (October 2018 and December 2020 by Solovyeva et al. [27] and our observations, respectively), although this time difference is rather small with respect to the variability timescales for these sources [5,12,19]. Therefore, we confirmed the previous classifications for three out of four sources, leaving us with 6 new B[e]SGs (including the reclassified Fe star and the candidate NGC7793-1) and 2 LBVc. The majority (∼72%) of our findings are genuine discoveries and, as such, contribute greatly to the pool of extragalactic B[e]SGs, in particular.

Separating the Two Classes with Photometry
The total numbers of B[e]SGs and LBVs (even including candidates) are definitely small. Combined with the uncertainty pertaining to their roles in stellar evolution theory (e.g., B[e]SGs are not predicted by any code) it is easy to grasp why we really need larger samples and from different galactic environments, to fully understand these sources. Photometric data are typically used to pinpoint interesting candidates. These kinds of diagnostics exist mainly for IR, due to the presence of dust around these objects.
Bonanos et al. [93] found the B[e]SG, LBVs and RSGs to be among the most luminous sources in the mid-IR, using a color-magnitude diagram (CMD) with a combination of near-IR (2MASS) and mid-IR (Spitzer) J-[3.6] and [3.6]- [4.5] for the massive stars in the Large Magellanic Cloud (with a similar work for the Small Magellanic Cloud presented in [94]). In the most recent census of B[e]SGs, Kraus [7] presented color-color diagrams (CCD) to highlight the separation between B[e]SGs and LBVs (see their Figure 5). Indeed, by using the 2MASS near-IR colors H − K and J − H and mid-IR WISE W2-W4 and W1-W2 the two classes are distinct. This is the result of the hot dust component in the B[e]SGs, (formed in the denser disk/ring-like CSE closer to the star) which intensifies the near-and mid-IR excesses, compared to the LBVs (which form dust further away as the wind mass-loss and/or outburst material dissipates). Therefore, the location of a source in these diagrams may be used to verify its nature. We attempted to replicate these aforementioned diagrams by adding the new sources. However, one strong limitation was the lack of data for our sample. For the mid-IR WISE [95] we found data for 5 (out of 11) sources (see Table 4). Using the data for 21 stars (excepting LHA 120-S 111) provided in [7] we plot, in Figure 7, the WISE colors for the MC sources and our 5 objects. We notice that, in general, the newly discovered sources are almost consistent with the loci of the MC sources, with the exception of NGC55-1. The new B[e]SG extend the W2-W4 color further to the red, while the LBVc NGC55-3 extended the W1-W2 color to the blue. Errors were plotted in the cases where they were available 10 . The errors provided for NGC55-1 were (numerically) small and placed it within the locus of LBV. However, caution should be taken with WISE photometry, as the resolution from W1 to W4 worsens significantly, and, combined with the distance of our galaxies, the photometric measurements could be strongly affected by confusion due to crowding (e.g., for both NGC 55 and NGC300 at ∼2 Mpc). Combined with the position (and the uncertainty) of the LBVc NGC55-3, we might also be looking at a systematic offset of these populations. Unfortunately, the points in this plot are too scarce to make a robust examination of how the different galactic environments (e.g., metallicity, extinction effects) affect the positions of these populations. We were unable to construct the J-H vs. H-K CCD because of the lack of 2MASS data for our sources (only for NGC3109-1 did data exist; 2MASS point source catalog; [96]), due to the shallowness of the survey and the distances of our target galaxies. However, we were able to acquire J photometry from the VHS DR5 for 5 of our sources (including NGC3109-1; [97]). Equipped with both J and [3.6] photometry we plot, in Figure 7, the equivalent CMD plot presented in [93], where the underlying MC objects were the same as in [7]. We notice excellent agreement of all new sources to their corresponding classes.
Once again we were hampered by a lack of data for our sample. We could remedy this using the complete data from Spitzer and Gaia surveys (missing NGC253-1 from our sample without Gaia data). In order to consider the MC sources, we used the Gaia DR3 [66,98] and Spitzer data from the SAGE survey [93,94]. This time, we only lost two targets (CPD-69 463 and LHA 120-S 83 without Spitzer data), but were still left with 19 sources.  [7]) and our sample (for 5 out of 11 sources with WISE data). In general, the separation also holds for the new sources, with the exception of NGC55-1 (see text for more). (Right) IR CMD combining near-IR J-band (available for only five of our sources) with Spitzer [3.6]. We notice that, in this case, the newly found sources are consistent with the positions of the MC sources.
In Figure 8 we present the optical (Gaia) CMD, plotting BP-RP vs. M G band. We notice the lack of any correlation in the optical. We included all our sample and the MC sources from [7] (except for two sources without a complete dataset in both Gaia and Spitzer surveys). (Right) The mid-IR (Spitzer) CMD using the IR color [3.6]- [4.5] vs. M [4.5] . In this case, there is a significant improvement in the separation between the two classes. The position of NGC7793-1 favors a B[e]SG nature (see text for more).
In Figure 8 we also present the mid-IR (Spitzer) CMD, plotting [3.6]- [4.5] vs. M [3.6] band. The separation between the two classes becomes more evident in this case. The presence of hotter dusty environments becomes more significant for B[e]SGs, as they looked redder than LBVs (with a [3.6]- [4.5] range between 0.5 to 0.65 mag). They also tend to be much more luminous in the [3.6] than the LBVs. We highlighted the position of NGC7793-1 in this plot. Although, from its spectrum alone, we could not determine a secure classification (due to issues with the obtained spectrum) it is located among the B[e]SGs of our sample and of the MCs. Therefore, we considered it a candidate B[e]SG. A future spectrum is needed to verify the existence of the [O I] λ6300 line, similar to the rest of the secure B[e]SGs in our sample.
We also tried to combine the optical and IR data in a CMD where we plot the [3.6]-[4.5] vs. M G magnitude (Figure 9). The result was actually similar to the previous IR CMD (as the x-axis did not change). In this case, the plot can be more helpful, as the LBVs are populating the upper left part of the plot. Therefore, very bright optical sources with IR color up to ∼0.5 mag were most probably LBVs, while sources with color > 0.5 mag would be B[e]SG (at almost any G magnitude).

Metallicity Dependence of Populations
In this section, we examine the populations of the two classes as a function of metallicity. For this, we plot the cumulative distribution function with metallicity ( Figure 10), considering all detected and known objects in our sample of galaxies. Namely, the numbers presented in Table 2, as well as the two LBVs in NGC55 [79], one in NGC 6822 and three in IC 10 [90], resulting in 7 B[e]SGs (including the NGC7793-1 candidate) and 10 LBVs in our sample of 12 galaxies.
We notice the presence of B[e]SGs at metallicity as low as ∼0.14 Z (WLM). The current work is the first to detect these sources at such low metallicities. The population of LBVs begins at ∼0.21 Z (NGC 3109), and then increases steadily as we move towards higher metallicities. B[e]SGs presents an important step (increase) around ∼0.4 Z . In total, the two populations do not look significantly different. We have to be cautious interpreting this figure, however, due to the low number of statistics and completeness issues, as, for example, depending on the angle under which we observe a galaxy, we may not be able to fully observe its stellar content (e.g., NGC 253).

Conclusions
In this work, we report the detection of 6 secure B[e]SGs, 1 candidate B[e]SG, and 4 LBV candidates sources, of which 6 B[e]SGs and 2 LBVs are new discoveries. They are based on spectroscopic and photometric diagnostics, supplemented with RVs that are consistent with their host galaxies. By inspecting the available IR (2MASS, WISE, Spitzer) and optical (Gaia) CMDs we find that the new sources are totally consistent with the loci of these populations from MCs. This adds further support regarding their natures. Building the cumulative distribution function of both populations with metallicity we notice the presence of B[e]SGs at environments with Z∼ 0.14 Z , which increases the pool of extragalactic B[e]SGs and, especially, at lower metallicities. This is particularly important in order to investigate (with increased samples) these phases of massive stars. Since B[e]SGs and LBVs are among the classes with the most important episodic and outburst activities they provide valuable information on the role of episodic mass loss and insights into stellar evolution in general.  Data Availability Statement: Photometry and 1D extracted spectra will become available through the VizieR/CDS catalog tool.
in collaboration with the Instituto de Astronomía of the Universidad Autónoma de México. OSIRIS was funded by GRANTECAN and the National Plan of Astronomy and Astrophysics of the Spanish Government. This work was based, in part, on observations made with the Spitzer Space Telescope, which is operated by the Jet Propulsion Laboratory, California Institute of Technology, under a contract with NASA. This work made use of data from the European Space Agency (ESA) mission Gaia (https: //www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC, https://www.cosmos.esa.int/web/gaia/dpac/consortium). Funding for the DPAC was provided by national institutions, in particular, the institutions participating in the Gaia Multilateral Agreement. . This publication made use of data products from the Two Micron All Sky Survey, which is a joint project of the University of Massachusetts and the Infrared Processing and Analysis Center/California Institute of Technology, funded by the National Aeronautics and Space Administration and the National Science Foundation. This publication used data products from the Wide-field Infrared Survey Explorer, which is a joint project of the University of California, Los Angeles, and the Jet Propulsion Laboratory/California Institute of Technology, funded by the National Aeronautics and Space Administration. This work used Astropy http://www.astropy.org: a community-developed core Python package and an ecosystem of tools and resources for astronomy [99][100][101], NumPy (https://numpy.org/; [102]), and matplotlib (https://matplotlib.org/; [103])

Conflicts of Interest:
The authors declare no conflict of interest.

Abbreviations
The following abbreviations are used in this manuscript: