The extreme colliding-wind system Apep: resolved imagery of the central binary and dust plume in the infrared Free

ABSTRACT

The recent discovery of a spectacular dust plume in the system 2XMM J160050.7–514245 (referred to as ‘Apep’) suggested a physical origin in a colliding-wind binary by way of the ‘Pinwheel’ mechanism. Observational data pointed to a hierarchical triple-star system, however, several extreme and unexpected physical properties seem to defy the established physics of such objects. Most notably, a stark discrepancy was found in the observed outflow speed of the gas as measured spectroscopically in the line-of-sight direction compared to the proper motion expansion of the dust in the sky plane. This enigmatic behaviour arises at the wind base within the central Wolf–Rayet binary: a system that has so far remained spatially unresolved. Here, we present an updated proper motion study deriving the expansion speed of Apep’s dust plume over a 2-year baseline that is four times slower than the spectroscopic wind speed, confirming and strengthening the previous finding. We also present the results from high angular resolution near-infrared imaging studies of the heart of the system, revealing a close binary with properties matching a Wolf–Rayet colliding-wind system. Based on these new observational constraints, an improved geometric model is presented yielding a close match to the data, constraining the orbital parameters of the Wolf–Rayet binary and lending further support to the anisotropic wind model.

1 INTRODUCTION

Wolf–Rayet (WR) stars embody the final stable phase of the most massive stars immediately before their evolution is terminated in a supernova explosion (Carroll & Ostlie 2006). They are responsible for some of the most extreme and energetic phenomena in stellar physics, launching fast and dense stellar winds that power high mass-loss rates (Crowther 2007). When found in binary systems with two hot wind-driving components, a colliding-wind binary (CWB) is formed, which may produce observational signatures from the radio to X-rays. Among the wealth of rare and exotic phenomenology associated with CWBs, perhaps the most unexpected is the production of copious amounts of warm dust (Williams, van der Hucht & The 1987; Zubko 1998; Lau et al. 2020b), which can occur when one binary component is a carbon-rich WR star (WC) whose wind provides favourable chemistry for dust nucleation. However, considerable controversy surrounds the mechanism by which some of the most dust-hostile stellar environments known (high temperatures and harsh UV radiation fields) are able to harbour protected dust nurseries that require high densities and relatively low temperatures for nucleation (Crowther 2007; Cherchneff 2015).

A major clue was presented with the discovery of the ‘Pinwheel nebulae’ (Tuthill, Monnier & Danchi 1999; Tuthill et al. 2008) that pointed to dust formation along the wind-shock interface, which is co-rotating with the binary orbit. Subsequent inflation as the dust is carried outwards in the spherically expanding wind results in structures that encode the small-scale orbital and shock physics in the orders-of-magnitude larger plume structures. Images of such systems can reveal the dynamics of the embedded CWB at scales otherwise difficult to resolve. Given the rarity of WR stars in the Galaxy and the brevity of their existence, the discovery and monitoring of these systems provide valuable insight into the evolution of massive stars (van der Hucht 2001). However, only a small handful of such dust plumes have been found to be large and bright enough for imaging with current technologies in the infrared (Monnier et al. 2007).

The recent discovery of ‘Apep’ (2XMM J160050.7–514245) by Callingham et al. (2019) presented a CWB with a spectacular spiral plume visible in the mid-infrared. Extending ∼12 arcsec across, Apep’s dust plume is two orders of magnitude larger than the stereotypical pinwheel nebulae WR 104 (Tuthill et al. 2008; Soulain et al. 2018) and WR 98a (Monnier, Tuthill & Danchi 1999) given their respective distances. Further analysis revealed the system to boast exceptional fluxes from the radio, infrared and X-rays, placing it among the brightest CWBs known and a rival to famous extreme CWB systems such as η Carinae.

Although near-infrared imagery resolved a double at the heart of the system consisting of a companion star north of a highly IR-luminous central component, the 0.7 arcsec separation of the visual double is an order of magnitude too wide to precipitate conditions required to create dust via the Pinwheel mechanism. Instead, the suggestion arose that the system must be a hierarchical triple star with a much closer WR binary driving the dust formation lying hidden within the central component itself.

The most enduring enigma to the underlying physics driving Apep was laid bare in the finding that the spectroscopically determined wind speed of |$3400\pm 200\,$| km s⁻¹ was larger than the apparent expansion speed of the dust plume by a factor of 6, assuming a distance of 2.4 kpc (Callingham et al. 2019). The model attempting to reconcile these conflicting data invoked a strongly anisotropic wind, never previously observed in CWBs, arising from the WR system, suggested to be in the form of a slow equatorial flow and a fast polar wind. If real, such a wind speed anisotropy could be generated by at least one rapidly rotating WR star in the system (Callingham et al. 2019). Such near-critical rotation, in turn, is the essential ingredient required to transmute the ultimate fate of the star as a supernova into a potential long-duration gamma-ray burst (LGRB) progenitor under the currently referenced collapsar model (Woosley 1993; Thompson 1994; MacFadyen & Woosley 1999; MacFadyen, Woosley & Heger 2001; Woosley & Heger 2006). If confirmed, this would constitute the first system of this type yet observed (Callingham et al. 2019).

The detailed and clean geometry of the spiral nebula encodes properties of the inner binary, but is not well reproduced by existing models. Certain features suggest a stage in its history when dust production was turned on and off, while others seem difficult to fit with any level of precision using standard recipes for the expected morphology. Although a plausible illustrative geometric model has been successfully sketched out (Callingham et al. 2019), it fails to reproduce features accurately, with extensive searches of the available degrees of freedom offering little promise. This difficulty seems likely to point to additional model complexity.

To address these open questions, we observed Apep with the European Southern Observatory’s (ESO) Very Large Telescope (VLT) at near- and mid-infrared wavelengths. Observations with the NACO and VISIR instruments are described in Section 2. The detection of the central binary is provided in Section 3 and the proper motion of the dust plume is studied in Section 4. These results are used to constrain a geometric model presented in Section 5, and implications for the distance and wind speed of the system are discussed. The findings of this study are summarized in Section 6.

2 OBSERVATIONS

2.1 NACO observations

Near-infrared imagery of the Apep system was obtained with the NACO instrument (Lenzen et al. 2003; Rousset et al. 2003) on the VLT at Paranal Observatory, with observational dates and basic parameters listed in Table 1. The 2016 epoch consisted of more orthodox (filled pupil) imagery, immediately splitting the central region into a 0.7 arcsec binary denoted by Callingham et al. (2019) as the ‘central engine’ and the ‘northern companion’.

Because of the likely role of the central engine in driving the dust plume, a second epoch of imagery was motivated spanning the J, H, Ks, L^′, and M^′ bands and employing the 7- and 9-hole aperture masks (Tuthill et al. 2010): an observing methodology demonstrated to recover information up to and beyond the formal diffraction limit of the telescope (Tuthill et al. 2000a,b, 2006; Monnier et al. 2007). These observations were executed in the final proposal period of NACO before it was fully decommissioned in October 2019. For each filter and mask combination, observations of Apep and a point spread function (PSF) reference star were interleaved for the purpose of calibrating the interferograms of Apep against the telescope-atmosphere transfer function.

2.1.1 Data reduction

We reduced NACO imagery using software developed within our group (github.com/bdawg/MaskingPipelineBN), which performed background sky subtraction, gain correction, iterative bad pixel interpolation, cosmic ray removal, and interferogram centring and windowing. Bad frames with outlying flux distributions were removed to produce stacked sequences of calibrated image frames. Sample reduced images are displayed in Fig. A1. Data acquired with the M^′ filter were of poor quality – the signals in individual frames were below the threshold needed to perform basic cleaning and were not processed further.

The next step performed by default in the pipeline is to centre and window the stellar image from the larger raw data frame, however, the presence of the relatively close 0.7 arcsec binary made it problematic to extract separately isolated interferograms for both components. We developed two bespoke methods to resolve this issue. One method simply scaled the usual data window function to be very tight, retaining the component of interest and suppressing the remainder (including the unwanted binary component). The second applied a multiplicative window to suppress the unwanted component, after which the data could be handled as normal by the default package. In practice, these two were found to produce similar results and there were no cases where significant structures were found in one but not the other reduction method. An exception to these successful strategies was encountered in L-band data in which diffraction spreads the light sufficiently so that the central engine and northern companion merge and cannot be readily separated. The L-band data were therefore not used for model fitting.

As the spatial structures of the central component and the northern companion are both open to investigation, we performed all three possible calibration combinations: (a) calibrating the central component against the separate-in-time PSF reference star, (b) calibrating the central component against the northern companion, and (c) calibrating the northern companion against the PSF reference star.

Following established practice for masking interferometry data (Tuthill et al. 2000a), we sampled Fourier spectra for each windowed data cube (the sampling pattern depending on the mask and wavelength) that were normalized to yield complex visibility data. We accumulated robust observables comprising the squared visibility and the closure phase (argument of the bispectrum) over the data cube, with the statistical diversity yielding the uncertainties on the derived quantities. In order to extract high-resolution information from the full-pupil (non-masked) data, we performed an identical analysis to these sets of data with pseudo-sampling corresponding to that of an 18-hole mask. Although full-frame speckle interferometry codes might offer superior processing for these data, past success with this strategy of pseudo-masking has delivered successful outcomes and allowed for direct comparison with all data passing through the same reduction pipeline.

We discarded badly calibrated source-calibrator pairs and outlying visibility data for individual baselines upon inspection of squared-visibility plots. Small statistical errors of the closure phase (estimated based on scatter across the cube) are likely an underestimate of systematic errors and can bias subsequent model fitting. We set closure phase uncertainties below the threshold of 2/3 of the median uncertainty to this threshold value to avoid subsequent overfitting to these data points.

2.2 VISIR observations

Mid-infrared imagery of Apep was acquired with the VISIR instrument on the VLT (Lagage et al. 2004). Observations taken in 2016 revealed the ∼ 12 arcsec circumstellar structure in the form of an elaborately detailed spiral plume, while proper motions registered with 2017 data yielded a dust expansion speed of |$570\pm 70\,$| km s⁻¹ at their distance of 2.4 kpc (Callingham et al. 2019). Data from the 2017 epoch suffered from a minor issue in the construction of the observing sequence files: the chop/nod throw setting was too small, so that custom analysis code was needed to disentangle overlapped regions yielding uncontaminated imagery of the plume.

Given the critical nature of the proper motion study in motivating novel physics, we requested a third epoch of VISIR imagery to eliminate any chance for systematic offsets due to the problems with the 2017 observing scripts. This third epoch of VISIR data was taken in 2018, doubling the time-baseline available to validate and refine the previously established proper motion estimates. As a secondary goal, data were also taken with VISIR’s Q3 filter to search for more distant, cooler dust. A log of all epochs of mid-infrared data taken is also given in Table 1.

The reduction of VISIR data followed guidelines provided by the European Southern Observatory’s VISIR manual (de Wit 2019) to perform chopping and nodding subtractions. We calibrated the relative fluxes into units of Jy based on VISIR standard observations of HD133550 (7.18 Jy at J8.9 Filter), HD 178345 (8.21 Jy at B11.7 Filter) and HD145897 (1.84 Jy at Q3 Filter) taken on the same night as Apep observations with the corresponding filters. Given the high level of precision required to extract the modest proper-motion displacements of plume structures, we used a distant faint star, 2MASS J16004953-5142506, about 10 arcsec away from Apep and just bright enough to appear in stacked data to validate the consistency of the plate scale and orientation over the 2-year interval. This yielded an upper limit on the magnitude of any deviations of plate scale of |$0.4{{\ \rm per\ cent}}$|⁠. Final stacked images obtained from the new 2018 data are given in Fig. 1 and are discussed in Section 4.

3 RESULTS: THE CENTRAL BINARY IN THE NEAR-INFRARED

As the engine driving the colliding-wind plume and the origin of the presumed new physics underlying the discrepant wind speeds, the most critical observational challenge was to open a clear window on to the central binary. This was accomplished by extracting spatial structures at the highest possible angular resolutions from the near-infrared data, however, this process is not straightforward. The region has a number of complex elements (photospheres and dust) that vary with wavelength, are only partially resolved, and are not always well constrained by varying data quality over the observational record. Competing models of increasing complexity for structure within the central component were confronted with these data and their plausibility discussed and evaluated in the sections next.

3.1 Close binary star models

3.1.1 Fitting to closure phase data

Starting with the simplest possible model, our binary star consists of two point sources and involves only three free parameters: binary separation, orientation, and flux ratio. We fitted this model using the program binary_grid (github.com/bdawg/MaskingPipelineBN) developed within our group for statistically robust identification of companions, particularly those at high contrast ratios. The most likely model is that which globally minimizes the χ² error between the model and closure phase data, found by sampling a grid covering the parameter space and subsequently optimizing with the Powell method (Powell 1964). Given the Fourier coverage of our data, we estimated detection limits at confidence intervals of 99 per cent and 99.9 per cent using a Monte Carlo method. The code injects a large number of noise-only (no binary) simulations into the algorithm to estimate the likelihood of false binary detections across the separation–contrast parameter space, which are used to define the corresponding confidence thresholds.

The best-fitting solutions for the central binary largely (with one exception) calibrated against the PSF reference star [case (a) in Section 2.1.1] are shown in Table 2. The most certain detections produced by binary_grid exceeded 99.9 per cent confidence (for the IB_2.24 filter), firmly establishing the reality of a close binary star within the central engine. Furthermore, recovered binary parameters from Table 2 are largely consistent across epochs and across filters despite some examples of more scattered values. For example, visual inspection of the 2019 epoch revealed that the J-band data were of compromised quality. Taken as a whole, the near-infrared data affirm the existence of a binary with the fainter component ∼47 mas (K band) to the West.

The alternate strategy [(b) in Section 2.1.1] of calibrating the central binary against the northern companion yielded consistent results in the H and K bands to those shown in Table 2. Outcomes were also found to be robust under different choices of the windowing strategy to separate the 0.7 arcsec binary (also discussed in Section 2.1.1). Both findings further consolidate confidence in the detection made in this study.

On the other hand, calibrating the northern companion against the PSF reference star [(c) in Section 2.1.1] yielded inconsistent, low-confidence binary parameters. This argues for the northern companion being an isolated single star at the scales of contrast and spatial resolution probed by this study.

3.2 Models including embedded dust

3.2.1 Fitting to visibility data – 1-component model

While previous fits relied solely on closure phase data (as is customary in companion detection with non-redundant masking, Kraus et al. 2017), the squared visibilities show a drop-off towards longer baselines at all position angles, requiring the existence of resolved circumstellar matter. The most natural scenario is local dust in the immediate vicinity of the binary. In this and the subsections that follow, emphasis in the modelling is given to the K-band data for which the quality is significantly higher resulting in more tightly constrained fits, and for which thermal flux from the dust becomes more pronounced.

We note in passing that the addition of dust also considerably complicates interpretation of the nature of the binary. The apparent fluxes may be influenced by local non-uniform opacity and/or thermal re-radiation from dust, so that the contrast ratios presented in the previous section may not correspond directly to the two stars. This might permit, for example, a more nearly equal binary in luminosity as anticipated for a WR + WR (Callingham et al. 2019, 2020) despite the contrast data in Table 2.

In order to account for circumstellar material radiating in the near-infrared, we invoked models including spatially resolved components, the simplest being a circular Gaussian flux distribution. Fits obtained with this model provided an estimate of the overall dimensions of the resolved component and are tabulated in column (a) of Table 3. The J-band data were not used for visibility-only fitting due to poorer phase stability delivered by AO at shorter wavelengths. Fitting to the visibilities in the H and K bands with a simple Gaussian model yielded a full-width at half maximum (FWHM) of approximately 36 ± 4 mas, so that the spatial extent of the dust is comparable to the separation of the binary resolved in the previous section.

3.2.2 Fitting to visibility data – 2-component models

In order to gradually add complexity to our model, a useful intermediate step is the combination of a resolved Gaussian dust component with an unresolved, embedded photosphere centred at the peak. Confronting such a model to the data generally resulted in fits where the majority of the flux arose from the dust. For example, the data with the highest signal-to-noise ratio (narrow K band) yielded a 2-component model with a best-fitting value of 73 per cent of the total flux attributed to the resolved Gaussian.

To indicate the impact of this on the interpretation, we adjust our earlier ∼6:1 K-band binary flux ratio now recognizing the brighter component flux should be split two ways: a stellar photosphere and resolved dust shell. Performing this decomposition, we find that the embedded point sources are left with an approximately equal flux ratio (∼14 per cent of the total flux in each): a result that can be more neatly brought into accord with expectations from spectroscopic evidence (Callingham et al. 2019, 2020).

3.2.3 Fitting to visibility and closure phase data – 2-component models

With an estimate of the relative contribution from the resolved component, we experimented with more complex multicomponent models built with the interferometric data fitting program, LITpro (Tallon-Bosc et al. 2008). We replicated the model consisting of two simple point sources using LITpro’s gradient descent algorithm, now fitting to both the visibility and closure phase simultaneously in the H and K bands. This work successfully reproduced the binary detections presented in Table 2.

Simultaneous fitting to visibility and closure phase data was repeated with the 9-hole J-band data (previously flagged as of marginal quality). Given best estimate starting parameters, the algorithm yielded a roughly equal binary with a flux ratio of 1.0 ± 0.4. This appears to agree with the local dust hypothesis, a consequence of which is that the dominant contribution to the luminosity should shift from dust to photospheres as the observing band moves from H/K to J band. Indeed, although the poor data quality resulted in unconstrained fits with the simpler Gaussian source model, the drop in the squared visibilities towards longer baselines in the J band appears much less pronounced than in the H and K bands, suggesting a much reduced dust fraction and a more unresolved morphology.

3.2.4 Challenges with higher degrees of freedom

Physical expectations imply a 3-component model with two embedded point sources and a third resolved dust component (approximated here as a circular Gaussian, but it may also exhibit spatial structure). A model in which all three components can vary in flux, position and size (for the dust) introduces at least 7 degrees of freedom; parameters are also expected to vary with observing band and (possibly) over time.

Given limited and variable data quality, experimentation with models boasting such significant numbers of floating parameters were generally found to diverge or yield unreliable outcomes. With a little perspective, this is hardly surprising. All structures of interest are at or below the commonly stated diffraction limit of the instrument and are modelled with highly covariant parameters.

To mitigate this problem, a conservative approach was taken, permitting only the minimum free model parameters and strictly delimiting possible options based on astrophysical plausibility and outcomes of earlier fits with simpler models. Reliance on underlying physics eases the problem of requiring sufficient data to uniquely constrain the model.

3.2.5 Fitting to visibility and closure phase data – multicomponent models

To constrain our physically motivated 3-component (two stars plus dust shell) model, we begin by adopting the 47 mas binary star separation and 274° orientation parameters from the K-band fit in Table 2. These 2-component closure-phase-only fits should be only weakly influenced by the dust, and largely constant with wavelength and epoch. If we further adopt the idea that we are most likely dealing with a near-equal WR–WR binary (flux ratio ∼1; Callingham et al. 2020), then we may use the 14 : 14 : 72 per cent flux partition established in Section 3.2.2. Having established reasonable parameters of constraint for our model, the degree of freedom that remains is the spatial extent of the local infrared-bright hot dust: a quantity that would otherwise be entangled with ambiguities of the correct division of flux between dust and photosphere, and also with the true geometry of the binary. Fitting only to the FWHM of the Gaussian, the results of the χ² minimization are shown in Column (b) of Table 3. The H-band 9-hole data could not be fitted by LITpro since the algorithm was unable to converge for poorly calibrated visibility data in this set.

Compared to the single component Gaussian model fit to visibility data only, the FWHM of the Gaussian in the 3-component model is slightly lower, which is likely due to the added luminous point source displaced from the centre of the source distribution. The apparent increase in Gaussian dust FWHM in the J band is not interpreted as the physical extent of the dust, but is likely due to the assumed flux partition beginning to break down at short wavelengths.

3.2.6 Best near-infrared 3-component model

In summary, the complex models discussed above demonstrate the likely existence of local dust extending over a spatial scale comparable to the binary separation. Synthesizing information from model fitting and constraints based on basic stellar astronomy, we arrive at a plausible best-case model derived predominantly from K-band data to comprise an approximately equal binary separated by 47 ± 6 mas embedded in a dust shell fit by a 35 ± 7 mas-FWHM Gaussian. The position angle of the binary was determined to be 274 ± 2° in April 2016 and 278 ± 3° in March 2019. The binary photospheres rise to dominate the total flux in J band, while conversely further in the infrared at L band and beyond the dust rapidly becomes the dominant contribution, however declining data quality unfortunately presents challenges for the modelling at both.

4 RESULTS: THE EXTENDED DUST PLUME IN THE MID-INFRARED

The second set of observational data listed in Table 1 comprises a campaign of mid-infrared imaging with the VISIR instrument conducted to confirm the observed proper motions as well as to improve upon the derived parameters for the outflow.

The displacement of the dust plume is modest even over the 2-year time baseline now available, so that custom image processing procedures and algorithms were required for precise registration of the small (of order 1 pixel yr⁻¹) motions betraying wind-driven inflation of the circumstellar dust structure. We implemented an edge-detection algorithm (christened ridge_crawler) that has the capability locating edge-like features at the sub-pixel level. Images of Apep’s dust plume were high-pass filtered to accentuate edge-like features and suppress regions of constant flux, effectively producing images of the skeleton of the spiral plume as shown in Fig. 2. We then applied the algorithm, ridge_crawler, described in detail in Appendix  B, to extract the positions of the dust plume’s edges, which were used to determine the expansion speed of the plume.

4.1 The expansion rate of the dust plume

The filtered skeleton of Apep highlights a number of prominent but disjoint edges as shown in Fig. 2, requiring each segment to be sampled and analysed individually. Table 4 summarizes the results returned by ridge_crawler calculated with a pixel scale of 45 mas pixel⁻¹ (de Wit 2019).

The first point to note is that the 2018 data confirm the dust to exhibit an anomalously low proper motion. All values recovered in Table 4, even the most extreme, are a factor of several too low to reconcile the angular rate of inflation on the sky with the line of sight wind speed of |$\sim 3400\,$| km s⁻¹ from spectroscopy at any plausible system distance, which was estimated to be ∼2.4 kpc (Callingham et al. 2019).

At a greater level of detail, the values exhibit scatter within a single epoch and between epochs larger than errors, even accounting for the intrinsically small level of the signal (less than VISIR’s diffraction limit of ∼ 0.2 arcsec in the J8.9 filter (de Wit 2019)). Turning firstly to the latter, differences in mean expansion rate between the three combinations of epochs are unlikely from any real behaviour of the plume (such as an acceleration). We attribute these to an artefact introduced by the custom image processing procedures used to recover images from the 2017 epoch (see Section 2). Given the consistent high-quality data for both the 2016 the 2018 epochs, the values derived from this combination are therefore considered most reliable.

On the other hand, some real diversity of speeds from ridges a to f within any single epoch is to be expected as different edges show differing apparent expansion rates projected on to the plane of the sky. While 2D images of an optically thin 3D structure do result in most of the prominent edges close to the sky-plane (the effect of ‘limb brightening’), departures are significant enough to cause measurable diversity as described in Appendix B3.

To systematically study such projection effects, we constructed a model of the expanding dust plume (discussed in Section 5.3) accounting for its 3D structure, and calibrated the apparent motion of the edge-like structures identified above against the true (physical) expansion speed. Details of this calibration procedure are described in Appendix B3. Our results suggest that ridges a, c, d, and f (group I) expand at approximately the true expansion rate of the system, while ridges b and e (group II) suffer apparent slowing due to projection effects. We therefore estimate the true expansion speed of the dust plume using the mean speed of group I ridges between the 2016 and 2018 epochs, obtaining a final value of 80 ± 11 mas yr⁻¹.

This new estimate corresponds to a linear expansion speed of approximately |$910\pm 120\,$| km s⁻¹ assuming the previously published distance of 2.4 kpc (Callingham et al. 2019). At this distance, the dust expansion speed is then approximately four times slower than the spectroscopically determined wind speed. With the inclusion of a 2-year time baseline, this analysis has therefore confirmed a consistently slow dust expansion. Since no current model or observation of CWBs can explain such a drastic discrepancy between the outflow speed derived from the proper motions and that from spectroscopy, this finding strengthens the argument for an alternative model that deviates from the traditional picture of CWBs. One such model is discussed in Section 5.5.

4.2 Spectral energy distribution

The multiple epochs of Apep imagery across a number of infrared spectral bands provide information on the spectral energy distribution (SED) of the system. Building on the SED presented by Callingham et al. (2019), we implemented custom procedures to extract the SEDs of individual components within Apep using NACO and VISIR data, namely the central binary, northern companion and extended dust plume.

With half of the NACO data volume targeting calibrator stars whose photometric properties are well determined, we extracted spectral flux densities of the central binary and northern companion by calibrating pixel values of background-subtracted, windowed sub-images of either component against those of the the corresponding calibrator star. NACO imagery also placed an upper bound of the flux contribution from the extended dust plume of below 5 per cent, which is not visible at these wavelengths at the SNR of the instrument.

We used VISIR imagery to estimate the relative flux contributions from the central binary using PSF-fitting photometry. Since the northern companion lies on the shoulder of the much stronger mid-infrared peak from the central region, we isolated its relative flux by first subtracting the azimuthal average of the main peak. We attributed the remaining background-subtracted flux to the dust plume. The final SED is presented in Fig. 3.

We note that Apep exists in a complex region with extended background emissions (as revealed by the Spitzer 8 μm image, Werner et al. 2004). The HiGAL 70-micron (Molinari et al. 2016), AKARI N60 (Doi et al. 2015), and MSX (Price et al. 2001) surveys had effective beam sizes larger than Apep’s imaged dust plume, and so are vulnerable to contamination from background emission or cooler Apep emission from earlier episodes of dust production. We omitted WISE data (Wright et al. 2010) from the SED due to the low SNR despite their use of profile-fitted magnitudes to correct for effects of saturation.

The data appear to suggest variation in the mid-infrared flux over time. The MSX (c. 1996) and AKARI (c. 2006) values are both larger than the VISIR flux by a factor of a few across the mid-infrared. For example, the N-band flux decreases from MSX to AKARI (δt ∼ 10 yr) and to VISIR (δt ∼ 20 yr), which is consistent with the interpretation of an expanding and cooling dust plume. The measurements at shorter wavelengths were consistent with NACO fluxes, suggesting that the underlying flux heating the dust did not significantly vary over the interval covered by the observations. Future work is required to explore the mid-infrared variability of the system in more detail and evaluate its consistency with variations in the plume geometry.

5 DISCUSSION

5.1 Detection of Apep’s central binary

The results in Section 3 establish that Apep’s ‘central engine’ (Callingham et al. 2019) consists of a close binary with a secondary component approximately 47 ± 6 mas (K-band mean) to the West of the primary component. This result rules out the possibility of the secondary component being a neutron star or a black hole – which could conceivably drive a spiral plume by gravitational reflex motion, analogous to that witnessed in some mass-losing red giant systems such as AFGL 3068 (Mauron & Huggins 2006) – or that of a single-star dust production mechanism driving the prominent spiral nebula. This conclusion is supported by the findings of Callingham et al. (2020), in which an analysis of the visible-infrared spectrum taken with XSHOOTER on the VLT found that the central engine consists of two WR stars of subtypes WC8 and WN4-6b, respectively, identifying Apep as the first double classical WR binary to be observationally confirmed.

Furthermore, Marcote et al. (in preparation) were able to image the wind-collision region of the central binary at radio frequencies with very-long-baseline interferometry using the Australian Long Baseline Array (LBA) in July 2018. In addition to observing the colliding-wind shock front at the expected location of Apep’s central binary, the authors derived a best-fitting binary position angle of 277° based on the geometry of the shock front. The best-fitting binary orientation angles of 274 ± 2 (April 2016) and 278 ± 3 (March 2019) that we derived in Section 3 based on the NACO H/K-band data are in close agreement with the results of Marcote et al. (in preparation) based on the LBA data and confirm the CWB interpretation.

Although the detection of a change in the binary position angle is only marginally significant given the errors, the direction of change is indicative of a counter-clockwise orbit is consistent with the expected direction given the geometry of the spiral plume and of the right magnitude, estimated to be 2 ± 1° over interval, given the expected binary period (Section 5.3).

The multicomponent model fitting results also suggest the existence of local dust on a spatial scale comparable to the binary, which complicates the interpretation of the fitted flux ratios between the two stellar photospheres. The existence of a third component is supported by the SED (Fig. 3) of the central binary, in which there appears to be a large contribution from a component at a temperature significantly lower than that expected of stellar photospheres.

Given the existence of local hot dust, the fitted contrast ratios of the double point source model presented in Table 2 is a reflection of the sum of the dust and primary stellar components’ fluxes compared to the secondary stellar component’s flux. The K-band absolute magnitudes of WC8 and WN4-6b stars are typically −5.3 ± 0.5 and −4.6 ± 0.7, respectively (Rate & Crowther 2020). Given the spectral types classified by Callingham et al. (2020) and the positional assignment of the two stars based on Marcote et al. (in preparation; assuming the WN component launches a larger wind momentum), we expect the flux ratio of the two photospheres to be between 0.2 and 1.6. To achieve a contrast ratio of 6 ± 2 in the K band, we require dust centred near the WC component to contribute approximately between 50 per cent and 80 per cent of the total flux. The estimate provided in Section 3.2.2, though indicative only, is consistent with this range expected from typical WR fluxes.

The dust emission appears to diminish in flux towards J band, though we interpret the best-fitting values with caution due to the marginal data quality of this subset of data. The existence of dust around the WC8 star also appears to be supported by the observation that the C iv 2.08 μm and C iii 2.11 μm lines from the WC8 in the SINFONI spectrum of the central binary (Callingham et al. 2019) appear to be more diluted than the He ii 2.189 μm line primarily associated with the WN star.

5.2 Dust properties

The total measured J8.9 and Q3 fluxes from Apep are F_J8.9 = 18.09 ± 1.81 Jy and F_Q3 = 38.75 ± 7.75 Jy, where |$\sim 8.4{{\ \rm per\ cent}}$| of the J8.9 flux arises from the central system. The emission from the central system in the Q3 filter is negligible (⁠|$\lesssim 5{{\ \rm per\ cent}}$|⁠) compared to that of the overall nebula. Reddening by interstellar extinction is corrected using the ISM extinction law derived by Mathis (1990) and the estimated extinction towards Apep of A_V = 11.4. The extinction correction factors for the J8.9 and Q3 fluxes are 1.48 and 1.26, respectively. From the dereddened fluxes, we derive an isothermal dust temperature of |$T_d = 227^{+20}_{-15}$| K and dust mass (Equation 2) of |$M_d = 1.9^{+1.0}_{-0.7}\times 10^{-5}$| M_⊙ in the spiral plume assuming a distance of 2.4 kpc towards Apep. Assuming this dust was formed in ∼38 yr, the dust production rate for Apep is |$\dot{M}_d \sim 5\times 10^{-7}$| M_⊙ yr⁻¹.

5.3 Geometric model

In order to develop a physical model of Apep, we first construct a geometric model for the dust plume. Callingham et al. (2019) presented a plausible geometric spiral model based on the standard Pinwheel mechanism (Tuthill et al. 2008), matching the gross features in the data. However, on close comparison, the locations of edges and sub-structures provided a relatively poor match when overlaid. Despite fairly extensive exploration of the parameter space, varying the spiral-generative parameters did not offer significant improvement.

Here, we present a new modelling effort aimed at creating a geometric model systematically linking the spiral nebula of Apep to the underlying binary. We expand upon previous generative models, adding eccentricity as an extra degree of freedom to the binary orbital parameters. Details of model construction are described in Appendix  C.

5.3.1 Fitted model

A geometric model fitted to the mid-infrared data is presented in Fig. 4, in which the structural skeleton of the new model closely reproduces that of Apep’s spiral plume extracted from the mid-infrared images.

The best-fitting parameters of the model are presented in Table 5, where the wind speed of 80 ± 11 mas yr⁻¹ was constrained by the proper motion study in Section 4. The model is also consistent with the near-infrared position angle of the binary discussed in Section 5.1, which we used to constrain the true anomaly at the time of observation.

Along with assumptions of Apep’s distance (discussed in Section 5.4) and near-infrared angular separation that together constrain the semimajor axis, the orbit of Apep’s central binary is almost uniquely constrained. The remaining two-fold degeneracy arises from the sign of the inclination which the geometric model (based on an optically thin surface) does not constrain. In either case, the low inclination of the binary implies a mostly polar view on to the spiral plume.

5.3.2 Eccentricity and episodic dust production

Most notably, a relatively high orbital eccentricity of 0.7 ± 0.1 was required to fit the model of Apep. There exist examples of persistent dust producing WR stars (Williams & van der Hucht 2015) in which a prehistory of tidal interactions leading to orbital circularisation has been invoked to explain the observed lack of dependence of dust production on orbital phase (Williams 2019). Circular orbits have been confirmed for CWBs such as WR 104 (e < 0.06; Tuthill et al. 2008) and WR 98a (Tuthill et al. 1999) based on the circular geometry of the Pinwheel spiral.

On the other hand, the archetype episodic dust producer WR 140 with a very high orbital eccentricity (e ≃ 0.9) (Fahed et al. 2011) does not form the classic spiral geometry, owing to the very brief duration of dust formation near periastron. Apep, like WR 140, is too wide for tidal circularization, even during past giant phases. Although we have no historical record of variations in infrared flux (which would require time baselines of ∼125 yr), the central binary of Apep seems a much closer analogue to the WR 140 system than to persistent dust producers such as WR 104.

Models for episodic dust producers such as WR 140 explain the periodicity of dust production in terms of the binary’s orbital phase, with periastron approaches creating wind collisions strong enough to mediate dust nucleation (Williams et al. 1990; Usov 1991). Our new model for Apep’s plume suggests that the production of dust is indeed centred near periastron, with the central true anomaly of the dust-producing zone offset by ∼18° after periastron (or ∼ 1 yr). The dust, which spans ∼264° over the orbit, was produced over a period of ∼38 yr between 1936 and 1974.

In light of the geometric model, we can clearly identify the first and final dust rings of the spiral plume. The ratio of their radii (∼6.5 and ∼12 arcsec) are consistent with that between the number of years ago that they were produced at the time of observation. Assuming episodic dust production, our model also predicts the next cycle of dust production to begin at around 2061, with the next periastron passage occurring at around 2071.

The distinct ring-like features corresponding to the turning on and off of dust production appears to support the threshold effect of dust production as seen in WR 140 when conditions of nucleation are briefly reached near periastron (Williams et al. 2009). Observationally, this threshold effect may be confirmed from the detection of more distant dust from a prior dust production cycle concentric to the prominent spiral plume of Apep. Expectations for the location of the previous coil of the plume are given in Fig. C1. Future studies may benefit from instruments such as ALMA to search for distant dust structures expected from the threshold effect.

Another observational confirmation of episodic dust production may come from the production of dust from a newer cycle close to the central binary. Based on the expected time frame of the upcoming episode of dust production, it is unlikely that hot, local dust as part of a new cycle of dust spiral via the Pinwheel mechanism is responsible for the bright mid-infrared peak at the central engine. Alternative dust production mechanisms, such as those leading to eclipses observed in WR 104 (Williams 2014), may instead be responsible.

Observations with instruments capable of delivering higher angular resolution (e.g. SPHERE or MATISSE on the VLT) may reveal the nature of the local dust and search for the potential upcoming cycle of colliding-wind mediated dust production.

5.3.3 Opening angle

We also note that the opening angle of 2θ_{w, IR} = 125 ± 10° derived from this geometric model of Apep is smaller than the value of 2θ_{w, radio} = 150 ± 12° derived by Marcote et al. (in preparation) based on the geometry of the radio shock front at the wind-collision region. In interpreting this apparent discrepancy, it is important to note that the two values, 2θ_{w, IR} and 2θ_{w, radio}, may be reflective of the shock structure at different orbital phases of this eccentric binary. Specifically, 2θ_{w, IR} is a measurement of the opening angle when the shock was radiative to allow for efficient cooling and hence dust production. Following dust production centred near periastron, the shock structure may have been gradually modulated as the binary moved further apart, potentially resulting in the shocked region flaring out from the contact discontinuity as the wind shock becomes adiabatic (Pittard & Dawson 2018).

According to theoretical relationships derived by Gayley (2009), a radiative shock with opening angle 2θ_{w, IR} derived in this study is consistent in terms of wind-momentum ratio, η, with an adiabatic shock with mixing having an opening angle 2θ_{w, radio} derived by Marcote et al. (in preparation). The radiative opening angle would then imply a wind momentum ratio of η_IR = 0.21 ± 0.07 based on the relationship derived by Gayley (2009, equation 9), a value lower than η_radio = 0.44 ± 0.10 (Marcote et al., in preparation) calculated using the adiabatic opening angle.

The orbital modulation may hence explain the broadening of the shock region reflected in the 2θ_w,radio measurement, and is possibly associated with the brightening of the 843-MHz radio flux reported by Callingham et al. (2019), which has so far lacked a satisfactory explanation. However, this suggestion requires detailed modelling in future studies to verify.

It is worth noting that the geometric model sketched above is intended to serve as an argument for self-consistency and plausibility of the underlying physical parameters. The results provide evidence that the morphology of Apep’s spiral can indeed be produced via the Pinwheel mechanism in a CWB, and that the underlying binary is likely in an eccentric orbit. The geometric model paves the way to a full radiative transfer model that may more accurately represent the variations in brightness across the spiral plume, adding a layer of physical modelling on to what appears to be a sound geometrical construct.

5.4 Enclosed mass and distance to Apep

Spectroscopic studies have estimated that the distance of Apep is |$2.4^{+0.2}_{-0.5}$| kpc (Callingham et al. 2019), which was more recently revised to |$2.0^{+0.4}_{-0.3}$| kpc (Callingham et al. 2020). If we were to assume the most recent distance estimate and the binary separation discussed in Section 3, we derive a physical WR binary separation of |$95^{+19}_{-14}$| au at the time of NACO observations (corrected for projection). Together with constraints provided by the geometric model, the orbit of the central binary is uniquely constrained with a semi-major axis of |$56^{+24}_{-17}$| au. Using a Keplerian model, we immediately estimate the enclosed mass of the WR binary to be |$11^{+13}_{-6}$| M_⊙.

It is understood that the mean mass of WC8 stars is around 18 M_⊙ (Sander et al. 2019), although variations over a range of about 7 M_⊙ (Crowther 2007) are possible. WN stars have a much larger mass range with extreme cases up to 80 M_⊙, but with the mean values typically under 20 M_⊙. Given the spectral classification determined by Callingham et al. (2020), we expect the enclosed binary mass of Apep to lie between 20 and 40 M_⊙.

The mass estimate relying on a distance of |$2.0^{+0.4}_{-0.3}$| kpc is therefore lower than the most likely mass range. However, adopting the earlier |$2.4^{+0.2}_{-0.5}$| kpc (Callingham et al. 2019) estimate would imply a semimajor axis of |$67^{+20}_{-23}$| au and an enclosed mass of |$19^{+12}_{-12}$| M_⊙, which exhibits a larger overlap with the expected mass range. This analysis appears to argue for a distance towards the upper end of the spectroscopically determined range. We hereon adopt the |$2.4^{+0.2}_{-0.5}$| kpc value for subsequent calculations.

Confirmation of the ∼2.4 kpc distance also argues for a common distance for the WR binary and the northern O-type supergiant in a hierarchical triple. Although chances of finding such field stars aligned within an arcsecond seems remote, the apparent lack of participation of the northern companion either in the radio or infrared has motivated debate on whether the stars are physically associated or not (Callingham et al. 2020).

5.5 The wind speed dichotomy

Adopting the |$2.4^{+0.2}_{-0.5}$| kpc system distance (Callingham et al. 2020), the linear expansion speed of the dust translates to |$910\, ^{+210}_{-290}\,$| km s⁻¹. Comparing this value to the spectroscopic wind speed, this study supports the original finding that we are indeed observing a dust expansion significantly slower than either spectroscopic wind speed of the component WR stars: |$3500\pm 100\,$| km s⁻¹ or |$2100\pm 200\,$| km s⁻¹, launched by the two stellar components of the WR binary (Callingham et al. 2020). The underlying conundrum presented by Apep in which the plane-of-sky (proper motion) expansion is starkly at odds with well-established line-of-sight (spectroscopic) windspeed is therefore reaffirmed.

Furthermore, we note that the Apep system seems unique in this behaviour: other Pinwheel systems including prototypes WR 140 (Williams et al. 1990, 2009), WR 104 (Tuthill et al. 1999; Harries et al. 2004; Tuthill et al. 2008; Soulain et al. 2018), and WR 112 (Lau et al. 2017, 2020a; previously suspected as displaying ‘stagnant shells’) all exhibit dust motions in accord with expectations from spectroscopy.

We therefore encourage ideas from the stellar wind community: this conundrum now seems shored up by firm data on the winds, the identification of the stellar components and the distances. So far the main contender remains the anisotropic wind model proposed by (Callingham et al. 2019), depicted in Fig. 5, where the WR binary is proposed to be capable of launching wind speeds varying across latitude. The central WR star drives a significantly slower wind in the equatorial direction than in the polar direction (which is probed by spectroscopy), resulting in dust formation embedded within the slower terminal wind in the equatorial plane.

This may indicate rapid stellar rotation – a phenomenon associated with the generation of anisotropic winds (Groh et al. 2010; Steffen et al. 2014; Gagnier et al. 2019). Identifying Apep as containing a rapidly rotating WR has broader astrophysical significance. Under the collapsar model, rapidly rotating WR stars are likely progenitors of LGRBs (Woosley 1993; Thompson 1994; MacFadyen & Woosley 1999; MacFadyen et al. 2001; Woosley & Heger 2006). The anisotropic wind model would make Apep as a potential candidate for an LGRB progenitor, the first such system to be discovered in the Galaxy.

Future work may also like to explore the alternative theory of magnetic confinement which may be capable of slowing magnetic equatorial winds (Ud-Doula 2003).

6 CONCLUSIONS

The analysis of near-infrared NACO imagery primarily in H/K-bands revealed the structure of the dust engine embedded in the centre of Apep’s Pinwheel nebula for the first time, delivering the detection of a binary with a 47 ± 6 mas apparent separation. The binary at the centre of Apep is consistent with spectroscopic expectations (Callingham et al. 2020), and the observed position angles of 274 ± 2° (April 2016) and 278 ± 3° (March 2019) were robust against the choice of model and are in close agreement with radio observations for the colliding-wind shock front (Marcote et al., in preparation). The change in orientation between the two epochs, although of only modest statistical significance, is in the right direction and of the right magnitude for the expected binary orbit. Further towards the thermal-infrared (from L band onwards), flux from local hot dust rose to dominate the energetics within the region of the central binary.

Our new geometric model incorporating eccentricity offers an improved fit to the details of the dust plume compared to previously published work, constraining the orbital parameters and distance to the CWB. The model favours an orbital period of 125 ± 20 yr, an eccentricity of 0.7 ± 0.1 and a shock cone full opening angle of 125 ± 10°, and suggests that the binary is presently near apastron. Assuming Apep is an episodic dust producer, the best-fitting orbital parameters predict that the upcoming episode of dust production will begin in ∼41 yr and continue for ∼38 yr.

Using three epochs of mid-infrared VISIR imagery spanning a 2-yr time baseline, improved proper motion studies extracted an updated 80 ± 11 mas yr⁻¹ expansion speed, translating to approximately |$910\, ^{+210}_{-290}\,$| km s⁻¹ at a distance of |$2.4^{+0.2}_{-0.5}$| kpc (Callingham et al. 2019) favoured by orbital constraints. Confirmation of the observational dichotomy between the slowly expanding dust plume despite spectroscopically confirmed fast radial winds (Callingham et al. 2020) provides support for the anisotropic wind model which puts forward Apep as the first potential LGRB progenitor candidate to be identified in the Galaxy.

The authors encourage further observational and theoretical progress on this fascinating system. New data with the ability to reveal presently hidden structural detail (both close to the central binary and beyond the main mid-infrared dust plume), such as SPHERE, MATISSE, or ALMA imaging, may be particularly useful, as is radiative transfer modelling of the plume and theoretical exploration of the colliding-wind physics with anisotropic mass-loss.

SUPPORTING INFORMATION

Apep.gif

Please note: Oxford University Press is not responsible for the content or functionality of any supporting materials supplied by the authors. Any queries (other than missing material) should be directed to the corresponding author for the article.

ACKNOWLEDGEMENTS

This work was performed in part under contract with the Jet Propulsion Laboratory funded by NASA through the Sagan Fellowship Program executed by the NASA Exoplanet Science Institute.

YH, PT, and AS acknowledge the traditional owners of the land, the Gadigal people of the Eora Nation, on which the University of Sydney is built and this work was carried out. JRC thanks the Nederlandse Organisatie voor Wetenschappelijk Onderzoek (NWO) for support via the Talent Programme Veni grant. BJSP acknowledges being on the traditional territory of the Lenape Nations and recognizes that Manhattan continues to be the home to many Algonkian peoples. We give blessings and thanks to the Lenape people and Lenape Nations in recognition that we are carrying out this work on their indigenous homelands. BM acknowledges support from the Spanish Ministerio de Economía y Competitividad (MINECO) under grant AYA2016-76012-C3-1-P and from the Spanish Ministerio de Ciencia e Innovación under grant PID2019-105510GB-C31.

This research used NASA’s Astrophysics Data System; LITpro, developed and supported by the Jean-Marie Mariotti Centre (Tallon-Bosc et al. 2008); the ipython package (Pérez & Granger 2007); scipy (Jones et al. 2001); numpy (Van Der Walt, Colbert & Varoquaux 2011); matplotlib (Hunter 2007); and astropy, a community-developed core python package for Astronomy (Astropy Collaboration 2013).

DATA AVAILABILITY

The data underlying this article are available on the ESO Archive under observing programmes 097.C-0679(A), 097.C-0679(B), 299.C-5032(A), 0101.C-0726(A), and 0102.C-0567(A).

REFERENCES

APPENDIX A: NEAR-INFRARED IMAGES

A series of cleaned near-infrared images of Apep observed with NACO are presented in Fig. A1. Both the central binary (centred) and northern companion appear in each frame, the relative flux between which can be seen to vary across the spectral bands.

APPENDIX B: EDGE DETECTION

B1 Displacement visualization

In order to register displacements of image components over time, specific features or structures must be identified and tracked. High-pass filtering is a useful tool in this context, enhancing features such as edges, essentially revealing the skeleton of Apep’s dust. Such underlying structural elements of the plume, illustrated in Fig. B2, can be accurately registered at each epoch. Even without further processing, the expansion of the plume may be directly confirmed in a simple difference image between epochs as given in Fig. B1. This reveals that the edges of the spiral plume from the older epoch of all pairs are always completely contained within those of the newer epoch, implying a consistent outwards motion of the plume over time. Although such confirmation is helpful, it is difficult to extract robust estimates of the rate of expansion from such difference images; an issue that motivated the more advanced procedures discussed next.

B2 Displacement extraction algorithm

Given the small dust displacements of Apep’s plume of order one pixel per year, we implemented a custom algorithm, ridge_crawler to locate edges at the sub-pixel level.

The algorithm is initialized with a user-specified location and direction, from which the it begins taking discrete steps along the crest of the ridge in the image data. At each iteration, ridge_crawler looks n pixels ahead in the direction interpolated from its previous path, takes the profile of a slice of pixels in the orthogonal direction, and fits a parabola through the three brightest pixels along that slice. ridge_crawler then calculates the location of the peak of the parabola, records it as the next sample along the ridge, and moves to this location ready for a new iteration. The action of this algorithm is depicted in Fig. B2.

The algorithm fits spline functions to the sample points of each ridge and densely samples the splines. It then projects rays from the centre of the image, intercepting the splines at both epochs, calculating the radial displacement in polar coordinates across the range of shared angles between each corresponding pair of splines. Sample code is available at github.com/yinuohan/Apep.

We note that the errors associated with the extracted dust expansion speed are not statistically independent. Both the time-varying PSF and the ridge_crawler algorithm introduce correlated errors into the displacement measurements for each point along a given ridge. To estimate the error of the mean displacement of each ridge, we calculated an adjusted SEM using the ‘true’ number of independent sample points (adjusted for the autocorrelation distance of each ridge).

B3 Deriving expansion speed from ridge displacements

The 2D geometry of the dust plume arises from the projection on to the sky plane of an optically thin 3D structure. Using the geometric model for the plume, we calibrated the apparent expansion speed of individual ridges against the physical outflow speed of the system. Since the ridge_crawler algorithm introduces correlated error (newer sample points along each ridge are not independent from previous ones), multiple independent measurements were conducted. Four pairs of images of Apep’s dust plume were generated using the geometric model separated by 646 d pairwise (consistent with the time between the 2016 and 2018 VISIR epochs) and a 1 yr offset between different pairs. Each edge-like structure identified in the mid-infrared images was also identified in the model images and sampled with ridge_crawler pairwise, the displacement of which were averaged across pairs. The results are presented in Table B1.

The results suggest that expansion speed of ridges a, c, d, and f are consistent with the true expansion rate of the system, and therefore these features arise from limb-brightened edges that lie in or close to the plane of the sky. On the other hand, ridges b and e are significantly slower, implying that dust producing these likely has some projection to the line of sight. We exploited this finding to inform the ridges chosen to estimate the physical expansion speed of the dust plume, as described in Section 4.

 

APPENDIX C: GEOMETRIC MODEL

C1 Model construction

We adopt a geometric model based upon the Pinwheel mechanism, where we assume that the production of dust occurs on the surface of a conical shock – the downstream shape that results from the collision of the two winds (the curved region at the nose of the shock is too small to be seen here). This is modelled by forming circular rings at the wind–wind stagnation point, which are inflated and carried linearly away from the geometric centre of the system at the terminal wind speed (w). The geometry of the wind–wind interface carrying the dust is specified by the cone opening angle (θ_w), a physical parameter determined by the wind momentum ratio. The present position angle of the binary is defined by a time offset (T₀) from periastron. The simulated orbit begins a specified number of orbital periods (⁠|$n_{\rm circ}$|⁠) ago in which dust production occurs over a specified range of orbital phase. The binary star evolves with an orbital period (P) and eccentricity (e), and the system is further specified by additional binary-star orbital elements associated largely with viewing angle, the position angle of the ascending node (Ω), argument of periastron (ω), and the inclination of the orbital plane (i), used to project the model into a final image. In essence, the model effectively creates a conical surface wrapped into a spiral by the orbit of the binary and viewed at user-specified angles in 3D space.

C2 Numerical optimization approaches

Model optimization with a numerical algorithm such as the Markov Chain Monte Carlo (MCMC) was attempted, but we were unable to find conditions where the code would converge. The greatest hindrance to this approach is the lack of an acceptable metric that can quantify correspondence between ridges in the image and model at arbitrary clocking angles, a problem made worse by the richness of structures that such a model can generate and the high-dimensional search space. Indeed, since it is the geometry and location of the spatial features that the model is attempting to create, pixel values in isolation carry little structural information when simply applying a χ² difference metric between the data and model-generated images.

In an attempt to address this challenge, we developed an algorithm to convert both the high-pass filtered data and model-generated images into binary pixel values, effectively turning them into ‘pencil sketches’, which was designed to preserve only the main skeleton of the plume. The pencil sketches were subsequently Gaussian-blurred (to provide smooth image gradients) before computing the χ² difference between the data and model, which can be fed into a numerical optimization algorithm. However, the structural outlines generated by the model remained many and complex and, together with the high-dimensional nature of the search space, the MCMC was unable to converge. Following these attempts, we instead manually explored the parameter space informed by our knowledge of the physical parameters of the binary.

C3 Visualization

Fig. C1 provides additional images to assist the interpretation of the geometric model and its implications. A pole-on view of the dust plume, the hypothetical geometry of the dust plume without the threshold effect and the location and geometry of a more distant dust cycle predicted by the threshold effect are presented. An additional animated figure showing the formation of the dust plume over time is available as an online supplementary material.

Author notes