© The Authors 2026

1. Introduction

Mass loss from massive stars governs the chemical enrichment of the interstellar medium. Post-red supergiants (post-RSGs) are prime laboratories in which to study the morphology, composition, and kinematics of recently ejected material, and to reconstruct mass-loss histories, since they can undergo very strong outflows (up to ∼ 10⁻⁴ M_⊙ yr⁻¹; Mauron & Josselin 2011). The class is small: only a few well-studied examples exist, such as IRC +10420 (Koumpia et al. 2022) and the ‘Fried Egg’ nebula (Lagadec et al. 2011; Koumpia et al. 2020). For a recent overview of post-RSGs and their circumstellar environments, see Gordon & Humphreys (2019); and for RSGs, their mass loss, and subsequent evolution, see van Loon (2025).

The complex nebula shapes observed around such objects can be caused by several factors, such as magnetic fields or rotation, but it is often thought that the presence of one or more companions is the main source of the angular momentum excess necessary for their shaping (Balick & Frank 2002). Observational studies and models have shown that a large proportion of stars – especially solar-type and massive stars – reside in binary or multiple systems (Whitworth & Lomax 2015; Fuhrmann et al. 2017; Merle 2024). Companions must have an impact on the mass-loss compositions and rate from massive stars. Studying the mass loss history of post-RSGs and the influence of binaries is key to understanding matter ejection from massive stars before the supernova explosion.

AFGL 4106 provides an excellent case study for investigating how binarity influences mass loss in massive stars. It was examined in the late 1990s by Molster et al. (1999) and van Loon et al. (1999), both of which reached similar conclusions: the system is a binary, consisting of an evolved star – likely a post-RSG – and a less evolved companion, embedded within a circumbinary nebula produced by one or both stars. The primary (star 1) is estimated to have a temperature of 7250 ± 250 K, while the companion (star 2) has a temperature of 3750 ± 250 K. The binary system was first resolved in 2006 using Hubble Space Telescope (HST) observations, leading to a separation measurement of 0.3″with a position angle (PA) of 270° (Bobrowsky et al. 2006).

The Spectro-Polarimetric High-contrast Exoplanet REsearch (SPHERE) on the Very Large Telescope provides an exquisite ∼20 milliarcsec angular resolution in the optical and ∼50 milliarcsec in the near-infrared, enabling us to better characterise the central stars of the system. SPHERE’s coronagraph and polarimetric capabilities enable a deep and detailed study of the nebula around the central binary system, allowing, for the first time, a study of the mass-loss history of a resolved binary system including a post-RSG star.

In this paper, we present a detailed analysis of the evolved massive binary system AFGL 4106; it is based on high-angular-resolution, high-contrast observations obtained with VLT/SPHERE in both direct and polarimetric imaging modes across multiple optical and near-infrared wavelengths. These data reveal the morphology and scattering properties of the surrounding dusty nebula, putting constraints on the mass-loss history of the system.

This paper is structured as follows. In Section 2, we describe the SPHERE observations and data reduction. Section 3 presents the morphological and photometric analysis, including the spectral energy distributions of each component; and in it we discuss the radiative-transfer modelling used to constrain the stellar and dust parameters. In Section 4, we interpret the results in the context of stellar evolution, mass-loss history, and nebular shaping mechanisms. Finally, our conclusions are summarised in Section 5.

2. Observations and data reductions

2.1. SPHERE observations

Observations were carried out with the extreme adaptive optics imager and polarimeter SPHERE (Beuzit et al. 2019) installed at the Unit Telescope 3 Nasmyth focus of the Very Large Telescope (VLT). Its adaptive optics module SAXO (Fusco et al. 2006) provides a correction from the atmosphere turbulence to reach the highest angular resolution possible. The main component of the adaptive optics system is the Shack-Hartmann wave-front sensor (WFS), with 41 × 41 actuators working at a frequency of 1380 Hz on the deformable mirror.

We obtained intensity images of AFGL 4106 using the SPHERE instruments ZIMPOL (optical, 0.5–0.9 μm) and IRDIS (near-infrared, 0.95–2.32 μm). To enhance contrast in the ZIMPOL data, we acquired polarimetric maps, as dust-scattered light is polarised, while stellar emission is not, revealing the circumstellar environment at flux levels several orders of magnitude fainter than the stars. For selected optical filters, we combined polarimetry with a 77.5 mas coronagraph to further suppress the primary star’s intensity and improve contrast. Additional data were obtained without a coronagraph to probe the innermost regions of the system. IRDIS observations were performed in classical imaging mode, since no coronagraph was large enough to mask both stars; we increased the exposure time and allowed the central region to saturate in some IRDIS frames to better detect the nebula. A summary of the SPHERE observations is provided in Table A.1.

Both ZIMPOL and IRDIS data were reduced using the SPHERE data-reduction centre (Delorme et al. 2017), which makes use of the Zurich ZIMPOL pipeline (Schmid et al. 2018) and the SPHERE pipeline (Pavlov et al. 2008) for the basic reduction steps. The SPHERE data centre makes use of custom-designed python routines to provide the ZIMPOL polarisation map. In polarisation mode, ZIMPOL measures the Stokes parameters Q and U; the polarisation maps are then computed using the following equations:

$$\begin{array}{c}\hfill p_L=\sqrt{q^2+u^2};I=\frac{I_Q+I_U}{2};\theta =arctan\left(\frac{q}{u}\right)\end{array}$$(1)

$$\begin{array}{c}\hfill \text{with:}q=\frac{Q_{+}+Q_{-}}{2I_Q};u=\frac{U_{+}+U_{-}}{2I_U},\end{array}$$(2)

and where p_L is the polarized intensity, I the total intensity, θ the polarization angle, q and u the normalized Stokes parameters, I_Q and I_U the total intensity measured in the Q and U cycles, respectively, and Q₊, Q₋, U₊ and U₋ the Stokes parameters measured in each half cycle.

2.2. Auxiliary photometry

The spectral energy distribution (SED) of the system is affected by the properties of the two stars and the interaction of their light with dust around them. To construct the SED of AFGL 4106, we retrieved archival multi-wavelength photometric data covering both the optical and infrared ranges. This broad spectral coverage allows us to characterise the contributions from both the stellar components and the surrounding dust. Additionally, we incorporated ISO-SWS spectra (Sloan et al. 2003) to complement the infrared portion of the SED. The spectra were integrated over ISO filter response curves to extract corresponding photometric points. The data we used are presented in Table A.2.

3. Analysis

3.1. The central stars

From Hubble Space Telescope (HST) imaging, the separation between the two central stars has been estimated as 0.3″(Bobrowsky et al. 2006) with a position angle of ∼270° for the secondary. However, the theoretical diffraction limit (λ/D) of HST at 1.644 μm is 0.141″. This implies that the separation of 0.3″is just above the Nyquist limit. With our SPHERE ZIMPOL images, we reached a theoretical diffraction limit of 0.014″in the V band Fig. 1, which is about one order of magnitude better than HST resolution. This allows us to obtain a better estimation of the separation. We performed a photocentre fit on each components in the V-band image in order to measure the separation. For this, we used the iteratively-weighted centre-of-gravity (IWCoG) technique (Vyas et al. 2009), where the photocentre is defined as

Fig. 1. ZIMPOL polarised intensity map of AFGL 4106 in filter V.

$$\begin{array}{c}\hfill ({\text{x}}_{\text{c}},{\text{y}}_{\text{c}})=\frac{{\sum }_{\mathit{ij}}{\text{X}}_{\mathit{ij}}{I}_{\mathit{ij}}{W}_{\mathit{ij}}}{{\sum }_{\mathit{ij}}{I}_{\mathit{ij}}{W}_{\mathit{ij}}}0<i,j<1023,\end{array}$$(3)

where (x_c, y_c) are the photocentre coordinates, X is the position of the pixel, I the value associated with the pixel, and W a Gaussian weight, calculated using a radial profile of the star. We finally find a separation of 0.272 ± 0.017 arcsec with a position angle of ∼90°. At a distance of 3.19${}_{-0.19}^{+0.45}$ kpc (see Sect. 4.1), this corresponds to a physical separation of 865 ± 102 au. Our measurement differs from that reported by Bobrowsky et al. (2006). This discrepancy arises from the choice of the reference star, which introduces a 180 degree° ambiguity in the position angle. They indeed used star 2 as a reference, while we used star 1 (the one in the east in all images).

Fig. 2. Three-colour image of AFGL 4106 in Band I_PRIM, J and H. ZIMPOL image is re-scaled to the IRDIS FoV. The colour and contrast are enhanced differently in the inner and outer parts of the nebula. The red arrow indicates the direction of the proper motion of the system.

3.2. The nebula

3.2.1. ZIMPOL optical-imaging polarimetry

With its narrow field of view (∼3.5″ × 3.5″) and high angular resolution (about 20 mas, Schmid et al. (2018)), observations with the ZIMPOL instrument reveal the morphology of the inner part of the nebula. To enhance the contrast, we used a coronagraph (seen as a black circle with spiders) and polarimetric measurements. The light from the star is unpolarised, while dust scatters this light. Images in polarised light enabled us to partially mask the central star and reveal the dust around it. Figure 3 presents the polarised-intensity image obtained with the ZIMPOL I_PRIM filter (central wavelength: 789.7 nm).

Fig. 3. ZIMPOL polarised-intensity map of AFGL 4106 in filter I_PRIM.

The central circular shape, with a diameter of ∼1 arcsec, is the adaptive optics correction halo and, therefore, not physical. The ZIMPOL polarised-intensity image Fig. 3 highlights two prominent cavities, which are referred to as holes 1 and 2. These cavities, each measuring approximately 2″ in length and 1″ in width, are positioned symmetrically on either side of a bright northeast-southwest axis.

In the polarised-intensity images, the two stars appear embedded in a non-uniform cavity of ∼2 arcsec in diameter, including the two voids described before. The inner edges of this cavity are clearly non-spherical, and a dichotomy is observed between the northwest rim, where four loops of ∼0.3 arcsec in diameter are observed, and the southeast displays five comet-like features pointing towards the centre. The external parts of the observed nebula are more uniform, albeit still clumpy.

3.2.2. IRDIS near-infrared imaging

The IRDIS intensity images reveal the outer parts of the nebula. As no coronagraphs were large enough to mask the central stars, the central stars were heavily saturated in some sequences to reveal the nebula.

The nebula has a diameter of approximately 8.6″. It exhibits a variety of intricate structures at multiple scales. Inside the nebula, the saturated stars appear has black dots. The two double lines crossing the stars represent the spider of M2. The ∼1.5 arcsec diameter circle is the AO correction halo. Two filaments are observed on the southeast and west sides of the nebula, about 2 arcsec from the centre, clearly visible in Fig. 4. Another one is present in the northeast, about 3 arcsec away from the centre. Finally one is present in the northwest about 2 arcsec away from the centre. The gaps between the filaments are not empty, but rather clumpy.

Fig. 4. IRDIS intensity map of AFGL 4106 in filter J.

On a larger scale, in the IRDIS images we observe sharp, triangle-shaped lobes exiting the central star radially. They are observed at various distances between 1.5 and 2.0 arcsec from the central star. They appear to be sharper towards the south and west of the nebula.

They are similar to the lobes observed on CPD-568032 (Chesneau et al. 2006) or in the starfish nebula (Sahai et al. 2005) and are also similar to the sporadic, clumpy ejecta imaged around VY CMa (Humphreys et al. 2021). Mostly visible in filter J Fig. 4 and continuum J Fig. B.7, these lobes are not visible at larger wavelengths. The most prominent structure is the southern one denoted as S lobe on Fig. 2, which is ∼2″ long. The other lobes are smaller, and all are ∼0.5″ long.

Throughout the nebula, we observe filamentary structures that likely arise from interactions between the circumstellar material and the interstellar medium (ISM). van Loon et al. (1999) reported similar filaments on larger scales and interpreted them as bow shocks produced by the motion of AFGL 4106 through the ISM. Using Gaia DR3 astrometry (Gaia Collaboration 2023), we derive a proper motion of μ_α = −6.68 ± 0.33 mas yr⁻¹ and μ_δ = −1.46 ± 0.35 mas yr⁻¹, corresponding to a total motion of 6.84 mas yr⁻¹ with a position angle of 257°. This direction is indicated by the red arrow in Fig. 2. However, the fact that the filaments are not spatially aligned with the proper-motion direction disfavours the bow-shock interpretation. Instead, these features are more likely the relics of past mass-ejection events from the system, subsequently shaped by the interaction with the surrounding ISM.

The iso-intensity map of the nebula displays variations in the shape of the ejecta Fig. 5a. Looking up at the eccentricity of these iso-intensity curve, $e=\sqrt{1-\frac{a^2}{b^2}}$, where a is the semi-major axis and b the semi-minor axis, and e seems to change with distance to the central star. We observe oscillations with a pseudo-period of ∼2.560 arcsec, given that we only have 1.5 periods. The gap observed between 1.8 and 2.2 arcseconds corresponds to a sharp variation of ellipticity, likely caused by a rotation of the elliptical contours by π/2, resulting in a swap between the semi-major and semi-minor axes. This suggest a change in the ejecta preferential direction across time.

Fig. 5. Zoomed-in view of intensity contours of AFGL 4106 nebula in J band is shown in panel a, along with the corresponding eccentricity of these contours in panel b. The top panels present the raw data, while the bottom panels show the residuals after subtracting a linear fit.

Assuming a typical stellar-wind velocity of 20–30 km/s (Molster et al. 1999) and a distance of 3.19 kpc (see discussion in Sect. 4.1), we estimate oscillation periods ranging from approximately 1300 to 1950 years, depending on the adopted wind speed.

3.3. Spectral energy distribution of both stars

Using an ISO-SWS spectrum covering the 2.36–45.3 μm range, Molster et al. (1999) previously estimated the temperatures of the system’s components through SED fitting, finding a temperature of 7250 ± 250 K for the primary and 3750 ± 250 K for the secondary, concluding that the binary system consists of an A–F and an M type star of similar luminosity.

Using the first resolved observations of both components, we performed aperture photometry on our SPHERE data to isolate the SED of each star in the system across various optical and near-infrared wavelengths. Flux ratios were derived from the SPHERE data by integrating the flux within circular apertures of radius 1.22 ⋅ λ/D, corresponding to the diffraction-limited core and ensuring that the full dispersed flux is captured. Among the 16 available filters, only five provide a suitable combination of PSF quality and unsaturated stellar images, allowing us to perform aperture photometry on both stars.

Since SPHERE data are not flux-calibrated, the total flux for each star was obtained from photometric auxiliary data. These fluxes were corrected from the interstellar reddening using the PyAstronomy unred routine (Czesla et al. 2019), making use of the Fitzpatrick (1999) parametrisation to estimate the wavelength dependence of interstellar extinction. This parametrisation uses spectral observation to derive the shape of the extinction curve, fitted to the following formula:

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill \frac{\text{E}(\lambda -\text{V})}{\text{E}(\text{B}-\text{V})}& ={\text{C}}_2+{\text{C}}_1{\lambda }^{-1}+\frac{{\text{C}}_3}{{[{\lambda }^{-1}-\frac{{\lambda }_0^{-2}}{{\lambda }^{-1}}]}^2+{\gamma }^2}\hfill \\ \hfill & +{\text{C}}_4[0.539{({\lambda }^{-1}-5.9)}^2+0.0564{({\lambda }^{-1}-5.9)}^3]\hfill \end{array},\end{array}$$(4)

where E(λ–V) is the extinction at wavelength λ, E(B–V) the extinction in the B and V bands, and

$$\begin{array}{c}\hfill \{\begin{array}{cc}\hfill {\text{C}}_1& =-0.824+4.717\times R^{-1}\hfill \\ \hfill {\text{C}}_2& =2.030+3.007\times {\text{C}}_1\hfill \\ \hfill {\text{C}}_3& =3.23\hfill \\ \hfill {\text{C}}_4& =0.41\hfill \\ \hfill {\lambda }_0^{-1}& =4.595\mathrm{\mu }{\text{m}}^{-1}\hfill \\ \hfill \gamma & =1.051\mathrm{\mu }{\text{m}}^{-1}\hfill \end{array}\end{array}$$(5)

for the Milky Way.

The interstellar extinction value derived by Molster et al. (1999), E(B−V) = 1.07 mag, was compared with the value derived by de Laverny (in preparation, private communication). This latter estimate was obtained using Gaia/GSP-Spec atmospheric parameters (Recio-Blanco et al. 2023) to compare observed and theoretical stellar colours in the Gaia bands and thus determine the total reddening along the line of sight. For AFGL 4106, GSP-Spec reports T_eff = 6019 K, log g = 1.75, and [M/H] = −0.57 dex, adopting the calibrations as a function of T_eff described in Recio-Blanco et al. (2023, 2024). With these parameters, the expected intrinsic colour is (B_P−R_P)₀ = 0.65 ± 0.013 mag, while the observed Gaia/DR3 colour is (B_P−R_P)_obs = 2.46 ± 0.01 mag, yielding (B_P−R_P) = 1.809 ± 0.013 mag.

It is worth noting that that GSP-Spec has analysed a Gaia/RVS spectrum that is actually a combined spectrum of the unresolved binary system, which may slightly affect the atmospheric parameter determination. Moreover, as AFGL 4106 is likely variable, the derived extinction is estimated from an average spectrum obtained over the ∼2.5 years of Gaia/DR3 observations. Using our adopted stellar temperature (T_eff = 6723 K, see Sect. 3.4), the expected intrinsic colour becomes (B_P−R_P)₀ ≃ 0.5 mag, leading to a slightly higher reddening, E(B_P−R_P) ≈ 1.95. This difference remains small and does not significantly affect the other derived parameters, such as luminosity, but should be mentioned.

To facilitate comparison with the extinction from Molster et al. (1999), we converted E(B_P−R_P) to E(B−V) using the extinction coefficients A(B_P)/A(V) and A(R_P)/A(V) from the ESA/Gaia photometric system documentation (Riello et al. 2021). Therefore, we have

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill E({\mathrm{B}}_{\mathrm{P}}-{\mathrm{R}}_{\mathrm{P}})& =({\mathrm{B}}_{\mathrm{P}}-{\mathrm{R}}_{\mathrm{P}})-{({\mathrm{B}}_{\mathrm{P}}-{\mathrm{R}}_{\mathrm{P}})}_0\hfill \\ \hfill & =A({\mathrm{B}}_{\mathrm{P}})-A({\mathrm{R}}_{\mathrm{P}})\hfill \\ \hfill & =(\frac{A({\mathrm{B}}_{\mathrm{P}})}{A(V)}-\frac{A({\mathrm{R}}_{\mathrm{P}})}{A(V)})\times R_V\times E(B-V)\hfill \\ \hfill & \simeq 1.392\times E(B-V)\hfill \end{array}.\end{array}$$(6)

Adopting the standard value of R_V = 3.1 for the Milky Way interstellar medium, the total extinction reported by Gaia corresponds to E(B − V) = 1.30 ± 0.01 mag. This value is consistent with that reported by Molster et al. (1999), considering that they found an extinction due to the circumstellar material of E(B − V) = 0.22 ± 0.05 mag. This leads to a total extinction in agreement with the one used by Molster et al. (1999). In the following analysis, we adopt the same interstellar reddening: E(B–V) = 1.07 mag. The measurements of the photometry are presented in Table 1 and plotted in Fig. 6.

Fig. 6. SED of both components of system, as calculated by extracting photometry from the two stars, resolved in our SPHERE images. Larger uncertainties on the CntJ filter are due to the strong influence of the surrounding nebula.

For the first time, we resolved the binary system in the core of AFGL 4106 at multiple wavelengths in the optical and near-infrared, allowing us to build SEDs for both objects. The five photometric data points clearly indicate that the primary star is hotter and more luminous than its companion. To have a better estimation of the stellar properties, we performed radiative-transfer modelling of the system and its surrounding dust shell, as described in the next section.

3.4. Star properties

Figure 6 give us an overview of the SED of each component of the system in five bands. We used the 1D radiative transfer code DUSTY (Version 4, Ivezic & Elitzur (1997)) to model the SED of the system, adopting the dust composition identified by Molster et al. (1999). The model aims to fit the photometric data obtained from VizieR¹ as well as photometric points constructed by convolving the ISO-SWS spectrum with the ISO filter profiles Table A.2. We used the standard MRN dust size distribution (Mathis et al. 1977) and opted for the solution for radiatively driven winds. The best fit is determined with a minimisation of the χ² using a Markov chain Monte Carlo method. We created synthetic fluxes by convolving the model spectrum with the filter profiles of each observational instrument. We then used the χ² defined by Yang et al. (2023):

$$\begin{array}{c}\hfill {\chi }^2=\frac{1}{\mathrm{N}-\mathrm{p}-1}{\displaystyle \sum \frac{{[1-f(\lambda ,\mathrm{model})/f(\lambda ,\mathrm{Obs})]}^2}{f(\lambda ,\mathrm{model})/f(\lambda ,\mathrm{Obs})}},\end{array}$$(7)

where f(λ) = F(λ)/F(K) is the flux at a specific wavelength, λ, normalised to the flux in the K band; N is the number of photometric data points (5 in the case of Table 1 and 47 in the case of Table A.2); and p is the number of free parameters (5 to 7). This modified χ² balances the contributions from both optical and infrared fluxes by effectively flattening the spectrum. In contrast, a standard χ² tends to underweight the infrared data due to the steep decline in flux at longer wavelengths of the Planck function.

We initially used the SEDs of each component, as shown in Fig. 6, to estimate the stellar temperatures through a first MCMC run. Based on these results, we then constrained the temperatures to narrow intervals centred on the initial estimates. Using the full photometric dataset available from VizieR, we subsequently refined the model to determine the properties of the circumstellar shell and the relative luminosities of the two components in a second MCMC run². The final results are presented in Fig. 7 and the model parameters in Table. 2.

Fig. 7. SED fitting using DUSTY version 4.0 at a distance of 3.19 kpc. The orange diamonds show the photometric measurements, and the black line is the best-fit model consisting of the attenuated flux (dashed light blue), the scattered flux (dot-dashed grey), and dust emission (dotted brown).

With such luminosities and temperatures, those stars have radii of R₁ = (209 ± 12) R_⊙ and R₂ = (567 ± 99) R_⊙ at 3.19 kpc. This indicates that these stars are both supergiants in two different stages of evolution. The primary seems to be the more massive, and the nebula was formed when it expelled its outer layers, leaving the star in the post-RSG phase. The companion is likely in a less evolved stage, with a radius similar to giant stars, and may be surrounded by a possible expanding shell caused by wind at the surface.

DUSTY computes both the inner radius of the dust shell and an estimate of the mass-loss rate. For our system, the inner radius is calculated as r_int = 6.46 × 10¹⁴ m, corresponding to 1.35 arcsec at a distance of 3.19 kpc. This radius is slightly less than the 1.6–1.9 arcsec derived by van Loon et al. (1999). To estimate the mass-loss rate, we applied the rescaling relation for dust-driven winds described by van Loon (2000):

$$\begin{array}{c}\hfill \dot{M}={\dot{M}}_{\mathrm{DUSTY}}{\left(\frac{L}{{10}^4{L}_{\odot }}\right)}^{3/4}{\left(\frac{r_{\mathit{gd}}}{200}\frac{{\rho }_d}{3}\right)}^{1/2},\end{array}$$(8)

where r_gd = 100 is the gas-to-dust ratio assumed by Molster et al. (1999) and ρ_d = 3.3 g cm⁻³ represents the bulk density of the dust. Using these parameters, we derive a total mass-loss rate of $\dot{M}\simeq 4\times {10}^{-3}{\left(\frac{d}{3.19\mathrm{kpc}}\right)}^{1.5}{M}_{\odot }{\mathrm{yr}}^{-1}.$

4. Discussion

4.1. Distance to AFGL 4106

The distance to AFGL 4106 is debated. Gaia’s parallax gives a distance of 0.94 ± 0.26 kpc (Gaia Collaboration 2016, 2023), but when applying the Bailer-Jones method (Bailer-Jones et al. 2021), using a probabilistic approach to estimate stellar distances and a prior constructed from a three-dimensional model of our Galaxy, we find a distance of 1.47${}_{-0.4}^{+0.59}$ kpc. Due to its proximity with the Carina nebula Molster et al. (1999) assumed that AFGL 4106 was at a distance of 2.61 kpc

To improve the distance estimation, we combined the interstellar extinction measurement from Molster et al. (1999) with the value derived from GSP-Spec data (de Laverny, in preparation, private communication), and compared them with the extinction maps from G-TOMO (Lallement et al. 2022; Vergely et al. 2022). G-TOMO delivers three-dimensional extinction maps within a 5 × 5 × 0.8 pc region centred on the Sun, using Gaia EDR3 and 2MASS photometry and parallaxes. This provides access to distance-dependent extinction curves in the solar neighbourhood, which can then be inverted: given the extinction, we can estimate a star’s distance. Taking an extinction value of E(B − V) = 1.07 ± 0.2 mag, we derive a distance of 3.19${}_{-0.19}^{+0.45}$ kpc.

As these are massive stars, their lifetime is particularly short – around 11 Myr for a 15 M_⊙ star – so they are necessarily located close to the cluster where they were formed. The two closest star forming regions are the Carina nebula and the open cluster Westerlund 2, which are, respectively, at a projected distance of ∼2.1° and ∼1.8° of AFGL 4106 and at a distance of 2.61 kpc and 4.13 kpc. In the Carina reference frame, the proper motion is oriented in the direction of the south, with an absolute value of 4 mas yr⁻¹ (Gaia Collaboration 2023). However, the proper motion of AFGL 4106 is not consistent with an origin in the Carina Nebula, as shown in Fig. 8. When re-referenced to the Carina reference frame, the proper motion vector does not point towards the nebula, indicating a different origin or kinematic history. In the Westerlund referential frame, the proper motion is also directed towards the south, with an absolute value of 4.5 mas yr⁻¹. However, in this frame, the motion appears to be moving directly away from the cluster.

Fig. 8. Relative proper motion of AFGL 4106 with the Carina (red arrow) and the open cluster Westerlund 2 (black arrow).

4.2. Radial profiles

The radial profiles of both stars were computed from ZIMPOL and IRDIS images, and the full width at half maximum (FWHM) was derived for each filter Fig. 9. The measured FWHM values were then compared to the diffraction limit, λ/D (black line), and to the IRDIS and ZIMPOL-modelled point spread function (PSF) generated using the Maoppy code (Fétick et al. 2019) (blue line).

Fig. 9. Radial profile of both stars of AFGL 4106 at different wavelengths.

In each filter, both stars appear resolved. However, despite the distance to AFGL 4106, the photosphere of the companion (R₂ = 567 × $\frac{\mathrm{d}}{3.19\text{kpc}}$ R_⊙ ∼0.82 mas) cannot be resolved with our observations. This indicates that the resolved structure is likely a circumstellar shell. While the exact size of this shell is challenging to determine, the FWHM at longer wavelengths tends to align with the modelled PSF. This suggests a shell size of approximately 65 mas. At a distance of 3.19kpc, this corresponds to an envelope size of roughly 208 au (∼80 R₂).

4.3. Shaping and mass

4.3.1. HR diagram

Based on the SED fitting results, we placed both components of the system on a Hertzsprung-Russell diagram Fig. 10. Their positions were then compared with the evolutionary tracks from Ekström et al. (2012), considering both rotating and non-rotating stellar models. At a distance of 3.19 kpc, both stars align well with the theoretical tracks.

Fig. 10. HR diagram of two components of AFGL 4106 at 3.19 kpc. Evolutionary tracks are from Ekström et al. (2012). Isoradius lines are overplotted with dotted lines.

The temperature and luminosity of the primary are consistent with those of a rotating star with an initial mass of approximately 15 M_⊙. According to the evolutionary tracks, such a star would be around 13.8 Myr old, with a current mass of 14.5 M_⊙. However, this inferred age is significantly older than the estimated age of the candidate star-forming region, which is only 2–3 Myr (Zeidler et al. 2015). This discrepancy may indicate either an incorrect assumption about the system’s origin or a deviation from standard evolutionary pathways due to binary interaction.

Similarly, the secondary’s temperature and luminosity are consistent with those of a non-rotating star with an initial mass of about 9 M_⊙, corresponding to an age of roughly 30 Myr and a current mass of 8.7 M_⊙. Once again, this age is incompatible with that of Westerlund 2 and, more importantly, it does not agree with the age inferred for the primary in the previous paragraph. This discrepancy further suggests that standard single-star evolutionary tracks may not adequately describe the system’s past.

We also placed both components on a Hertzsprung-Russell diagram Fig. 11 alongside luminous blue variables, yellow hypergiants, and RSGs from Oudmaijer et al. (2009).

Fig. 11. HR diagram of two components of AFGL 4106 at 3.19 kpc. Stars from Oudmaijer et al. (2009) are overplotted. The dotted line corresponds to the Humphreys-Davidson limit.

The secondary falls clearly within the RSG region, consistent with its classification. In contrast, the primary is located just below the so-called yellow void. It is too hot to be classified as an RSG and too faint to be considered a yellow hypergiant, placing it instead among the yellow supergiants (YSGs). Intermediate-mass stars can experience a ‘blue loop’ during their evolution, temporarily evolving from the RSG phase to warmer YSG or even blue-supergiant (BSG) phases. However, the primary’s mass (M > 12 M_⊙) is likely too high for this mechanism to occur under standard evolutionary models. The system’s binarity is a factor that may alter its evolutionary pathway. However, the binary is too wide to have affected the evolution of the primary. For stars to truly affect each other they should be close, at most within a few astronomical units (Sana et al. 2012). This leaves the possibility that AFGL 4106 experienced interaction with a third, close companion, or even a past merger, open, which might help explain a blue-loop-like excursion or its unusually low post-RSG luminosity.

4.3.2. Nebular shaping: Influence of the companion

We see in Sect. 4.1 that the nebula shows signs of oscillation (Fig. 5). These oscillations might come from the variation of the orbital speed of a companion, alongside the orbit curve (El Mellah et al. 2020). When the ratio η = v_wind/v_orb between the wind speed and orbital velocity is low, the ejecta is constrained in the orbital plane. Thus, depending on the angle of observation, the shape of this ejecta might be elliptical, whereas when η is larger, the ejecta is spherical, so the nebula will appear spherical for all lines of sight.

If the observed variation in ellipticity is attributed to the orbital motion of the binary, by applying Kepler’s third law with an oscillation period of approximately 1300–1950 years (depending on the assumed wind velocity of 20–30 km s⁻¹, Sect. 3.2.2), the resulting total system mass is estimated to be in the range of 175–400 M_⊙. However, this estimate is inconsistent with the stellar parameters derived in the previous section. Moreover, such high total mass is improbable, as the formation of a binary system comprising two stars each exceeding 80 M_⊙ is considered unlikely based on current star formation models.

The presence of a third object to explain the ellipticity oscillations is plausible but remains unsupported observationally: no such companion is detected in either the imaging or the spectroscopic data. Assuming an unresolved companion at the resolution limit of our observations (20 mas), Kepler’s law yields a total mass of only 0.06–0.15 M_⊙, far from the ∼15 M_⊙ estimated in Sect. 4.3.1.

van Loon et al. (1999) reported radial-velocity variations with a period of 4500 days, interpreted as the orbital motion of the known binary and implying stellar masses ≤15 sin⁻³i M_⊙, which is consistent with our results. If, instead, the radial velocity oscillation were attributed to an orbit involving a third companion located at 20 mas, the resulting total mass would be ∼1700 M_⊙, which is unphysical. Matching the mass estimated in Sect. 4.3.1 would require a separation smaller than ∼4.1 mas (i.e. ∼13 au at 3.19 kpc).

Such a close companion could, in principle, exist and might account for the radial-velocity variations reported by van Loon et al. (1999). However, it cannot explain the large-scale oscillations observed in the nebula, which would require a much longer orbital period. Therefore, while the presence of a close companion cannot be ruled out, it is unlikely to be responsible for the observed nebular morphology.

The shaping of the nebula is likely the result of a complex dynamic driven by, first, the influence of one or many companions and, second, by a complex mechanism of ejection. It is also important to note that, the semi minor axis is almost aligned with the proper motion of the system. The interaction with the interstellar medium (ISM) may be responsible for the peculiar morphology of the nebula, particularly its large-scale elliptical shape with an eccentricity of approximately 0.5, as the stellar wind is decelerated by the ISM through the formation of a ‘shock front’.

5. Conclusions

The main results of this study can be summarised as follows:

• AFGL 4106 is a binary system composed of a post-RSG and an RSG with a separation of 0.272 ± 0.017 arcsec and temperatures T₁ = 6723 ± 196 K and T₂ = 3394 ± 264 K.

• The primary star appears slightly brighter than its companion, with a bolometric-luminosity ratio of L₁/L₂ = 2.09, primarily attributed to its temperature being approximately two times higher.

• Both stars have high bolometric luminosities, with L₁ = (7.9 ± 0.18) ${(\frac{\mathrm{d}}{3.19\text{kpc}})}^2$ ×10⁴ L_⊙ and L₂ = (3.8 ± 0.11) ${(\frac{\mathrm{d}}{3.19\text{kpc}})}^2$ ×10⁴ L_⊙ with d in kpc.

• The radii R₁ = (209 ± 12)$(\frac{\mathrm{d}}{3.19\text{kpc}})$ R_⊙ and R₂ = (567 ± 99)$(\frac{\mathrm{d}}{3.19\text{kpc}})$ R_⊙ reveal that the stars are in two close, but different, stages of evolution: the primary in a post-RSG phase and the companion in a RSG phase.

• The stellar evolution model gives an estimation of 15 M_⊙ for the primary and 9 M_⊙ for the secondary.

• The proximity to the open cluster Westerlund 2 and the relative movement of both systems suggest that the star could have originated from the cluster. This gives an upper bound on the age of the stars ≲2–3 Myr.

• The binary is surrounded by a shell of ∼80 R₂.

• The complex morphology of the nebula reveals a non-isotropic and non-continuous matter ejection from the system, as well as a strong interaction with the ISM.

This study of AFGL 4106 will enable us to constrain the dust and evolutionary models in this type of system, to better understand the formation and evolution of such nebulae. Further calculation and hydrodynamic models, combined with observations of the dynamics of the nebula, will help us to study the physics of such systems and the final step of star evolution.

Acknowledgments

This work was supported by the French National Research Agency (ANR) through the PEPPER project (ANR-20-CE31-0002). I would like to thank the referee for their constructive comments, which helped improve the quality of this paper.

References

Appendix A: Observation and Data

Appendix B: SPHERE Images

Fig. B.1. ZIMPOL Polarized intensity map of AFGL 4106 in filter N_R.

Fig. B.2. IRDIS intensity map of AFGL 4106 in filter BrG.

Fig. B.3. IRDIS intensity map of AFGL 4106 in filter PaB.

Fig. B.4. IRDIS intensity map of AFGL 4106 in filter H.

Fig. B.5. IRDIS intensity map of AFGL 4106 in filter H2.

Fig. B.6. IRDIS intensity map of AFGL 4106 in filter HeI.

Fig. B.7. IRDIS intensity map of AFGL 4106 in filter CntJ.

Fig. B.8. IRDIS intensity map of AFGL 4106 in filter CntK1.

Fig. B.9. IRDIS intensity map of AFGL 4106 in filter CntK2.

Fig. B.10. IRDIS intensity map of AFGL 4106 in filter CO.

All Tables

All Figures

	Fig. 1. ZIMPOL polarised intensity map of AFGL 4106 in filter V.
In the text

	Fig. 2. Three-colour image of AFGL 4106 in Band I_PRIM, J and H. ZIMPOL image is re-scaled to the IRDIS FoV. The colour and contrast are enhanced differently in the inner and outer parts of the nebula. The red arrow indicates the direction of the proper motion of the system.
In the text

	Fig. 3. ZIMPOL polarised-intensity map of AFGL 4106 in filter I_PRIM.
In the text

	Fig. 4. IRDIS intensity map of AFGL 4106 in filter J.
In the text

	Fig. 5. Zoomed-in view of intensity contours of AFGL 4106 nebula in J band is shown in panel a, along with the corresponding eccentricity of these contours in panel b. The top panels present the raw data, while the bottom panels show the residuals after subtracting a linear fit.
In the text

	Fig. 6. SED of both components of system, as calculated by extracting photometry from the two stars, resolved in our SPHERE images. Larger uncertainties on the CntJ filter are due to the strong influence of the surrounding nebula.
In the text

	Fig. 7. SED fitting using DUSTY version 4.0 at a distance of 3.19 kpc. The orange diamonds show the photometric measurements, and the black line is the best-fit model consisting of the attenuated flux (dashed light blue), the scattered flux (dot-dashed grey), and dust emission (dotted brown).
In the text

	Fig. 8. Relative proper motion of AFGL 4106 with the Carina (red arrow) and the open cluster Westerlund 2 (black arrow).
In the text

	Fig. 9. Radial profile of both stars of AFGL 4106 at different wavelengths.
In the text

	Fig. 10. HR diagram of two components of AFGL 4106 at 3.19 kpc. Evolutionary tracks are from Ekström et al. (2012). Isoradius lines are overplotted with dotted lines.
In the text

	Fig. 11. HR diagram of two components of AFGL 4106 at 3.19 kpc. Stars from Oudmaijer et al. (2009) are overplotted. The dotted line corresponds to the Humphreys-Davidson limit.
In the text

	Fig. B.1. ZIMPOL Polarized intensity map of AFGL 4106 in filter N_R.
In the text

	Fig. B.2. IRDIS intensity map of AFGL 4106 in filter BrG.
In the text

	Fig. B.3. IRDIS intensity map of AFGL 4106 in filter PaB.
In the text

	Fig. B.4. IRDIS intensity map of AFGL 4106 in filter H.
In the text

	Fig. B.5. IRDIS intensity map of AFGL 4106 in filter H2.
In the text

	Fig. B.6. IRDIS intensity map of AFGL 4106 in filter HeI.
In the text

	Fig. B.7. IRDIS intensity map of AFGL 4106 in filter CntJ.
In the text

	Fig. B.8. IRDIS intensity map of AFGL 4106 in filter CntK1.
In the text

	Fig. B.9. IRDIS intensity map of AFGL 4106 in filter CntK2.
In the text

	Fig. B.10. IRDIS intensity map of AFGL 4106 in filter CO.
In the text