AU-scale radio imaging of the wind collision region in the brightest and most luminous non-thermal colliding wind binary Apep
B. Marcote, J. R. Callingham, M. De Becker, P. G. Edwards, Y. Han, R. Schulz, J. Stevens, P. G. Tuthill
MMon. Not. R. Astron. Soc. , 1–9 (2020) Printed 14 December 2020 (MN L A TEX style file v2.2)
AU-scale radio imaging of the wind collision region in the brightestand most luminous non-thermal colliding wind binary Apep
B. Marcote ,(cid:63) , J. R. Callingham , , M. De Becker , P. G. Edwards , Y. Han , R. Schulz ,J. Stevens , P. G. Tuthill Joint Institute for VLBI ERIC, Oude Hoogeveensedijk 4, 7991 PD Dwingeloo, The Netherlands Leiden Observatory, Leiden University, PO Box 9513, 2300 RA, Leiden, The Netherlands ASTRON, Netherlands Institute for Radio Astronomy, Oude Hoogeveensedijk 4, Dwingeloo, 7991 PD, The Netherlands Space sciences, Technologies and Astrophysics Research (STAR) Institute, University of Liège, Quartier Agora, 19c, Allée du 6 Août, B5c, B-4000 Sart Tilman, Belgium CSIRO Astronomy and Space Science, Australia Telescope National Facility, PO Box 76, Epping, NSW 1710, Australia Sydney Institute for Astronomy (SIfA), School of Physics, The University of Sydney, NSW 2006, Australia
Accepted 2020 December 11 / Received 2020 December 5
ABSTRACT
The recently discovered colliding-wind binary (CWB) Apep has been shown to emit lumi-nously from radio to X-rays, with the emission driven by a binary composed of two Wolf-Rayet (WR) stars of one carbon-sequence (WC8) and one nitrogen-sequence (WN4–6b).Mid-infrared imaging revealed a giant spiral dust plume that is reminiscent of a pinwheelnebula but with additional features that suggest Apep is a unique system. We have conductedobservations with the Australian Long Baseline Array to resolve Apep’s radio emission onmilliarcsecond scales, allowing us to relate the geometry of the wind-collision region to thatof the spiral plume. The observed radio emission shows a bow-shaped structure, confirmingits origin as a wind-collision region. The shape and orientation of this region is consistentwith being originated by the two stars and with being likely dominated by the stronger windof the WN4–6b star. This shape allowed us to provide a rough estimation of the opening an-gle of ∼ ◦ assuming ideal conditions. The orientation and opening angle of the emissionalso confirms it as the basis for the spiral dust plume. We also provide estimations for thetwo stars in the system to milliarcsecond precision. The observed radio emission, one orderof magnitude brighter and more luminous than any other known non-thermal radio-emittingCWB, confirms it is produced by an extremely powerful wind collision. Such a powerfulwind-collision region is consistent with Apep being a binary composed of two WR stars, sofar the first unambiguously confirmed system of its kind. Key words: radiation mechanisms: non-thermal – binaries: close – stars: individual (Apep)– radio continuum: stars – techniques: interferometric
A significant fraction of massive stars are found in binary or highermultiplicity systems (Sana et al. 2014), implying the existence of alarge reservoir of systems where the powerful winds of these mas-sive stars can collide and produce strong shocks. Such systems areclassified as colliding-wind binaries (CWBs), and the region wherethe two stellar winds collide is often known as wind-collision re-gion (WCR). WCRs are extremely e ffi cient environments to accel-erate particles up to relativistic energies (Eichler & Usov 1993),producing emission from radio to gamma-rays. WCRs are excel-lent laboratories to study particle acceleration and non-thermal pro-cesses since they operate in a unique parameter space to other astro-physical shocks. While the processes operating in WCRs are equiv- (cid:63) E-mail: [email protected] . alent to the ones in supernova remnants or interstellar bow-shocks,WCRs allow studies at higher mass, photon, and magnetic energydensities, while keeping the regions simple enough to be well de-scribed by current models and theory (e.g. Pittard & Dougherty2006; Reimer et al. 2006).The subset of CWBs known to display non-thermal emis-sion are referred to as Particle-Accelerating Colliding-Wind Bina-ries (PACWBs). The most complete census of PACWBs containsaround 40 systems (De Becker & Raucq 2013; De Becker et al.2017), with a large fraction of them involving evolved massive starscharacterized by strong, high kinetic power stellar winds. Amongthese PACWBs, many contain a Wolf-Rayet (WR) star, which arethe end point in the life of massive O-type stars, and are charac-terized by fast line-driven stellar winds ( (cid:38) km s − ) and signif-icantly higher mass-loss rates ( (cid:38) − M (cid:12) yr − ) than their progen-itor form. Therefore, the WCRs created by these types of stars are © a r X i v : . [ a s t r o - ph . S R ] D ec B. Marcote et al. expected to be more powerful than the stereotypical systems com-posed of OB-type stars.However, as lifetimes of WR stars ( ∼ yr; Meynet &Maeder 2005) are significantly shorter than their hydrogen-burningprogenitor phases, most of the known powerful PACWBs are com-posed of a single WR star and an OB-type companion. It is ex-pected only a few binary systems composed of two WR stars existin our Galaxy. Only a potential system of this kind has been re-ported so far: WR 48a (Williams et al. 2012; Zhekov et al. 2014),although such a classification is contentious since the spectra ofWR 48a does not possess emission lines indicative of two WR stars.The source 2XMM J160050.7 − ≈ . ∼
100 yr) CWB(Callingham et al. 2019). One peculiarity of Apep inconsistent withstandard CWB physics is that the measured expansion velocity ofthe spiral dust plume is at nearly an order of magnitude slower thanthat expected from the spectroscopically measured stellar winds.Callingham et al. (2019) suggested that one possible explanationfor this mismatch of measured velocities is that one of the WR starsin Apep is rapidly-rotating, producing a slow and dense equatorialwind. Such a conclusion would make Apep a potential Galacticprogenitor to long-duration gamma-ray bursts (e.g. Cantiello et al.2007).Callingham et al. (2019) predicted the existence of an unre-solved massive binary at the centre of Apep, and resolved a po-tential companion O-type supergiant ≈ . v ∞ , WC = ±
200 km s − and v ∞ , WN = ±
100 km s − , respectively (Callingham et al. 2020).Additionally, through the use of aperture masks with the NACOcamera on the VLT, Han et al. (2020) have constrained the separa-tion between the two WR stars to be consistently D = ± = ± ◦ (2016 April 28) and 278 ± ◦ (2019 March 21–24). By model fitting the dust structure, a distanceof ∼ . h m . s . s . s − ◦ D ec li n a ti on ( ) h m . s . s . s Right Ascension (2000)45 . . . . − ◦ . D ec li n a ti on ( ) − ) Figure 1.
Top panel: Mid-infrared 8 . µ m image of Apep displaying thedust pattern resembling the pinwheel nebulae but with additional intricateand overall larger structure (Callingham et al. 2019). The central white starrepresents the position of the central engine, while the o ff set green trian-gle represents the position of a northern O-type supergiant. Bottom panel:Radio image of the central engine obtained by the LBA at 13 cm on 2018July 18 showing the WCR. Contours start at three times the rms noise levelof 38 µ Jy beam − , and increase by factors of √
2. Only three contours areshown for clarity. The proper-motion corrected
Gaia
DR2 position is shownby the white star, representing the convolution of the position of the WC8and WN4–6b stars. The black crosses identify the location of the positive-flux clean components that were generated during cleaning (see Sect. 2).Their sizes are proportional to the flux density for each component. Thesynthesized beam is 11 . × . , PA = ◦ , and is shown at the bottom-left corner. For reference, 10 and 0 .
02 arcsec represent ≈
25 000 and 50 au,respectively, at the assumed source distance of ∼ . binary systems and their individual stars (see e.g. Dougherty et al.2005; Benaglia et al. 2015).To understand the nature of Apep and its radio emission, weconducted VLBI radio observations with the Southern HemisphereLong Baseline Array (LBA; Edwards & Phillips 2015) that allowedus to resolve the emission associated with the wind-collision re- © , 1–9 adio imaging the wind collision region of Apep gion. In Section 2 we describe these observations and the data re-duction. Section 3 details the obtained results and the modelling ofthe radio-emitting region. The implications of these results for un-derstanding the dynamics of Apep are discussed in Section 4. Wesummarise the conclusions of this paper in Section 5. Apep was observed with the LBA at 13 cm (2.28 GHz central fre-quency) on 2018 July 18 from 05:00 to 17:00 UTC (project codeV565). Ten stations participated in this observation: the AustralianTelescope Compact Array (ATCA), Ceduna, Hartebeesthoek, Ho-bart, Katherine, Mopra, Parkes, Tidbindilla, Warkworth, and Yarra-gadee. The data were recorded with a total bandwidth of 64 MHzand divided during correlation with the DiFX software correlator(Deller et al. 2011) into four subbands of 32 channels each. Moststations recorded full polarization, except Parkes, Tidbindilla andWarkworth, that only recorded right circular polarization.The source PKS B1622 −
297 (J1626 − −
638 was used forphasing-up ATCA. PMN J1603 − − . ◦ away from Apep) was used as phase calibrator in aphase-referencing cycle of 3.5 min on target and 1.5 min on thecalibrator. As a result, Apep was observed for a total of ≈ . (Greisen 2003) andD ifmap (Shepherd et al. 1994) following standard procedures. A-priori amplitude calibration was performed using the known gaincurves and system temperature measurements when recorded oneach station during the observation. Nominal System EquivalentFlux Density (SEFD) values were used for Ceduna, Hobart, Kather-ine, Tidbindilla, Warkwork, and Yarragadee. We manually removedbad data: missing polarizations, bad subbands, or times and fre-quencies a ff ected by radio frequency interference. We first cor-rected for the instrumental delays and bandpass calibration usingJ1626 − clean algorithm (Clark 1980). Apep is detected on milliarcsecond scales from the LBA data asa bow-shaped radio source with a peak brightness of 20 . ± .
04 mJy beam − and a total flux density of 58 . ± . ∼ ×
30 mas with a centroid, inferred from a Gaussian fit-ting to the radio image, located at the (J2000) position α WCR = h m . s ± . , δ WCR = − ◦ (cid:48) . (cid:48)(cid:48) ± . .
16 and 0 . α and δ , respectively), the uncertaintiesin the absolute International Celestial Reference Frame position ofthe phase calibrator (0 .
26 mas; Beasley et al. 2002; Gordon et al.2016), and the estimated uncertainties from the phase-referencingtechnique (0 .
06 and 0 .
19 mas; Pradel et al. 2006) added in quadra-ture. No additional sources of significant ( > σ ) radio emission are The Astronomical Image Processing System (AIPS) is a software pack-age produced and maintained by the National Radio Astronomy Observa-tory (NRAO). reported in an 6 × area centred on Apep above the rmsnoise level of 40 µ Jy.The obtained position for the radio emission lies close to theposition provided by the
Gaia data release 2 (DR2) Catalog (GaiaCollaboration et al. 2016, 2018) for the central binary (with sourceID 5981635832593607040) after accounting for proper motion (seebottom panel of Fig. 1). The separation between both positions ( ≈ < σ o ff set after accounting for uncertainties.In any case, we note that the Gaia position is composed of opticalemission arising from both stars of the central binary, which remainunresolved in their data. The radio and optical positions are thusexpected to be close but not fully agree, and it is unclear how thebinary motion of the stars impacts the precision of the
Gaia
DR2position .The reported radio flux density from the LBA data is roughlya factor of two lower ( ∼
60 against ∼
120 mJy) than the flux den-sity measured from ATCA data (Callingham et al. 2019), where thesource remains unresolved on ∼ . ∼
20% due to the intrinsic calibration procedures related toVLBI arrays, the observed di ff erence is too large to be explainedsolely by standard amplitude calibration uncertainties. The miss-ing flux is also unlikely to be related to intrinsic source variabil-ity. While the radio emission is known to be variable (Callinghamet al. 2019), such variability is only observed on very long (severalyear) timescales. However, we note that Apep is significantly re-solved by the LBA data, which is sensitive to milliarcsecond scales.This is clear from both the image plane, where the source shows asize significantly larger than the synthesized beam by roughly upto a factor of five (see Fig. 1), and from the ( u , v )-plane, where thelongest baselines of the LBA did not detect any significant emis-sion. We note that the LBA lacks short baselines (the shortest base-lines are 114 and 207 km, from ATCA–Mopra and Mopra–Parkes,respectively). Therefore, it is expected that a significant fraction ofthe radio emission is resolved out in the LBA data. The compact-ness of the emission in the ATCA data guarantees that all the radioemission reported for Apep belongs to the WCR (Callingham et al.2019).It is important to note that the missing flux density likely rep-resents only a tiny fraction of the total peak brightness in the LBAimage, given the small synthesized beam relative to the source size.We can consider that the missing ∼ × ∼ − , which represents only ∼
10% of the observedpeak brightness. In a more realistic scenario where at least part ofsuch missed emission arises from a more extended region, the as-sociated peak brightness would be even lower. We thus concludethat the observed radio region is a trustworthy representation of thepredominant location of radio emission in Apep.Additionally, the lack of compact radio emission outside theWCR (above the rms noise level of 40 µ Jy) also implies that no The parallax reported in the
Gaia
DR2 Catalog is labelled as not reli-able as it can be seen by the large values of the astrometric goodness offit and astrometric excess noise. Callingham et al. (2019, Sect. 1.2.1 in thesupplementary information) provided a possible explanation for the largeuncertainties for this system: the elongated pixels from
Gaia only resolvethe emission from the central binary and the northern O-type star for someorientations. This can thus bias the parallax measurements. ©000
Gaia only resolvethe emission from the central binary and the northern O-type star for someorientations. This can thus bias the parallax measurements. ©000 , 1–9 B. Marcote et al. significant interaction occurs between the central binary and thenorthern O-type supergiant. We also clarify that the individual stel-lar winds of each of the three stars are not detected. Their individ-ual thermal emission is estimated to lie below the 0.1-mJy level at13 cm, using the formalism proposed by Wright & Barlow (1975).At a distance of about 2.5 kpc, such winds would be more likelydetected at shorter wavelengths, considering the thermal nature ofthe emission.
The observed bow-shaped radio emission is representative of aWCR. However, the reported radio flux densities make Apepbrightest and most luminous radio-bright CWB discovered to date.The presence of two WR stars, both exhibiting powerful winds,may introduce a significant di ff erence in the properties involved atthe WCR with respect to the more common CWBs comprising oneWR star and an OB star. For example, the emission could be bol-stered by a significantly higher kinetic power budget for the winds,and / or stronger magnetic fields in the shock than typical. A largerenergy budget and magnetic field strength could go some way toexplaining why Apep displays synchrotron radio emission one or-der of magnitude brighter and more luminous than other CWBs.We note that the observed synchrotron emission is expected to beproportional to the magnetic field as B / and to the electron den-sity (Longair 2011), where the latter is related to ˙ M v − (Cantó et al.1996). The significantly slower winds and higher mass-loss rates ofthe WC star could then bolster the radio emission.Furthermore, the separation between the two stars may be op-timal for luminous non-thermal radio emission: being close enoughto produce a powerful shock, but far enough to avoid free-free ab-sorption at gigahertz-frequencies. The size of the radio photosphere(radius of a sphere at which the optical depth is equal to one) ofeach individual wind can be estimated according to the approachadopted for instance by De Becker et al. (2019). This quantity al-lows for a rough estimate of the capability of a given stellar wind toattenuate radio emission due to free-free absorption. At 13 cm, oneobtains values shorter than 10 au. Assuming the WCR is locatedroughly midway between the two stars, the projected separationbetween each star and the stagnation point at ∼ . ∼
59 au. The synchrotron emitting region is thus clearly away fromthe denser parts of the winds and is thus not a ff ected by free-freeabsorption. This is consistent with the absence of signatures of ab-sorption have been reported to date at gigahertz frequencies (Call-ingham et al. 2019), and the observed radio morphology, that wouldbe expected to be centred on the stagnation point of the shock, tooccur at some intermediate position between the stars.For comparison, we note that η Carinae, the most powerfulCWB, only shows thermal radio emission due to its closer binarysemi-major axis of ≈ . η Carinae’s WCR is likely to be com-pletely free-free absorbed due to the aforementioned conditions(De Becker & Raucq 2013).
The central engine in Apep is known to host a WC8 and WN4–6b star (Callingham et al. 2020). The separation between the twostars has been measured by near-infrared NACO data, taken in 2016April 28 and 2019 March 21–24, to be consistently D = ± = ± ± ◦ , respectively (Han et al. 2020). However, the absolute positions of the two starsremained unknown from these data. We note that the orbital pe-riod of the central binary is estimated to be ∼
100 yr, and thus wedo not expect significant changes ( ∼ ± ◦ . This value for the PA is consistentwith the orientation observed in the curved radio emission of Apep(see Fig. 1). We can thus confirm that the observed radio emissionarises from the interaction of only these two stars.A Gaussian fitting to the region revealed that while the radioemission is clearly elongated in the North-South direction (with asemi major axis — half width at half maximum — of ∼
18 mas),its extension in the direction between the two stars (roughly East-West) is comparable to the size of the synthesized beam ( ∼ (cid:46)
12 mas (or (cid:46)
30 au at the assumed distance of ∼ . A full understanding of the observed WCR is only possible by con-ducting a full magnetohydrodynamic and radiative-transfer simula-tion of the collision between the two stellar winds, that go beyondthe purpose of this manuscript. However, we note that even a sim-plistic scenario assuming two ideal (spherically symmetric) stellarwinds can be typically taken as a first (rough) approach to charac-terize the observed broad picture of the radio emission, as success-fully demonstrated in other CWBs (see e.g. Blomme et al. 2010;Benaglia et al. 2015), and test its consistency with respect to theproperties derived from the IR data (Callingham et al. 2019, 2020;Han et al. 2020).Following this path, and following the details provided in Ap-pendix A, one could compare the observed bow-shaped radio emis-sion to the ideal contact discontinuity (CD) shape expected undersuch scenario. Figure 2 shows this ideal CD with respect to the ob-served emission, suggesting that the observed curvature can still beunderstood under ideal conditions – and for the given resolution.The ideal CD shape provides both a rough estimation of thewind momentum rate ratio – under the assumption of sphericallysymmetric and homogeneous winds – of η = . ± .
08 (see equa-tion A3) and of the putative positions of the two stars (see Ap-pendix A for details) at the epoch of the LBA observations (2018July 18): α WC = h m . s ± ,δ WC = − ◦ (cid:48) . (cid:48)(cid:48) ± . , (1)for the WC8 star, and α WN = h m . s ± ,δ WN = − ◦ (cid:48) . (cid:48)(cid:48) ± . , (2)for the WN4–6b star. We note that the quoted uncertainties take intoaccount both the statistical uncertainties of the radio data and the © , 1–9 adio imaging the wind collision region of Apep . s . s . s . s h m . s Right Ascension (2000)45 . . . . − ◦ . D ec li n a ti on ( ) WC WN
Figure 2.
Ideal contact discontinuity (CD) compared to the shape ofApep’s radio emission. Contours start at three times the rms noise levelof 38 µ Jy beam − , and increase by factors of √
2, and represent the sameimage presented in the bottom panel of Fig. 1. The solid blue line and theshadowed light blue region represent the ideal CD shape for the mean valueand 1- σ confidence interval of η = . ± .
08, respectively. The predictedpositions for the two stars (WC and WN) from the given CD are representedby the blue crosses, where uncertainties on position (relative to the positionof the stagnation point) represent the 1- σ confidence interval. The synthe-sized beam is 11 . × . , PA = ◦ , and is shown at the bottom-leftcorner. For reference, 0 .
02 arcsec represent 50 au at the assumed sourcedistance of ∼ . curvature on the ideal CD shape, and the systematic uncertaintiesfor the absolute position (not shown in Fig. 2).The ideal CD shape shown in Fig. 2 would also imply an open-ing angle of 2 θ w = ± ◦ , which is in fact consistent with thevalue of ∼ ◦ estimated from optical / near-infrared spectroscopicobservations (Callingham et al. 2020), and slightly higher than thevalue of ∼ ◦ determined from the model of the dust expansion(Callingham et al. 2019; Han et al. 2020). We note that, in any case,we would expect the opening angle from the WCR to deviate fromthe later one for several reasons, not just because of the ideal con-ditions assumed in the WCR. First, it could be expected that theopening angles of the shocks have not reached yet their asymp-totic values at the region where the radio emission is produced (seee.g. Pittard & Dougherty 2006), but such value is reached at thescales where the dust structure is observed. Secondly, while theradio WCR is created and dominated by the direct collision of thetwo stellar winds, the dust plume is more sensitive to the significantclumpiness and turbulence of the winds of the WR stars (see e.g.Crowther 2007) and of the region with the mixed, shocked, winds.Furthermore, a larger opening angle for the WCR than that derivedfrom the dust plume could actually be expected due to the condi-tions necessary for dust formation to occur (Tuthill et al. 2008). Insummary, only a small portion of the winds are expected to collideat or near the stagnation point (the radio-emitting WCR). Most ofthe wind interactions take place at larger distances along the spi-ral plume (Tuthill et al. 2008), which already deviate significantlyfrom any ideal wind-interaction scenario. This would explain whythe typically observed WCR, including this one in Apep, can beconsistently explained by such ideal scenarios with no significantdeviations, while the full modelling of the system (including thedust plume) would require a more elaborated scenario. Finally, we note that the assumed ideal CD – assuming idealand spherical winds – would imply a mass-loss rate ratio of˙ M WC / ˙ M WN = . ± .
15 (considering the measured terminal line-of-sight wind velocities of the two WR stars from spectroscopic ob-servations, Callingham et al. 2020, of v ∞ , WC = ±
200 km s − and v ∞ , WN = ±
100 km s − ). This value is actually con-sistent with the typical mass loss rates expected for these kindsof stars: Galactic WC8 stars show an average mass loss rate of˙ M WC8 ∼ − . M (cid:12) yr − (Sander et al. 2019), while WN4–6 starsexhibit mass loss rates ranging 10 − . to 10 − . M (cid:12) yr − (Hamannet al. 2019). We can then adopt an average value of ˙ M WN4–6b ∼ − . M (cid:12) yr − , which would imply a ratio of ∼ .
63, consistentwith the aforementioned one.However, the expansion of the dust plume has been measuredto be significantly slower ( ∼
600 km s − ) than the ∼ − the dust is expected to inherit from the collision of the winds withthe aforementioned mass loss rates and terminal wind speeds (Call-ingham et al. 2020). One possible solution to this discrepancy isthat the wind of at least one of the stars, likely the WC8, is highlyasymmetric, with a putative dense, slow equatorial wind potentiallyproduced by a rapid rotation of the star (Callingham et al. 2019,2020; Han et al. 2020). In this case, one could expect its stellarwind to be as slow as (cid:46) − to explain the dust expan-sion. In such a case, that would imply a mass-loss rate ratio of˙ M WC8 / ˙ M WN4–6b ∼ .
5. Assuming that the WC8 star still exhibitsa typical mass-loss rate, we estimate ˙ M WN4–6b ∼ − . M (cid:12) yr − ,which still lies within the observed range of mass-loss rates forthese types of stars (Hamann et al. 2019). Assuming, on the otherhand, an average mass-loss rate for the WN4–6b star, we estimate˙ M WC8 ∼ − . M (cid:12) yr − , which would imply a value lower than theones observed for WC8 stars (Sander et al. 2019). However, we re-mark that this value would be an estimation for the mass-loss rateof the WC8 star only valid at equatorial latitudes, as we are underthe assumption of non-spherically-symmetric winds.The two extreme cases (symmetric or highly-asymmetricwinds) provide then consistent results with the observed curvatureof the WCR. The existence of a more likely scenario where atleast one of the winds may exhibit a inhomogeneous profile wouldimply that the wind momentum rate ratio estimation must beredefined by taking into account such profile (not fully understoodto date) and the inclinations of the two star spins with respect tothe orbital plane.To summarize, we have imaged the radio emission of Apep,whose structure is consistent with a WCR produced by the WC8and WN4–6b stars, and thus directly resolves the nature of Apepas a CWB. The orientation, position, and opening angle of the ob-served bow-shaped structure are consistent with the positions de-rived from the NACO (Han et al. 2020) and X-SHOOTER data(Callingham et al. 2020). Additionally, even a scenario assumingtwo spherical winds with ideal conditions provide a consistent de-scription of the system.Following VLBI observations of the source with a cadenceof years would allow us to reconstruct the orbital motion of thetwo stars and the evolution of the WCR. These data would unam-biguously and accurately constrain the orbit of Apep. These data,together with the spectral information from hundreds of megahertzto tens of gigahertz, would allow us to fully characterize the radio-emitting region; constraining not only the dynamic properties of thewinds but also the particle density and magnetic field at the WCR.In particular, it would be interesting whether, or at what phase ofthe orbit, the WCR becomes fully free-free absorbed. Such a mea- ©000
5. Assuming that the WC8 star still exhibitsa typical mass-loss rate, we estimate ˙ M WN4–6b ∼ − . M (cid:12) yr − ,which still lies within the observed range of mass-loss rates forthese types of stars (Hamann et al. 2019). Assuming, on the otherhand, an average mass-loss rate for the WN4–6b star, we estimate˙ M WC8 ∼ − . M (cid:12) yr − , which would imply a value lower than theones observed for WC8 stars (Sander et al. 2019). However, we re-mark that this value would be an estimation for the mass-loss rateof the WC8 star only valid at equatorial latitudes, as we are underthe assumption of non-spherically-symmetric winds.The two extreme cases (symmetric or highly-asymmetricwinds) provide then consistent results with the observed curvatureof the WCR. The existence of a more likely scenario where atleast one of the winds may exhibit a inhomogeneous profile wouldimply that the wind momentum rate ratio estimation must beredefined by taking into account such profile (not fully understoodto date) and the inclinations of the two star spins with respect tothe orbital plane.To summarize, we have imaged the radio emission of Apep,whose structure is consistent with a WCR produced by the WC8and WN4–6b stars, and thus directly resolves the nature of Apepas a CWB. The orientation, position, and opening angle of the ob-served bow-shaped structure are consistent with the positions de-rived from the NACO (Han et al. 2020) and X-SHOOTER data(Callingham et al. 2020). Additionally, even a scenario assumingtwo spherical winds with ideal conditions provide a consistent de-scription of the system.Following VLBI observations of the source with a cadenceof years would allow us to reconstruct the orbital motion of thetwo stars and the evolution of the WCR. These data would unam-biguously and accurately constrain the orbit of Apep. These data,together with the spectral information from hundreds of megahertzto tens of gigahertz, would allow us to fully characterize the radio-emitting region; constraining not only the dynamic properties of thewinds but also the particle density and magnetic field at the WCR.In particular, it would be interesting whether, or at what phase ofthe orbit, the WCR becomes fully free-free absorbed. Such a mea- ©000 , 1–9 B. Marcote et al. surement would provide another independent estimate of the mass-loss rate present in the system (Dougherty et al. 2003). In general,we would expect these additional parameters would help reveal thefundamental di ff erences in Apep with respect to other CWBs toexplain why radio emission is so bright. The LBA data presented here allowed us to resolve the radio emis-sion of Apep, revealing the presence of a strong wind collision re-gion produced by the WC8 and WN4–6b stars. The observed ori-entation and opening angle of such bow-shock region are consis-tent with the dust spiral structure observed by Callingham et al.(2019) and modelled by Han et al. (2020). We confirm that the fullstructure (radio-emitting WCR and dust pinwheel nebula) are pro-duced by these two WR stars in the same shock and no additionalinteractions are required from third components, like the nearbyO-type supergiant of the system. Apep is the first known particle-accelerating colliding-wind binary unambiguously composed oftwo WR stars, and is the brightest and most luminous PACWB byan order of magnitude (De Becker & Raucq 2013).Apep thus belongs to a select group of binaries composed oftwo WR stars. As the WR is a short-lived stage in the life of stars,it is expected that this kind of binaries are rare. Apep is thus a valu-able object that allows us to study in detail these extreme wind in-teractions and the evolution of these types of stars. Given that thesesystems are potential progenitors of long gamma-ray bursts (e.g.Cantiello et al. 2007), the resulting studies would provide valuabledata to constrain the evolution and dynamics of these systems priorto cataclysm.To finalize, we have shown how VLBI observations have beenonce again a successful approach to unveil the nature of CWBs,and the only approach able to resolve the radio emission arisingfrom the WCR. By studying the morphology of this region, we havebeen able to constrain the properties of the system like the wind-momentum rate ratio and the opening angle of the shock and, forthe first time, indirectly estimate the positions of the two stars in thesky. Multi-epoch observations, combining very-high-resolution andmulti-frequency data, spread along a decade ( ∼
10% of the orbit)should allow us to trace a significant variation of these propertiesto better characterize Apep and the evolution of the two stars.
ACKNOWLEDGMENTS
The authors thank the anonymous referees for the insightful com-ments that have significantly improved the manuscript. We alsothank P. A. Crowther, I. Negueruela, and P. Williams for use-ful discussions, and C. Day for helping during the observations.The Long Baseline Array is part of the Australia Telescope Na-tional Facility which is funded by the Australian Government foroperation as a National Facility managed by CSIRO. This paperincludes archived data obtained through the Australia TelescopeOnline Archive ( http://atoa.atnf.csiro.au ). This work wassupported by resources provided by the Pawsey SupercomputingCentre with funding from the Australian Government and the Gov-ernment of Western Australia. This work has made use of datafrom the European Space Agency (ESA) mission
Gaia ( ), processed by the Gaia
Data Pro-cessing and Analysis Consortium (DPAC, ). Funding for the DPAC has been provided by national institutions, in particular the insti-tutions participating in the
Gaia
Multilateral Agreement. BM ac-knowledges support from the Spanish Ministerio de Economía yCompetitividad (MINECO) under grant AYA2016-76012-C3-1-Pand from the Spanish Ministerio de Ciencia e Innovación undergrants PID2019-105510GB-C31 and CEX2019-000918-M of IC-CUB (Unidad de Excelencia “María de Maeztu” 2020-2023). JRCthanks the Nederlandse Organisatie voor Wetenschappelijk Onder-zoek (NWO) for support via the Talent Programme Veni grant. YHacknowledges the traditional owners of the land, the Gadigal peo-ple of the Eora Nation, on which the University of Sydney is builtand some of this work was carried out. This research made useof APL py , an open-source plotting package for P ython (Robitaille& Bressert 2012); A stropy , a community-developed core P ython package for Astronomy (Astropy Collaboration et al. 2013); andM atplotlib (Hunter 2007). This research made use of NASA’s As-trophysics Data System. REFERENCES
Astropy Collaboration et al., 2013, A&A, 558, A33Beasley A. J., Gordon D., Peck A. B., Petrov L., MacMillan D. S.,Fomalont E. B., Ma C., 2002, ApJS, 141, 13Benaglia P., Marcote B., Moldon J., Nelan E., De Becker M.,Dougherty S. M., Koribalski B. S., 2015, A&A, 579, A99Blomme R., De Becker M., Volpi D., Rauw G., 2010, A&A, 519,A111Callingham J. R., Tuthill P. G., Pope B. J. S., Williams P. M.,Crowther P. A., Edwards M., Norris B., Kedziora-Chudczer L.,2019, Nature Astronomy, 3, 82Callingham J. R., Crowther P. A., Williams P. M., Tuthill P. G.,Han Y., Pope B. J. S., Marcote B., 2020, MNRAS, 495, 3323Cantiello M., Yoon S. C., Langer N., Livio M., 2007, A&A, 465,L29Cantó J., Raga A. C., Wilkin F. P., 1996, ApJ, 469, 729Clark B. G., 1980, A&A, 89, 377Crowther P. A., 2007, ARAA, 45, 177D’Agostini G., 2003, Reports on Progress in Physics, 66, 1383De Becker M., Raucq F., 2013, A&A, 558, A28De Becker M., Benaglia P., Romero G. E., Peri C. S., 2017, A&A,600, A47De Becker M., Isequilla N. L., Benaglia P., 2019, A&A, 623,A163Deller A. T., et al., 2011, PASP, 123, 275Dougherty S. M., Pittard J. M., Kasian L., Coker R. F., WilliamsP. M., Lloyd H. M., 2003, A&A, 409, 217Dougherty S. M., Beasley A. J., Claussen M. J., Zauderer B. A.,Bolingbroke N. J., 2005, ApJ, 623, 447Edwards P. G., Phillips C., 2015, Publication of Korean Astro-nomical Society, 30, 659Eichler D., Usov V., 1993, ApJ, 402, 271Gaia Collaboration et al., 2016, A&A, 595, A2Gaia Collaboration et al., 2018, A&A, 616, A1Gordon D., et al., 2016, ApJ, 151, 154Greisen E. W., 2003, in Heck A., ed., Astrophysics and SpaceScience Library Vol. 285, Information Handling in Astronomy -Historical Vistas. p. 109, doi:10.1007 / © , 1–9 adio imaging the wind collision region of Apep Longair M. S., 2011, High Energy Astrophysics, third edn. Cam-bridge University Press, CambridgeMadura T. I., Gull T. R., Owocki S. P., Groh J. H., Okazaki A. T.,Russell C. M. P., 2012, MNRAS, 420, 2064Meynet G., Maeder A., 2005, A&A, 429, 581Pearson K. F. R. S., 1900, The London, Edinburgh, and DublinPhilosophical Magazine and Journal of Science, 50, 157Pittard J. M., Dougherty S. M., 2006, MNRAS, 372, 801Pradel N., Charlot P., Lestrade J.-F., 2006, A&A, 452, 1099Reimer A., Pohl M., Reimer O., 2006, ApJ, 644, 1118Robitaille T., Bressert E., 2012, APLpy: Astronomical Plot-ting Library in Python, Astrophysics Source Code Library(ascl:1208.017)Sana H., et al., 2014, ApJS, 215, 15Sander A. A. C., Hamann W. R., Todt H., Hainich R., Shenar T.,Ramachandran V., Oskinova L. M., 2019, A&A, 621, A92Shepherd M. C., Pearson T. J., Taylor G. B., 1994, in Bulletin ofthe American Astronomical Society. p. 987Stevens I. R., Blondin J. M., Pollock A. M. T., 1992, ApJ, 386,265Tuthill P. G., Monnier J. D., Danchi W. C., 1999, Nature, 398, 487Tuthill P. G., Monnier J. D., Lawrance N., Danchi W. C., OwockiS. P., Gayley K. G., 2008, ApJ, 675, 698Williams P. M., et al., 2009, MNRAS, 395, 1749Williams P. M., van der Hucht K. A., van Wyk F., Marang F.,Whitelock P. A., Bouchet P., Setia Gunawan D. Y. A., 2012, MN-RAS, 420, 2526Wright A. E., Barlow M. J., 1975, MNRAS, 170, 41Zhekov S. A., Tomov T., Gawronski M. P., Georgiev L. N.,Borissova J., Kurtev R., Gagné M., Hajduk M., 2014, MNRAS,445, 1663
APPENDIX A: CONTACT DISCONTINUITY FROM TWOSPHERICALLY-SYMMETRIC WINDS
As presented in the main paper, the bow-shaped radio emissionobserved in Apep is fully consistent with emission arising from aWCR. Given that the shock is expected to be adiabatic due to thescales of the system (Stevens et al. 1992), and that the two shockfronts remain unresolved in the LBA data (see Sect. 4.1), we wouldexpect the contact discontinuity (CD) to be a plausible proxy todescribe the morphology of the radio emission.In this appendix we provide the detailed description of a CDin the ideal case of a binary system where the two stellar winds arespherically-symmetric and homogeneous, following the descrip-tion done by Cantó et al. (1996). We note that this scenario does notproperly characterise the environment of Apep, where an asymmet-ric wind may be expected, but its simplicity allowed us to providesome meaningful estimations for the system with a minimal num-ber of assumptions.
A1 Derivation of the CD
Under the assumption of spherically-symmetric stellar winds andadiabatic shocks, if these two shocks remain close enough so theyare not resolved, then they would be expected to encompass the CD.Following Cantó et al. (1996), the CD can be then characterizedanalytically as the interaction front where the ram pressure of theideal stellar winds of the two stars (WC and WN in the case ofApep) are balanced (see Fig. A1): ρ WC v , ⊥ = ρ WN v , ⊥ (where ρ WC WN rR WC R WN θ WC θ WN D xyα δ PA ξ Figure A1.
Schematic illustration of the contact discontinuity (CD) and thenotation used in this study. The two stars are represented by the blue circlesand the thick blue line represents the WCR, where the ram pressure of thetwo stellar winds (thin blue arrows) is balanced. and v ⊥ are the density and the wind velocity normal to the front,respectively). It has been shown (following Cantó et al. 1996) thatthis front can be parametrized as: r ( θ WC , θ WN ) = D sin θ WN sin − ( θ WC + θ WN ) , (A1)where r is the separation of the front to the WC star for given θ WC and θ WN angles, and D is the separation between the two stars. SeeFig. A1 for a representation of this front. The two angles are relatedby the expression θ WN tan − θ WN = + η (cid:16) θ WC tan − θ WC − (cid:17) , (A2)where η is the wind momentum rate ratio defined as: η ≡ ˙ M WC v WC ˙ M WN v WN = (cid:32) R WC R WN (cid:33) , (A3)where ˙ M , v , and R are the mass-loss rate, stellar wind velocity, anddistance to the stagnation point relative to the WC and WN stars,respectively. We note that η is independent of the orbital inclina-tion, and thus can be directly measured for any system once thepositions of the two stars and the WCR are measured. We remarkthat under the described scenario, the two stellar winds are assumedto be spherically-symmetric, and thus this value of η is assumed toremain constant.The asymptotic angle θ WC , ∞ of the bow shock (correspondingto r → ∞ ) can be found from equation (A2) given that both anglesmust verify θ WC , ∞ + θ WN , ∞ = π : θ WC , ∞ − tan θ WC , ∞ = π (1 − η ) − , (A4)which is related to the shock full opening angle 2 θ w = π − θ WC , ∞ ).This angle can be estimated from the morphology typically ob-served in resolved radio-emitting WCRs, and can then be used toestimate the positions of the two stars. Both the separation betweenthe stars and the stagnation point can be directly obtained by fol- ©000
Schematic illustration of the contact discontinuity (CD) and thenotation used in this study. The two stars are represented by the blue circlesand the thick blue line represents the WCR, where the ram pressure of thetwo stellar winds (thin blue arrows) is balanced. and v ⊥ are the density and the wind velocity normal to the front,respectively). It has been shown (following Cantó et al. 1996) thatthis front can be parametrized as: r ( θ WC , θ WN ) = D sin θ WN sin − ( θ WC + θ WN ) , (A1)where r is the separation of the front to the WC star for given θ WC and θ WN angles, and D is the separation between the two stars. SeeFig. A1 for a representation of this front. The two angles are relatedby the expression θ WN tan − θ WN = + η (cid:16) θ WC tan − θ WC − (cid:17) , (A2)where η is the wind momentum rate ratio defined as: η ≡ ˙ M WC v WC ˙ M WN v WN = (cid:32) R WC R WN (cid:33) , (A3)where ˙ M , v , and R are the mass-loss rate, stellar wind velocity, anddistance to the stagnation point relative to the WC and WN stars,respectively. We note that η is independent of the orbital inclina-tion, and thus can be directly measured for any system once thepositions of the two stars and the WCR are measured. We remarkthat under the described scenario, the two stellar winds are assumedto be spherically-symmetric, and thus this value of η is assumed toremain constant.The asymptotic angle θ WC , ∞ of the bow shock (correspondingto r → ∞ ) can be found from equation (A2) given that both anglesmust verify θ WC , ∞ + θ WN , ∞ = π : θ WC , ∞ − tan θ WC , ∞ = π (1 − η ) − , (A4)which is related to the shock full opening angle 2 θ w = π − θ WC , ∞ ).This angle can be estimated from the morphology typically ob-served in resolved radio-emitting WCRs, and can then be used toestimate the positions of the two stars. Both the separation betweenthe stars and the stagnation point can be directly obtained by fol- ©000 , 1–9 B. Marcote et al. lowing equations (A1) and (A2) to be: R WC = η / D + η / , R WN = D + η / , (A5)where R WC and R WN are defined as in Fig. A1. A2 Placing the CD in the plane of the sky
Equations (A1) and (A2) define the CD region for the referencesystem shown in Fig. A1, where the origin of coordinates is placedat the center of the WC star and the x axis in the direction of the WNstar. For systems where the orbital plane lies almost in the plane ofthe sky (as in the case of Apep, with an inclination of ∼ ± ◦ ; Hanet al. 2020), the radio emission would be seen as a conical structurewhere the front envelope can be parametrized by the following CDcurve (in Cartesian coordinates): x CD = r cos θ WC ,y CD = r sin θ WC ; (A6)and, by using the relations in equation (A1): x CD = D (cid:104) + sin( θ WN − θ WC ) sin − ( θ WC + θ WN ) (cid:105) ,y CD = D (cid:104) cos( θ WN − θ WC ) sin − ( θ WC + θ WN ) − tan − ( θ WC + θ WN ) (cid:105) . (A7)In absence of significant absorption, as in the case of Apep, theradio emission is expected to be maximal at the stagnation point,and quickly decay for larger values of θ . Under this scenario, andfor a system like Apep (with an orbit on the plane of the sky), it canthen be expected that the given envelope would gather the bulk ofthe radio emission following projection reasoning (e.g. such enve-lope would collect the highest values of the column density).A generalized form for equation (A7) is however needed giventhat the system can have a random orientation and is placed at somesky coordinates ( α, δ ). Furthermore, the positions of the stars were a-priori unknown. Only the position of the radio-emitting WCRcould be directly determined from the LBA data. We thereforechose a better reference system where the origin of coordinates isplaced at the stagnation point of the CD (which is expected to co-incide with the peak of the reported radio emission, α WCR , δ
WCR ),and rotated by an angle ξ . We can then recover the CD curve in skycoordinates: (cid:32) α CD cos δ WCR δ CD (cid:33) = (cid:32) cos ξ − sin ξ sin ξ cos ξ (cid:33) · (cid:32) x CD − R WC y CD (cid:33) + (cid:32) α WCR cos δ WCR δ WCR (cid:33) , (A8)where α CD , δ CD are functions of θ WC and θ WN , but it can reducedto only θ WC by numerically solving equation (A2). We note thatthis transformations assume that the declination can be consideredconstant for the full CD region (i.e. the system subtends a smallangle in the sky). As it can be seen, the final curve only dependson the following parameters: the position of the stagnation point inthe sky ( α WCR , δ
WCR ), the position angle of the system ( ξ , whichis related to the PA of the system mentioned in the main text by ξ = π − PA), the separation between the two stars ( D ), and thewind-momentum rate ratio ( η ) . The code used to obtain the CD curve is part of the
Binaries pack-age that can be publicly found at https://github.com/bmarcote/binaries under GPLv3 license.
A3 Comparison with WCR VLBI images
As mentioned in Sect. 3, typical radio images are described as acollection of point-like components (so-called clean components;see Fig. 1). Each component contains information of its position( α i , δ i ) and flux ( S i ). One can then use the positions of all cleancomponents with positive flux to fit the expected CD curve.To compare the expected ideal CD with respect to the ob-served morphology of the radio-emitting WCR in Apep we firstfixed the parameters D and PA (as both are provided by Han et al.2020), and the position of the stagnation point to match the cen-troid of the radio image (see Sect. 3). Only the eta parameter wasthen free. To compare which geometry best explained the observedcurvature, we computed di ff erent curves for 500 di ff erent values of η ranging 0.2–0.8, each of them covering angles of | θ WC | (cid:54) ◦ .These limit values were chosen by a manual inspection of Fig. 2.From preliminary trials we guaranteed that the given range for η contained any reasonable value that could reproduce the observedcurvature, and the aforementioned values of θ WC were the ones cov-ering the region with significant radio emission (i.e. where all cleancomponents were located; see Fig. 2).The determination of the most plausible values of η was per-formed by computing the χ values from the cumulative separa-tions of each clean component position ( α i , δ i ) to the CD curve( α CD , δ CD ). Given that the CD curve (for a given η ) is defined asa discrete set of positions ( α i CD , δ i CD ), the χ values were computedas χ ( η ) = (cid:88) j min i (cid:16) ( α j − α i CD ) cos δ WCR (cid:17) | α i CD − α WCR | cos δ WCR + (cid:16) δ j − δ i CD (cid:17) | δ i CD − δ WCR | . (A9)That is, we estimated the angular separation between each cleancomponent and the closest point to the CD curve. The resulting χ is then determined following the Pearson generalized form (Pear-son 1900), where the separations between the data and the expectedpositions are weighted by the expected positions, which are relativeto the stagnation point.We then proceeded to compute the χ probability density func-tion (p.d.f.) to obtain the most probable value of η , which can beapproximated under Gaussian assumptions to be f ∝ exp( − χ / η = (cid:90) η f ( η ) d η,σ η = (cid:90) ( η − ¯ η ) f ( η ) d η. (A10)As mentioned in the text, we obtained a value of η = . ± . ff erences in the final results withrespect to the aforementioned approach and we recovered the samevalues of η .Finally, even when the obtained results explain accurately theobserved emission and produce consistent physical parameters forApep, it is worth to remark the few caveats that underlie our model,specially when applying it to the case of Apep. The main ones areobviously related to the assumed ideal conditions and ideal CD © , 1–9 adio imaging the wind collision region of Apep shape. We note that the described model assumes the radio emis-sion to follow the ideal CD shape. The radio emission is expected tobe confined between two shock fronts, and we assume both of themto be close enough to be well described by the CD (as mentionedin Sect. 3). This approximation has been widely used in previousstudies of CWBs, showing consistent results (see e.g. Blomme et al.2010; Benaglia et al. 2015). Furthermore, Fig. 2 shows that the ra-dio emission fits nicely the predicted CD curve, and thus we canbe confident that such approximation is also accurate in the case ofApep for the given resolution that we achieved with the LBA data.We note that the strong radio emission from Apep (roughly one or-der of magnitude brighter than any other known non-thermal radioemitting CWB), makes the system an exceptional case as it allowsus to obtain reliable VLBI images with a large signal-to-noise, min-imizing the typical concerns on these kinds of analyses (Pittard &Dougherty 2006).The derived opening angle, on the other hand, may be lessreliable. The radio emission is likely to be produced at the regionwhere the opening angles of the shocks have not reached yet theirasymptotic values (see e.g. Pittard & Dougherty 2006). And finally,the compared mass-loss rates are actually sensitive to the stellarwind velocities and the η estimation. In addition, the possibilitythat a system like Apep shows anisotropic winds for at least one ofthe stars would produce – in this model – an average value of η thatsmears such di ff erences.In any case, we note that in our model we did not considerneither the inclination of the system nor the brightness profile of theWCR. Whereas the former one is not expected to have a significante ff ect on the η parameter (Apep is an almost face-on system, Hanet al. 2020, and the bow-shape curvature does not vary significantlyfor low inclination systems), the brightness profile could allow us toestimate the physical conditions at the WCR. However, this wouldnot have any e ff ect on the estimated wind momentum rate ratio andit would imply a significant number of additional assumptions andideal conditions. We thus consider that such analysis is beyond thescope of this paper and a more realistic, (magneto-)hydrodynamicalmodels, would be always preferable. ©000