The SPHERE view of Wolf-Rayet 104

Context. WR104 is an emblematic dusty Wolf-Rayet star and the prototypical member of a subgroup hosting spirals that are mainly observable with high-angular resolution techniques. Previous aperture masking observations showed that WR104 is likely an interacting binary star at the end of its life. However, several aspects of the system are still unknown. This includes the opening angle of the spiral, the dust formation locus, and the link between the central binary star and a candidate companion star detected with the Hubble Space Telescope (HST) at 1". Aims. Our aim was to directly image the dusty spiral or"pinwheel"structure around WR104 for the first time and determine its physical properties at large spatial scales. We also wanted to address the characteristics of the candidate companion detected by the HST. Methods. For this purpose, we used SPHERE and VISIR at the Very Large Telescope to respectively image the system in the near-and mid-infrared. Both instruments furnished an excellent view of the system at the highest angular resolution a single, ground-based telescope can provide. Based on these direct images, we then used analytical and radiative transfer models to determine several physical properties of the system. Results. Employing a different technique than previously used, our new images have allowed us to confirm the presence of the dust pinwheel around the central star. We have also detected up to 5 revolutions of the spiral pattern of WR104 in the K-band for the first time. The circumstellar dust extends up to 2 arcsec from the central binary star in the N-band, corresponding to the past 20 years of mass loss. Moreover, we found no clear evidence of a shadow of the first spiral coil onto the subsequent ones, which likely points to a dusty environment less massive than inferred in previous studies. We have also confirmed that the stellar candidate companion previously detected by the HST is gravitationally bound to WR104 and herein provide information about its nature and orbital elements.


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 is 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. lines; and 3.) [WR] 1 -the nuclei of some planetary nebulae, i.e. low-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 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 14 C or 26 Al; Tatischeff et al. 2010).
WR stars are classified according to their emission spectra, which depend on where they fall on the 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 interstellar medium and have extremely high dust-formation rates, with values up toṀ = 10 −6 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 , WR104 , WR118 (Millour et al. 2009b), WR48a, WR112, WR137, WR140 (Marchenko & Moffat 2007;Lau et al. 2017), and the Quintuplet Cluster near the Galactic Center (Tuthill et al. 2006)). These dusty spirals in such massive binary systems are interpreted as the consequence of a collision between the WR wind with 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 , Tuthill et al. 2006, Tuthill et al. 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 . 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 1 the square brackets are used to destinguish these from their massive counterparts 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 twenty years 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.
Thanks to 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, etc.), 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 organised as follows. In Section 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 Section 4 by an analytical model of the dust emission in the spiral, complemented by a preliminary qualitative radiative transfer model. In Section 5, we interpret and discuss our results and finally, in Section 6, we conclude and put this work in perspective.

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 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 YJ bands (0.95-1.35 µm) with a spectral resolution of R = 50, (Claudi et al. 2008).
At most epochs, for all filters, we observed WR104 together with a non-resolved star (either Sgr 4, HD163955 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 3 = 0.40 ± 0.03 mas for 4Sgr and θ UD =   (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 1 or 2 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 multiwavelength study (from 1 µm to 12.4 µm) is shown in Figure A.1.

Data reduction and astrometry
We reduced the SPHERE/IRDIS data using an internallydeveloped Python program, 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 in two parts. Then, all frames are shifted, added, and averaged using dithering information and cross-correlation 4 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 the images to correct for the residual hot pixels.
The IFS data were reduced using the SPHERE Data Reduction and Handling (DRH) automated pipeline (Pavlov et al. 2008) at the SPHERE Data Center 5 (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, crosstalk, 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: (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: 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 of July 2016.

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 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) 6 . 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 socalled 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 −7 . 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 abovementioned 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 −7 ) 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 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)

VISIR observation
VISIR is a mid-infrared imager and spectrograph installed at the Cassegrain focus of Unit Telescope 3 of the VLT (Lagage 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". Note that 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 Section 5.3). We also cannot confirm the spiral shape far from the star because of the relatively lower spatial resolution of the VISIR instrument compared to SPHERE (250 mas 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 data set, especially regarding the PSF reference star, needs to be treated with caution.

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 −2 ) 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 1.08 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 powerlaw 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 χ 2 red 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). Table 3. Power law fits of the average radial profiles. The q parameter corresponds to the power-law index fitted in the selected zone of the radial profile (see text for details). The radial profiles are presented in Fig.D

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 Section 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 spiral represented by the grey, dashed line in the lowerright panel of Fig. 1. Using a cubic interpolation on the pixel coordinates, we compute the flux value 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 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.

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

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 the PSF 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 about the reflective nature of the dust emission. Nevertheless, due to the lack of angular resolution of IFS data, we cannot distinguish if this HeI emission came from the Wolf-Rayet 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. Fig. 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.  We also detect a difference in the position of the photocenter 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 9 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. . Different colour images using Y, HeI, and J bands (blue, green and purple parts of Fig. 6). The photocenters 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 where the two spectra have been extracted in Fig. 6: from the star (inner region) and the pinwheel (between the two circles).

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 HST in 1998 and reported in Wallace et al. (2002). According to Wallace et al. (2002), "the additional component [...] is likely to be physically 9 atmospheric dispersion corrector 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 detection rather than astrometry, we retrieve their HST data (3 filters: 336, 439, 555nm) and recompute the astrometry. Our IRDIS-SPHERE data and these HST data were taken 18 years apart. We fit 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 values d at all epochs and filters is provided in Table 4. The uncertainties given here are the statistical errors scaled to the χ 2 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 et al. 2016) provides a new and precise position measurement of WR104. The 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 Figure 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 Density (SED) of the stars with PoWR 10 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 = 120000 L , see Table 6). Finally, we compute the flux density using our estimation of the distance (D = 2.58 ± 0.12 kpc, see Section 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 use a Levenberg-Marquardt (LM) minimization of the χ 2 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 of approximately 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. 11 A precise determination of this parameter is beyond the scope of the determination of the flux ratio. Fig. 9 shows the χ 2 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 χ 2 map subtracted by the minimum χ 2 value (χ 2 red = 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 (this study), we are able to provide a lower limit on the effective temperature of the companion star such that T B ≥ 33000 K (see Table 5). To further constrain the temperature of this companion, we would require additional UV data so as to account for the peak of the emission of the star.  Fig. 9. Reduced χ 2 maps for the flux ratio fitting. Here, we present ∆χ 2 r = χ 2 r −χ 2 r,min and show the respective 1σ and 2σ confidence intervals with red and purple lines.

Simple Archimedean spiral
The simplest approach to describe the circumstellar environment of WR104 is to use an Archimedean spiral (also called arithmetic 11 We set log(g) to 4.5 to compute the Kurucz models with an effective temperature between 15,000 and 45,000K. Given the high luminosity of companion B compared to the inner binary (factor 4 in luminosity), the companion is likely to be another massive star which is consistent with a high value of log(g). spiral 12 ) to follow the dust distribution around the central binary star . 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. 13 As such, we consider that the dust expansion velocity is set by the WR wind speed, so V dust V ∞,WR , where V ∞ is the terminal wind 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).  Table 3, taking L OB = 80000L and T OB = 30000 K. 3 They use a 1/2 luminosity ratio between WR and O, based on an assumed V band ratio of 1/2 from Crowther (1997).
Reported terminal velocities of WC9 stars range from 1220 km/s (Howarth & Schmutz 1992;Crowther 2007) up to 12 The Archimedean spiral has the property that any ray from the origin intersects successive coils of the spiral with a constant separation distance. 13 Phenomenon occurring in the interaction within the massive binary system, whereby the WR wind is decelerated by the weaker O wind close to the star (Gayley et al. 1997). 1600 km/s (Sander et al. 2012). Taking the lower value of 1220 km·s −1 , we find a distance of D = 2.58 ± 0.12 kpc, which is in good agreement with previous measurements and make WR104 a potential member of the Sgr OB1 stellar association (Lundstrom & Stenholm 1984). Taking the most recent value of 1600 km·s −1 (Sander et al. 2012), we obtain a distance of D = 3.38 ± 0.15 kpc, that would place WR104 behind the Sgr OB1 association.
Future measurements of distances with Gaia, expected with the data releases DR2 (Gaia Collaboration et al. (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.

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. 14 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 and O 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 the WCZ (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: (4) Note that 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 ratio star , arbitrarily set at a wavelength of 1 µm. 15 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): We fixed the number of spiral turns to 3.5 to be able to correctly fit the steps and the sky orientation defined by θ 0 (the spiral orientation at the reference date of April 1998). 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.

Dust formation zone
Step Opening angle 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 nondeconvolved images. The fitting is a least squares calculation based on the pixel values: where O i is the i-th 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 χ 2 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 2000 = 0.45, χ 2 = 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 IRAS catalog (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 θ 0 is given at the reference date of April 1998 (MJD = 57512). Our two epochs allow us to confirm the 241.5 day orbital period of the pinwheel.

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 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 wellestablished 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 near infrared (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, with each one corresponding to a physical hypothesis about the dust nucleation: -Model 1 consists of uniformly thick rings of linearlyincreasing 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 interface far 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 O component 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 overdensity 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 overand 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 16 WC9 component, and a dust-to-gas ratio ξ: 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 χ 2 criterion 17 for comparison. 16 According to the previously reported values of the mass-loss rates (Table 6), the mass-loss of the OB component is negligible compared to the WR mass-loss. The resulting total dust mass differs only by 0.6% when neglecting the OB component. 17 We assume that the uncertainties in the data are dominated by the speckle noise. The errors are determined by computing the unbiased standard deviation for each pixel (Eq. 2). 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. Yet, the WR104 images in the J band are definitely resolved with a significant fraction of dust emission (Fig. C.1, 2nd row at 1.21 µm). For this reason, we do not include the J band in our comparison.
Considering the H and K band data only, the resulting χ 2 red computation for all considered dust-to-gas ratios is shown in Fig.12. It appears that Model 1 (thick rings) provides the best match (χ 2 red = 1.31), with ξ = 1%. Fig. 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 χ 2 red = 2.15 for Model 2 and χ 2 red = 2.25 for Model 3.

Curvilinear profile and shadowing effect
In section 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 (χ 2 red = 1.24). Our model was based on a power-law decrease of the temperature with a constant index (q 2000 = 0.45 ± 0.03) and a dust sublimation temperature of 2000 K. This Article number, page 11 of 20 A&A proofs: manuscript no. aanda-RT   Fig. 13. The 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. 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 10. 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 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 χ 2 red for the models with and without the shadowing effect. The model with a gap seems to be less representative (χ 2 gap = 4.3 compared to χ 2 nogap = 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 MATISSE spectro-interferometer or the planned Extremely Large Telescope (ELT)).

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 value of the different χ 2 red do not allow us to clearly discriminate between the three models of the dust geometry. Fig. 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 data set. 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 Spectral Energy Distribution (SED), which is indicative of circumstellar material. The reference star also presents a very low dynamic range with a poor signal-to-noise ratio, 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-togas ratio. All models favour an intermediate dust-to-gas ratio, between 1 and 10%. This translates into a masses of the dust rings between 5.3 × 10 −8 M and 5.3 × 10 −7 M . For comparison, Harries et al. 2004 set the mass-loss rate ofṀ = 3 × 10 −5 M / yr in their model and determined a dust-creation rate ofṀ = 8 × 10 −7 M / yr, which yields a dust-to-gas ratio of ξ Harries = 2.7%. This leads to a ring mass of 5.3 × 10 −7 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.

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, i.e. 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. 18 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 = 20M 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 = 20M . According to the WR star's spectral type as a WC9 (Sander et al. 2012), we set its mass to M WR = 10M . We treat the central binary star (M WR+O = 30M ) as a single star and compute the orbital period of the companion B star as ≈ 17000 ± 1300 days.
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 fifteenth outer dust coil, produced 10 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, etc.) 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 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. 2017. 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.

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 et al. 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 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 catalog. 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 give 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 Wolf-Rayet 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 the case 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.

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 years. Further, 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 highdynamic 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 signal-to-noise ratio) and a better reference star in the J and K bands (4Sgr 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.     In addition to Section 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.

Appendix B: Detailed Log of observations
As presented in Section 3.4, we used different atmosphere models to determine the effective temperature of the third component of WR104, i.e. companion B. Fig. 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 Sec. 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.
The infrared excess of WR104 compared to the companion indicates the presence of relatively hot dust around the central binary star, i.e. 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 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).