© ESO 2018

1 Introduction

The study of massive stars is important for many aspects of stellar evolution and cosmic enrichment. For example, a significant fraction of massive stars are part of in binary systems (Sana et al. 2013). The stellar components of these binaries are good candidates to generate gravitational waves by the merging of the resulting black holes or neutron stars.

Wolf-Rayet (WR) stars are hot stars with broad emission lines. These lines originate from a strong and optically thick wind, which is driven by radiation pressure and reaches a velocity of up to a thousand kilometers per second. This supersonic wind of WR stars exerts a major influence on the immediate surroundings of the star.

WR stars come in three groups: 1) cWR: classical, He-burning, H-poor remnants of massive stars; 2) WNLh: very luminous and massive main sequence late (cooler) stars with a hydrogen-enriched wind and a spectrum dominated by N III–V and He I–II lines; and 3) [WR]¹: the nuclei of some planetary nebulæ, that is, low to medium mass stars on their way to becoming white dwarfs (WDs). In this study, we focus on cWR and the final stage of massive stellar evolution.

Massive WR stars exist at the limit of exploding as a supernova of type SNIbc and are good candidates for long gamma-ray bursts (Dessart et al. 2017). WR stars also play a dominant role in the enrichment of the local interstellar medium (ISM) with their important mass-loss and their supernovae explosions. They are particularly one of the contributors of key chemical elements for planet formation (such as ¹⁴C or ²⁶ Al; Tatischeff et al. 2010).

WR stars are classified according to their emission spectra, which depends on their WR evolutionary stage. Those of class WN exhibit many emission lines of nitrogen and helium, while class WC stars manifest emission lines of carbon, oxygen, and helium. The standard accepted sequence of evolution places WC stars at the very end of the WR stage, before the very short oxygen-rich WO stage and the final supernova (Groh et al. 2014). The chronological evolution of WC stars is then classified into sub-stages from WC4 (youngest) to WC9 (oldest).

In addition, several WC8 and the majority of WC9 stars exhibit a strong infrared excess, which is reminiscent of hot and warm dust produced by the central source (Allen et al. 1972). These dust-producing WRs contribute to the enrichment of the ISM and have extremely high dust-formation rates, with values up to Ṁ = 10⁻⁶M_⊙ yr⁻¹ (Harries et al. 2004) corresponding to approximately 2% of the standard accepted wind mass-loss of late WC stars (Puls et al. 2008).

Several WR stars in binary systems with hot OB-type stars rich in oxygen and hydrogen are believed to harbour spiral, dusty structures called “pinwheel nebulae” (e.g. WR98a, Monnier et al. 1999; WR104, Tuthill et al. 1999; WR118, Millour et al. 2009b; WR48a, WR112, WR137, WR140, Marchenko & Moffat 2007, Lau et al. 2017; and the Quintuplet Cluster near the Galactic Centre, Tuthill et al. 2006). These dusty spirals in such massive binary systems are interpreted as the consequence of a collision between the WR wind and that of the hot, massive companion star, which can generate a violent shock interaction at the interface of the two stellar winds (Pittard 2009; Lamberts et al. 2012). This wind collision zone (WCZ) between the WR and the OB components offers ideal conditions to reach critical densities and enable dust nucleation, especially where the mixing between the carbon-rich wind of the WR star and H-rich wind of the OB star becomes significant (Hendrix et al. 2016). This newly-produced dust created within the shocked region then follows a ballistic trajectory forming a spiral pattern described in detail in the literature (Monnier et al. 1999; Tuthill et al. 2006, 2008; Millour et al. 2009b; Lau et al. 2017). This pattern follows an Archimedian spiral in the case of a circular orbit (as for WR104),or is more complex (producing arcs) when the inner binary orbit is elliptical (as can be seen e.g. in WR140 or WR48a).

Among permanent dust-makers, WR104 stands out as a nearly face-on (inclination angle i ≤ 16°) pinwheel nebula that was detected with the aperture-masking technique on the Keck telescope (Tuthill et al. 1999). After its discovery, the system was further investigated using the same aperture-masking technique, showing the rotation of the spiral and providing some details about the spiral properties (Tuthill et al. 2008). A radiative transfer model was developed by Harries et al. (2004) to explain the apparent flux saturation of the inner region of the spiral, detected by Tuthill et al. (1999). This saturation was interpreted as the presence of optically thick dust in the inner region and a dust formation region located at about 10 mas from the central binary star.

After 20 yr of observations based on reconstructed images from the aperture-masking technique and Fourier plane sampling, we present the first direct images of the WR104 system obtained with the Spectro-Polarimetric High-contrast Exoplanet Research (SPHERE) instrument at the Very Large Telescope (VLT) in Chile. The advantage of direct imaging is the reliability of the object’s brightness distribution.

Using the new SPHERE images with high dynamic range, our aims are to: 1) constrain and confirm the general characteristics of the system (spiral step, or radial separation between successive spiral coils, orientation, and so on), 2) provide an unprecedented view and physical description of the large scale of the pinwheel, up to 30 revolutions. Having access to these large spatial scales is essential to be able to define a set of physical models describing the dust distribution along the spiral arms and thus constrain the history of dust production over the last decade.

Our paper in organized as follows. In Sect. 2, we describe the SPHERE and VLT Imager and Spectrometer for mid-InfraRed (VISIR) observations and the related data reduction processes. Section 3 presents a first direct analysis of the images, followed in Sect. 4 by an analytical model of the dust emission in the spiral, complemented by a preliminary qualitative radiative transfer model. In Sect. 5, we interpret and discuss our results and finally, in Sect. 6, we conclude and put this work into perspective.

2 Observations and data processing

2.1 SPHERE observations

SPHERE is a high-performance imaging instrument equipped with an Extreme Adaptive Optics system (XAO) installed at the Nasmyth focus of Unit Telescope 3 of the Very Large Telescope (Beuzit et al. 2008). Our observations were obtained between April and July 2016 in several filters across the J, H, and K near-infrared (NIR) bands with the IRDIS infrared camera of SPHERE, (Infra-Red Dual-beam Imager and Spectrograph; Dohlen et al. 2008; Langlois et al. 2014). As detailed in Table 1, we took images in many narrow-band filters; a more detailed log of the data is presented in Table B.1. Concurrent with our dual-band IRDIS observations, we collected data using the Integral Field Spectrometer (IFS) of SPHERE covering the Y J bands (0.95–1.35 μm) with a spectral resolution of R = 50 (Claudi et al. 2008).

At mostepochs, for all filters, we observed WR104 together with a non-resolved star (either Sgr 4, HD 163955 of spectral type B9V, or TYC 6295-803-1), which we used to estimate the point spread function (PSF). We checked that these stars are unresolved by IRDIS (λ∕D≈ 30–35 mas, with λ the wavelength of incoming radiation and D the lens diameter, while the diameters are θ_UD ² = 0.40 ± 0.03 mas for 4Sgrand θ_UD = 0.40 ± 0.03 for TYC 6295-803-1) using the Searchcal/JMMC catalogue (Bonneau et al. 2006).

The observations required 5 h of telescope time and were taken without the SPHERE coronagraph. Due to the extreme brightness of WR104 in the H and K bands, an additional neutral-density filter was necessary in conjunction with most filters (with filter ND2, this yielded an attenuation of a factor 100 for the H and K bands, and with filter ND1, an attenuation of a factor of 10 for the J band).

All images were taken using the dithering technique, consisting of shifting the target on the detector between different integrations (a 4 × 4 pattern was used with a shift of one or two pixels). A total of 288 images were recorded for WR104 (with 4 × 4× 3 images for each filter) in classical imaging mode, and 96 images in dual-band-imaging (DBI mode), in which two different filters are used for two different parts of the camera.

The spectral response of the filters used for our multi-wavelength study (from 1 to 12.4 μm) is shown in Fig. A.1.

2.1.1 Data reduction and astrometry

We reduced the SPHERE/IRDIS data using an internally-developed Python programme, consisting of standard calibrations (dark current subtraction, flat field correction) and cosmetic processes (bad-pixel correction). The bad-pixel map is generated using the dark images, and a median filter is applied to be compared with the raw dark image. For classical imaging, using the same filter, IRDIS detector frames are first split into two parts. Then, all frames are shifted, added, and averaged using dithering information and cross-correlation³ in both parts of the camera. In DBI mode, two different filters are used in both parts of the detector and processed separately. A second cosmetic step is then applied to the images to correct for the residual hot pixels.

The IFSdata were reduced using the SPHERE Data Reduction and Handling (DRH) automated pipeline (Pavlov et al. 2008) at the SPHERE Data Centre⁴ (DC; Delorme et al. 2017). We applied basic corrections for bad pixels, dark current, and flat field and complemented the DRH pipeline with additional steps to improve the wavelength calibration, cross-talk, and bad pixel correction (Mesa et al. 2015).

Absolute orientation and pixel scales of the images were calculated using the parallactic angle (PARANG) and the absolute calibration provided by the SPHERE consortium. The current respective estimates of the pixel scale and true north (TN) angle for IRDIS are 12.255 ± 0.009 mas pixel⁻¹ and TN = −1.75 ± 0.08° (Maire et al. 2016). Our images were obtained using the pupil-stabilized mode, and therefore they were also derotated from the zero point angle of the instrument pupil (PUPIL_offset). The measured values are stable around an average value of −135.99 ± 0.11°. For IFS, the provided calibrated values are a pixel scale of 7.46 ± 0.02 mas pixel⁻¹ and a north angle of TN = −102.18 ± 0.13° (Maire et al. 2016). Finally, the orientation angle α_corr was calculated for each wavelength to align north upward and east to the left, according to: $${\alpha }_{\text{corr}}=-({\text{PARANG+TN+PUPIL}}_{\text{offset}}).$$(1)

The uncertainties of SPHERE data are computed taking the speckle noise, photon noise, and detector noise into account. We determine the errors for each pixel using the unbiased standard deviation calculation on all the individual frames (almost 300 frames for each filter). We specifically follow Eq. (2) for random variables following a normal distribution: $$S=k_n{\left[\frac{1}{n-1}{\sum }_i{(X_i-M)}^2\right]}^{1/2},$$(2) $$k_n=\sqrt{\frac{n-1}{2}}\frac{\Gamma \left(\frac{n-1}{2}\right)}{\Gamma \left(\frac{n}{2}\right)},n\ge 2,$$(3)

where, n is the frame number, M the mean value of the cube (the final image in this case), X_i the individual pixel value, k_n the corrector factor, and Γ(x) the gamma function (Euler integral of the second kind).

Figure C.1 presents all the reduced images of WR104 and the associated calibrators. Accounting for the orbital period of the system of 241.5 days, the images were all rotated to share the same orbital phase as that of 21 July 2016.

2.1.2 Image deconvolution

SPHERE is an instrument designed to directly detect planets around stars other than our own Sun. As a consequence, the instrument is designed to provide a very high Strehl ratio (≥ 90% in band K; Zurlo et al. 2014), meaning that deconvolution is an easier process than with previous Adaptive Optics (AO) equipped instruments. We deconvolve the images of WR104 using the associated PSF and the Scaled Gradient Projection (SGP) algorithm (Bonettini et al. 2009), implemented in the IDL programming language within the AIRY image reconstruction package (Correia et al. 2002), developed within the CAOS Problem-Solving Environment (Carbillet & La Camera 2018)⁵. The IFS data were deconvolved with a standard Richardson–Lucy algorithm (RL) implemented in Python (Lucy 1974). Considering the Poisson nature of the noise, the deconvolution problem is based either on the minimization of the Kullback–Leibler functional or the Csiszár I-divergence (Csiszár 1991) to which we can add a penalty term weighted by a real (positive) number, called a regularization parameter.

SGP is a gradient method that permits acceleration of the deconvolution process. RL can be seen as a particular case of the scaled gradient method, with a constant step-length. In contrast, SGP is an optimization method based on an adaptive strategy for the step-length parameter and is more efficient (from a computational point-of-view) than RL. In all cases, the algorithms are iterative and must be stopped before the noise amplification exceeds the fit of the data. In fact, RL-based methods need an early stoppage of the iterations to avoid the so-called checkerboard effect in the case of diffuse objects (Bertero et al. 2009).

In the case of RL, to determine the best compromise between the deconvolution efficiency and artefact creation, we tested a range of numbers of iterations between 10 and 300. We found that 60 iterations avoid any artefact creation and increase of the noise. In the case of SGP, we stopped the iterations when the objective function (sum of the Kullback–Leibler function and the regularization term) was approximately constant, given a tolerance of 10⁻⁷. In addition to the main SGP deconvolution process, we also applied a pre-processing of the reduced PSF and acquired images. During pre-processing, we first set the minimum flux of each image’s data to zero, then normalized the PSF to a unit integral, and added a small constant to the images in order to avoid divisions by zero during the deconvolution process. After this, the background value was evaluated and then considered within the algorithm (see Bertero et al. 2009).

Moreover, in the case of SGP, and to avoid the above-mentioned checkerboard effect when it appears (which was always the case except for the K-band data), we used a second-order Tikhonov regularization defined by the discrete Laplacian of the unknown object. The choice of the associated regularization parameter (β) was chosen in each specific data case as a trade-off between the checkerboard effect and a regularized reconstruction of the observed spiral shape, which was too smooth. The goal was to obtain the best possible fit of the spiral shape. In practice, the SGP deconvolution process of the data leads to:

• no regularization necessary in the K band (114 iterations to reach the fixed tolerance of 10⁻⁷) and in the CO filter (169 iterations);

• β= 0.01 (94 iterations) in the H band;

• β= 0.05 (188 iterations) in the FeII filter;

• β= 0.1 (82 iterations) in the J band.

It is clear here that shortening the wavelength (and hence roughly lowering the data’s Strehl ratio) implies increasing the regularization needed for the reconstruction, which is what is logically expected.

In Fig. 1, we show the deconvolution of WR104 in the J, H, and K bands as well as the FeII filter. The deconvolved images exhibit a clear spiral pattern, confirming for the first time with direct imaging the pinwheel nature of WR104 (aperture masking is considered here an indirect imaging technique, since it involves interferometric techniques and image reconstruction). More specifically, in the K-band deconvolved images, five turns of the spiral are detected. On the other hand, the raw images show flux in up to ten turns of the spiral (see Sect. 3.2).

Fig. 1 SGPdeconvolution in J, H, and K bands, as well as the FeII filter. Contour levels are 1, 5, 10, 20 and 50% of the peak. All images are cophased to the second epoch at 21 July 2016. The images are represented with a power normalization scale of 0.3. The field of view is 980 mas. We can clearly see the spiral pattern and the first revolution in all reconstructed images.

2.2 VISIR observation

VISIR is a mid-infrared imager and spectrograph installed at the Cassegrain focus of Unit Telescope 3 of the VLT (Lagage et al. 2004). Our observations were obtained during 1 h of Science Verification Time (SVT) in March 2016 with the new Raytheon AQUARIUS 1024 × 1024 pixel arrays and the new coronagraphic mode. We used the Annular Groove Phase Mask (AGPM-12.4 μm; Delacroix et al. 2012) and the 4-Quadrant Phase Mask (4QPM-10.5 μm; Kerber et al. 2014). Chopping and nodding were used to remove the sky background and thermal background. The chopping was parallel to the nodding, with a chop throw of 10′′ . With the chopping, we obtain images alternatively with and without the coronagraph centred on the target. We observed the PSF calibrator star HD169916 (spectral type K0III) unresolved by VISIR (λ∕D≈ 300 mas, while θ_UD = 3.90 ± 0.21 mas; Richichi et al. 2005) in the N band. The same star was also used as photometric calibrator with the corresponding absolute flux tabulated in Cohen et al. (1999).

Data reduction

The VISIR data were reduced with the official ESO-pipeline, which cleans the data from bad pixels. Nodding images are created by averaging the images in the two positions of the chopper. Since the images were obtained with a coronagraph with chopping in parallel mode, the final image contains a column of three images aligned along the north/south axis. The first and last images are two negative images without the coronograph, while the central image is with the coronagraph and has twice the integration time of the other two.

The corresponding images, in classical mode on top and coronagraphic mode at the bottom, are shown in Fig. 2 at 10.5 μm and Fig. 3 at 12.4 μm. WR104 is clearly more extended than the PSF reference star. This is confirmed by the radial profile of Fig. 4. We cannot distinguish any relevant structure around WR104, especially at the location of its second companion star (see Sect. 5.3). We also cannot confirm the spiral shape far from the star because of the relatively lower spatial resolution ofthe VISIR instrument compared to SPHERE (250 vs. 30–35 mas). Nevertheless, dust is present far from the star, at least up to 2 arcsec. This implies that a fraction of the dust can survive after many orbits (30 orbits = 2′′ ). However, we are unable to measure this fraction with the data at hand.

We note the low quality of the PSF reference star obtained with the AGPM coronagraph. The PSF seems to be elongated in the NE–SW direction. This issue may be due to a very low wind speed, favouring dome turbulence, and/or to the fact that the observations were made at the very end of the night, during twilight, after a relatively bad night (seeing ≥1.2 arcsec). Therefore, the AGPM dataset, especially regarding the PSF reference star, needs to be treated with caution.

Fig. 2 QPM coronagraph at 10.5 μm. Left panels: WR104, right panels: PSF reference star HD169916. Top panels: classical imaging. Bottom panels: coronograph mode. All images are displayed with a power normalization scale of 0.3. Contours at 2, 5, 10, and 20% of the maximum are represented. The theoretical resolution limit by the Rayleigh criterion of 1.22λ∕Dis also shown as a white circle in the upper left panel.

Fig. 3 AGPM coronagraph at 12.4 μm. Same description as for Fig. 2.

Fig. 4 Solidlines: mean radial profile of the non-deconvolved SPHERE and VISIR images of WR104. All profiles are shifted for readability reasons. Dotted lines: radial profiles of the PSF reference star. The diamonds indicate the theoretical angular resolution of the 8 m telescope at the given wavelength. We also show the position of the companion B star of WR104 (vertical line) for information (see Sect. 5.3).

3 Data analysis

3.1 Radial profile

The radiation intensity profiles at different wavelengths can tell us more about the type of physics behind the observed emission. Considering a smooth spherical shell with radius r, a power-law intensity profile with a power index of −2 (i.e. I(r) ∝ r⁻²) can be associated to optically thin emission, while a different power index indicates a more complex process (thermal emission, back warming, etc.; Rybicki & Lightman 1986).

We compute the azimuthally-averaged radial profiles on our raw (non-deconvolved) images for all available filters between 10.8 and 12.4 μm. All profiles are normalized to 1 at the peak and shown in Fig. 4 for comparison.

We compute the power-law index of the intensity decrement for all filters (the q parameter in Table 3). In practice, to avoid any problems of sub-resolution, we select a zone of the profile that is unaffected by the limited resolution of each instrument at the wavelength considered (between 40 mas in the J band and 400 mas in the N band). We only consider the region above the noise floor (around 700 mas for SPHERE and 3 arcsec for VISIR) and vary the fitting region a few tens of mas around these limits to reliably estimate the uncertainties of the power-law index. The resulting indices are presented in Table 3. Since the uncertainty on the radial profile varies widely across the different filters, the chi-squared values ${\chi }_{\text{red}}^{\text{2}}$ are not completely representative of goodness of fit. In Fig. D.1, we present the results of the different fits with the associated uncertainties in the data. The figure shows the differences between the filters and the clear deviation from a power law in the J band (CntJ and HeI filters).

3.2 Curvilinear profile of the pinwheel

Given the complexity of the dust distribution, we follow Tuthill et al. (2006) and compute the flux as a function of angular displacement along the spiral. With our data, we are able to detect the flux along the spiral for up to 15 coils (compared to 2 in Tuthill et al. 2006), as can be seen in Fig. 5.

First, we adjust the deconvolved images to our phenomenological model, further described in Sect. 4.2. This first fit is then refined with a comparison to the raw (non-deconvolved) image. We then use the fitted origin, orientation, and step of the spiral to compute the curvilinear flux over 15 turns in the J, H, and K bands and extract the flux at the expected position of the spiralrepresented by the grey, dashed line in the lower-right panel of Fig. 1. Using a cubic interpolation on the pixel coordinates, we compute the flux value in 10-degree increments of the spiral’s azimuthal coordinates.

The resulting flux exhibits a similar shape in all filters, with the best signal-to-noise ratio (S/N) in the K band. The flux decreases monotonically with radius over the first two to three coils, without any clear breaks. The ups and downs observed in the profile (especially after the first coil) are not real and seem to come from the convolution process with the PSF. Similar to Tuthill et al. (2006), we also observe a systematic offset between the start of the spiral and the brightest pixel. This maximum, occurring at an azimuth of 90° (J and H bands) and 110° (K band), is consistent with the offset of 90° in azimuth (or 12 mas in distance),attributed by Tuthill et al. (2006) to the dust formation location locus.

Fig. 5 Image flux along the fitted spiral (J, H, and K bands for the non-deconvolved images). Black dotted lines correspond to the first 15 turns of the spiral. We also note the presence of artefacts due to the K-band filter.

3.3 Spatially resolved spectroscopy with IFS

The two spectra of WR104 presented in Fig. 6 are extracted from the IFS data cube and calibrated using the associated PSF calibrator. We integrate the flux over all 39 spectral channels to give a raw spectrum both for thePSF reference star and WR104. We extracted one spectrum in the inner region (representing the star), and the other in the outer region (representing the pinwheel); these two regions of interest are shown in Fig. 7.

In the two resulting spectra, we can recognize only one emission line close to 1.08 μm, which we identify as a helium line. The same spectral feature appears in both spectra, which is an indication concerning the reflective nature of the dust emission. Nevertheless, due to the lack of angular resolution of the IFS data, we cannot distinguish if this HeI emission came from the WR star itself or from the hot-shocked plasma into the WCZ (both are contained in the inner region of the image).

The power of IFS data derives from their spatio-spectral nature. We use the spectrum from Fig. 6 and the IFS data cube to compute three narrow-band images around 1.01 μm (Y-band image, blue region in Fig. 6), 1.09 μm (HeI-band image, green region in Fig. 6), and 1.20 μm (J-band image, purple region in Fig. 6). We normalize all the images by the number of corresponding spectral channels in order to compare their fluxes. Figure C.2 shows the resulting images and the associated reference star for the PSF. Furthermore, we subtract the three resulting images to compute the colour index images (J-Y, He-Y and He-J) presented in Fig. 7. We attribute the main part of the continuum emission to the circumstellar dust and the HeI emission to the central binary star. Therefore, the different colours highlight the different components of the system:

• In the (J-Y) image, the dust (red) is separated from the star (blue) and highlights the spiral structure of the system.

• In the (He-Y) image, most of the observed flux comes from the dust and reveals the dusty environment of the star.

• In the (He-J) image, the main part of the observed flux comes from the star and reveals the position of the star.

These colour images easily reveal the complex dusty environment of WR104, especially in the case of the WR star, where the HeI line emission dominates this part of the spectrum (i.e. 0.95–1.35 μm for the IFS-SPHERE instrument).

We alsodetect a difference in the position of the photocentre in the images for the Y and J bands, offset by ≈ 10 mas, with a position angle of θ_offset = −140°. This angle does not seem to be correlated with the parallactic angle ( − 104°). This offset might be due to a problem of atmospheric dispersion correction. However, the SPHERE ADC ⁶ was designed to allow <2 mas displacement with wavelength, and the PSF data exhibit an offset of only ≈ 1.8 mas. Therefore, we conclude that this offset is likely to be real.

Fig. 6 Spectra of inner (blue) and outer (orange) regions of WR104 in the near infrared (NIR; the spectral line at 1.08 μm corresponds to HeI). The corresponding images, integrated over the coloured spectral ranges, are shown in Fig. C.2.

Fig. 7 Different colour images using Y,HeI, and J bands (blue, green, and purple parts of Fig. 6). The photocentres of the Y (black cross) and J bands (blue cross) are represented as well as the best fit of the spiral, shown by a grey dashed line. The dashed green circles both represent regions from which the two spectra have been extracted in Fig. 6: from the star (inner region) and the pinwheel (between the two circles).

3.4 The companion B of WR104

Our J-band images of WR104 exhibit the presence of a companion at a distance of approximately 1′′ from the central binary star. This companion star, hereafter referred to as “companion B”, was first detected with the Hubble Space Telescope (HST) in 1998 and reported in Wallace et al. (2002). According to Wallace et al. (2002), “the additional component [...] is likely to be physically related to the WR star [...]”. We now have means of verifying this thanks to the new extremely precise proper motions of stars provided by the Gaia mission (Zacharias et al. 2017).

Because Wallace et al. (2002) focused on companion detec- tion rather than astrometry, we retrieve their HST data (three filters: 336, 439, 555 nm) and recompute the astrometry. Our IRDIS-SPHERE data and these HST data were taken 18 yr apart. We fit two-dimensional (2D) Gaussian functions with a standard Levenberg–Marquardt algorithm to represent the central WR+O binary star and the companion B. The resulting vector of separation at all epochs and filters are provided in Table 4. The uncertainties given here are the statistical errors scaled to the χ² values (using the covariance matrix). We find a small linear motion of the companion B toward WR104 of 4 ± 1 mas between the two epochs and an angular motion compatible with zero.

The first Gaia data release DR1 (Gaia Collaboration 2016) provides a new and precise position measurement of WR104. The US Naval Observatory CCD Astrograph Catalog (UCAC5) uses a combination of Gaia DR1 measurements and the Naval Observatory Merged Astrometric Dataset (NOMAD) to compute precise proper motions of millions of stars, including WR104. We retrieved the following UCAC5 proper motions for WR104: pm_RA = −3.1 ± 1.0 mas yr⁻¹; pm_dec = − 2.4 ± 1.0 mas yr⁻¹ (Zacharias et al. 2017). If companion B were a background star (i.e. with a negligible proper motion), WR104 should have moved relative to it by 68 ± 17 mas between the HST and SPHERE observations: such a large motion would have been detectable. This increase is not detected, so companion B is unlikely to be a background star.

In Fig. 8, we show the expected motion of companion B relative to WR104 under the assumption that it is a background star (green dashed line with the position at different dates). In the figure inset, we show the measured position of companion B relative to WR104 made with SPHERE in 2016 (blue cross) and the HST in 1999 (orange crosses). This demonstrates a common proper motion (rate and angle) between WR104 and companion B.

That common proper motion provides an argument in favour of a gravitational link between companion B and the central binary WR+O star, as explained in Moe & Di Stefano (2017). Therefore, assuming that both stars are at the same distance, we identify companion B and WR104 as a common proper motion binary.

We also consider the effective temperature of the companion B star (defined as T_B ). To do so, we start by measuring the flux ratio between the WR+O central star and the companion star in all the available filters (see Fig. E.2). We adjust the flux ratio using different atmosphere models to represent the spectral energy distribution (SED) of the stars with PoWR⁷ models for the WR star(Sander et al. 2015) and Kurucz models for the two companion stars (OB-type inner binary and companion B; Castelli & Kurucz 2004). The WR+O binary is surrounded by dust (the pinwheel nebula), which is represented by a blackbody with temperature T_dust. We assume a WC9-subtype for the WR star to obtain the corresponding SED with the PoWR models (Sander et al. 2012). To select the appropriate Kurucz atmosphere model, we assume the spectral type of the inner OB-type star to be B0.5V as reported by Rosslowe & Crowther (2015). We scale the SED of WR+O to match the total luminosity of the binary (L_WR+O = 120 000 L_⊙; see Table 6). Finally, we compute the flux density using our estimation of the distance (D = 2.58 ± 0.12 kpc, see Sect. 4.1).

We scale the flux density of the pinwheel with the solid angle seen from Earth Ω _dust. Then, we compute the flux density of companion B using another scale parameter r_B , representing the radius of the star in units of solar radii r_⊙.

We usea Levenberg–Marquardt (LM) minimization of the χ² to find the best solution of the flux ratio between companion B and the WR+O+dust component (namely WR104). The fit provides the uncertainties on the scale parameter Ω_dust and the star B radius r_B using the covariance matrix. We used a constant, relative uncertainty on the flux ratio of 3% (corresponding to an uncertainty of the flux ofapproximately 2% found with the SPHERE data). Then, we perform a scan of the parameter space {T_B, T_dust} to provide a more robust 1σ uncertainty. For this study, we set the log(g) parameter in the Kurucz model of companion B to 4.5⁸. A precise determination of this parameter is beyond the scope of the determination of the flux ratio. Figure 9 shows the χ² map for a range of effective temperatures of the companion B (T_B between 15 000 and 45 000 K) and for the dust temperature (T_dust between 1000 and 3000 K). The figure contains the χ² map subtracted by the minimum χ² value (${\chi }_{\text{red}}^2=2.3)$ as well as the corresponding 1 and 2σ confidence intervals.

With the wavelengths used in this study (U, B, and V bands from Wallace et al. (2002) and J, H, and K bands from this study, we are able to provide a lower limit on the effective temperature of the companion star such that T_B ≥ 33 000 K (see Table 5). To further constrain the temperature of this companion, we would require additional ultraviolet (uv) data so as to account for the peak of the emission of the star.

Fig. 8 Measured positions of companion B (CompB) with respect to the position of WR104, set to the epoch of the SPHERE observation. The predicted motion of a hypothetical background star relative to WR104 is presented with a green dashed line. The measured positions of companion B are shown in blue (SPHERE in 2016) and orange (HST in 1999) crosses in the inset. The widths of the crosses indicate the uncertainties of the measured positions.

Fig. 9 Reduced χ2maps for the flux ratio fitting. Here, we present $\Delta {\chi }_r^2={\chi }_r^2-{\chi }_{\text{r},\text{min}}^2$and show the respective 1 and 2σ confidence intervals with red and purple lines.

4 Modelling the dust emission in the spiral

4.1 Simple Archimedean spiral

The simplest approach to describe the circumstellar environment of WR104 is to use an Archimedean spiral (also called arithmetic spiral⁹ ) to follow the dust distribution around the central binary star (Tuthill et al. 1999). The Archimedean spiral takes as hypothesis that the dust is created and accelerated up to a constant velocity and then follows a ballistic trajectory (linear radial motion).

Based on the determined period of P = 241.5 ± 0.5 days, the best fit of the spiral on our images provides an angular speed of v_ang = 0.273 ± 0.013 mas day⁻¹, which corresponds to the radial outward motion of the dust. It yields a separation of the individual turns of the spiral (“steps”) of 0.273 mas day⁻¹ × 241.5 mas = 66 mas.

Previous studies showed that the central binary system has a circular orbit (Tuthill et al. 2008). As a consequence, the dust nucleation locus follows the orbital motion of the system. Therefore, the spiral pattern is directly linked to the orbital period of the system and the speed of the expanding dust. Tuthill et al. (2008) demonstrated that the wind velocity at large radii (at least at the dust formation locus) is strongly dominated by the WR wind. The latter overwhelms the O star wind and radiative braking¹⁰ . As such, weconsider that the dust expansion velocity is set by the WR wind speed, so V _dust ≃ V _∞,WR, where V _∞ is the terminalwind velocity, with only minor adjustments.

This means that if we know the physical speed of the dust over an orbital period, we can directly determine the distance of the system. Conversely, a precise measurement of the distance can yield a unique measurement of the WR-dominated wind speed. In Table 6, we provide all the physical parameters of WR104 found either in the literature (Crowther 1997; Monnier et al. 2002) or in tables assuming typical values for the WR and O stars (Sander et al. 2012; Fierro et al. 2015). These parameters are the temperature (T), the stellar mass (M), the terminal wind velocity (V _∞), the gas mass-loss (Ṁ), and the total luminosity (L).

Reported terminal velocities of WC9 stars range from 1220 km s⁻¹ (Howarth & Schmutz 1992; Crowther 2007) up to 1600 km s⁻¹ (Sander et al. 2012). Taking the lower value of 1220 km s⁻¹, we find a distance of D = 2.58 ± 0.12 kpc, which is in good agreement with previous measurements and makes WR104 a potential member of the Sgr OB1 stellar association (Lundstrom & Stenholm 1984). Taking the most recent value of 1600 km s⁻¹ (Sander et al. 2012), we obtain a distance of D = 3.38 ± 0.15 kpc, which would place WR104 behind the Sgr OB1 association.

Future measurements of distances with Gaia, expected with the data releases DR2 (Gaia Collaboration 2018, see Sect. 5.4) and DR3, will provide a very good estimate of the still unconstrained terminal wind velocity of the WR star, in the case when these WR stars are surrounded by a resolved pinwheel nebula.

4.2 Phenomenological model

To probe more of the system’s physical parameters, we update the model from Millour et al. (2009a), available in the public distribution of fitOmatic¹¹. This updated version relies on physical parameters of the object such as the temperature laws, the dust-formation regions, and the opening angle of the shock. As mentioned in the Introduction, it is the shock caused by the collision between the winds from the WR andO stars that results in the dusty spiral. This model provides a reasonable match to the flux of the inner part of the spiral.

The model is composed of the pinwheel and an inner binary system with a separation d_bin fixed to 1 mas (2.6 au at 2.6 kpc). The pinwheel is constructed with a succession of rings growing linearly with the distance from the central binary and positioned following an Archimedean spiral pattern with a given step. This spiral model is motivated by the hypothesis that the central binary is on a circular orbit and that the dust nucleation occurs at the interface of the shocked winds of theWCZ (Lamberts et al. 2012; Hendrix et al. 2016).

The flux contributions of the two components are calculated at different wavelengths using blackbodies. We chose 45 000 K for the WR star, 30 000 K for the O star (Table 6) and assumed that each ring of the pinwheel emits as a blackbody at a temperature defined as a decreasing function of the ring’s distance to the central system. The beginning of the pinwheel is defined by the parameter r_nuc corresponding to the dust nucleation location. The temperature T follows a power law with the nucleation temperature (T_nuc at the r_nuc distance) and a power-law index q: $$T(r)=T_{\text{nuc}}\times {\left(\frac{r}{r_{\text{nuc}}}\right)}^{-q}.$$(4)

We consider here that the nucleation temperature T_nuc is identical to the dust-sublimation temperature T_sub. The binary star contribution to the total flux is calculated using the scale parameter C_B∕P, arbitrarily set at a wavelength of 1 μm¹². We set the semi-opening angle of the spiral to α = 17.5° as determined by the momentum flux ratio η of the winds and the empirical relation from Eichler & Usov (1993): $$\eta =\frac{{\dot{M}}_{OB}{v}_{\infty /OB}}{{\dot{M}}_{WR}{v}_{\infty /WR}},$$(5) $$\alpha =2.1\left(1-\frac{{\eta }^{4/5}}{4}\right){\eta }^{1/3},\text{ for \hspace{0.17em}}{\text{10}}^{\text{-4}}\le \eta \le \text{1}.$$(6)

We fixed the number of spiral turns to 3.5 to be able to correctly fit the steps and the sky orientation defined by θ₀ (the spiral orientation at the reference date of May 2016). We also set the inclination of the orbital plan to zero (i =0°). An illustration of this phenomenological model and its main parameters is presented in Fig. 10.

We use a two-stage model-fitting approach: first we obtain rough values of the parameters using the deconvolved K-band image only, and then we refine the parameters with the non-deconvolved images. The fitting is a least squares calculation based on the pixel values $${\chi }^2={\sum }_i\frac{{(O_i-M_i)}^2}{e_i^2},$$(7)

where O_i is the ith pixel’s intensity, e_i its associated error, and M_i the prediction from the model. We use a Levenberg–Marquardt (LM) minimization of the χ² to find the best solution. To ensure that we have not converged on a local minimum and that our fit is robust, we also scan the parameter space of the relevant parameters (step of the spiral, orientation, x0, y0).

Table 7 presents the best-fit parameters of our phenomenological model. Taking into account the resolution of the VLT in the K band (λ∕D = 60 mas) and the pixel size of IRDIS (12.25 mas), we are able to set a sub-resolution constraint on the step and orientation of the spiral, as well as the dust-formation location (d_shock). By assuming two different sublimation temperatures (T_sub = 1500 K as used by Harries et al. 2004; T_sub = 2000 K for carbon from Kobayashi et al. 2011), we also derive a new measurement of temperature along the spiral.

We find better agreement for the higher sublimation temperature of 2000 K (q₂₀₀₀ = 0.45, χ² = 1.24). This is consistent with the hypothesis that the dust is composed of amorphous carbon. Indeed, since WR104 is a WC star, it produces large amounts of carbon. We assume that the chemistry of dust nucleation is based on carbon and not silicates. By fitting the SED of different WR stars of the Infrared Astronomical Satellite (IRAS) catalogue (such as WR104), Zubko (1998) showed that amorphous carbon seems to be more representative of the dust in the dusty WR star.

The uncertainties are given with a 1σ confidence level, and each parameter is constrained one by one around a given range of values (e.g. 0.3 < q < 0.5). To compare the orientation obtained with SPHERE with those from Tuthill et al. (2008), the phase θ₀ is given at the reference date of May 2016 (MJD = 57512). Our two epochs allow us to confirm the 241.5 day orbital period of the pinwheel.

Fig. 10 Phenomenological model used in this work. The important parameters are the dust formation zone (purple), opening-angle (green), and step of the spiral (blue). Projection angles on the sky are the same as in Millour et al. (2009b). Left panel: spiral as seen from the Earth (face-on,i = 0°) and right panel: model seen in its equatorial plane (edge-on; i = 90°).

4.3 Radiative transfer modelling

Although a complete radiative transfer model of WR104 is beyond the scope of this paper (see the study by Harries et al. 2004), a qualitative comparison based on radiative transfer is necessary to further interpret the new WR104 data presented here.

To do so, we compare our images with a three-dimensional (3D) axisymetric radiative transfer model based on a geometric model made of consecutive rings. These rings represent the consecutive coils of the dusty Archimedean spiral at a given azimuth. The corresponding 3D density grid is then fed into the well-established publicly available code RADMC3D (Dullemond et al. 2012), which computes the self-consistency of temperature distributions and produces mock images.

RADMC3D uses a Monte-Carlo method and launches photon packets from the central star into the dust density grid in order to compute the dust temperature at each point of the grid. The scattering source function is then computed at each wavelength through an additional Monte-Carlo run. We assumed an isotropic scattering. A last step uses ray tracing to compute the mock images.

The NIR continuum emission we are modelling is dominated by the dust. Since we do not aim for a detailed fit of the spectral features already done elsewhere, our approach is to assume realistic standard opacity laws. For the spiral, we expect a sub-micron grain size due to the rapid growth of the dust nuclei (Zubko 1998), so we assume a pure amorphous carbon composition with a grain size of 0.1 μm. In this model, the total dust mass is degenerate with the grain size, and results for the total mass are valid only for the grain size considered.

Three dust density models are presented in Fig. 11, each one corresponding to a physical hypothesis about the dust nucleation.

• Model 1 consists of uniformly thick rings of linearly increasing height H separated by the spiral step and an opening angle α set by the WCZ opening angle. The theoretical basis is that the dust is formed downstream from the reconfinement shock behind the O star, where the dust density can reach very high values for low temperatures (Lamberts et al. 2012). The dust then follows a ballistic trajectory, and its size grows linearly with the distance, which yields a uniform density in the rings.

• Model 2 consists of hollow rings, with the width of the walls set by w (expressed as a fraction of the total ring size H). Such a model is associated with the hypothesis of dust nucleation at the interface of the shock. The wind collision interfacefar from the central binary star offers ideal conditions to reach critical densities. This makes dust nucleation possible especially where the mixing between the carbon-rich wind of the WR component and the H-rich wind of the Ocomponent becomes significant (Hendrix et al. 2016). As in Model 1, the dust rings will then grow linearly to create successive coils with a constant fraction of over-density in the walls of the rings.

• Model 3 also consists of hollow rings. This time the first ring also has a hole of size h(expressed as a fraction of the ring height) in both walls facing the central stellar source. This setup simulates variations or breaks in the dust density, which could stem from the turbulence inside the shocked region. In this case, dust formation is not completely uniform. Lamberts et al. (2012) and Hendrix et al. (2016) discussed the possibility of hydrodynamical instabilities at work in the WCZ, which would then create a non-uniform density distribution. For instance, thin shell instabilities occur when cooling is important so that the shocked zone narrows to a thin layer, which is easily perturbed (Vishniac 1994). Such instabilities provoke strong distortions of the whole colliding region (Pittard 2009; van Marle et al. 2011). Therefore, if the dust is formed in this turbulent region, the resulting density distribution will be non-uniform with over- and under-densities.

For the three types of models, the masses M of the rings are independent of distance and calculated using the period of the spiral P, the mass-loss rate Ṁ of the dominant¹³ WC9 component, and a dust-to-gas ratio ξ: $$M=\xi \dot{M}\frac{P}{365.25}.$$(8)

We compare our model intensity profiles convolved with the PSF to the observed radial intensity profiles in the J, H, and K bands at a given azimuth, which is set such that the first ring is positioned at 30 mas from the central binary star. To avoid any perturbation caused by the bad quality of the PSF reference star (low dynamic range), we make the comparison between 30 and 300 mas.

We test a range of dust-to-gas ratios (from ξ = 0.1% to ξ = 50%) on the three models, with a fixed value of w = 10% and h = 30%. For each model, we compute a χ² criterion¹⁴ for comparison.

For all models considered, we cannot satisfactorily reproduce the radial profile obtained with SPHERE. Particularly in the J band, the emission from the dusty rings, including thermal emission and scattering, is negligible relative to the star. However, the WR104 images in the J band are definitely resolved with a significant fraction of dust emission (Fig. C.1; second row at 1.21 μm). For this reason, we do not include the Jband in our comparison.

Considering the H- and K-band data only, the resulting ${\chi }_{\text{red}}^2$ computation for all considered dust-to-gas ratios is shown in Fig. 12. It appears that Model 1 (thick rings) provides the best match (${\chi }_{\text{red}}^2=1.31$), with ξ = 1%. Figure 13 shows the corresponding observed and modelled radial profiles in the H band. As shown in Fig. 12, the agreement of Models 2 and 3 with the H and K band data is not as good. Nevertheless, the best match is obtained for ξ = 5% with ${\chi }_{\text{red}}^2=2.15$ for Model 2 and ${\chi }_{\text{red}}^2=2.25$ for Model 3.

Fig. 11 Different dust density distributions considered in our radiative transfer models. From top to bottom (edge-on view of the rings model delineated by dashed lines): full rings, hollow rings, and hollow rings with a central hole. Further details can be found in the text.

Fig. 12 Comparison of the ${\chi }_{\text{red}}^2$ for the three models considered: Model 1 (solid blue line), Model 2 (dashed orange line), and Model 3 (dotted green line).

Fig. 13 Comparison between the best model (Model 1; ξ = 1%) and the H-band data. We represent the radial profile for the model (green), reference star (dashed grey), model convolved with the PSF reference star (red), and the data (blue). The comparison points are represented with red crosses.

5 Discussion

5.1 Curvilinear profile and shadowing effect

In Sect. 4.2, we presented the different results obtained from the fit of our phenomenological model. We showed that our model is in good agreement with the K-band non-deconvolved data (${\chi }_{\text{red}}^2=1.24)$. Our model was based on a power-law decrease of the temperature with a constant index (q₂₀₀₀ = 0.45 ± 0.03) and a dust sublimation temperature of 2000 K. This value is very close to the theoretical power-law index of 0.5 corresponding to an optically thin dust distribution (Dullemond & Monnier 2010). The continuous decrease in temperature implies a continuous curvilinear intensity profile along the spiral. This is in contradiction with the rapid drop in temperature, and correspondingly in intensity, detected after the first spiral turn by Harries et al. (2004). They showed that an optically thick first turn could explain such a temperature drop by shadowing the rest of the spiral. With complementary data, Tuthill et al. (2008) derived a curvilinear intensity profile, which also exhibited an important drop after the first coil.

On this basis, we further investigate the possibility of shadowing effects in this system. For this purpose, we create an artificial gap factor in the spiral intensity profile of our best phenomenological model. This gap represents the shadow cast by the first turn of the spiral onto the second and mimics a significantly optically thick first turn. Following Tuthill et al. (2008), we set this gap factor to ten. We compare the K-band data with the curvilinear profile of the model with and without gap over three revolutions.

In Fig. 14, we show both the models with an optically thin first turn (solid grey line) and with an optically thick first turn (dashed grey line). Then we convolve these models with the PSF reference star (orange lines) and compute the ${\chi }_{\text{red}}^2$ for the models with and without the shadowing effect. The model with a gap seems to be less representative (${\chi }_{\text{gap}}^2=4.3$ compared to${\chi }_{\text{nogap}}^2=1.5$).

This conclusion is in apparent contradiction with the results of Harries et al. (2004). They considered an optically thick inner part of the spiral that would produce steps in the curvilinear profile; such steps were detected afterwards by Tuthill et al. (2008). The previous studies were based on image restoration techniques that are known to produce artefacts in the flux’s spatial distribution. Here, we based our analysis on direct imaging data with extreme adaptive optics and had a Strehl ratio of 0.87, which also has its flaws. Both image restoration and direct imaging were used at the limit of angular resolution. In our case, the PSF convolution significantly smooths possible intensity variations in the curvilinear profile. Images with a better resolution will be needed to definitely reach a conclusion on this aspect (e.g. images from the European Southern Observatory’s new Multi-AperTure mid-Infrared SpectroScopic Experiment (MATISSE) spectro-interferometer or the planned Extremely Large Telescope (ELT)).

Fig. 14 Comparison of the curvilinear flux in the K band over three orbits. We represent the best model with the gap (dashedline) and without it (solid line) after the first turn. There is no apparent evidence of the gap in the data (${\chi }_{\text{gap}}^2=4.3$,${\chi }_{\text{nogap}}^2=1.5$).

5.2 Radiative transfer and dust mass

Our first conclusion from the radiative transfer model is that in the J band, we overestimate the stellar contribution with respect to the observations. This suggests the presence of circumstellar dust along the line of sight, which could obscure the star mostly in the J band, and also very likely in the V band. At this point of the study, we cannot infer a geometry of this dust and can only consider uniform attenuation. This conclusion is supported by the fact that the WR104 system has regular eclipses in the optical (up to 3 mag; Williams 2014).

Considering the H- and K-band data only, the comparable values of thedifferent ${\chi }_{\text{red}}^2$ do not allowus to clearly discriminate between the three models of the dust geometry. Figure 13 shows the strong impact of the convolution of the PSF on the modelled intensity profile. All differences between the three models are significantly attenuated by the convolution process and cannot be satisfactorily constrained by our dataset. This is partly due to the relatively bad quality of our reference star, especially in the K band. It appears that the chosen reference star presents an infrared excess in its SED, which is indicative of circumstellar material. The reference star also presents a very low dynamic range with a poor S/N, which is probably due to the unoptimized choice of neutral density and exposure time during the observation.

Nevertheless, a trend emerges from the study of the dust-to-gas ratio. All models favour an intermediate dust-to-gas ratio, between 1% and 10%. This translates into masses of the dust rings between 5.3 × 10⁻⁸ M_⊙ and 5.3 × 10⁻⁷ M_⊙. For comparison, Harries et al. (2004) set the mass-loss rate of Ṁ = 3 × 10⁻⁵ M_⊙ yr⁻¹ in their model and determined a dust-creation rate of Ṁ = 8 × 10⁻⁷ M_⊙ yr⁻¹, which yields a dust-to-gas ratio of ξ_Harries = 2.7%. This leads to a ring mass of 5.3 × 10⁻⁷ M_⊙. Our results are compatible with these values but seem to suggest a less massive dusty environment than inferred by Harries et al. (2004). This raises the question of the efficiency of dust nucleation processes around WR stars.

5.3 The third star

The third “companion B” star located at ≈1′′ from the central binary is confirmed to be a hot star gravitationally linked to the inner binary. This confirms that WR104 is at least a hierarchical triple massive system.

Assuming that companion B is orbiting in the same plane as the central binary star, that is, in the plane of the sky (the inclination of WR104 is close to 0), we can estimate the physical separation between companion B and WR104. Taking into account the distance of the system D = 2.58 ± 0.12 kpc (see Sect. 4.1) and assuming a circular orbit (e = 0), the semi-major axis would be a ≈ 2480 ± 120 au¹⁵ .

To estimate the orbital period of this triple system, we can make assumptions on the stellar masses of the two O star components (companion B and the OB companion of the WR star). The luminosity and the measured effective temperature (Crowther 1997) yield an inferred a mass of M_OB = 20 M_⊙ and a corresponding age of approximately 7 Myr along the main sequence (CMFGEN models; Fierro et al. 2015). According to Wallace et al. (2002) and our study of the effective temperature of companion B (Sect. 3.4), the companion is very likely another massive OB type star. Therefore, we also set the mass of the companion B at M_OB = 20 M_⊙. According to the WR star’s spectral type as a WC9 (Sander et al. 2012), we set its mass to M_WR = 10 M_⊙. We treat the central binary star (M_WR+O = 30 M_⊙) as a single star and compute the orbital period of the companion B star as $$P=\sqrt{\frac{4{\pi }^2{a}^3}{G(M_{\text{WR+O}}+M_{\text{OB}})}}\approx 17000\pm 1300\text{\hspace{0.17em} days}.$$(9)

The orbital plane of the central binary star is close to the plane of the sky (≤ 16°). The projected separation between the central binary and companion B is 977 mas, and the distance between two spiral steps is measured to be 66 mas. If, as we assumed at the beginning of this section, the tertiary companion B star orbits around WR104 in the same orbital plane as the central binary, we expect it to eventually encounter the 15th outer dust coil, produced ten years ago. Such an encounter might be seen either with hot dust around companion B or a wealth of other phenomena (X-rays, bow shock)

We do not detect hot dust (infrared excess) around companion B. This means that it may lie outside of the plane of the binary. Based on the estimate of the opening angle of the wind collision zone of 35° (see Eq. (5)), we find that “B” may be located beyond 35°.

A misaligned set of orbital planes in triple systems has already been observed in other massive stars systems, like Algol (Baron et al. 2012) and σ Ori (Schaefer et al. 2016). Recent N-body numerical simulations report that the non-alignment of triple systems is very common, especially for large projected separations (>1000 AU; Tokovinin 2017). The formation of multiple massive stellar systems is still not fully understood. Multiple isolated systems are probably the result of the collapse of massive cores (Krumholz et al. 2009). The gravitational instabilities occurring during the formation cause protostellar disk fragmentation, allowing massive multiple systems to form and to overcome the uv radiation pressure barrier (Wolfire & Cassinelli 1987; Motte et al. 2018). Some hypotheses are favoured to explain the misalignment of triple systems: possible accretion of gas with randomly aligned angular momentum at the epoch of star formation, or dynamical processes such as the eccentric Kozai mechanism (Antognini & Thompson 2016).

Hierarchical triple systems are possible candidates to be at the origin of black holes or neutron star mergers (Silsbee & Tremaine 2017). Therefore, the WR104 system could be a progenitor of future compact object merger and gravitational wave emission.

5.4 Gaia DR2 distance and wind velocity

As discussed in Sect. 4, an interesting characteristic of the pinwheel nebula is the possibility to obtain a distance estimation using some hypotheses about the wind velocity. Now, if we are able to find the distance another way (parallaxes for instance), we can retrieve the dust velocity of the pinwheel and the associated wind speed. We used the last distance determination provided by the Gaia DR2 (Gaia Collaboration 2018). With the Gaia parallax (plx_G = 0.2431 ± 0.0988 mas), we can naively compute the new distance of the WR104 system as D_Gaia = 1∕plx_G = 4.11 ± 1.67 kpc, which would place WR104 further than expected, but is still consistent with the previous values. Therefore, if we compute the corresponding dust speed, we find V _dust = 1945 ± 795 km s^-1 which is compatible but higher than previously reported.

Bailer-Jones et al. (2018) used a more sophisticated method to compute the geometrical distance using the Gaia DR2 catalogue. The non-linearity of the transformation and the asymmetry of the resulting probability distribution does not allow a simple inversion of the Gaia parallaxes. They used a Bayesian procedure and a galactic model to give a more accurate estimate of the distance with a 68% confidence interval. This method gives a distance for WR104 of D_Bailer = 3.64 + 1.92∕−1.02 kpc, which is still higher than expected but nevertheless compatible with our distance estimate (D = 2.58 ± 0.12 kpc; Sect. 4.1). This latest distance measurement allows us to infer the associated dust velocity as V _dust = 1721 ± 911 km s⁻¹. We used the upper limit of the distance uncertainty (1.92 kpc) to propagate the error on the speed. This velocity is higher than the terminal wind speed of the WR component, usually used as the dust velocity at the dust nucleation locus (Tuthill et al. 2008). The relatively high uncertainty of the Gaia DR2 prevents us from drawing a conclusion about the exact dust velocity. Nevertheless, the trend is towards a higher speed than expected, which could be explained by dust nucleation closer to the star or more acceleration by the wind of the OB star. This raises the question about the still unconstrained dust formation locus as well as the velocity of dust launch.

Nevertheless, some warnings were reported concerning the Gaia DR2 parallax measurements, in particular in cases where the source is bright (G < 13; Lindegren et al. 2018), extended, distant (D > 100 pc), and lacking a Gaia measurement of radial velocity. Since WR104 fulfils all these criteria, we are cautious with this naive distance estimate.

6 Conclusions

With direct imaging using the SPHERE and VISIR instruments on the VLT, we have confirmed the spiral structure of the WR104 system for the first time. We probed the extension of the dust spiral over 15 revolutions with SPHERE and 30 with VISIR, with an unprecedented dynamic range. This corresponds to the history of mass loss in the last 20 yr. Furthermore, we determined a step of the spiral of 66 ± 3 mas, leading to a distance of the system of 2.58 ± 0.12 kpc, thus refining the previous estimate by Tuthill et al. (2008).

Based on the IFS data, and the model fitting of IRDIS data, we confirmed the presence of the dust formation zone approximately 12 mas away from the central binary. The determination of the dust formation region is done at the performance limit of the SPHERE instrument and needs to be confirmed with the higher angular resolution capabilities of interferometric instruments. In the future, the next generation of interferometric instruments, such as MATISSE (Lopez et al. 2014; Matter et al. 2016), will be able to reveal the spiral shape of WR104 at its peak-emission (L-band) and with unprecedented spatial resolution (i.e. 4 mas).

We also confirm that the third component of the system, named “companion B”, discovered with the HST (Wallace et al. 2002), could be gravitationally bound to the inner binary. This companion likely orbits in a different plane compared to the central binary and provides information on the formation of the triple system (competitive or core accretion).

We created an updated version of the phenomenological model presented in Millour et al. (2009b) and the curvilinear profile of the pinwheel to show that the shadowing effect previously reported does not seem to occur around WR104. This conclusion needs to be confirmed with a new K-band high-dynamic range measurement, which would allow the system to be followed during a longer time lapse.

The radiative transfer models used in this work seems to favour a dust-to-gas ratio between 1% and 10%. Unfortunately, our dataset does not allow us to distinguish between the different hypotheses of the process of dust nucleation. This highlights the need to acquire a new SPHERE dataset with larger bandwidth filters (to improve the S/N) and a better reference star in the J and K bands (4 Sgr appears to be a good candidate considering the better quality of the H band).

Furthermore, to make conclusions about the flow dynamics as well as the dust formation process and its exact location, we need more sophisticated models, including coupling hydrodynamical models with the radiative transfer processes. This accurate determination is only achievable with a thorough comparison between high-level models and new interferometric data.

Acknowledgements

This work uses the SPHERE Data Centre, jointly operated by OSUG/IPAG (Grenoble), PYTHEAS/LAM/CeSAM (Marseille), OCA/Lagrange (Nice), and Observatoire de Paris/LESIA (Paris) and supported by a grant from Labex OSUG@2020 (Investissements d’Avenir – ANR10 LABX56). We acknowledge support from the CNRS, University of Nice Sophia Antipolis, and the Observatoire de la Côte d’Azur. We used NASA’s Astrophysics Data System Bibliographic Services. This research made use of Astropy¹⁶ , a community-developed core Python package for astronomy (Astropy Collaboration 2013). Support for A.L. was provided by an Alfred P. Sloan Research Fellowship, NASA ATP Grant NNX14AH35G, and NSF Collaborative Research Grant 1715847 and CAREER grant 1455342.

Appendix A Filter transmissions

Fig. A.1 Transmission curves of all filters used in the present work. The exact transmission curves of the newly installed coronagraphic filters of VISIR are not yet available, so we present the bandwidth communicated by the ESO.

Appendix B Detailed log of observations

Appendix C All SPHERE images

Fig. C.1 All reduced SPHERE images with the corresponding PSF calibrators. All images are rotated to be phased with the last SPHERE observation epoch (21 July 2016). The true orientation is shown on all panels.

Fig. C.2 IFS images in the different parts of the spectra (see Fig. 6). Top panels: blue part of Fig. 6 (Y band, λ = 1.01 ± 0.05 μm). Middle panels: green part of Fig. 6 (He line, λ = 1.09 ± 0.03 μm). Bottom panels: purple part of Fig. 6 (J band, λ = 1.20 ± 0.07 μm). We represent the reduced image (left), the PSF calibrator (middle), and the deconvolution (right). The FWHM of the PSF is also represented by a white circle in all the panels of the left column.

Appendix D Radial profiles and uncertainties

Fig. D.1 Results of a power-law fit on the radial profiles. Top panels: two good fits with a representative estimation of the uncertainty in the data (in light blue). Middle panels: good fit but with data of poor quality. Bottom panels: deviation of the radialprofile from a power law in the J band.

Appendix E Complementary model fitting results

In addition to Sect. 4.2, we present a comparison of the best-fit phenomenological model with the K-band image. The step and orientation of the spiral is represented in Fig. E.1 as a dashed red line. The K-band data are well-fitted by the model, with residuals within 1σ compared to the standard deviation.

Fig. E.1 Comparison of the best-fit model and the reduced image in K band. The residuals are represented in terms of σ (standard deviation of the image).

As presented in Sect. 3.4, we used different atmosphere models to determine the effective temperature of the third component of WR104, companion B. Figure E.2 shows the results of the best fit of the flux ratio between the central binary star and the companion star. We present the resulting flux density at the distance of WR104 (D = 2.58 ± 0.12 kpc; see Sect. 4.1) of the four components in Fig. E.3. The total flux density is comprised of two Kurucz atmosphere models representing the companion B and the OB-type inner companion: a Potsdam Wolf-Rayet (PoWR) model representing the WR star and a blackbody representing the dust around the central binary.

Fig. E.2 Flux ratio between the central binary star surrounded by dust and the companion B star. Our best fit is represented in orange (Tdust = 2200 K, TB = 45 000 K).

Fig. E.3 Best fits of the flux density from: 1) two Kurucz atmosphere models representing the companion B (green dashed line) and the OB-type inner companion (orange line), 2) one PoWR model representing the WC9 star (blue line), and 3) one blackbody model representing the dusty pinwheel (red dashed line). The flux ratio is obtained by comparison between theWR+OB flux density (grey line) and the companion B. We also show the total flux density (black dashed line) and position of the different filters used (vertical black dotted lines).

The infrared excess of WR104 compared to the companion indicates the presence of relatively hot dust around the central binary star, that is, the pinwheel nebula. Given the relative flatness of the flux ratio in the B and V bands, the HST data provide important information about the stellar temperature. This indicates that the companion B star has an effective temperature comparable to the OB star of the inner binary (which dominates the total flux), and is probably an OB type star itself, thereby confirming the findings of Wallace et al. (2002).

References

All Tables

All Figures

	Fig. 1 SGPdeconvolution in J, H, and K bands, as well as the FeII filter. Contour levels are 1, 5, 10, 20 and 50% of the peak. All images are cophased to the second epoch at 21 July 2016. The images are represented with a power normalization scale of 0.3. The field of view is 980 mas. We can clearly see the spiral pattern and the first revolution in all reconstructed images.
In the text

Fig. 2 QPM coronagraph at 10.5 μm. Left panels: WR104, right panels: PSF reference star HD169916. Top panels: classical imaging. Bottom panels: coronograph mode. All images are displayed with a power normalization scale of 0.3. Contours at 2, 5, 10, and 20% of the maximum are represented. The theoretical resolution limit by the Rayleigh criterion of 1.22λ∕Dis also shown as a white circle in the upper left panel.

In the text

	Fig. 3 AGPM coronagraph at 12.4 μm. Same description as for Fig. 2.
In the text

Fig. 4 Solidlines: mean radial profile of the non-deconvolved SPHERE and VISIR images of WR104. All profiles are shifted for readability reasons. Dotted lines: radial profiles of the PSF reference star. The diamonds indicate the theoretical angular resolution of the 8 m telescope at the given wavelength. We also show the position of the companion B star of WR104 (vertical line) for information (see Sect. 5.3).

In the text

	Fig. 5 Image flux along the fitted spiral (J, H, and K bands for the non-deconvolved images). Black dotted lines correspond to the first 15 turns of the spiral. We also note the presence of artefacts due to the K-band filter.
In the text

	Fig. 6 Spectra of inner (blue) and outer (orange) regions of WR104 in the near infrared (NIR; the spectral line at 1.08 μm corresponds to HeI). The corresponding images, integrated over the coloured spectral ranges, are shown in Fig. C.2.
In the text

Fig. 7 Different colour images using Y,HeI, and J bands (blue, green, and purple parts of Fig. 6). The photocentres of the Y (black cross) and J bands (blue cross) are represented as well as the best fit of the spiral, shown by a grey dashed line. The dashed green circles both represent regions from which the two spectra have been extracted in Fig. 6: from the star (inner region) and the pinwheel (between the two circles).

In the text

Fig. 8 Measured positions of companion B (CompB) with respect to the position of WR104, set to the epoch of the SPHERE observation. The predicted motion of a hypothetical background star relative to WR104 is presented with a green dashed line. The measured positions of companion B are shown in blue (SPHERE in 2016) and orange (HST in 1999) crosses in the inset. The widths of the crosses indicate the uncertainties of the measured positions.

In the text

	Fig. 9 Reduced χ2maps for the flux ratio fitting. Here, we present $\Delta {\chi }_r^2={\chi }_r^2-{\chi }_{\text{r},\text{min}}^2$and show the respective 1 and 2σ confidence intervals with red and purple lines.
In the text

	Fig. 10 Phenomenological model used in this work. The important parameters are the dust formation zone (purple), opening-angle (green), and step of the spiral (blue). Projection angles on the sky are the same as in Millour et al. (2009b). Left panel: spiral as seen from the Earth (face-on,i = 0°) and right panel: model seen in its equatorial plane (edge-on; i = 90°).
In the text

	Fig. 11 Different dust density distributions considered in our radiative transfer models. From top to bottom (edge-on view of the rings model delineated by dashed lines): full rings, hollow rings, and hollow rings with a central hole. Further details can be found in the text.
In the text

	Fig. 12 Comparison of the ${\chi }_{\text{red}}^2$ for the three models considered: Model 1 (solid blue line), Model 2 (dashed orange line), and Model 3 (dotted green line).
In the text

	Fig. 13 Comparison between the best model (Model 1; ξ = 1%) and the H-band data. We represent the radial profile for the model (green), reference star (dashed grey), model convolved with the PSF reference star (red), and the data (blue). The comparison points are represented with red crosses.
In the text

	Fig. 14 Comparison of the curvilinear flux in the K band over three orbits. We represent the best model with the gap (dashedline) and without it (solid line) after the first turn. There is no apparent evidence of the gap in the data (${\chi }_{\text{gap}}^2=4.3$,${\chi }_{\text{nogap}}^2=1.5$).
In the text

	Fig. A.1 Transmission curves of all filters used in the present work. The exact transmission curves of the newly installed coronagraphic filters of VISIR are not yet available, so we present the bandwidth communicated by the ESO.
In the text

	Fig. C.1 All reduced SPHERE images with the corresponding PSF calibrators. All images are rotated to be phased with the last SPHERE observation epoch (21 July 2016). The true orientation is shown on all panels.
In the text

Fig. C.2 IFS images in the different parts of the spectra (see Fig. 6). Top panels: blue part of Fig. 6 (Y band, λ = 1.01 ± 0.05 μm). Middle panels: green part of Fig. 6 (He line, λ = 1.09 ± 0.03 μm). Bottom panels: purple part of Fig. 6 (J band, λ = 1.20 ± 0.07 μm). We represent the reduced image (left), the PSF calibrator (middle), and the deconvolution (right). The FWHM of the PSF is also represented by a white circle in all the panels of the left column.

In the text

	Fig. D.1 Results of a power-law fit on the radial profiles. Top panels: two good fits with a representative estimation of the uncertainty in the data (in light blue). Middle panels: good fit but with data of poor quality. Bottom panels: deviation of the radialprofile from a power law in the J band.
In the text

	Fig. E.1 Comparison of the best-fit model and the reduced image in K band. The residuals are represented in terms of σ (standard deviation of the image).
In the text

	Fig. E.2 Flux ratio between the central binary star surrounded by dust and the companion B star. Our best fit is represented in orange (Tdust = 2200 K, TB = 45 000 K).
In the text

Fig. E.3 Best fits of the flux density from: 1) two Kurucz atmosphere models representing the companion B (green dashed line) and the OB-type inner companion (orange line), 2) one PoWR model representing the WC9 star (blue line), and 3) one blackbody model representing the dusty pinwheel (red dashed line). The flux ratio is obtained by comparison between theWR+OB flux density (grey line) and the companion B. We also show the total flux density (black dashed line) and position of the different filters used (vertical black dotted lines).

In the text