© The Authors 2025

1 Introduction

The late evolutionary stages of high-mass stars affect the shape and evolution of the Milky Way. As these massive stars evolve through a series of transitional phases to their endpoint as core-collapse supernovae (Langer et al. 1994), their intense UV fields and conspicuous mass-loss deeply transform the dynamics, structure, and chemical composition of their environs. Among these stages, a short interlude of ~10⁴ yr known as the luminous blue variable phase (hereafter LBV) stands out as a main contributor in terms of its radiative and mechanical output (Humphreys & Davidson 1994) through high mass-loss rates ($\dot{M}$ up to 10⁻⁴ M_⊙ yr⁻¹) and occasional eruptions (shedding up to some ~10 M_⊙ on extremely short timescales). These episodes of enhanced mass loss that may rival the energetic output of supernovae (Smith 2013) leave an indelible footprint in the outskirts of these stars and give rise to complex and heterogeneous circumstellar nebulae of dust and gas (Weis 2001, 2011).

This circumstellar material has been probed by diverse observational strategies. Infrared continuum and spectroscopy revealed the existence of different circumstellar dust populations around some stars, traced by successive mass-loss events (Jiménez-Esteban et al. 2010; Vamvatira-Nakou et al. 2013, 2015; Buemi et al. 2017), and yielded valuable insights into the chemistry and dust lifecycle around these sources (Gull et al. 2020). Continuum emission at centimetre wavelengths traced the distribution of ionised gas around the stars and detected their thermal stellar winds, which enabled accurate estimates of their mass-loss budget (Duncan & White 2002; Umana et al. 2010; Ingallinera et al. 2016) and, under some assumptions, their present-day mass-loss rate (Umana et al. 2005, 2012). Molecular spectroscopy studies of the millimetre window revealed molecular structures that trace swept-up material and previous events of mass ejection. This component is critical to our understanding of the physical processes that drive the evolution of the circumstellar material (i.e. energy and mass-loss balance) (Rizzo et al. 2008; Bordiu et al. 2019, 2021) and its chemical complexity (Rizzo et al. 2014; Morris et al. 2020; Bordiu et al. 2022; Rizzo et al. 2023).

The millimetre region of the electromagnetic spectrum is, a priori, the region in which cold dust and ionised gas can contribute significantly through continuum emission. Observations in this range are therefore needed to characterise and constrain the relative importance of each component. Despite its potential relevance, there has been limited research on the (sub) millimetre continuum emission of LBV stars. In the Galaxy, we only found observational studies of G79.29+0.46 (Higgs et al. 1994), Wd1-243 (Fenech et al. 2018), η Car (Bordiu et al. 2019), and AG Car (Bordiu et al. 2021). In the first two cases, there are only flux measurements at the central source, while in AG Car, it was detected the central source and weak emission from the surrounding nebula. In η Car, the millimetre continuum maps revealed the stellar wind and the extended ejecta from the 1890 outburst.

A handful of extragalactic LBV stars in the Large Magellanic Cloud have also been observed: S61, RMC 127, and RMC 143 (Agliozzo et al. 2017a,b, 2019, respectively), and all of them are known to harbour dense circumstellar nebulae that are visible at infrared or centimetre wavelengths. In the case of S61, neither the central star nor the shell were detected; RMC 127 appeared as a point-like source, but the nebula was not detected; and RMC 143 was clearly detected, surrounded by a massive resolved nebula of dust and ionised gas, different in morphology from the circumstellar material observed at longer wavelengths.

The disparity of results in such a small source sample underlines the relevance of the (sub) millimetre characterisation as a valuable tool in the study of these objects. The millimetre regime is the range in which thermal dust and bremsstrahlung emission compete, so that observations in this band are fundamental for separating their contribution to the overall mass-loss budget.

To mitigate this lack of sampling and gain a better understanding of the processes governing emission in this part of the spectrum, we used the NIKA2 continuum camera at the IRAM 30 m telescope and carried out a pilot study of the 1 and 2 mm continuum emission of five Galactic LBV stars (see Table 1). We list the sources below.

• HD 168607. Located at a distance of ~1.8 kpc (Gaia Collaboration 2022), this source is a confirmed LBV star of spectral type B9Ia+ (Clark et al. 2012) that shows no signs of circumstellar material (Hutsemekers et al. 1994).

• HD 168625. Projected in the sky at ~1 arcmin and slightly closer to HD 168607, this star is a supergiant of spectral type B8Ia+ (Clark et al. 2012). It is enshrouded in an hourglass-shaped nebula that was likely produced by binary interaction (Martayan et al. 2016; Mahy et al. 2016). It is also recognised as a nearly perfect analogue of the circumstellar environment of SN1987A (Smith 2007). Because substantial variability is not fully demonstrated (Clark et al. 2012), it remains a candidate LBV. It is one of the less luminous stars of its type (Smith et al. 2004).

• [GKF2010]MN87 (hereafter MN87). This source was first discovered through its associated bright and asymmetric shell with a radius of 28 arcsec in the MIPSGAL survey (Mizuno et al. 2010). It was promptly identified as as a massive star based on infrared spectroscopy (Wachter et al. 2010; Gvaramadze et al. 2010). It was later proposed to be a LBV candidate on the basis of high-resolution radio observations (Ingallinera et al. 2014, 2016). The central source of the nebula is compatible with a thermal stellar wind and a clumpy ionised shell. Its distance is unknown. At an assumed radius of 0.3–0.6 pc that is typical of other LBV shells, the most likely range of distances is 2–4 kpc.

• [GKF2010]MN101 (hereafter MN101). This source was also discovered in MIPSGAL as a nebula with a diameter of ~20 arcsec and was proposed to be a LBV candidate by NIR spectroscopy of the central source (Wachter et al. 2011; Flagey et al. 2014). Ingallinera et al. (2014, 2016) studied its radio emission and revealed a central source with a canonical thermal wind and an ionised nebula co-spatial with the infrared, characterised by a negative spectral index. Finally, Bordiu et al. (2019) detected a CO torus that is wrapped around the infrared nebula, with a low [¹²CO/¹³CO] ratio that is indicative of CNO-processed material. The distance to this source is uncertain; the most likely range of geometric distances estimated from Gaia parallaxes (Bailer-Jones et al. 2021) is 3.0–6.2 kpc.

• G79.29+0.46. This source is an archetypal LBV candidate (Waters et al. 1996) surrounded by a series of nested infrared shells (Voors et al. 2000; Jiménez-Esteban et al. 2010) associated with successive mass-loss episodes. The innermost shell has an ionised counterpart detected in radio-continuum observations (Umana et al. 2011; Agliozzo et al. 2014), and a molecular counterpart traced by CO emission (Rizzo et al. 2008) that partially interacts with the nearby infrared dark cloud (IRDC) G79.3+0.3 (Palau et al. 2014). Some warm hotspots of ammonia have also been reported in and within the shell (Rizzo et al. 2014). The distance to G79.29+0.46 is 1.75 kpc (Gaia Collaboration 2022).

The continuum maps delivered by the NIKA2 camera have fields of view (FOV) of several arcminutes, enabling the simultaneous study of the stars and their surroundings, with the potential to reveal multiple co-existing emission mechanisms. In Sect. 2, we describe the observations. In Sect. 3 we present the main results, including the overall maps, zoom-ins around the stars, and a study of the millimetre spectral indices. In Sect. 4, we complement the information by constructing spectral energy distribution (SED) diagrams from the optical to radio wavelengths. In Sect. 5, we identify and model the components in the SEDs, which allowed us to draw some general ideas about the dominant emission mechanisms. We conclude with a general discussion of the survey and the individual sources in Sect. 6.

2 Observations

The observations were made with the NIKA2 continuum camera (Perotto et al. 2020) available at the IRAM 30 m radio telescope during three consecutive runs in October 2017 (project 044-17; PI R. Rizzo). The weather was stable during all the sessions, with measured opacities at 225 GHz between 0.2 and 0.3, which is equivalent to a precipitable water vapour between 3 and 5 mm.

For the four regions of interest, we performed a total of 126 individual on-the-fly maps with a duration of three to five minutes each. The scanning angles were 0,45, and 90 degrees to avoid stripping effects in the reconstruction of the final maps. An additional 15% were added in size to the definition of the maps to account for a better reconstruction of the background noise. The pointing was regularly checked, with corrections always below 4 arcsec. The FOV for all sources but G79.29+0.46 was 6.5 arcmin. In G79.29+0.46, we observed with 8.6 arcmin of FOV in order to include the previously known infrared shells and the nearby IRDC.

The NIKA2 camera has three arrays (two at 260 GHz and one at 150 GHz) filled with thousands of kinetic inductance detectors (KIDs) that provide an accurate combination of the responses from off-source detectors to reconstruct the noise template that is to be subtracted from the data. To produce the final images, we used a data-analysis pipeline based on MOPSIC (Zylka 2013) and developed by the NIKA2 collaboration prior to the first public release of PIIC¹. In order to preserve large-scale features, we used an iterative procedure. In the first iteration, the software accounted for the absolute calibration of the flux and the opacity correction and merged all collected maps. As output of this first iteration, we had a first attempt of the total map. This was used to estimate the S/N on the map and extrapolate the region with maximum brightness. The extrapolated structure of the source was used as a mask that was introduced in a second round of iteration. At the end of this latter, we obtained a convergence map that resulted from the difference between the two previous maps. When the brightness of the convergence map exceeded a well-defined threshold in the central region, the process continued and a third iteration started. This method is based on minimising the residuals on the final map and allowed us to recover the largest scales in case of extended emission.

The rms map values of the four fields are depicted in Fig. 1. The values were computed as averages over concentric rings with a width of 2 arcsec. At 1.15 mm (260 GHz), the noise map is rather uniform up to 120–140 arcsec of the centre, while at 2 mm (150 GHz), it extends well beyond 200 arcsec.

Finally, we smoothed the maps with a Hanning 3 × 3 kernel to improve the S/N. The resulting angular resolution (half-power beam width; HPBW) after this step was 12.5 and 18.5 arcsec at 1.15 and 2 mm, respectively. With the smoothing, the rms was reduced by a factor of about 0.36 in all cases.

Fig. 1 Noise map of the observed fields. The rms noise was computed in concentric rings with a width of 2 arcsec. The noise in the 1.15 mm maps is plotted in blue, and the noise corresponding to the 2 mm maps is plotted in red. The field names (abridged) are indicated near the 1 mm curves.

3 NIKA2 results

3.1 Emission from the stars and their surroundings

The final processed maps are presented in Fig. 2. Overall, we note emission from the stars, circumstellar structures, and other likely unrelated features. A close-up of the emission at the positions of the stars is depicted in Fig. 3, and the zoom at 1.15 mm and 2 mm is depicted in the top and bottom rows, respectively. To facilitate a visual inspection of the point-like nature of the sources, we highlight a contour corresponding to 50% of the continuum peak, together with the HPBW at each frequency.

HD 168625 and HD 168607 are detected in both bands, but more clearly at 2 mm. The entire field is relatively devoid of significant emission, except around the stars. In the case of HD 168625, the source appears to be slightly resolved. At 1.15 mm, it shows a faint extended emission plateau that peaks towards the south-west and might be linked to the complex hourglass-shaped nebula that enshrouds the star at infrared wavelengths (Smith 2007). At 2 mm, however, the emission is rather circular, and the maximum agrees well positionally with the star. On the other hand, HD 168607, which lacks reported circumstellar material, is unresolved at 2 mm and very probably also at 1.15 mm, although the low S/N prevents a more detailed analysis.

MN87 is not detected in either band. The field is dominated by a series of filaments that are brighter at 1.15 mm. These structures are well correlated with the dust bands seen in the images at 250 and 500 μm from Herschel/SPIRE (Molinari et al. 2010).

The field of MN101 is dominated by a bright emission region towards the north that likely corresponds to a molecular cloud that is not related to the star. Towards the centre of the map, two components are barely resolved: A central compact emission that matches the position and size of the infrared nebula, possibly embedded in a plateau, and a more diffuse component located about 30 arcsec to the south-east, corresponding to the brightest region of the CO structure reported by Bordiu et al. (2019). While the brightness of the two components is similar at 1.15 mm, the central source dominates at 2 mm.

In the field of G79.29+0.46, the bright IRDC stands out. It extends from the south-east to the west. The central source is well detected in the two bands, with notably brighter emission at 2 mm. A shell-like structure surrounding the star from the north-east is also detected faintly at 1 mm and clearly at 2 mm. This structure has an approximate radius of 100 arcsec, which means that it is most likely associated with the innermost infrared dust shell, as is visible in Spitzer images (Jiménez-Esteban et al. 2010). This millimetre shell is bounded by the CO shells reported by Rizzo et al. (2008).

Table 2 lists the flux densities of the central sources at 1.15 and 2 mm along with their associated uncertainties, centroids, and sizes. They were all derived from Gaussian fitting. The last column also shows the spectral indices, which are discussed in the next section.

Fig. 2 Sky maps of the millimetre continuum emission in the four fields observed with NIKA2. The maps at 1.15 mm are shown at the top row, and the 2 mm maps are shown at the bottom. The field that contains HD 168607 and HD 168625, labelled HD168*, is centred on HD 168625. The colour scales (indicated at the right of each row) are in Jy beam−1. HPBWs are drawn near the bottom left corners of the HD168* images. The central areas indicate the zoomed region displayed in Fig. 3.

3.2 Millimetre spectral indices

We computed the spectral indices of the central sources between the two NIKA2 frequencies (α_mm hereafter) following the convention S_ν ∝ ν^α. Despite the low S/N at 1 mm, it is evident overall that all the detected sources except MN101 exhibit values close to a thermal stellar wind (α_mm ~ 0.6), which is indeed expected for the strong winds driven by LBV stars.

On the other hand, we built the maps of α_mm for all the observed fields after convolving them to a common beam of 20 arcsec. We only used pixels with brightness above a 2σ threshold in the convolved maps. The resulting maps are shown in Fig. 4. They reveal a diverse range of features that are likely due to multiple emission mechanisms.

HD 168625 shows a somewhat negative index, consistent with an ionised H II region. As a consequence of the convolution, the emission at 1 mm is severely diluted and at the centre displays the contribution of the central source plus that of an additional component, which is probably close-in young stellar ejecta. In contrast, HD 168607 displays slightly positive values of α_mm that are more compatible with a thermal stellar wind. This is in line with the values reported in Table 2.

In the case of MN87, the dominant emission arises primarily from filaments that cross the whole field. These features have high values of α_mm, which is typical of thermal dust emission. Considering their morphology as well, we suggest that the millimetre emission in this field very likely arises from foreground or background Galactic clouds.

MN101 presents a more complex scenario in which the central region depicts negative values of α_mm. Again, the NIKA2 beam plays a critical role: The size of the infrared nebula is comparable to the beam, so that the negative spectral index likely captures the combined contribution of the stellar wind and the nebular material, which is known to show non-thermal emission in the centimetre regime (Ingallinera et al. 2016). The nebula is embedded in a free-free plateau and is surrounded by thermal dust (in yellow). This agrees well with well-known molecular structures (Bordiu et al. 2019). The highest spectral indices in this enshrouding structure are located south-west of the star, in coincidence with the brightest part of the CO cloud.

G79.29+0.46 is perhaps the most puzzling object. The central region exhibits a rather negative α_mm, in contrast with the positive value derived for the central source in Table 2 (i.e. not convolved). While the fitted value agrees relatively well with the α = 0.6 derived by Higgs et al. (1994) at 0.8, 1, and 1.2 mm, the negative value observed in the map is a consequence of the severe dilution of the 1 mm central point source (around 0.39).

G79.29+0.46 is known to host a set of concentric shells at angular distances of 100–200 arcsec from the exciting star (Rizzo et al. 2008; Jiménez-Esteban et al. 2010; Agliozzo et al. 2014). The inner shell is clearly visible in Fig. 4, where it shows a remarkable stratification at its eastern side around position (100″, −20″), with three distinct components. Farther away, there is a transition from negative values, which might indicate shocked material, to flat values corresponding to the ionised shell, and finally, to positive values associated with the thermal dust surrounding the shell. A secondary, less pronounced stratification is evident in the south-west region, where the shell interacts with a nearby IRDC (Rizzo et al. 2008; Palau et al. 2014), and where we measure the highest values of α_mm (red in Fig. 4).

Fig. 3 Zoomed sky maps around the central sources. The regions correspond to the regions indicated in Fig. 2. To span the full dynamic range, the colour scale is different for each of the maps. The contours correspond to 30, 50 (highlighted in white), and 70% of the maximum value. HPBWs are indicated near the bottom of the HD168* images. The contours at 50% are highlighted to facilitate a visual comparison with the beam and give an idea about the point-like nature of the central sources.

4 SED analysis

4.1 Building

Aiming to study the relative contribution of the different physical processes that dominate the emission of our targets, we constructed their SEDs along virtually the whole electromagnetic spectrum. We were interested in the study of the circumstellar gas and dust and therefore focused the analysis on the infrared, millimetre, and radio domains.

For each target star, we proceeded in a number of sequential steps that we describe below.

• We constructed a Python script to retrieve all the photometry data available in VizieR. The script was validated after comparing its outputs with the VizieR photometry tool². The structure of the script outputs is similar to that of VizieR, with the addition of the wavelength in microns.

• We initially ran the script with a search radius of 3 arcsec.

• In order to check possible data that were lost beyond 3 arcsec (typically, in the infrared to centimetre wavelengths), we performed a second search by changing the search radius to 30 arcsec.

• We removed duplicated information, which gave more priority to original data and to data that included flux uncertainties, instead of compiling other surveys and data without quoted uncertainties.

• In the VizieR database³, we checked all the literature information about the targets. This step was performed to incorporate possible data from catalogues that were not included in the photometry tool.

• We used the Simbad service⁴ to check and incorporate data from observational papers that were not included in the VizieR database.

Table A.1 (in Appendix A.1) lists all the datasets we used to build the SEDs. In addition to the catalogue names, the VizieR table number, the wavelength range, and the bibliographic reference are included. The Appendix A.2 enumerates a series of additional information about each of the stars that we noted in our tailored search for multi-wavelength photometry.

Fig. 4 Spectral index maps between the two NIKA2 frequencies. The plotted field and marks are the same as in Fig. 2. For each map, all pixels with brightness below 2σ were masked.

4.2 Clarification on particular datasets

Gaia. With the Gaia DR3 photometry data, we additionally checked the consistency of the G-band flux with respect to the BP and RP fluxes. This effect is a consequence of different background estimates for the G and the BP and RP bands⁵. While the G flux is determined by fitting a profile function, the BP and RP fluxes are determined by aperture photometry. The different approach may be critical in our case, where the targets are often very red due to surrounding warm dust and, at the same time, affected by contamination of other (bluer) stars due to their low Galactic latitude. The consistency is measured by the so-called excess factor, defined as the sum of fluxes in the BP and RP bands, divided by the flux in the G band; this factor is slightly above one in most cases. However, we noted very high excess factors in all the sources (1.41, 1.38, 2.14, 2.17, and 1.83 for the sources ordered as in Table 1). This fact might be interpreted by an incorrect measurement of the background in G (due to the inclusion of nebular material), rather than by contamination in BP and RP bands. In consequence, we did not consider the G-band fluxes for the SED fitting.

GLIMPSE. Its photometry is potentially affected by extended nebular emission, mainly due to the additional contribution of line emission and PAHs. A priori, this contamination may affect MN101, MN87, and HD 168625, which are the sources in which the nebular material is not completely detached from the star. Upon inspection of the fluxes, we confirm that MN87 and MN101 are not affected, as their nebular emission is only detected at 24 μm (Mizuno et al. 2010). In contrast, in HD 168625, the VizieR photometry tool only reports flux at 5.8 μm, while the catalogue itself gives data at 3.6, 4.5, and 5.8 μm (all of them are well below the fluxes from other catalogues, however). Consequently, we removed these data points from the HD 168625 SED.

Radio. The radio data we used to build the SEDs come from targeted high-resolution observations and all-sky surveys conducted at various radio facilities. Consequently, the angular resolution can vary by up to a factor of ~10. While the highest-resolution data points measure the flux of the central stellar source directly, the lowest-resolution data points may include emission from the star and the circumstellar material, except when the star is completely detached from the nebula. When we were unable to separate the two components, the data points contained the sum of the flux densities corresponding to the stellar wind and the surrounding nebula (modelled as bremsstrahlung, as explained in Sect. 5.1).

Fig. 5 Spectral energy distribution of HD 168607 in logarithmic scale. The catalogues are indicated at the top right. The corresponding error bars of all points are plotted. The NIKA2 points are shown with blue stars.

4.3 Resulting SEDs

The resulting SEDs are depicted in Figs. 5–9. All of them are highly reddened, as their visual or near-IR peaks are in the range from 1 to 3 μm. This behaviour is indeed expected, considering the long distances of the stars, their location in the Galactic plane, and their nature as dust producers.

Roughly, we distinguish two types of SEDs. The first type is represented by HD 168607 (Fig. 5) and G79.29+0.46 (Fig. 9), where the fluxes decrease monotonically from IR to millimetre to centimetre wavelengths. In both cases, the NIKA2 fluxes are significantly higher than the radio fluxes (up to one order of magnitude). No resolved ejecta structures fall within the NIKA2 beam in the two sources (see Fig. 3), so that the measured flux is most likely the combined contribution of the stellar wind and close-in unresolved dust.

The second type of SEDs is that of HD 168625 (Fig. 6) and MN101 (Fig. 8). The two sources are immersed in conspicuous gaseous and dusty nebulae that contribute to the emission observed with NIKA2. Both SEDs are rather flat in the millimetre-centimetre range, with NIKA2 fluxes comparable to centimetre wavelength fluxes (slightly below in HD 168625 and slightly above in MN101), of the order of 10 mJy. In the infrared, both stars show a strong bump that peaks at ~20–30 μm, which probably indicates a large amount of warm circumstellar dust.

MN87 is a puzzling object. It exhibits strong infrared emission without a well-defined peak. This behaviour is not observed in the other four sources (for which we applied the same processing) and is likely related to the different instrument apertures. In MN87, this effect seems to be a critical factor due to widespread Galactic contamination, as indicated by the dominance of SPIRE dust structures in the NIKA2 maps (see Sect. 3.1). At centimetre wavelengths, the hot star (more properly, its stellar wind) is clearly detected, with flux densities comparable to those of MN101 (Ingallinera et al. 2016).

Fig. 6 Spectral energy distribution of HD 168625.

Fig. 7 Spectral energy distribution of [GKF2010]MN87.

5 Multi-wavelength modelling

5.1 Fitting methodology

The environs of LBV stars are strongly affected by their strong UV radiation and powerful winds. Under adequate physical conditions, LBVs are also efficient dust producers. The continuum emission in these environments is dominated by a combination of thermal dust radiation and free-free emission produced by stellar winds and ionised nebular gas. Non-thermal emission may also be present under special circumstances, such as colliding winds in binary stars, as in HD 168625 (Martayan et al. 2016; Mahy et al. 2016), after mass ejection events, or in the presence of shocks.

The SEDs we compiled and curated in the previous sections demonstrate that our targets depict all the above components with diverse significance at different wavelength ranges. Even more, the millimetre range is key to constraining the contribution of the warm or cold dust, and at the same time, the highly ionised stellar wind. We therefore proceeded to model the observed SEDs with the aim to evaluate the relative importance of these components. The models are intended to be as simple as possible, with the minimum number of parameters capable of providing robust estimates of the different physical components. Because the SEDs always arise from a single position and lack spectroscopic information, we did not attempt to provide fine details about the space distribution and velocity field of our targets. In some cases, we assumed reasonable physical quantities after determining their influence in the flux predictions. For each modelled component, Appendix B details these aspects in depth.

The first step was to correct the stellar fluxes by interstellar extinction according to the Bayestar19 3D dust model (Green et al. 2019). We assumed the distances quoted in Table 1, adopting 3 and 4 kpc for MN87 and MN101, respectively. We used the Python package dustmaps (Green 2018), which provides extinction values for grizy and JHK bands, and interpolated the available photometry by a cubic spline. After this initial correction, we noted still significant remaining reddening in all the sources. This is indeed expected considering the presumed presence of circumstellar dust close to the extended star atmospheres. A follow-up determination of the total extinction was subsequently performed for wavelengths smaller than 4.6 μm (for details see Green et al. 2019), which allowed us to separate the interstellar and the circumstellar extinction at each optical or near-IR band. For this step, we considered effective temperatures (T_eff) from 8000 to 20 000 K, which is the range corresponding to all the known LBVs (Smith et al. 2004).

We modelled the extinction-corrected SEDs as a sum of independent flux contributors. We included four types of components: (1) a central star, modelled as a black body; (2) one or more circumstellar dust clouds, conceived as black or grey bodies at different temperatures; (3) stellar wind, modelled as a homogeneous and isotropic mass loss at a constant rate; and (4) free-free emission arising from a homogeneous ionised nebula surrounding the hot star. Appendix B describes all these components in detail, including the free parameters and assumptions involved.

It is challenging to model a SED from the optical to the radio bands (i.e. more than six decades). The spectral regions have diverse coverage, aperture, and sampling; and all the flux measurements are also potentially affected by intrinsic variability of the sources, geometry, and composition, among other factors.

We followed an iterative fitting procedure to ensure that each component was constrained progressively. First, we divided each SED into four spectral regions, bounded by the wavelengths in which the contribution of two components was comparable. In general terms, the four regions are (1) optical to a few microns, dominated by the contribution of the stellar photosphere; (2) near-infrared up to some tens of microns, with contributions mainly from the star and a hot dust component; (3) mid-infrared up to the NIKA2 points, dominated by warm- or cold-dust components and free-free emission; and (4) millimetre and centimetre wavelengths, mostly dominated by ionised gas and/or the stellar wind. We adopted a stepwise approach, starting from the optical region and progressively incorporating additional components, as described below.

• We started by fitting the optical part of the SED by the stellar photosphere alone. The parameter space representative of LBV stars (T_eff from 8000 to 20 000 K and stellar radius R_star from 30 to 300 R_⊙) was employed as a kick-off exploration and was later refined through χ² minimisation.

• Then, we proceeded to fit the near-infrared region. The star component was combined with an additional greybody representing a hot-dust component. The stellar parameters were kept within a narrow range around the previously determined parameters, while the dust parameters were probed in wide ranges (R_dust from 0.5 to 20 au, (T_dust from 500 to 3000 K, and β from 0 to 2). We repeated the χ² minimisation to locate the new best-fit parameters for the star and the hot-dust components.

• As a next step, we refined the fit over the two previous ranges to ensure consistency between the two components and further constraining the parameter space around previously found minima, thus obtaining a new set of best-fitting parameters.

• The above steps were repeated by iteratively incorporating the remaining components at the other wavelength ranges, until the whole SED was fitted and the best-fit parameters for all the components were determined.

Fig. 8 Spectral energy distribution of [GKF2010]MN101.

Fig. 9 Spectral energy distribution of G79.29+0.46.

5.2 Fitting results

The NIKA2 upper limits of MN87 are not sufficient to constrain all the SED components and keep a physical meaning. On one hand, the strong infrared emission of this source requires a very high emissivity index (β ≫ 2) to obtain the necessary drop at NIKA2 frequencies. On the other hand, the ~10 mJy level of the radio emission requires a deeply negative value of α. In order to find some explanation, both parts of the spectrum must be taken into account; it is known that non-thermal processes are important in this source (Ingallinera et al. 2016), and the mid-infrared SED points may be heavily contaminated by background emission.

For the other stars, Table 3 summarises the best-fit results. It includes the fundamental stellar parameters (T_eff and radius R_star), the interstellar and circumstellar absorption in g band, the parameters of the dust components (dust temperature T_dust, the characteristic radius R_dust, and the total dust mass M_dust) and, depending on the case, the current mass-loss rate $\dot{M}$ of the stellar wind, or the parameters of the free-free model (electron density n_e and size of the ionised region R_ff). Furthermore, Fig. 10 displays the four best-fit models, together with the individual components and the de-reddened photometry. The four SEDs are plotted at the same scale to facilitate visual comparison.

It is noteworthy that the two stars with the highest T_eff, HD 168607 and G79.29+0.46, are also the two whose SEDs were fitted by including a stellar wind component. In contrast, HD 168625 and MN101 were fitted by the low-end of the usual values of T_eff, and both were fitted by free-free components. The circumstellar extinction is high in all cases, being HD 168607 the less reddened of the sample. The warm- or cold-dust masses are remarkably similar, of the order of 10⁻³ M_⊙. The hot dust masses, however, are three orders of magnitude lower, except for the much higher M_dust of HD 168625, which is in line with the huge IR emission from its dust nebula. The mass-loss rates inferred from the stellar winds of HD 168607 and G79.29+0.46 are of the order of a few 10⁻⁶ M_⊙ yr⁻¹ (scaled to v_∞ = 100 km s⁻¹; see Appendix B for details).

In the four cases, the relevance of the NIKA2 data is evident. At 1–2 mm, it is clearly noted that the observed flux results as the combined contribution from the ionised stellar wind (or ionised nebula) and the cold dust. The NIKA2 fluxes are vital to constrain the relative contribution of the two components arising from the far-infrared and from centimetre wavelengths.

Below, we compare our fitting results with those available in the literature.

HD 168607. The best-fitting model predicts an effective temperature of 17 000 K, which is higher than the value of ~10⁴ K determined by Clark et al. (2012). The luminosity resulting from our fitting is 10^5.88±0.19 L_⊙, which agrees well with the estimates provided in the same work, and inferred from a visual inspection to one of its figures. A hot-dust component lies very close to the star (T_dust = 1100 K) and causes most of the emission in the near-infrared, and a more extended but still unresolved cold-dust component (35 K) that causes the bump observed around 100μm. The radio part of the SED is well fitted by a stellar wind. Leitherer et al. (1995), based on 8.64 GHz continuum, determined $\dot{M}$ = 2.3 ± 0.5 × 10⁻⁶ M_⊙ yr⁻¹, using 2.2 kpc and 140 km s⁻¹ for the distance and for v_∞, respectively. Their value of M corrected for the distance of 1.84 kpc shifts to 1.5 ± 0.3 × 10⁻⁶ M_⊙ yr⁻¹. On the other hand, when we correct our $\dot{M}$ considering v_∞ = 140 km s⁻¹, we obtain 1.4 ± 0.4 × 10⁻⁶ M_⊙ yr⁻¹, which agrees perfectly with Leitherer et al. (1995).

HD 168625. According to Table 3, the luminosity derived from our model is 10^4.93±0.27 L_⊙. The strong infrared emission was fitted by a combination of a blackbody and a greybodies. The blackbody corresponds to a compact dust component of T_dust = 500 K (below the 1000–2000 K of the other targets by far), while the greybody represents a more extended warm-dust component (T_dust = 140 K). This is the only source in which a greybody component was required to fit the data because it was the only way to reconcile the far-infrared emission with the NIKA2 millimetre measurements.

The best-fit emissivity index, β = 1.0 ± 0.2, indicates the dust to some extent. Theoretical models (Draine & Lee 1984) predict β values between 1 and 2, with the higher values in the diffuse ISM (Planck Collaboration 2014) and the lower values in dense cores (Forbrich et al. 2014) and circumstellar disks (Friesen et al. 2018). Physical conditions in the outskirts of evolved massive stars differ markedly from those in the ISM, which are characterised by diffuse radiation fields, low temperatures, and low densities. The close environs of these hot stars are bathed by strong UV radiation fields, but at the same time, self-shielding is possible due to the continuous mass injection through the dense stellar winds.

Therefore, the value of β derived from the fitting is indicative of large grains that formed through continuous accretion, favoured by the current photospheric temperature. Kochanek (2011) proposed that episodes of enhanced mass loss might be produced during the cold phase of LBV stars and create appropriate conditions (in terms of density and self-shielding) for the grain growth. This scenario is further supported by the study of the dust around the LBV HR Car, performed by Buemi et al. (2017), which revealed a radial gradient in β, with lower values near the star and an increase at larger distances, where smaller grains mix with the ISM.

In order to proceed to a fair comparison of our estimates with other published parameters, it must be considered that the assumed distance to HD 168625 may significantly differ from one paper to the next, with obvious systematic deviations of the quoted values. We therefore re-scaled the published parameters to 1.53 kpc (Table 1) when necessary.

O’Hara et al. (2003) imaged the nebula at several near- and mid-IR wavelengths, assuming a distance of 2.8 kpc. After scaling their distance to ours, their fundamental stellar parameters closely agree with our estimates: T_eff = 14 000 K and R_star = 82 M_⊙. They also obtained a dust mass of ~7.5 × 10⁻⁴ M_⊙, again after correcting to our distance. This is somewhat lower than the values shown in Table 3, but still within the uncertainties. Even so, the difference may indicate a real value of κ₀ that is slightly higher than we assumed in this work, probably closer to 0.9 cm² g⁻¹.

The stellar luminosity computed by us is substantially lower than that of Mahy et al. (2016), which is 10^5.58±0.11 L_⊙ also assuming a distance of 2.8 kpc. This difference cannot be explained in terms of uncertainties. After correcting their value to a distance of 1.53 kpc, however, we obtained ~10^5.05 L_⊙, which is fully compatible with our luminosity estimates. In the next section, we discuss the implications of this revised lower luminosity.

Our warm-dust temperature is in accordance with several studies. O’Hara et al. (2003) found dramatic spatial changes in the physical dust properties in the nebula, with a representative value of T_dust = 120 K. Strong temperature variations were also found by Umana et al. (2010), with values between 110 and 210 K (one of the figures impressively show a large dust plateau with a characteristic temperature of ~120 K). Blommaert et al. (2014) included this source in a sample of evolved stars observed with Herschel/PACS at the 69 μm olivine band, finding a strong signature of the band and determining a dust temperature range of 50–150 K. Furthermore, Arneson et al. (2018) made a sophisticated dust model mainly based on data from SOFIA and determined T_dust = 70 ± 40 K and M_dust = 2.5 ± 0.1 × 10⁻³ M_⊙. To summarise, the values of the cold-dust component quoted in this work agree remarkably well with the above works and underline the validity and reliability of the fitting.

The stellar wind model of Umana et al. (2010) is highly appropriate to separate the star contribution from the total flux (i.e. star plus nebula) at centimetre wavelengths, but it is not able to account for the (relatively low) NIKA2 fluxes. Furthermore, a bright and highly ionised nebula is known to be excited by HD 168625. We therefore proceeded to model the ionised gas by a free-free component, which (in combination with the greybody) successfully explained the SED at far-infrared, millimetre, and centimetre wavelengths.

MN101. To the best of our knowledge, there are no estimates of the fundamental stellar parameters of MN101. Our best-fit T_eff is exactly equal to that of HD 168625, while R_star is slightly smaller. These parameters are compatible with a low-luminosity relatively cold LBV (Smith et al. 2004). The strong infrared emission that peaks in the 20–30 μm range is well modelled by an au-size hot-dust component at 1800 K, plus a colder one at 160 K that extends to 100 au. Considering the similarity between MN101 and HD 168625, we modelled the centimetre wavelength emission using a simple free-free model, which reconciled the different contributions in the vicinity of the star. In Fig. 10, we note that the thermal 160 K dust and the free-free emission from the ionised gas provide comparable fluxes to the total millimetre emission observed by NIKA2.

We emphasise that this fitting constitutes just a first (and rather crude) approximation. In any case, it provides the first estimates of the fundamental star parameters and some global values of its closest circumstellar material (the good coincidence of our parameters in HD 168625 and those reported in the literature supports our approach). Other different mechanisms may clearly be at play in the stellar surroundings (e.g., strong stellar winds and non-thermal emission; Ingallinera et al. 2016), but they are not significant enough compared with the 160 K dust plus the ionised region associated with the nebula in the centimetre and millimetre wavelength domains.

G79.29+0.46. The best-fit stellar parameters of this source are similar to those of HD 168607. We computed a luminosity of 10^5.78±0.21 L_⊙. This luminosity and the fitted T_eff (18 000 K) place this target in the middle of the S Doradus instability strip of the HR diagram (Smith et al. 2004). This luminosity is also close to that quoted by Agliozzo et al. (2014), 10^5.4 L_⊙ in a study that focused on fitting the SED in the visual and near-infrared.

Waters et al. (1996), based on infrared spectroscopy, determined a mass-loss rate of ~10⁻⁶ M_⊙ yr⁻¹ (for a measured wind velocity of 94 km s⁻¹), which excellently agrees with our values. The wind velocity was also measured by Voors et al. (2000), who determined a value of 110 km s⁻¹; therefore, the $\dot{M}$ quoted in Table 3 seems close to the real one, considering that it is scaled to 100 km s⁻¹ of the wind velocity. It is worth noting that in Voors et al. (2000), the necessary stellar input to explain the infrared brightness distribution also closely agrees with our estimates: T_eff = 20 000 K and R_star = 76 R_⊙ (equivalent to a luminosity of 10^5.92 L_⊙).

Kraemer et al. (2010) also modelled the infrared part of the central source and determined a dust temperature of 1800 K. This is not far from our results, but these authors only considered data up to 24 μm and did not attempt to include data at larger wavelengths. In contrast, Agliozzo et al. (2014) fitted the whole infrared range by summing five greybodies and adopting β = 2 for all of them; the temperature of these components ranges from 40 to 1200 K, which is consistent with our two (hot and cold) dust components.

Fig. 10 Best-fit models. The target names are indicated in the top left corners of each panel. Individual components (star, hot dust, warm or cold dust, and stellar wind or Bremsstrahlung) are indicated by the blue, orange, red, and green lines, respectively. The combination of the models (best fit) are shown in bold black lines. The grey area encloses the combined uncertainties of the fitting. The magenta dots with error bars depict the de-reddened photometry. The four charts are depicted with the same wavelength and flux density scales to facilitate the comparison among the targets.

6 Concluding remarks

We detected the central sources in four out of the five surveyed LBVs at both millimetre bands. We also found extended emission that reaches up to some 1–2 pc from the stars, although not all the features seem to be clearly associated with the targets. This is a step forward in the field because the millimetre continuum has only been observed in a handful of massive stars (Leitherer & Robert 1991; Higgs et al. 1994; Agliozzo et al. 2017b, 2019; Fenech et al. 2018).

With the 30 m beam size at the NIKA2 frequencies, our data can resolve structures that are typically larger than ~10⁴ au at the distances of our targets. We distinguished several features in the millimetre maps: (1) unresolved sources, corresponding to the stellar winds, probably enshrouded by ionised gas and dust; (2) a plateau of extended emission, particularly noted in HD 168625; (3) a totally detached shell in G79.29+0.46; (4) some Galactic clouds that are not clearly linked to the stars; and (5) the IRDC in the G79.29+0.46 field.

These features may be characterised not only by their morphology, but also by their spectral indices (α_mm). For the point sources, except for MN101, α_mm is close to 0.6, which is the expected result for a canonical stellar wind. This trend becomes different when α_mm is computed for the extended emission. The Galactic clouds in the fields of MN87 and MN101 and the IRDC in the field of G79.29+0.46 have the highest values of α_mm (well above 2.5), which is typical of thermal dust emission. At the opposite end, there are the negative values found in the close environs of the central sources, where the stellar wind contribution seems diluted by more extended nebular material. In the field of G79.29+0.46, we noted two areas with clear gradients of α_mm, from negative to positive values farther away from the position of the star; we propose an explanation of this behaviour by some sort of low-velocity shocks acting on one of the wellknown shells (Rizzo et al. 2008), and also in an interaction region with the IRDC (Palau et al. 2014).

In order to provide a better understanding of the relative importance of the excitation mechanisms, we complemented our NIKA2 data by an exhaustive search of archives, surveys, and dedicated observations available in the Virtual Observatory. For each target, we constructed six-decade SEDs with information from the visual to the radio regimes.

After we curated the data, we distinguished two types of SEDs that are determined by the relative flux density at millimetre wavelengths compared to their centimetre counterparts. The main feature of the first group (to which HD 168607 and G79.29+0.46 belong) is that the millimetre flux density is significantly higher (one order of magnitude) than the centimetre flux. In turn, the second group is characterised by comparable flux densities at millimetre and centimetre wavelengths. In the five cases we observed, the SEDs show that the sources are highly reddened and depict strong near- and mid-infrared emission.

We also proceeded with a (somewhat generic) modelling of the SEDs, with the objective of identifying different physical components and their relative contribution. The approach we used kept the number of components and parameters as small as possible. In this way, we were able to confirm or predict the existence of the modelled components, without including sophisticated information about structure, composition, or any other properties.

The fitting was satisfactory overall, considering that the data spanned six decades of the electromagnetic spectrum, the diversity of the instruments included, and the minimal number of free parameters employed. Notably, it allowed us to unveil unresolved hot-dust components at only four to six times the stellar radius in all cases but HD 168625, where the 500 K gas may extend up to 50–60 R_star. It is known that these close-in structures are found in other massive evolved stars, such as B[e] supergiants, which are surrounded by circumstellar disks Zickgraf et al. (1985). In the context of LBV stars, the detection of au-sized hot-dust structures (500 to 2000 K) around active and quiescent sources may be explained by a dynamics where the ejecta expelled during a past mass-loss event falls back to the star, as proposed by some hydrodynamical simulations (Owocki 2013; Owocki et al. 2019). This fallback would produce dense rotating disks in which dust grains can survive against the hot radiation from the central star and grow efficiently, helped by self-shielding. Very high resolution IR and millimetre observations are needed to test this plausible but highly speculative interpretation.

In HD 168607, the mass-loss rate was better constrained by the combination of a better knowledge of its distance, the use of a measured wind velocity, and a fitting that considered the contribution of more extended cold dust.

The revised luminosity of HD 168625 would place it in the lowest part of the LBV zone in the HR diagram (Smith et al. 2004). This star has long been regarded as a twin of Sk-69 202, the blue supergiant identified as the progenitor of SN1987A (Walborn et al. 1987), owing to the striking similarities of their circumstellar material (Smith 2007). With the revised parameters from this work, the resemblance of HD 168625 and Sk-69 202 is even stronger: HD 168625 is nearly at the same position in the HR diagram as it was Sk-69 202 before the SN explosion (Smartt et al. 2002, see Fig. 6), that is, with nearly identical physical properties.

The non-detection of MN87 is intriguing. The upper limits for the central source at both NIKA2 frequencies are well below the expected fluxes from thermal cold dust or ionised gas. While a non-thermal contribution might help to explain these low fluxes, its inclusion is challenging considering the significant IR emission and the thermal spectral index of the central source at centimetre wavelengths (Ingallinera et al. 2016).

We presented the first estimates of the fundamental parameters of the central source of MN101. These results are compatible with a LBV close to its coldest value (T_eff = 11 000 K). Its SED is very similar to that of HD 168625, although its hot-dust component is the most compact (1–2 au) and hottest of the sample (T_dust = 1800 K).

In G79.29+0.46, the sub-millimetre, millimetre, and centimetre-wavelength flux (and consequently, the energy input from the hot star) remained stationary during at least 22 years (see Higgs et al. 1994; Crossley et al. 2007; Umana et al. 2010, and our Appendix A.2). In the absence of a systematic photometric monitoring of G79.29+0.46, the apparent steadiness of the stellar wind and the dust properties strongly supports the candidate LBV status of this source (Voors et al. 2000).

It is remarkable that the stars that share a given type of SED also display similar T_eff. Specifically, the two stars with a monotonically decrease from infrared to centimetre ranges (HD 168607 and G79.29+0.46) have the highest T_eff, and their millimetre and centimetre SEDs are well-fitted only by stellar winds, while the other sources show the opposite behaviour (low T_eff and free-free components). There is no clear explanation of this effect, although the sample is admittedly small for any conclusion.

This pilot study demonstrates the importance of studying the millimetre wavelength continuum emission in this type of object. The NIKA2 frequencies fall between the far-infrared and the radio domains, which are overall constrained by cold dust thermal emission and Bremsstrahlung mechanism, respectively. Therefore, the NIKA2 frequencies are very well suited to understand the relative contribution of these two mechanisms in detail, with the consequent determination of important parameters for the energy that is injected by the LBVs in their surrounding circumstellar material. Likewise, it suggests the existence of different families of SEDs that are closely tied to the stellar parameters. Follow-up works addressing a larger and statistically significant sample of LBV stars can shed further light on these aspects.

Acknowledgements

We would like to thank the IRAM staff for their kind and professional support during the observations and the NIKA2 core team for providing the data analysis software that was used to reconstruct the maps shown in the article. J.R.R. acknowledges support by PID2022-137779OB-C41 funded by MCIN/AEI/10.13039/501100011033 by “ERDF A way of making Europe”. This research has made use of the Spanish Virtual Observatory (https://svo.cab.inta-csic.es) project funded by MCIU/AEI/10.13039/501100011033/ through grant PID2023-146210NB-I00.

References

Appendix A Spectral energy distributions

Appendix A.1 Data origin

In Table A.1 we detail basic information about the catalogues, surveys, and observations used to build the SEDs.

Appendix A.2 Data notes on individual sources

In this section we report some particular notes which complement the information presented in the SEDs.

HD168607 (Figure 5)

• Pan-STARRS photometry shows significant dispersion (up to a factor of 4 in g) so it was discarded.

• The WISE band 4 measurement is heavily contaminated by nearby extended emission not related to the star. As no reliable flux could be determined, the measurement was discarded.

• We found point like emission at 70 and 100 μm with Herschel/PACS. The nearby HD168625 is a target of a project (observation IDs from 1342217763 to 1342217766), and the observed images include by chance the position of HD168607. We extracted a field around HD168607 from the original images and fitted a two-dimensional Gaussian source. The results are shown in Fig. A.1 and in Table A.2. At 160 μm the star is no longer detected.

• This source is catalogued in the GLOSTAR survey at 1.4 GHz (Yang et al. 2023), but not in CORNISH, MAGPIS, or THOR. In the search for references not available in VizieR, we found a survey made with ATCA at 8.54 GHz (Leitherer et al. 1995), which was added to the SED.

Fig. A.1 Herschel/PACS images around HD168607 at 70 μm (left) and 100 μm (right). Equatorial coordinates are relative to the star position. The continuum emission at these wavelengths is reported for the first time in this work.

HD168625 (Figure 6)

• The source is included in the IRAS Point Source Catalogue, although not detected at 100 μm.

• In the PACS images at 70 and 100 μm the star is not detected. In Hi-Gal, the nebula is listed as two different sources (the brightest parts of the nebula); we consequently summed up the two fluxes, and assumed uncertainties of 10%.

• A complete study in the infrared was done by Arneson et al. (2018), who included their SOFIA measurements together with 2MASS, AKARI, WISE, MSX, Spitzer and Herschel data. We cross-checked the spectral fluxes in the paper with our ones obtained from VizieR, and added the SOFIA and Herschel/PACS data to our SED.

• The WISE band 4 flux is clearly overestimated, up to three times the fluxes of MIPS, MSX, IRAS and AKARI in nearby wavelengths. Such effect is probably due to the large aperture of WISE at 22 μm, which is critical in presence of the strong IR nebula. We therefore remove this point from the SED.

• The source is not included in the THOR catalogue, but we downloaded the images and detected some extended emission. We averaged three frequencies (1690, 1820 and 1950 MHz) to get a reliable source detection, and incorporated the average to the SED. This source was also observed by Leitherer et al. (1995) at 8.54 GHz using ATCA, and we therefore also included those measured values.

• Umana et al. (2010) observed this target at high resolution with the VLA at 8.54 GHz, detecting the central point source (0.42 mJy) while resolving out most of the extended nebula; we therefore decided to take this measurement into account for the fitting, as explained below. The VLASS survey also detected this point source at 2–4 GHz, with a surrounding ring of 5–10 arcsec.

• After a close inspection to the data, we noted that the THOR flux is higher than the adjacent ones. We explored whether this is due to an incorrect flux measurement⁶, or due to the possibility of reaching the turnover frequency close to the VLASS frequency (3 GHz). To do so, we downloaded and analysed the individual images available in the THOR Image Server⁷. These images correspond to epochs 1.2 and 2.2, but only in the first one the source is detected. Our measured flux is compatible with that provided in the archive, so we confirm its validity and used the tabulated value and its uncertainty.

MN101 (Figure 8)

• The Gaia G flux is anomalously high with respect to the other Gaia bands, and significantly above measurements of other surveys. This effect may be due to the presence of circumstellar dust which contaminates the star photometry, as explained in Sect. 4.2.

• This source is included in most of the radio surveys searched. The flux densities from THOR were quite discrepant with respect to others, so we recomputed them using aperture photometry. Our values are more consistent with the other surveys, and we therefore decided to incorporate them to the SED, instead of the catalogued ones.

• The nebular feature was observed with the VLA (Ingallinera et al. 2014) at 1.4 and 5 GHz, with flux densities of 10.2 ± 0.8 and 10.5 ± 0.1 mJy, respectively. Considering that the angular extent of the source (~20 arcsec) is comparable to the NIKA2 beam, we preferred these values for the fitting. Higher resolution 5 GHz observations by Ingallinera et al. (2016) revealed a central point source (7.0 ± 0.2 mJy), with a spectral index compatible with a thermal stellar wind and an ionised circumstellar material characterised by a non-thermal spectral index.

G79.29+0.46 (Figure 9)

• The source is highly reddened. G79.29+0.46 has been observed in three different epochs by SDSS; the measurements differ up to a factor of two, probably due to stellar variability. The PAN-STARRS survey, made some years after the SDSS data release, shows a similar behaviour. In fact, G79.29+0.46 is included in the catalogues of Heinze et al. (2018) and Chen et al. (2020) of variable stars. We decided to keep different points for SDSS and PAN-STARRS to reflect the source variability.

• Agliozzo et al. (2014, referred to as Agl14 in the SED) observed this source at 8.46, 4.96, and 1.4 GHz using the VLA interferometer, and made a SED study in the infrared domain using observations of Herschel, Spitzer, ISO and WISE. We cross-checked our values in common (e.g., from WISE and 2MASS) and completed our SED by incorporating their values from Spitzer and Herschel.

• VLASS does not include this source in its catalogue, although it is noted some point-like emission in the images corresponding to the three epochs observed. We therefore computed the integrated fluxes and averaged them with a inverse-square weighting of the flux uncertainties.

• The NVAS archive also provides images in the C-, X-, K-, and Q-bands. We computed the integrated fluxes in these images as well, and incorporated them to the SED. The flux uncertainties in C- and X-bands are notoriously large in comparison to those of the targetted observations of Agliozzo et al. (2014) at similar frequencies; we therefore decided to keep these NVAS data in the SED, but remove them as input for the fitting.

• Higgs et al. (1994) made a multi-wavelength study of the central source and the surrounding shells, including data from IRAS and targetted observations of several instruments from 0.8 mm to 5 GHz. The fluxes measured at all these frequencies are consistent with our NIKA2 values and the above mentioned radio surveys (NVAS, VLASS, Agl14), which demonstrates that G79.29+0.46 remained stationary at centimetre- and millimetre-wavelengths at least during the last 20 years.

Appendix B Modelling details

Central star

The contribution of the star was modelled as a black body, parametrised by the effective temperature T_eff, the stellar radius R_star, and the distance d to the star. In this way, the flux density is: $$S_{\nu }=\mathrm{\Omega }{B}_{\nu }(T_{\text{eff}})=\mathrm{\Omega }\frac{2h{\nu }^3/c^2}{\exp \left(h\nu /k{T}_{\text{eff }}\right)-1}$$(B.1)

with Ω being the solid angle subtended by the source, Ω = π (R_star/d)², and B_ν(T_eff) the Planck function.

The simplicity of this model is adequate for the purpose of this article. The uncertainties introduced by this simplification is by far below those due to the variability of the stars and the heterogeneity of the collected visual and near-infrared data. This fitting is sufficient to provide good approximations of T_eff and R_star, and to constrain the neighbour range of wavelengths, where the hot dust is dominant.

Circumstellar dust

The emission arising from this component was modelled by modified black bodies (sometimes referred to as grey bodies), parametrised by a dust temperature T_dust, a characteristic radius of the emitting region R_dust (expected much larger than R_star), and a dust opacity law (Ossenkopf & Henning 1994): $${\kappa }_{\nu }={\kappa }_0{\left(\nu /{\nu }_0\right)}^{\beta }$$(B.2)

In the above equation, κ_ν is the dust mass opacity in units of cm² g⁻¹, κ₀ is the dust opacity at the reference frequency ν₀ (here chosen at 230 GHz) and β is the dust emissivity index. The measured flux density of the modified black body results: $$S_{\nu }=\sigma {\kappa }_{\nu }{\mathrm{B}}_{\nu }(T_{\text{dust}}){\mathrm{\Omega }}_{\text{cloud }}$$(B.3)

where σ is the dust mass column density along the line of sight (in units of g cm⁻²) and Ω_cloud the solid angle subtended by the cloud of radius R_dust. If the total dust mass of the cloud is M_dust, we can roughly compute $$\begin{array}{c}\sigma =\frac{M_{\text{dust }}}{\pi R_{\text{dust }}^2}\\ {\kappa }_{\nu }={\kappa }_{{\nu }_0}(\nu /{\nu }_0{)}^{\beta }={\kappa }_{{\nu }_0}({\lambda }_0/\lambda {)}^{\beta }\\ \mathrm{\Omega }=\pi (R_{\text{dust }}/d{)}^2\end{array}$$

and the flux density results: $$S_{\nu }=M_{\text{dust}}{\kappa }_0{\lambda }_0^{\beta }{\lambda }^{-\beta }{\mathrm{B}}_{\nu }(T_{\text{dust}})/d^2$$(B.4)

For low to moderate gas densities, dust opacities of MRN grain (Mathis et al. 1977) with and without ice mantles ranges from 0.3 to 0.9 cm² g⁻¹. So far, we estimated the dust masses in Table 3 by assuming κ₀ = 0.5 cm² g⁻¹, which yields uncertainties due to unknown κ₀ within a factor of two. In this way, the can fit the density flux arising from a dust cloud in terms of β, M_dust, and T_dust.

The parameter β encapsulates key dust properties, such as composition, porosity, and grain size distribution. Therefore, single-β grey bodies are just a simplification without a specific physical meaning (Jones et al. 2013). Even so, they provide useful information about the dust behaviour. In our modelling we considered β as a free parameter, although it was indeed necessary only for the cold component of HD168625, where the best-fit was found for β = 1 ± 0.2. For the other cases, the best-fit value of β was zero, i.e., all other components are well-fitted as ideal black bodies.

When the best-fit corresponds to a black body, we can additionally derive R_dust because the density flux is also: $$S_{\nu }={\mathrm{B}}_{\nu }(T_{\text{dust}})\mathrm{\Omega }={\mathrm{B}}_{\nu }(T_{\text{dust}})\pi {\left(R_{\text{dust}}/d\right)}^2$$(B.5)

In practice, we fitted the SED through the expression: $$S_{\nu }=\frac{1}{(\lambda /1\mu \mathrm{m}{)}^{\beta }}{\mathrm{B}}_{\nu }(T_{\mathrm{d}\mathrm{u}\mathrm{s}\mathrm{t}})\pi (\hat{R}/d{)}^2$$(B.6)

where $\hat{R}$ is considered as a fiducial radius. By equating Eqs. B.4 and B.6 we can derive, after some algebra, a relationship between M_dust and $\hat{R}$: $$M_{\mathrm{d}\mathrm{u}\mathrm{s}\mathrm{t}}=\frac{\pi {\hat{R}}^2}{{\kappa }_0{\left({\lambda }_0/1\mu \mathrm{m}\right)}^{\beta }}$$(B.7)

When β ≠ 0, $\hat{R}$ is just a proxy to derive M_dust, without a real physical meaning. On the other hand, when the best-fit corresponds to a black body $(\beta =0)\hat{R}=R_{\text{dust}}$. Due to the uncertainties and assumptions involved, the masses must be considered just as a rough approximation.

Stellar wind

This component was modelled by a stationary, isotropic, and isothermal stellar wind, expanding at constant velocity (Panagia & Felli 1975; Wright & Barlow 1975). In the optically thick regime, the flux density of such a wind at the frequency ν is: $$\begin{array}{rl}{S}_{\nu }(\text{thick})& =5.12{\left(\frac{\nu }{10\mathrm{G}\mathrm{H}\mathrm{z}}\right)}^{0.6}{\left(\frac{T_{\mathrm{e}}}{{10}^4\mathrm{K}}\right)}^{0.1}{\left(\frac{\dot{M}}{{10}^{-5}{\mathrm{M}}_{\odot }/\mathrm{y}\mathrm{r}}\right)}^{4/3}\\ & \times {\left(\frac{\mu }{1.2}\right)}^{-4/3}{\left(\frac{{\nu }_{\mathrm{\infty }}}{{10}^3\mathrm{k}\mathrm{m}/\mathrm{s}}\right)}^{-4/3}(\overline{Z}{)}^{2/3}{\left(\frac{d}{\mathrm{k}\mathrm{p}\mathrm{c}}\right)}^{-2}\mathrm{m}\mathrm{J}\mathrm{y}\end{array}$$(B.8)

where T_e is the electronic temperature, $\dot{M}$ the mass-loss rate, μ the mean molecular weight (assumed 1.2), v_∞ the expansion velocity of the wind, and $\overline{Z}$ the average ionic charge (assumed 1). The optically thick assumption holds up to a certain turnover frequency ν_c, equal to $${\nu }_c=4.7\times {10}^{42}\frac{{\gamma }^{1/2}{g}_{\mathrm{f}\mathrm{f}}^{1/2}\overline{Z}}{T_{\mathrm{e}}^{3/4}{R}_{\mathrm{c}}^{3/2}}\frac{\dot{M}}{\mu v_{\mathrm{\infty }}}\mathrm{H}\mathrm{z}$$(B.9)

with γ being the number of electrons per ion and g_ff the free-free Gaunt factor, approximated as Leitherer & Robert (1991): $$g_{\mathrm{f}\mathrm{f}}=9.77(1+0.13\log \frac{T_{\mathrm{e}}^{3/2}}{\overline{Z}\nu })$$(B.10)

and R_c denoting the inner boundary radius of the ionised envelope (here assumed equal to R_star). Above ν_c the spectrum begins to flatten, approaching that of an optically thin H II region with α = −0.1. In this regime, the flux density results: $$S_{\nu }(\text{thin})=S_{{\nu }_c}(\nu /{\nu }_c{)}^{-0.1}$$(B.11)

We may note in Eq. (B.8) the weak dependence of the flux density with T_e. In any case, we tested different models by varying T_e from 10³ to 10⁵ K and noted that the changes in the flux were always smaller than the uncertainties of the observed fluxes. We therefore adopted T_e = 10⁴ K to avoid an unnecessary free parameter. Finally, there is a degeneracy in $\dot{M}$ and v_∞ because the flux density varies with ($\dot{M}$/v_∞)^4/3. For simplicity, we adopted v_∞ = 100 km s⁻¹ and leave $\dot{M}$ as the only free parameter.

There is a circular dependence between ν_c and g_ff because ν_c is frequency-dependent through g_ff. It is therefore necessary to adopt a fixed value of g_ff. The Gaunt factor g_ff varies roughly from 2 to 6 within its validity range (from centimetre to millimetre wavelengths), while ν_c falls in the mid-IR range. Consequently, we adopted a characteristic g_ff = 4.08, corresponding to a frequency of 30 GHz, for the calculation of ν_c. This causes an uncertainty in the turnover frequency up to a factor of 2, which has no significant impact on the millimetre/centimetre fluxes.

Bremsstrahlung

For this component, we modelled the free-free emission from an ionised, static, and isothermal H II region characterised by a radius R_ff. In the optically thin regime (Mezger & Henderson 1967), the flux density behaves as: $$S_{\nu }(\text{thin})=2k{T}_{\mathrm{e}}{\left(\frac{\nu }{c}\right)}^2\tau {\mathrm{\Omega }}_{\mathrm{f}\mathrm{f}}$$(B.12)

Ω_ff = π (R_ff/d)² is the solid angle of the ionised region and τ is the optical depth, given by $$\tau =3.28\times {10}^{-7}{\left(\frac{T_{\mathrm{e}}}{{10}^4\mathrm{K}}\right)}^{-1.35}{\left(\frac{\nu }{\mathrm{G}\mathrm{H}\mathrm{z}}\right)}^{-2.1}\left(\frac{EM}{\mathrm{p}\mathrm{c}{\mathrm{c}\mathrm{m}}^{-6}}\right)$$(B.13)

where EM represents the emission measure of the ionised region.

Solving Eq. (B.13) for τ = 1, it is possible to compute the turnover frequency ν_c. At frequencies below ν_c, the H II region enters in the optically thick regime, the spectrum steepens, and it follows the Rayleigh-Jeans approximation (S_ν ∝ ν²).

All Tables

All Figures

	Fig. 1 Noise map of the observed fields. The rms noise was computed in concentric rings with a width of 2 arcsec. The noise in the 1.15 mm maps is plotted in blue, and the noise corresponding to the 2 mm maps is plotted in red. The field names (abridged) are indicated near the 1 mm curves.
In the text

Fig. 2 Sky maps of the millimetre continuum emission in the four fields observed with NIKA2. The maps at 1.15 mm are shown at the top row, and the 2 mm maps are shown at the bottom. The field that contains HD 168607 and HD 168625, labelled HD168*, is centred on HD 168625. The colour scales (indicated at the right of each row) are in Jy beam−1. HPBWs are drawn near the bottom left corners of the HD168* images. The central areas indicate the zoomed region displayed in Fig. 3.

In the text

Fig. 3 Zoomed sky maps around the central sources. The regions correspond to the regions indicated in Fig. 2. To span the full dynamic range, the colour scale is different for each of the maps. The contours correspond to 30, 50 (highlighted in white), and 70% of the maximum value. HPBWs are indicated near the bottom of the HD168* images. The contours at 50% are highlighted to facilitate a visual comparison with the beam and give an idea about the point-like nature of the central sources.

In the text

	Fig. 4 Spectral index maps between the two NIKA2 frequencies. The plotted field and marks are the same as in Fig. 2. For each map, all pixels with brightness below 2σ were masked.
In the text

	Fig. 5 Spectral energy distribution of HD 168607 in logarithmic scale. The catalogues are indicated at the top right. The corresponding error bars of all points are plotted. The NIKA2 points are shown with blue stars.
In the text

	Fig. 6 Spectral energy distribution of HD 168625.
In the text

	Fig. 7 Spectral energy distribution of [GKF2010]MN87.
In the text

	Fig. 8 Spectral energy distribution of [GKF2010]MN101.
In the text

	Fig. 9 Spectral energy distribution of G79.29+0.46.
In the text

Fig. 10 Best-fit models. The target names are indicated in the top left corners of each panel. Individual components (star, hot dust, warm or cold dust, and stellar wind or Bremsstrahlung) are indicated by the blue, orange, red, and green lines, respectively. The combination of the models (best fit) are shown in bold black lines. The grey area encloses the combined uncertainties of the fitting. The magenta dots with error bars depict the de-reddened photometry. The four charts are depicted with the same wavelength and flux density scales to facilitate the comparison among the targets.

In the text

	Fig. A.1 Herschel/PACS images around HD168607 at 70 μm (left) and 100 μm (right). Equatorial coordinates are relative to the star position. The continuum emission at these wavelengths is reported for the first time in this work.
In the text