<figure><img src="image_1.png" /></figure>

# Revisiting the Impact of Dust Production from Carbon-rich Wolf–Rayet Binaries

Ryan M. Lau¹, J. J. Eldridge², Matthew J. Hankins³, Astrid Lamberts⁴, Itsuki Sakon⁵, and Peredur M. Williams⁶

¹Institute of Space & Astronautical Science, Japan Aerospace Exploration Agency, 3-1-1 Yoshinodai, Chuo-ku, Sagamihara, Kanagawa 252-5210, Japan
ryanlau@ir.isas.jaxa.jp

²Department of Physics, University of Auckland, Private Bag 92019, Auckland 1010, New Zealand

³Division of Physics, Mathematics, and Astronomy, California Institute of Technology, Pasadena, CA 91125, USA

⁴Université Côte d'Azur, Observatoire de la Côte d'Azur, CNRS, Laboratoire Lagrange, Laboratoire ARTEMIS, France

⁵Department of Astronomy, School of Science, University of Tokyo, 7-3-1 Hongo, Bunkyo-ku, Tokyo 113-0033, Japan

⁶Institute for Astronomy, University of Edinburgh, Royal Observatory, Edinburgh EH9 3HJ, UK

Received 2019 November 25; revised 2020 May 18; accepted 2020 June 12; published 2020 July 24

## Abstract

We present a dust spectral energy distribution (SED) and binary stellar population analysis revisiting the dust production rates (DPRs) in the winds of carbon-rich Wolf–Rayet (WC) binaries and their impact on galactic dust budgets. DustEM SED models of 19 Galactic WC “dustars” reveal DPRs of $\dot{M}_d \sim 10^{-10}-10^{-6} M_\odot$ yr$^{-1}$ and carbon dust condensation fractions, $\chi_C$, between 0.002% and 40%. A large (0.1–1.0 $\mu$m) dust grain size composition is favored for efficient dustars where $\chi_C \gtrsim 1$. Results for dustars with known orbital periods verify a power-law relation between $\chi_C$, orbital period, WC mass-loss rate, and wind velocity consistent with predictions from theoretical models of dust formation in colliding-wind binaries. We incorporated dust production into Binary Population and Spectral Synthesis (BPASS) models to analyze dust production rates from WC dustars, asymptotic giant branch stars (AGBs), red supergiants (RSGs), and core-collapse supernovae (SNe). BPASS models assuming constant star formation (SF) and a coeval $10^6 M_\odot$ stellar population were performed at low, Large Magellanic Cloud (LMC)-like, and solar metallicities ($Z = 0.001, 0.008, \text{and } 0.020$). Both constant SF and coeval models indicate that SNe are net dust destroyers at all metallicities. Constant SF models at LMC-like metallicities show that AGB stars slightly outproduce WC binaries and RSGs by factors of 2–3, whereas at solar metallicities WC binaries are the dominant source of dust for ~60 Myr until the onset of AGBs, which match the dust input of WC binaries. Coeval population models show that, for “bursty” SF, AGB stars dominate dust production at late times ($t \gtrsim 70$ Myr).

Unified Astronomy Thesaurus concepts: Circumstellar dust (236); Wolf-Rayet stars (1806); Infrared astronomy (786); WC stars (1793); Interstellar dust (836)

## 1. Introduction

Dust formation can, surprisingly, occur in the hostile circumstellar environment of Wolf–Rayet (WR) stars, which are descendants of massive O-stars characterized by fast winds ($\gtrsim 1000$ km s$^{-1}$), hot photospheres ($T_* \gtrsim 40,000$ K), and high luminosities ($L_* \sim 10^5 L_\odot$; Gehrz & Hackwell 1974; Williams et al. 1987; Crowther 2007). All of the known dust-forming WR stars, hereafter referred to as “dustars” (Marchenko & Moffat 2007), are of the carbon-rich (WC) subtype. These WC dustars typically exhibit late spectral subtypes (WC7–WC9), which are characterized by relatively cool WR photospheres ($T_* \sim 40,000–70,000$ K; Sander et al. 2019).

In order to explain the observed dust formation in the hostile environment of WC dustars, the binary nature is believed to play a key role: strong winds from the WC star collide with weaker winds from an OB-star companion and create dense regions in the wake of the companion’s orbit that allow for dust to condense (Williams et al. 1990; Usov 1991; Cherchneff 2015). The density enhancements in the wind collision region facilitate rapid radiative cooling that leads to dust formation in addition to shielding newly formed dust grains from UV photons emitted by the central binary. Unambiguous evidence of dust formation regulated by binary influence was presented in high spatial resolution near-IR images of the WC9 binary WR104 by Tuthill et al. (1999), revealing a remarkable dust “pinwheel.”

The binary orbital parameters, such as semimajor axis/orbital period and eccentricity, modulate the dust formation in

WC dustars. For example, systems with low eccentricity and short orbital periods (~yr) like WR104 exhibit persistent dust formation in continuous pinwheel plumes resembling an Archimedean spiral that change in position angle as the WC star and its binary companion move in their orbit (Tuthill et al. 2008). Alternatively, WC dustars with high orbital eccentricity and longer orbital periods ($\gtrsim 5$ yr) like WR140 exhibit periodic dust formation, where the onset of dust formation corresponds to periapse when the WC star and its companion are near their closest orbital separation (Williams et al. 2009a).

Despite the fact that WC dustars can produce copious amounts of dust ($\dot{M}_d \sim 10^{-8}-10^{-6} M_\odot$ yr$^{-1}$; Williams et al. 1987; Hankins et al. 2016; Hendrix et al. 2016) they have been commonly overlooked as significant sources of dust in the ISM of galaxies in the local and early universe. The input from WC dustars compared to leading dust producers like asymptotic giant branch (AGB; Boyer et al. 2012) stars and core-collapse supernovae (SNe; Dwek & Cherchneff 2011) is thought to be low, given their relative rarity as products of only the most massive O-stars. Furthermore, at low metallicities, it is difficult to form a WR star in the context of single star evolution (Conti 1975) due to a lack of metals that drive mass loss and the expulsion of the hydrogen envelope.

However, the influence of their binary nature on the formation of WR stars and their dust formation is not well-understood, especially given observations that reveal a majority ($\gtrsim 70\%$) of massive stars are in binaries that will eventually interact with their
companion over their lifetimes (Sana et al. 2012). Such binary
interaction enables the formation of WR stars through mechan-
isms such as envelope stripping via Roche-Lobe overflow
(Mauerhan et al. 2015; De Marco & Izzard 2017). This provides
formation channels of WC binaries even at low metallicities.
Revisiting the impact of WC dustars is motivated by the need for
additional dust input sources, since SNe may in fact be net dust
destroyers due to the shocks they drive into the surrounding ISM
(Temim et al. 2015).

There are, however, major challenges in addressing the dust
contribution from WC dustars in the ISM of galaxies across
cosmic time. The measured dust production rates of WC
dustars show a wide range of conflicting values in literature
throughout the past several decades. This is due in part to
uncertainties in distance estimates and the different techniques
used to model dust emission, which is also difficult given the
complex dust morphology observed around WC dustars.
Another challenge is that binary stellar evolution tracks are
complex and much more difficult to model than single star
evolution.

In this paper, we study the dust production from WC dustars
by conducting a consistent dust spectral energy distribution
(SED) analysis utilizing archival IR photometry and spectrosc-
opy of a sample of 19 Galactic dustars with distance estimates
from Gaia DR2 (Bailer-Jones et al. 2018; Rate &
Crowther 2020). Based on the results of our SED analysis,
we incorporate dust production in Binary Population and
Spectral Synthesis (BPASS; Eldridge et al. 2017; Stanway &
Eldridge 2018) models to investigate the dust input from WC
dustars relative to AGB stars, red supergiants (RSGs), and
core-collapse SNe at difference metallicities representative of
different epochs in cosmic time. We also model the dust input
in the well-studied Large Magellanic Cloud (LMC), where
AGB stars are claimed to be a significant source of dust and SN
are likely net dust destroyers (Riebel et al. 2012; Temim et al.
2015). This work presents the first implementation of BPASS
to investigate dust production from binary stellar populations.

The paper is organized as follows. First, we describe the
sample selection and the IR photometric and spectroscopic data
sets and extinction correction utilized for the SED modeling
(Section 2). In Section 3, we detail our dust SED modeling
approach with DustEM, present the results of our analysis with
descriptions of selected WC dustars, and compare our result
with previous studies and the theoretical dust formation model
by Usov (1991). In Section 4, we describe the BPASS dust
models and present the results of constant star formation
and coeval population models in environments corresponding
to low (Z = 0.001), LMC-like (Z = 0.008), and solar (Z =
0.020) metallicities. We also provide a comparison between
BPASS model results and the LMC, highlighting the
importance of star formation history in modeling the dust
input from massive stars. Finally, we discuss the astrophysical
implications of our study on the dust production from WC
dustars (Section 5).

## 2. WC Dustar Sample and Observations

2.1. Dustar Sample Selection

The goal of this analysis is to fit dust SED models to the
mid-IR emission from WC dustars, incorporating distance
information from Gaia DR2 (Rate & Crowther 2020) to derive
the dust “shell” radii, mass, production rates, and condensation

fraction. Dust SED models require spectroscopy or well-
sampled photometry between ~2 and 30 µm, where the
thermal emission from circumstellar dust dominates over the
photosphere and free-free emission from ionized winds in
the WC system (e.g., Williams et al. 1987). Beyond ~30 µm,
background emission from cooler ISM dust may start to
contaminate the IR SED. Dustar distances are adopted from the
recent study by Rate & Crowther (2020), who perform a
Bayesian analysis on Gaia DR2 parallaxes to 383 Galactic WR
stars with priors based on H II regions and dust extinction.

Since a non-negligible fraction of known WC dustars exhibit
variability on ≲1 yr timescales (Williams 2019), the known
variables in this sample are restricted to those that have
contemporaneous mid-IR observations or mid-IR spectroscopy
out to at least ~20 µm taken in a single epoch. The dustars
with no observed variability in this sample have ~10–20 µm
spectroscopic or photometric coverage with sufficiently high
spatial resolution at longer wavelengths (FWHM ≲ 6′′) to
distinguish possible contamination from surrounding IR
sources or ISM emission.

Out of the 88 currently known WC dustars in V 1.23 (2019
July) of the Galactic Wolf Rayet Catalogue (Rosslowe &
Crowther 2015),<fnref n="7" /> our sample consists of 19 dustars (shown in
Table 1) that meet the above criteria.

## 2.2. Observations and Archival Data

2.2.1. Mid-IR Imaging of WR48a with VLT/VISIR

Mid-IR imaging observations of WR48a (PID: 097.D-0707
(A); PI—Lau) were performed on the Very Large Telescope
(VLT) at the ESO Paranal observatory using the VLT
spectrometer and imager for the mid-infrared (VISIR; Lagage
et al. 2004) at the Cassegrain focus of UT3 on 2016 June 15.
Images of WR48a presented in this work were taken with the
NeII_2 filter (λ<sub>c</sub> = 13.04 µm, Δλ = 0.22) and obtained using
chopping and nodding to remove the sky and telescope thermal
background emission.

Given the size of the field of view with VISIR of
38" × 38" in comparison to the extent of the dust emission
from WR48a (~5"), an on-detector chop-nod configuration was
used with 20" chop and nod amplitudes. The total integration
time on WR48a was 15 minutes. Raw images were accessed
and downloaded from the ESO Science Archive Facility and
processed using the Modest Image Analysis and Reduction
(MIRA) software written by Terry Herter (see Herter et al.
2013).

The NeII_2 imaging observations of WR48a achieved near-
diffraction-limited imaging with a measured Gaussian full
width at half maximum (FWHM) of 0.45.

2.2.2. Space-based Mid-IR Spectroscopy and Photometry from ISO,
Spitzer, and WISE

ISO/SWS. Five WC dustars in this sample have archival
2.2–40 µm medium-resolution (R ≈ 250–600) spectra taken
by the Short Wavelength Spectrometer (SWS; de Graauw et al.
1996) on the Infrared Space Observatory (ISO; Kessler et al.
1996) that were obtained in the WRSTARS program (PI—van
der Hucht; van der Hucht et al. 1996). Archival ISO/SWS
spectra of these five dustars (WR48a, WR70, WR98a, WR104,
and WR118) were downloaded from the database of SWS

<footnote n="7">⁷ http://pacrowther staff.shef.ac.uk/WRcat/</footnote>
<table><caption>**Table 1**

WC Dustar Table</caption><thead><tr><th>WR Num</th><th>References</th><th>Alt Name</th><th>R.A.</th><th>Decl.</th><th>Spec. Type</th><th>Spec. References</th></tr></thead><tbody><tr><td>19</td><td>VII</td><td>LS 3</td><td>10 18 05.02</td><td>-58 16 25.90</td><td>WC5d</td><td>CDB98, W09b</td></tr><tr><td>48b</td><td>VII</td><td>SMSNPL 8</td><td>13 11 27.45</td><td>-63 46 00.80</td><td>WC9d</td><td>SMS99, VII</td></tr><tr><td>48a</td><td>VII</td><td>D83 1</td><td>13 12 39.65</td><td>-62 42 55.80</td><td>WC8vd+WN8</td><td>Z14</td></tr><tr><td>53</td><td>VII</td><td>HD 117297</td><td>13 30 53.26</td><td>-62 04 51.80</td><td>WC8d</td><td>WH00</td></tr><tr><td>59</td><td>VII</td><td>GSC 9004-3553</td><td>13 49 32.66</td><td>-61 31 42.20</td><td>WC9d</td><td>HH88, VII</td></tr><tr><td>70</td><td>VII</td><td>HIP 75863</td><td>15 29 44.70</td><td>-58 34 51.20</td><td>WC9vd+B0I</td><td>N95, W13a</td></tr><tr><td>95</td><td>VII</td><td>He3-1434</td><td>17 36 19.76</td><td>-33 26 10.90</td><td>WC9d</td><td>VII</td></tr><tr><td>96</td><td>VII</td><td>LSS 4265</td><td>17 36 24.20</td><td>-32 54 29.00</td><td>WC9d</td><td>HH88, VII</td></tr><tr><td>98a</td><td>VII</td><td>IRAS 17380-3031</td><td>17 41 12.90</td><td>-30 32 29.00</td><td>WC8-9vd</td><td>W95</td></tr><tr><td>103</td><td>VII</td><td>HIP 88287</td><td>18 01 43.14</td><td>-32 42 55.20</td><td>WC9d</td><td>VII</td></tr><tr><td>104</td><td>VII</td><td>Ve2-45</td><td>18 02 04.07</td><td>-23 37 41.20</td><td>WC9d+B0.5V</td><td>WH00</td></tr><tr><td>106</td><td>VII</td><td>IC14-8</td><td>18 04 43.66</td><td>-21 09 30.70</td><td>WC9d</td><td>VII</td></tr><tr><td>118</td><td>VII</td><td>GL 2179</td><td>18 31 42.30</td><td>-09 59 15.00</td><td>WC9d</td><td>CV90, VII</td></tr><tr><td>119</td><td>VII</td><td>The 2</td><td>18 39 17.91</td><td>-10 05 31.10</td><td>WC9d</td><td>WH00</td></tr><tr><td>121</td><td>VII</td><td>AS 320</td><td>18 44 13.15</td><td>-03 47 57.80</td><td>WC9d</td><td>VII</td></tr><tr><td>125</td><td>VII</td><td>V378 Vul</td><td>19 28 15.61</td><td>+19 33 21.4</td><td>WC7ed+O9III</td><td>W94, W19</td></tr><tr><td>137</td><td>VII</td><td>HIP 99769</td><td>20 14 31.77</td><td>+36 39 39.60</td><td>WC7pd+O9</td><td>SS90, W01</td></tr><tr><td>140</td><td>VII</td><td>HIP 100287</td><td>20 20 27.98</td><td>+43 51 16.30</td><td>WC7pd+O5</td><td>W90, W09a</td></tr><tr><td>124-22</td><td>KSF15</td><td>1695-2B7</td><td>19 27 17.98</td><td>16 05 24.6</td><td>WC9d</td><td>KSF15, W19</td></tr></tbody></table>

Note. Name, coordinates, and spectral subtype of the 19 WC dustars in our sample obtained from V 1.23 of the Galactic Wolf–Rayet Catalogue. The references correspond to the following: VII—van der Hucht (2001); KSF15—Kanarek et al. (2015); RC18—Rosslowe & Crowther (2018); CDB98—Crowther et al. (1998); W09b—Williams et al. (2009b); SMS99—Shara et al. (1999); Z14—Zhekov et al. (2014); WH00—Williams & van der Hucht (2000); HH88—van der Hucht et al. (1988); N95—Niemela (1995); W13a—Williams et al. (2013b); W95—Williams et al. (1995); CV90—Conti & Vacca (1990); SS90—Smith et al. (1990); W01—Williams et al. (2001); W90—Williams et al. (1990); W94—Williams et al. (1994); W09a—Williams et al. (2009a); W19—Williams (2019).

spectra processed and hosted by Sloan et al. (2003).<fnref n="8" /> The “sws” files, which are the final versions of the SWS spectra corrected for segment discontinuities and overlaps, were used for the SED analysis.

In order to verify the quality of the SWS data and identify possible normalization issues between SWS bands, the “pws” files that show the spectra prior segment-to-segment normalization were inspected. None of the five dustars with ISO/ SWS data exhibit segment-to-segment discontinuities larger than ~10% below 27.5 µm, which is the wavelength where the 3E band of the long-wavelength section starts. Due to the flux discontinuities and larger flux uncertainties beyond 27.5 µm, only the 2.2–27.5 µm data were used for the SED fitting. The spectra of WR48a, WR98a, WR104, and WR118 were smoothed by a median filter with a 51-element kernel, and the spectrum of WR70, which had an “sws” file with wavelengths sampled a factor of four times higher than the other four sources, was smoothed with a 201-element kernel.

**Spitzer/IRS.** Six WC dustars have archival mid-IR spectra taken by the Infrared Spectrograph (IRS; Houck et al. 2004) on the Spitzer Space Telescope (Werner et al. 2004) across multiple programs throughout the Spitzer “cold mission” (2004–2009). Hi-Res (R ~ 600) Spitzer/IRS spectra of WR140 (PID—124; PI—Gehrz), WR19, WR103, and WR53 (all three from PID 199; PI—Houck) were downloaded from the Combined Atlas of Sources with Spitzer IRS Spectra (CASSIS; Lebouteiller et al. 2015).<fnref n="9" /> The spectrum of WR19 was published in the WC dust chemistry study by Marchenko & Moffat (2017), and the spectrum of WR103 was published by Crowther et al. (2006) in their UV to mid-IR spectral analysis comparing WC9 and [WC9] stars. Ardila et al. (2010)

also included the spectra of WR53 and WR103 in their atlas of Spitzer stellar spectra. Broad emission line features from WR103, WR140, and WR19 were masked, and the spectra of WR140 and WR19 with lower signal-to-noise ratios were smoothed by a median filter with an 11-element kernel.

In order to verify the photometric accuracy of the long-wavelength (>20 µm) IRS spectra of WR140, the flux was checked against 24 µm photometry from the Multiband Imaging Photometer for Spitzer (MIPS; Rieke et al. 2004) taken in the same program and within ~1 month of the IRS observations (PID—124; PI—Gehrz). The 24 µm Spitzer/ MIPS photometry of WR140 extracted with a 35"-radius aperture and a 40"-50" background annulus is $F_{24,MIPS} = 1.07 \pm 0.05$.<fnref n="10" /> This is consistent with the IRS spectrum at 24 µm ($F_{24,IRS} \approx 0.8 \pm 0.4$ Jy).

The dustars WR48a and WR98a had coverage from both Spitzer/IRS (PID—40285; PI—Waters) and ISO/SWS, but the ISO spectra were adopted in this work for SED modeling due to its extended wavelength coverage and the sufficiently high signal-to-noise detection.

**Spitzer/IRAC and MIPS.** For the nonvariable dustars, mid-IR Spitzer photometry was adopted from the Galactic Legacy Infrared Mid-Plane Survey Extraordinaire (GLIMPSE; Churchwell et al. 2009) taken by the Infrared Array Camera (IRAC; Fazio et al. 2004) at 3.6, 4.5, 5.8, and 8.0 µm in addition to 24 µm photometry from the MIPS Galactic Plane Survey (MIPSGAL; Carey et al. 2009).

**WISE.** For dustars without any significant contamination or confusion from background emission or nearby sources with ≲10", four-band mid-IR photometry (3.4, 4.6, 12.1, and

<footnote n="8">⁸ https://users.physics.unc.edu/~gcsloan/library/swsatlas/atlas.html</footnote>

<footnote n="9">⁹ https://cassis.sirtf.com/</footnote>

<footnote n="10">¹⁰ An aperture correction factor of 1.08 was included in this flux calculation, in accordance with Table 4.13 of the MIPS Instrument Handbook for the aperture and background annulus radii used.</footnote>
<table><caption>**Table 2**

WC Dustar SED Model Details</caption><thead><tr><th>WR Num</th><th>r<sub>1</sub> Range (au)</th><th>r<sub>1</sub> Intervals</th><th>IR Archival Data</th><th>References</th></tr></thead><tbody><tr><td>19</td><td>(100, 10000)</td><td>50</td><td>Spitzer/IRS</td><td>MM17</td></tr><tr><td>48b</td><td>(50, 3000)</td><td>30</td><td>2MASS, WISE, Spitzer/IRAC+MIPS</td><td>Cu03, Cu13, Ch09, Ca09</td></tr><tr><td>48a</td><td>(10, 3000)</td><td>25</td><td>ISO/SWS</td><td>vdH96</td></tr><tr><td>53</td><td>(50, 3000)</td><td>30</td><td>2MASS, Spitzer/IRAC+IRS</td><td>Cu03, Ch09, Ar10</td></tr><tr><td>59</td><td>(50, 3000)</td><td>30</td><td>2MASS, WISE, Spitzer/IRAC</td><td>Cu03, Cu13, Ca09</td></tr><tr><td>70</td><td>(50, 3000)</td><td>30</td><td>ISO/SWS</td><td>S03</td></tr><tr><td>95</td><td>(10, 1000)</td><td>30</td><td>2MASS, ESO-3.8 m, Spitzer/MIPS</td><td>Cu03, W87, Ca09</td></tr><tr><td>96</td><td>(50, 3000)</td><td>30</td><td>2MASS, WISE, Spitzer/IRAC+MIPS</td><td>Cu03, Cu13, Ch09, Ca09</td></tr><tr><td>98a</td><td>(20, 3000)</td><td>25</td><td>ISO/SWS</td><td>vdH96</td></tr><tr><td>103</td><td>(10, 1000)</td><td>30</td><td>2MASS, ESO-3.8 m, Spitzer/IRS</td><td>Cu03, W87, Ar10</td></tr><tr><td>104</td><td>(20, 3000)</td><td>25</td><td>ISO/SWS</td><td>vdH96</td></tr><tr><td>106</td><td>(10, 1000)</td><td>30</td><td>2MASS, ESO-3.8 m, Spitzer/MIPS</td><td>Cu03, W87, Ca09</td></tr><tr><td>118</td><td>(10, 1000)</td><td>25</td><td>ISO/SWS</td><td>vdH96</td></tr><tr><td>119</td><td>(50, 1000)</td><td>30</td><td>2MASS, ESO-3.8 m, WISE</td><td>Cu03, W87, Cu13</td></tr><tr><td>121</td><td>(10, 1000)</td><td>30</td><td>2MASS, WISE, Spitzer/IRAC+MIPS</td><td>Cu03, Cu13, Ch09, Ca09</td></tr><tr><td>125</td><td>(50, 1000)</td><td>30</td><td>UKIRT</td><td>W94</td></tr><tr><td>124-22</td><td>(20, 3000)</td><td>30</td><td>2MASS, WISE, Spitzer/IRAC+MIPS</td><td>Cu03, Cu13, Ch09, Ca09</td></tr><tr><td>137</td><td>(50, 1000)</td><td>30</td><td>UKIRT</td><td>W01</td></tr><tr><td>140</td><td>(100, 3000)</td><td>30</td><td>TCS, UKIRT, Spitzer/IRS<sup>a</sup></td><td>W09a</td></tr></tbody></table>

**Notes.** The range of $r_1$ values (in au) and number of logarithmic intervals are shown for each dustar model. For the two-component dust models, the $r_2$ fitting parameters is the same as that of $r_1$, and $f(=L_{IR,1}/L_{IR,2})$ is searched in 20 logarithmic intervals between (0.1, 100). Archival IR photometry taken by 2MASS, WISE, and Spitzer were adopted from the MIPSGAL 24 µm point source catalog (Gutermuth & Heyer 2015). The following references are associated with the archival data of each dustar and are indicated in the rightmost column: MM17—Marchenko & Moffat (2017); Cu03—Cutri et al. (2003); Cu13—Cutri (2013); Ch09—Churchwell et al. (2009); Ca09—Carey et al. (2009); vdH96—van der Hucht et al. (1996); Ar10—Ardila et al. (2010); S03—Sloan et al. (2003); W87—Williams et al. (1987); W94—Williams et al. (1994); W01—Williams et al. (2001); W09a—Williams et al. (2009a).

22.0 µm) from the Wide-field Infrared Survey Explorer (WISE; Wright et al. 2010) in the ALLWISE program (Cutri 2013) were adopted. WISE W3 (12.1 µm) photometry is notably valuable since it bridges the wavelength gap between the Spitzer/IRAC Ch4 (8.0 µm) and MIPS 24 µm measurements.

The Spitzer, WISE, and ground-based 2MASS photometry provided in the MIPSGAL 24 µm point source catalog (Gutermuth & Heyer 2015) were primarily utilized for the nonvariable dustars without ISO spectroscopy.

### 2.2.3. Ground-based IR Photometry

Given their brightness in the mid-IR (≥Jy), WC dustars have been targeted by ground-based observatories in the IR for many decades (e.g., Williams et al. 1987). In the near-IR, JHK fluxes from 2MASS (Skrutskie et al. 2006) were incorporated in the SED modeling for the following dustars without ISO/SWS coverage: WR48b, WR53, WR59, WR95, WR96, WR103, WR106, WR119, WR121, and WR124-22. Mid-IR L-, M-, N-, and Q-band photometry published by Williams et al. (1987) taken by the ESO 3.6 m and the United Kingdom Infra-Red Telescope (UKIRT) was used for WR103, WR106, and WR119.

For WR137, one of the known periodic dust-makers, the SED is composed of contemporaneous *J*, *H*, *K*, *L*, *M*, *N*, and *Q* photometry taken in 1985.47 from UKIRT and published by Williams et al. (2001). Near-contemporaneous JHKLMNQ photometry taken by UKIRT in 1993.5 was also adopted for the SED of the periodic dust-maker WR125 (Williams et al. 1994).

In the unique case of the periodic dust-maker WR140, ground-based mid-IR observations with UKIRT by Williams et al. (2009a)

were linked with the space-based spectroscopic observations by Spitzer/IRS. Although near-IR photometry of WR140 was not obtained around the same time as the mid-IR UKIRT and IRS observations, JHK measurements were taken during the previous dust formation epoch at an identical orbital phase ($\phi = 0.43$; Williams et al. 2009a). Importantly, the epoch-to-epoch stability of the near-IR light curve is consistent to within 0.1 mag. The semi-contemporaneous and phase-consistent coverage of WR140 from 1.2 to 35 µm therefore enabled the dust SED analysis of this system.

Table 2 provides a summary of the archival data set used for each of the 19 dustars.

### 2.2.4. IR Extinction Correction

Many of the WC dustars in this sample are heavily extinguished ($A_V \gtrsim 5$) along the line of sight by the ISM, which is apparent from deep silicate absorption in their IR SEDs at 9.7 µm (e.g., Chiar & Tielens 2001). This silicate-dominated extinction is interstellar in nature because the circumstellar dust formed by the WC dustars is carbon-rich (Roche & Aitken 1984). In order to correct for interstellar extinction, the dustar SEDs were dereddened with the “Local ISM” extinction curve from Chiar & Tielens (2006) and adopting the visual extinction, $A_v$ (=1.1$A_V$), derived by Rate & Crowther (2020)<fnref n="11" /> except for WR98a due to its large distance and extinction uncertainties. The visual extinction toward WR98a is instead based on the value from van der Hucht (2001). The v filter ($\lambda_C = 5160 \AA$) is on the narrowband system introduced by Smith (1968) specifically to study WR stars. The adopted $A_v$ for each dustar is shown in Table 3. The

<footnote><sup>a</sup> Spitzer/IRS spectroscopy of WR140 is unpublished.</footnote>

<footnote n="11">¹¹ $A_v$ is referred to as $A_v^{WR}$ in Rate & Crowther (2020).</footnote>
<table><caption>Table 3
Adopted Dustar Parameters</caption><thead><tr><th>WR</th><th>L<sub>*</sub><br>(L<sub>☉</sub>)</th><th>References</th><th>T<sub>*</sub><br>(K)</th><th>Log(R<sub>t</sub>)<br>(R<sub>☉</sub>)</th><th>v<sub>exp</sub><br>(km s<sup>-1</sup>)</th><th>References</th><th>A<sub>v</sub></th><th>M&#775;<br>(M<sub>☉</sub> yr<sup>-1</sup>)</th><th>References</th><th>d<br>(kpc)</th><th>References</th><th>P<sub>orb</sub><br>(yr)</th><th>References</th></tr></thead><tbody><tr><td>19</td><td>4 × 10<sup>5</sup></td><td>S19</td><td>79400</td><td>0.3</td><td>2780</td><td>S19</td><td>5.19</td><td>4.1 × 10<sup>-5</sup></td><td>S19</td><td>4.33<sup>+0.78</sup><sub>-0.58</sub></td><td>RC20</td><td>10.1</td><td>W09b</td></tr><tr><td>48a</td><td>4 × 10<sup>5</sup></td><td>W12</td><td>50000</td><td>0.9</td><td>1000</td><td>this work</td><td>8.28</td><td>1.7 × 10<sup>-4</sup></td><td>Z14<sup>a</sup></td><td>2.27<sup>+0.92</sup><sub>-0.57</sub></td><td>RC20</td><td>32.5</td><td>W12</td></tr><tr><td>48b</td><td>2.5 × 10<sup>5</sup></td><td>S19</td><td>40000</td><td>1</td><td>1390</td><td>S19</td><td>5.73</td><td>2.2 × 10<sup>-5</sup></td><td>S19</td><td>5.12<sup>+1.25</sup><sub>-0.92</sub></td><td>RC20</td><td></td><td></td></tr><tr><td>53*</td><td>3.3 × 10<sup>5</sup></td><td>S19</td><td>50000</td><td>0.9</td><td>1800</td><td>S12</td><td>3.25</td><td>2.2 × 10<sup>-5</sup></td><td>S19</td><td>4.14<sup>+0.74</sup><sub>-0.56</sub></td><td>RC20</td><td></td><td></td></tr><tr><td>59*</td><td>5.8 × 10<sup>5</sup></td><td>S19</td><td>40000</td><td>1</td><td>1300</td><td>S12</td><td>6.43</td><td>3.3 × 10<sup>-5</sup></td><td>S19</td><td>3.57<sup>+0.69</sup><sub>-0.51</sub></td><td>RC20</td><td></td><td></td></tr><tr><td>70</td><td>7 × 10<sup>5</sup></td><td>W13a</td><td>40000</td><td>1</td><td>1390</td><td>S19</td><td>5.20</td><td>2.2 × 10<sup>-5</sup></td><td>S19</td><td>3.01<sup>+0.44</sup><sub>-0.34</sub></td><td>RC20</td><td>~2.8</td><td>W13a</td></tr><tr><td>95*</td><td>1.7 × 10<sup>5</sup></td><td>S19</td><td>45000</td><td>0.9</td><td>1900</td><td>S12</td><td>6.63</td><td>1.9 × 10<sup>-5</sup></td><td>S19</td><td>2.07<sup>+0.43</sup><sub>-0.31</sub></td><td>RC20</td><td></td><td></td></tr><tr><td>96*</td><td>2.5 × 10<sup>5</sup></td><td>S19</td><td>40000</td><td>1</td><td>1390</td><td>S19</td><td>5.48</td><td>2.2 × 10<sup>-5</sup></td><td>S19</td><td>2.64<sup>+0.58</sup><sub>-0.43</sub></td><td>RC20</td><td></td><td></td></tr><tr><td>98a</td><td>1.5 × 10<sup>5</sup></td><td>H16</td><td>45000</td><td>0.9</td><td>900</td><td>W95</td><td>13.79</td><td>2.2 × 10<sup>-5</sup></td><td>S19</td><td>1.9<sup>+0.58</sup><sub>-0.35</sub></td><td>M99</td><td>1.54</td><td>M99</td></tr><tr><td>103*</td><td>3.2 × 10<sup>5</sup></td><td>S19</td><td>45000</td><td>0.8</td><td>1190</td><td>S12</td><td>1.40</td><td>2.8 × 10<sup>-5</sup></td><td>S19</td><td>3.46<sup>+1.28</sup><sub>-0.77</sub></td><td>RC20</td><td></td><td></td></tr><tr><td>104</td><td>2.5 × 10<sup>5</sup></td><td>S19,M07</td><td>40000</td><td>1</td><td>1220</td><td>HS92</td><td>6.67</td><td>3 × 10<sup>-5</sup></td><td>C97</td><td>2.58<sup>+0.12</sup><sub>-0.12</sub></td><td>So18</td><td>0.66</td><td>T08</td></tr><tr><td>106*</td><td>1.7 × 10<sup>5</sup></td><td>S19</td><td>45000</td><td>0.8</td><td>1100</td><td>S12</td><td>4.61</td><td>1.6 × 10<sup>-5</sup></td><td>S19</td><td>3.07<sup>+0.56</sup><sub>-0.43</sub></td><td>RC20</td><td></td><td></td></tr><tr><td>118</td><td>2.5 × 10<sup>5</sup></td><td>S19</td><td>40000</td><td>1</td><td>1390</td><td>S19</td><td>13.68</td><td>2.2 × 10<sup>-5</sup></td><td>S19</td><td>2.49<sup>+0.78</sup><sub>-0.68</sub></td><td>RC20</td><td></td><td></td></tr><tr><td>119*</td><td>0.5 × 10<sup>5</sup></td><td>S19</td><td>45000</td><td>0.8</td><td>1300</td><td>S12</td><td>3.91</td><td>0.7 × 10<sup>-5</sup></td><td>S19</td><td>3.22<sup>+1.24</sup><sub>-0.73</sub></td><td>RC20</td><td></td><td></td></tr><tr><td>121*</td><td>1.4 × 10<sup>5</sup></td><td>S19</td><td>45000</td><td>0.8</td><td>1100</td><td>S12</td><td>5.31</td><td>1.4 × 10<sup>-5</sup></td><td>S19</td><td>2.23<sup>+0.30</sup><sub>-0.24</sub></td><td>RC20</td><td></td><td></td></tr><tr><td>125*</td><td>1.6 × 10<sup>5</sup></td><td>S12</td><td>55000</td><td>1.1</td><td>2000</td><td>S12</td><td>6.48</td><td>2.7 × 10<sup>-5</sup></td><td>S19</td><td>3.36<sup>+0.99</sup><sub>-0.65</sub></td><td>RC20</td><td>28.3</td><td>W19</td></tr><tr><td>137*</td><td>5.4 × 10<sup>5</sup></td><td>S12</td><td>56000</td><td>1</td><td>2000</td><td>S12</td><td>1.70</td><td>2.7 × 10<sup>-5</sup></td><td>S19</td><td>2.10<sup>+0.18</sup><sub>-0.16</sub></td><td>RC20</td><td>13.05</td><td>L05</td></tr><tr><td>140</td><td>1 × 10<sup>6</sup></td><td>W09a</td><td>56000</td><td>1</td><td>~2400</td><td>W09a</td><td>2.21</td><td>2 × 10<sup>-5</sup></td><td>Su15</td><td>1.64<sup>+0.11</sup><sub>-0.09</sub></td><td>RC20</td><td>7.94</td><td>M11</td></tr><tr><td>124-22</td><td>2.5 × 10<sup>5</sup></td><td>S19</td><td>40000</td><td>1</td><td>1390</td><td>S19</td><td>14.77<sup>b</sup></td><td>2.2 × 10<sup>-5</sup></td><td>S19</td><td>1.91<sup>+1.07</sup><sub>-0.72</sub>†</td><td>RC20</td><td></td><td></td></tr></tbody></table>

**Notes.** PoWR models are used for the radiation field of the dust heating source with the following parameters: the stellar heating source luminosity, $L_*$, the effecting temperature of the heating source, $T_*$, the “transformed radius” of the WC star (Sander et al. 2012, 2019), $R_t$ in Log units of $R_☉$. Additional model parameters are the wind velocity/dust expansion velocity, $v_{exp}$, the visual extinction toward the dustar, $A_v$, the mass-loss rate of the WC star, $\dot{M}$, the distance toward the dustar, $d$, and the dustar orbital period, $P_{orb}$. All of the visual extinctions were adopted from Rate & Crowther (2020), except for WR98a, whose extinction value is from van der Hucht (2001). WR names marked with apostrophes indicate systems whose stellar spectra were explicitly modeled with the PoWR grid by Sander et al. (2012), which provides $L_*$, $T_*$, $R_t$, $v_{exp}$, and $\dot{M}$. The † symbol indicates distance values that were flagged for astrometric excess noise >1 mas and large parallax uncertainty by RC20. The provided references correspond to the following: S19—Sander et al. (2019); W12—Williams et al. (2012); Z14—Zhekov et al. (2014); W13a—Williams et al. (2013b); H16—Hendrix et al. (2016); M07—Monnier et al. (2007); So18—Soulain et al. (2018); S12—Sander et al. (2012); W95—Williams et al. (1995); HS92—Howarth & Schmutz (1992); W09b—Williams et al. (2009b); W19—Williams (2019); C97—Crowther (1997); Su15—Sugawara et al. (2015); RC20—Rate & Crowther (2020); M99—Monnier et al. (1999); T08—Tuthill et al. (2008); L05—Lefèvre et al. (2005); M11—Monnier et al. (2011).

a WR48a mass-loss rate from Z14 was adjusted for the revised RC20 distance.

b Derived from $A_{Ks} = 1.58$ (RC20) and assuming $A_v = (0.107)^{-1}A_{Ks}$.
selection of the local ISM extinction curve by Chiar & Tielens (2006) was motivated by their use of WC dustar spectra to define the shape of the 1.25–25 µm extinction. Since several of the dustar spectra extended beyond 25 µm, a consistent $\lambda^{-2}$ power law was used in this work to extrapolate the extinction curve.

# 3. WC Dustar SED Modeling Results and Analysis

## 3.1. WC Dustar SED Modeling Overview

The complex “pinwheel” morphology of the spatially resolved circumstellar dust observed around several well-studied dustars such as WR104 (Tuthill et al. 1999; Soulain et al. 2018) and WR98a (Monnier et al. 1999) presents a complicated scenario for radiative transfer modeling. For example, Hendrix et al. (2016) performed radiative transfer dust models on 3D hydrodynamic simulations fitting the observed emission and morphology of the dust plume around WR98a. Alternatively, Williams et al. (1987) had fit simple analytic dust shell models to dustar SEDs assuming their dust emission is optically thin. Zubko (1998) provided an updated SED analysis using a detailed theoretical model of dust shells that accounts for grain formation physics and dynamics on a sample of Galactic WC dustars with photometry from Williams et al. (1987). Although the Zubko (1998) SED models were performed before the circumstellar emission from WC dustars were spatially resolved, the derived dust production rates and carbon condensation fraction provide valuable benchmarks for comparison.

We performed dust SED fits to the 19 WC dustars in our sample with the tool DustEM (Compiègne et al. 2011) with single or double dust component models. DustEM is a numerical tool that computes the dust emission in the optically thin limit heated by an input radiation field with no radiative transfer. Single dust component models were initially attempted for all dustars, and double components models were used for SEDs with unsatisfactory single-component fits. The two different SED modeling methods are described as follows.

1. Single-component dustar models are assumed to have a single geometrically thin dust shell centered on the WC +OB binary.

2. Double-component dustars are modeled with two geometrically thin dust ring components with different radii from the central binary, utilizing the technique described below in Section 3.2 in greater detail.

Note that, in the single-component model, the morphology of the dust shell (e.g., torus versus spherical shell) does not affect the dust emission spectrum.

In addition to circumstellar dust, the IR excess from WC dustars can originate from free-free emission in their ionized winds (Cohen et al. 1975). The IR excess due to free-free emission can be characterized by a power-law $F_{\text{ff}} \propto \lambda^{-0.96}$ (Morris et al. 1993) and may therefore dominate at shorter IR wavelengths such as the J-band ($\lambda = 1.25$ µm). In order to remove contamination from free-free emission in the SEDs, we subtract this free-free emission power-law model normalized at the J-band flux for dustars, where the J-band flux is $\geq 10\%$ of the mid-IR flux peak associated with circumstellar dust. In the unique case of WR140, where the emission from the central system has been spatially resolved from extended circumstellar dust, an interpolation is performed between the IR photometry

<figure><img src="image_2.png" /><figcaption>**Figure 1.** Top: Two-component dust model cross-section schematic. Bottom: Face-on “pinwheel” dust morphology schematic of continuous WC dustars overlaid with the approximated position of the two dust components in this model.</figcaption></figure>

of the stellar wind continuum measured by Williams et al. (2009a) to characterize and subtract this component.

## 3.2. Two-component Dust Model

In this section, we describe the framework of the two-component dust model. The “pinwheel” dust morphology from these systems is approximated as two concentric dust rings. Figure 1 presents a schematic illustration of the different components in this model and the comparison to the pinwheel dust, where the first dust component corresponds to the first dust spiral “coil” and the second component corresponds to the second continuation of this coil. In these models, the inner dust ring attenuates the radiation impinging on the outer ring and is described as follows.

The incident stellar radiative flux on dust components 1 and 2 at distances $r_1$ and $r_2$ (where $r_1 < r_2$) from the central heating source, respectively, can be described as

$$F_1 = \frac{L_*}{4\pi r_1^2} \quad \text{and} \quad F_2 = \frac{L_* e^{-\tau}}{4\pi r_2^2}, \qquad (1)$$

where $\tau$ is the optical depth through component 1. Deriving the radiative heating of component 1 is therefore straightforward since it only requires information of the heating source radiation field and its distance to the heating source. However, the radiative heating of component 2 requires information on the attenuated stellar flux through component 1. This attenuation can be related to the total IR luminosity radiated by
component 1, $L_{IR,1}$:

$$L_{IR,1} = L_{*}(1 - e^{-\tau})\Omega, \quad (2)$$

where $\Omega$ is the geometric coverage fraction of components 1 and 2 around the central heating source and it is assumed both components have the same coverage fraction (Figure 1). From Equation (2), it follows that $e^{-\tau} = 1 - \frac{L_{IR,1}}{L_{*}\Omega}$, and the incident flux on component 2 can then be expressed as

$$F_2 = \frac{1}{4\pi r_2^2} \left( L_* - \frac{L_{IR,1}}{\Omega} \right), \quad (3)$$

We assume that 100% of the stellar flux impinging on components 1 and 2 is absorbed and reradiated into the IR by the two components. This assumption is motivated by spatially resolved observations that reveal the IR emission from WC dustars are dominated by emitting regions within the first two pinwheel arcs (e.g., Tuthill et al. 2008). This implies that

$$\Omega = \frac{L_{IR,1} + L_{IR,2}}{L_{*}}, \quad (4)$$

from which it follows from Equation (3) that $F_2$ can be re-expressed as

$$F_2 = \frac{L_*}{4\pi r_2^2} \left( 1 - \frac{1}{1 + f^{-1}} \right), \quad (5)$$

where $f \equiv L_{IR,1}/L_{IR,2}$. With Equations (1) and (5), we have therefore arrived at a simplified two-component pseudo-radiative transfer model with $L_*$, $f$, $r_1$, and $r_2$ as the model parameters.

In a test for limit cases, it is apparent that when component 1 absorbs all of the stellar flux ($f \to \infty$), $F_2 = 0$. Conversely, when component 1 does not absorb any of the stellar flux ($f \to 0$), $F_2 = \frac{L_*}{4\pi r_2^2}$.

This two-component dust model is used for four dustars: WR48a, WR98a, WR104, and WR118.

### 3.3. Adopted and Derived SED Model Parameters

In our DustEM SED modeling, we input the luminosity and radiation field of the heating source, the distance to the dust components, dust composition, and the dust grain size distribution as model inputs. For the two-component models, we also input the IR luminosity ratio of the two dust components, $f$ (see Equation (5)).

We adopt the appropriate radiation field that is representative of the spectral subtype of each WC dustar from the Potsdam Wolf–Rayet Star (PoWR) model atmospheres (Gräfener et al. 2002; Sander et al. 2012, 2019). The radiation field is then scaled by the luminosity of the heating source and the distance to the dust components assuming the incident flux decreases as $F \propto r^{-2}$ (Equation (1)). While the radiation field and luminosity are adopted from previous studies, we performed a logarithmic grid search for the dust component distances and (for the two-component cases) the IR luminosity ratio of the two dust components. The grid search range and intervals for the different free parameters of each dustar model are provided in Table 2. Table 2 also lists the observatories from which the IR archival data were obtained for each dustar in the sample.

The dust grain size distribution in WC dustars is still disputed and is one of the issues addressed by this work. Theoretical studies on WC dustars winds by Zubko (1998)

predict that dust grains only grow up to sizes of $a \sim 0.01$ µm, which has been corroborated by dust SED analysis of observed IR emission from the dustar WR140 (Williams et al. 2009a). However, IR dust emission and extinction studies of WC dustars by different groups suggest the presence of large $a \gtrsim 0.5$ µm dust grains (Chiar & Tielens 2001; Marchenko et al. 2002; Rajagopal et al. 2007).

In order to address this grain size discrepancy, we test both small and large grain size distributions for the WC dustars in our sample that have been spatially resolved. Since heat capacity increases with increasing grain size, large grains need to be at a closer distance to the radiative heating source than small grains do in order to exhibit a similar temperature and spectral shape. Therefore, we can compare which grain size distribution in our SED models provides dust shell distances consistent with the observed circumstellar dust morphology constraints. The small (large) grain size distribution includes dust grains ranging from $a = 0.01-0.1$ µm (0.1–1.0 µm) with a number density distribution proportional to $n(a) \propto a^{-3}$. We note that the grain size has a minimal effect on the derived dust mass because this quantity is fit from the longer-wavelength IR emission ($\gtrsim 10$ µm) in the SED where the dust opacity is not as sensitive as the shorter-wavelength IR emission (e.g., Harries et al. 2004) with regard to grain size.

We assume that the dust grains are composed purely of amorphous carbon, which is consistent with their featureless IR spectra and the C-rich environment in the vicinity of the WC star (Williams et al. 1987; van der Hucht et al. 1996). We adopt the DustEM amorphous carbon (“amCBEx”) grains with refractive indices derived by Zubko et al. (1996) and Compiègne et al. (2011) and assume a bulk density of $\rho_b = 2.0$ g cm$^{-3}$.

By performing a least-squares fit of the DustEM model to the WC dustar SEDs, we derive the dust mass ($M_d$), dust temperature ($T_d$), and distances between the dust components and central system ($r$). For the continuous dust-makers with a dust expansion velocity of $v_{exp}$, the dust production rate can then be approximated by

$$(\text{continuous}) \dot{M}_d \sim \frac{M_d v_{exp}}{r}, \quad (6)$$

where $v_{exp}$ is assumed to be comparable to the stellar wind velocity $v_*$ or is derived directly from dust proper motion measurements in multi-epoch imaging.

For the periodic, noncontinuous dust-makers with a recurring dust formation timescale $P$ associated with the orbital period of the central binary (e.g., Williams et al. 2009a; Monnier et al. 2011), we approximate the dust production rate by

$$(\text{periodic}) \dot{M}_d \sim \frac{M_d}{P}. \quad (7)$$

Equation (7) effectively calculates the dust production rate averaged over the orbital period, whereas Equation (6) is an “instantaneous” dust production rate. For periodic, noncontinuous dust-makers, the instantaneous dust production rate calculation will vary depending on the phase the observations. The orbital period-averaged dust production rate calculation (Equation (7)) is therefore applied to the periodic, noncontinuous dustars WR 140, WR 137, WR125, and WR 19, while Equation (6) is applied to the other 15 dustars. Note that, for continuous dust-makers, the dust production rate averaged over
an orbital period will be the same as the instantaneous dust production rate. In the two-component dust models, we take a conservative approach and estimate the dust production rates based only on the dust mass and the distance to the first dust component, because cooler and more extended dust may contribute to the dust mass determined for the second component.

Finally, the fraction of carbon in the WC winds that condense into dust, $\chi_C$, is estimated based on the dust production rate, the adopted mass-loss rate for the WC wind, and an assumption that the WC wind is composed of 40% carbon by mass (Sander et al. 2019). The adopted WC dustar properties are summarized in Table 3.

### 3.4. WC Dustar SED Modeling Results

In this section, we present the results of our DustEM SED models for the 19 WC dustars in our sample. First, as a verification of our modeling approach, we highlight the results on the well-studied periodic dust-maker WR140, the continuous “pinwheel” dust-maker WR104, and the 32.5 yr orbital period dustar WR48a. Next, we present and describe the results on the remaining IR-luminous ($L_{IR}/L_* \gtrsim 0.1$) dustars, a selection of the single-component dustars, and additional systems with angular size constraints. Last, we discuss the overall results on dust formation and grain size properties, compare results from previous studies, and test our results against the theoretical model of dust formation in colliding winds by Usov (1991). The results from the DustEM SED models are summarized in Table 4.

#### 3.4.1. The Archetypal Periodic Dust-maker WR140

WR140 is a WC7+O5 system with a well-defined 7.94 yr orbital period, a high orbital eccentricity of $e = 0.90$, and a total luminosity of $L_* \sim 1 \times 10^6 L_\odot$ (Williams et al. 2009a; Monnier et al. 2011). The distance to WR140 measured by Gaia parallax is $d = 1640_{-90}^{+110}$ pc (Rate & Crowther 2020), which is consistent with previous distance estimates (Dougherty et al. 2005; Monnier et al. 2011).

Dust formation in WR140 only occurs when the binary system is at periapse. The proper motion of the resulting dust “arc” formed during periapse has been characterized by multi-epoch, spatially resolved mid-IR imaging as the arc propagates away from the central binary (Figure 2, Left; Williams et al. 2009a). These images allows us to verify our SED model-derived separation distances between the dust component and central heating source. For the line-of-sight extinction, we adopt $A_v = 2.21$ (Rate & Crowther 2020).

We perform single-component DustEM models to the stellar wind-subtracted 5–37 µm Spitzer/IRS spectrum taken on 2004 May, UKIRT L- and M-band imaging taken on 2004 June, and JHK imaging taken on 1996 August where WR140 was at an orbital phase ($\phi = 0.43$) nearly identical to that of the UKIRT and Spitzer observations (Figure 2, right). This orbital phase corresponds to 3.4 yr after IR maximum and periastron passage. The small (large) dust grain distribution SED fit provides a separation distance between the dust component and central heating source of $r = 1170_{-340}^{+310}$ au ($360_{-100}^{+160}$ au). Spatially resolved mid-IR imaging observations of WR140 at this orbital phase show that the extended dust arc is 700–900 mas from the central source (Williams et al. 2009a), which is closely consistent with the separation distance derived

from our small grain dust distribution model, 1170 au $\approx 700$ mas (Figure 2, left). This small grain distribution is also consistent with the dust emission analysis by Williams et al. (2009a), who deduced a characteristic grain size of 0.01 µm.

From the small grain dust SED, we derive a total dust mass of $M_d = (6.44_{-3.30}^{+3.84}) \times 10^{-9} M_\odot$, which implies a dust production rate of $\dot{M}_d = (8.1_{-4.2}^{+4.8}) \times 10^{-10} M_\odot \text{ yr}^{-1}$ (Equation (7)), given the periodic 7.94 yr dust formation timescales in WR140. SED fits from ground-based observations by Williams et al. (2009a) indicate a total dust mass of less than $2 \times 10^{-8} M_\odot$ at a phase of 0.56 around WR140, which consistently constrains our derived dust mass. Adopting a total mass-loss rate of $\dot{M} \approx 2 \times 10^{-5} M_\odot \text{ yr}^{-1}$ derived by X-ray monitoring observations of shocked gas in the wind collision region between the WC and O stars (Sugawara et al. 2015) and assuming a 40% carbon composition by mass, the carbon dust condensation fraction is $\chi_C \approx 0.01\%$.

#### 3.4.2. The Continuous “Pinwheel” Dust-maker WR104

WR104 is the prototypical, continuous dust-maker that exhibits a nearly face-on pinwheel (Figure 3, left) with an orbital period of 242 days (Tuthill et al. 1999, 2008). The central binary is composed of a WC9+OB star and has a Gaia-derived distance of $d = 2740_{-550}^{+720}$ pc. This is consistent with the well-defined distance determined from the location of its possible stellar association host, Sgr OB1, and observations of the proper motion of its circumstellar dust ($d = 2580 \pm 120$ pc; Soulain et al. 2018). Given the tighter constraints from their high spatial resolution imaging, we adopt the Soulain et al. (2018) distance for the dust SED model. The line-of-sight extinction toward WR104 is $A_v = 6.67$ (Rate & Crowther 2020). We adopt a total system luminosity of $L_* = 2.5 \times 10^5 L_\odot$, consistent with Monnier et al. (2007) and the mean WC9 luminosity derived by Sander et al. (2019).

Since single-component models fail to fit the ISO/SWS spectrum of WR104 between 2 and 27.5 µm, we apply the double-component SED models (Figure 3, Right) and derive separation distances between the central binary and dust components of $\eta = 240_{-80}^{+60}$ au ($70_{-10}^{+20}$ au) and $r_2 = 1300_{-740}^{+1130}$ au ($460_{-90}^{+110}$ au) for the small (large) dust distribution models. High spatial resolution IR imaging by Tuthill et al. (2008) shows that the distance between the central binary and the onset of the spiral dust plume where the IR emission peaks is $\sim 15$ mas (40 au). Figure 3 (left) shows the approximate location of the first dust component for small and large grain models compared to the observed dust morphology. Given the closer agreement of the first dust component, we favor the large grain distribution for WR104.

From our large grain dust SED model, we derive dust masses of $M_{d1} = (1.19_{-0.44}^{+0.70}) \times 10^{-6} M_\odot$ and $M_{d2} = (6.06_{-2.06}^{+3.09}) \times 10^{-5} M_\odot$ for the two dust components. Adopting a dust expansion velocity of $v_{exp} = 1220$ km s$^{-1}$ (Howarth & Schmutz 1992), the dust production rate based on component 1 is $\dot{M}_d = (4.39_{-0.97}^{+1.27}) \times 10^{-6} M_\odot \text{ yr}^{-1}$.

This dust production rate is a factor of $\sim 5$ higher than the $\dot{M}_d = 8 \pm 1 \times 10^{-7} M_\odot \text{ yr}^{-1}$ value derived by Harries et al. (2004) from 3D dust radiative transfer models. However, we can reconcile this discrepancy based on the following points:

1. Harries et al. (2004) adopt a distance to WR104 of 1600 pc as opposed to our value of 2580 pc, which accounts
<table><caption>Table 4
WC Dustar SED Model Results</caption><thead><tr><th>WR</th><th>r<sub>1</sub> (au)</th><th>r<sub>2</sub> (au)</th><th>T<sub>d1</sub> (K)</th><th>T<sub>d2</sub> (K)</th><th>d<sub>unc</sub> Factor</th><th>L<sub>IR</sub> (L<sub>☉</sub>)</th><th>L<sub>IR</sub>/L<sub>☉</sub> (%)</th><th>L<sub>IR,1</sub>/L<sub>IR,2</sub></th><th>M<sub>d1</sub> (M<sub>☉</sub>)</th><th>M<sub>d2</sub> (M<sub>☉</sub>)</th><th>Mdot (M<sub>☉</sub> yr<sup>-1</sup>)</th><th>χ<sub>C</sub> (%)</th></tr></thead><tbody><tr><td>19 (SG)</td><td>6870<sup>+3130</sup><sub>-3630</sub></td><td></td><td>250<sup>+80</sup><sub>-30</sub></td><td></td><td>(0.75, 1.39)</td><td>(5.13<sup>+1.48</sup><sub>-0.43</sub>)e + 0</td><td>0.00128<sup>+0.000037</sup><sub>-0.00011</sub></td><td></td><td>(3.98<sup>+3.76</sup><sub>-2.84</sub>)e-8</td><td></td><td>(3.94<sup>+3.72</sup><sub>-2.81</sub>)e-9</td><td>0.024<sup>+0.023</sup><sub>-0.017</sub></td></tr><tr><td>48a (SG)</td><td>350<sup>+100</sup><sub>-130</sub></td><td>2370<sup>+640</sup><sub>-900</sub></td><td>770<sup>+150</sup><sub>-70</sub></td><td>310<sup>+90</sup><sub>-40</sub></td><td>(0.56, 1.97)</td><td>(1.05<sup>+0.07</sup><sub>-0.01</sub>)e + 4</td><td>2.63<sup>+0.18</sup><sub>-0.03</sub></td><td>1.83<sup>+0.80</sup><sub>-0.95</sub></td><td>(1.42<sup>+1.12</sup><sub>-0.99</sub>)e-7</td><td>(9.82<sup>+5.84</sup><sub>-5.93</sub>)e-6</td><td>(8.46<sup>+3.48</sup><sub>-4.38</sub>)e-8</td><td>0.12<sup>+0.05</sup><sub>-0.06</sub></td></tr><tr><td>48b</td><td>180<sup>+90</sup><sub>-60</sub></td><td></td><td>910<sup>+150</sup><sub>-130</sub></td><td></td><td>(0.67, 1.55)</td><td>(3.43<sup>+0.66</sup><sub>-0.87</sub>)e + 2</td><td>0.14<sup>+0.03</sup><sub>-0.03</sub></td><td></td><td>(2.89<sup>+1.96</sup><sub>-1.36</sub>)e-9</td><td></td><td>(4.71<sup>+0.56</sup><sub>-0.96</sub>)e-9</td><td>0.053<sup>+0.006</sup><sub>-0.011</sub></td></tr><tr><td>53</td><td>240<sup>+30</sup><sub>-30</sub></td><td></td><td>860<sup>+40</sup><sub>-40</sub></td><td></td><td>(0.75, 1.39)</td><td>(7.27<sup>+0.34</sup><sub>-0.29</sub>)e + 2</td><td>0.22<sup>+0.01</sup><sub>-0.01</sub></td><td></td><td>(8.43<sup>+1.81</sup><sub>-1.67</sub>)e-9</td><td></td><td>(1.33<sup>+0.11</sup><sub>-0.11</sub>)e-8</td><td>0.152<sup>+0.12</sup><sub>-0.013</sub></td></tr><tr><td>59</td><td>210<sup>+30</sup><sub>-30</sub></td><td></td><td>1010<sup>+60</sup><sub>-50</sub></td><td></td><td>(0.73, 1.42)</td><td>(7.84<sup>+0.96</sup><sub>-0.79</sub>)e + 2</td><td>0.14<sup>+0.02</sup><sub>-0.01</sub></td><td></td><td>(3.91<sup>+0.68</sup><sub>-0.69</sub>)e-9</td><td></td><td>(5.10<sup>+0.14</sup><sub>-0.19</sub>)e-9</td><td>0.039<sup>+0.001</sup><sub>-0.001</sub></td></tr><tr><td>70</td><td>550<sup>+80</sup><sub>-70</sub></td><td></td><td>730<sup>+40</sup><sub>-40</sub></td><td></td><td>(0.79, 1.31)</td><td>(2.36<sup>+0.07</sup><sub>-0.07</sub>)e + 3</td><td>0.34<sup>+0.01</sup><sub>-0.01</sub></td><td></td><td>(6.60<sup>+2.40</sup><sub>-1.78</sub>)e-8</td><td></td><td>(3.51<sup>+0.65</sup><sub>-0.56</sub>)e-8</td><td>0.399<sup>+0.074</sup><sub>-0.063</sub></td></tr><tr><td>95 (LG)</td><td>40<sup>+70</sup><sub>-20</sub></td><td></td><td>1110<sup>+300</sup><sub>-340</sub></td><td></td><td>(0.72, 1.46)</td><td>(1.72<sup>+0.48</sup><sub>-0.96</sub>)e + 3</td><td>1.01<sup>+0.28</sup><sub>-0.57</sub></td><td></td><td>(0.93<sup>+2.18</sup><sub>-0.63</sub>)e-8</td><td></td><td>(9.36<sup>+2.00</sup><sub>-3.36</sub>)e-8</td><td>1.232<sup>+0.263</sup><sub>-0.443</sub></td></tr><tr><td>96</td><td>180<sup>+130</sup><sub>-60</sub></td><td></td><td>910<sup>+150</sup><sub>-170</sub></td><td></td><td>(0.70, 1.49)</td><td>(1.08<sup>+0.29</sup><sub>-0.45</sub>)e + 3</td><td>0.43<sup>+0.12</sup><sub>-0.18</sub></td><td></td><td>(9.09<sup>+6.67</sup><sub>-3.97</sub>)e-9</td><td></td><td>(1.48<sup>+0.03</sup><sub>-0.23</sub>)e-8</td><td>0.168<sup>+0.003</sup><sub>-0.026</sub></td></tr><tr><td>98a (LG)</td><td>60<sup>+10</sup><sub>-10</sub></td><td>370<sup>+190</sup><sub>-130</sub></td><td>960<sup>+70</sup><sub>-70</sub></td><td>340<sup>+80</sup><sub>-80</sub></td><td>(0.67, 1.70)</td><td>(1.74<sup>+0.05</sup><sub>-0.07</sub>)e + 4</td><td>11.63<sup>+0.32</sup><sub>-0.46</sub></td><td>5.46<sup>+5.83</sup><sub>-1.66</sub></td><td>(1.82<sup>+1.08</sup><sub>-0.68</sub>)e-7</td><td>(9.25<sup>+11.73</sup><sub>-5.29</sub>)e-6</td><td>(6.10<sup>+1.77</sup><sub>-1.38</sub>)e-7</td><td>6.93<sup>+2.02</sup><sub>-1.57</sub></td></tr><tr><td>103</td><td>240<sup>+40</sup><sub>-40</sub></td><td></td><td>860<sup>+50</sup><sub>-50</sub></td><td></td><td>(0.60, 1.88)</td><td>(4.92<sup>+0.38</sup><sub>-0.33</sub>)e + 2</td><td>0.15<sup>+0.01</sup><sub>-0.01</sub></td><td></td><td>(5.83<sup>+1.63</sup><sub>-1.29</sub>)e-9</td><td></td><td>(6.10<sup>+0.56</sup><sub>-0.52</sub>)e-9</td><td>0.054<sup>+0.005</sup><sub>-0.005</sub></td></tr><tr><td>104 (LG)</td><td>70<sup>+20</sup><sub>-10</sub></td><td>460<sup>+110</sup><sub>-90</sub></td><td>970<sup>+80</sup><sub>-70</sub></td><td>340<sup>+50</sup><sub>-40</sub></td><td>(0.91, 1.10)</td><td>(1.23<sup>+0.04</sup><sub>-0.04</sub>)e + 5</td><td>49.01<sup>+1.71</sup><sub>-1.48</sub></td><td>5.46<sup>+2.39</sup><sub>-1.66</sub></td><td>(1.19<sup>+0.70</sup><sub>-0.44</sub>)e-6</td><td>(6.06<sup>+3.09</sup><sub>-2.06</sub>)e-5</td><td>(4.39<sup>+1.27</sup><sub>-0.97</sub>)e-6</td><td>36.57<sup>+1.055</sup><sub>-8.07</sub></td></tr><tr><td>106 (LG)</td><td>50<sup>+120</sup><sub>-20</sub></td><td></td><td>1020<sup>+200</sup><sub>-370</sub></td><td></td><td>(0.74, 1.40)</td><td>(7.40<sup>+0.71</sup><sub>-4.36</sub>)e + 3</td><td>4.35<sup>+0.42</sup><sub>-2.58</sub></td><td></td><td>(6.39<sup>+23.69</sup><sub>-3.87</sub>)e-8</td><td></td><td>(2.97<sup>+1.14</sup><sub>-1.01</sub>)e-7</td><td>4.635<sup>+1.780</sup><sub>-1.586</sub></td></tr><tr><td>118 (LG)</td><td>30<sup>+10</sup><sub>-10</sub></td><td>260<sup>+120</sup><sub>-80</sub></td><td>1370<sup>+190</sup><sub>-170</sub></td><td>420<sup>+90</sup><sub>-80</sub></td><td>(0.53, 1.72)</td><td>(4.12<sup>+0.34</sup><sub>-0.26</sub>)e + 4</td><td>16.48<sup>+1.37</sup><sub>-1.04</sub></td><td>5.46<sup>+2.39</sup><sub>-1.66</sub></td><td>(5.58<sup>+6.24</sup><sub>-2.95</sub>)e-8</td><td>(6.61<sup>+6.99</sup><sub>-3.42</sub>)e-6</td><td>(6.27<sup>+2.78</sup><sub>-1.94</sub>)e-7</td><td>7.13<sup>+3.16</sup><sub>-2.20</sub></td></tr><tr><td>119</td><td>70<sup>+60</sup><sub>-30</sub></td><td></td><td>960<sup>+220</sup><sub>-200</sub></td><td></td><td>(0.60, 1.92)</td><td>(3.93<sup>+1.45</sup><sub>-1.19</sub>)e + 2</td><td>0.79<sup>+0.29</sup><sub>-0.24</sub></td><td></td><td>(2.50<sup>+3.52</sup><sub>-1.38</sub>)e-9</td><td></td><td>(9.80<sup>+2.90</sup><sub>-2.14</sub>)e-9</td><td>0.350<sup>+0.104</sup><sub>-0.076</sub></td></tr><tr><td>121</td><td>130<sup>+150</sup><sub>-50</sub></td><td></td><td>930<sup>+180</sup><sub>-240</sub></td><td></td><td>(0.80, 1.29)</td><td>(1.04<sup>+0.09</sup><sub>-0.54</sub>)e + 3</td><td>0.74<sup>+0.07</sup><sub>-0.38</sub></td><td></td><td>(7.92<sup>+9.90</sup><sub>-4.73</sub>)e-9</td><td></td><td>(1.79<sup>+0.13</sup><sub>-0.62</sub>)e-8</td><td>0.319<sup>+0.022</sup><sub>-0.110</sub></td></tr><tr><td>125 (SG)</td><td>490<sup>+170</sup><sub>-100</sub></td><td></td><td>570<sup>+50</sup><sub>-60</sub></td><td></td><td>(0.65, 1.68)</td><td>(9.78<sup>+0.32</sup><sub>-1.11</sub>)e + 2</td><td>0.61148<sup>+0.02005</sup><sub>-0.06943</sub></td><td></td><td>(1.00<sup>+0.61</sup><sub>-0.35</sub>)e-7</td><td></td><td>(3.54<sup>+2.15</sup><sub>-1.22</sub>)e-9</td><td>0.033<sup>+0.020</sup><sub>-0.011</sub></td></tr><tr><td>137 (SG)</td><td>660<sup>+240</sup><sub>-220</sub></td><td></td><td>630<sup>+110</sup><sub>-70</sub></td><td></td><td>(0.85, 1.18)</td><td>(2.02<sup>+0.37</sup><sub>-0.35</sub>)e + 2</td><td>0.03704<sup>+0.00686</sup><sub>-0.00652</sub></td><td></td><td>(1.20<sup>+0.64</sup><sub>-0.57</sub>)e-8</td><td>(9.21<sup>+4.90</sup><sub>-4.36</sub>)e-10</td><td>0.009<sup>+0.005</sup><sub>-0.004</sub></td><td></td></tr><tr><td>140 (SG)</td><td>1170<sup>+310</sup><sub>-340</sub></td><td></td><td>570<sup>+80</sup><sub>-50</sub></td><td></td><td>(0.89, 1.14)</td><td>(6.40<sup>+0.02</sup><sub>-0.19</sub>)e + 1</td><td>0.00640<sup>+0.00002</sup><sub>-0.00019</sub></td><td></td><td>(6.44<sup>+3.84</sup><sub>-3.30</sub>)e-9</td><td>(8.11<sup>+4.83</sup><sub>-4.15</sub>)e-10</td><td>0.010<sup>+0.006</sup><sub>-0.005</sub></td><td></td></tr><tr><td>124-22</td><td>190<sup>+80</sup><sub>-60</sub></td><td></td><td>900<sup>+140</sup><sub>-110</sub></td><td></td><td>(0.39, 2.43)</td><td>(1.03<sup>+0.10</sup><sub>-0.20</sub>)e + 1</td><td>0.00412<sup>+0.00049</sup><sub>-0.00071</sub></td><td></td><td>(9.67<sup>+6.49</sup><sub>-4.60</sub>)e-11</td><td>(1.49<sup>+0.26</sup><sub>-0.35</sub>)e-10</td><td>0.0017<sup>+0.0003</sup><sub>-0.0004</sub></td><td></td></tr></tbody></table>

Note. Best-fit DustEM SED models provide the radius of component 1, r₁, the dust temperature of component 1, T<sub>d1</sub>, the total IR luminosity of emitting dust, L<sub>IR</sub>, the fraction of stellar luminosity reradiated in the IR by dust, L<sub>IR</sub>/L<sub>*</sub>, the dust mass in component 1, M<sub>d1</sub>, the dust production rate, Mdot, and the mass fraction of carbon in the WC wind condensed into dust, χ<sub>C</sub>. For the two-component dust models, the radius (r<sub>2</sub>), temperature (T<sub>d2</sub>), dust mass (M<sub>d2</sub>), and luminosity ratio L<sub>IR,1</sub>/L<sub>IR,2</sub> for component 2 are fit. The T<sub>d</sub> values correspond to the dust temperature of the smallest grain component in the grain size distribution: a = 0.01 µm and 0.1 µm for the small and large distributions, respectively. The ±1σ uncertainties from the model fitting are shown for each derived parameter. The ±1σ uncertainties in the adopted distance are provided as multiplicative factors in the “d<sub>unc</sub> Factor” column which applies to L<sub>*</sub>, L<sub>IR</sub>/L<sub>*</sub>, M<sub>d1</sub>, M<sub>d2</sub>, Mdot, and χ<sub>C</sub>. Note that r₁, r₂, T<sub>d1</sub>, and T<sub>d2</sub> are independent of the adopted dustar distance. The WR names with “(LG)” or “(SG)” indicate models where the large or small grain size distribution were favored based on spatial size constraints. All other dustars SED models used the small grain size distribution.
<figure><img src="image_3.png" /><figcaption>**Figure 2.** Left: Mid-IR 12.5 µm image of WR140 taken by Gemini/Michelle (Williams et al. 2009a) in 2003 December, overlaid with circles centered on WR140 representing the fitted distances of emitting dust for the large (solid) and small (dashed) grain size distribution models. The small grain dust model best matches the observed location of the dust plume. Right: Best-fit *DustEM* SED model of WR140 using the small grain size distribution fit to UKIRT L- and M-band photometry taken in 2004 June (Williams et al. 2009a) and Spitzer/IRS spectroscopy taken in 2004 May. Gray bars correspond to the 1σ flux uncertainty.</figcaption></figure>

<figure><img src="image_4.png" /></figure>

<figure><img src="image_5.png" /><figcaption>**Figure 3.** Left: Near-IR reconstructed interferometric image of WR104 taken with the CH4 ($\lambda_c$ = 2.27 µm) filter on the Keck I Telescope reproduced from Figure 1 in Tuthill et al. (2008), centered on the origin of their best-fit Archimedean spiral model (dotted line) and overlaid with red circles representing the location of the fitted dust component 1 for the large (solid) and small (dashed) grain models. Contour levels are 0.4%, 1%, 2%, 5%, 10%, and 50% of the peak emission. The large grain dust model best matches the observed location of the centralized peak dust emission. Right: Best-fit two-component *DustEM* SED model of WR104 using the large grain size distribution fit to ISO/SWS spectroscopy. Gray bars correspond to the 1σ flux uncertainty.</figcaption></figure>

<figure><img src="image_6.png" /></figure>

for a factor of (2580/1600)$^2$ ≈ 2.6 higher for our derived dust mass.

2. The large grain size distribution we adopt versus the small grain distribution used by Harries et al. (2004) implies a closer/shorter distance of the dust component from the central binary—and thus a higher dust production rate (see Equation (6)).

When applying a small grain dust model and adopting a distance of 1600 pc we arrive at a consistent dust production rate of $\dot{M}_d \approx 6 \times 10^{-7} M_\odot \text{yr}^{-1}$.

The dust production rate from our large grain SED model shows that WR104 is producing the most dust out of the all the dustars in our sample and is likely the highest dust-maker among all known WC dustars (Table 4). Assuming a WC mass-loss rate of $3 \times 10^{-5} M_\odot \text{yr}^{-1}$ (Crowther 1997; Harries et al. 2004)

and a 40% carbon mass fraction, we infer a carbon dust condensation fraction of $\chi_C \approx 40\%$. This carbon dust condensation efficiency is notably higher than the ~few percent derived by Harries et al. (2004) and implies that a substantial fraction of the available carbon in the WC wind forms dust.

### 3.4.3. The Longest-period (32.5 yr) Dustar WR48a

WR48a is a continuous but variable dust-maker hosting a WC8 star with a proposed O-star companion. WR48a is the longest-period dustar known, with an IR light curve-derived period of 32.5 yr, and also exhibits episodes of brief dust formation (Williams et al. 2012). The total luminosity of the system is $L_* \sim 4 \times 10^5 L_\odot$, and its Gaia-derived distance of $d = 2270_{-570}^{+920} \text{ pc}$ is slightly less than previous distance estimates based on its suggested association with galactic
<figure><img src="image_7.png" /><figcaption>**Figure 4.** Top: Mid-IR 12.3 µm image of WR48a taken by Gemini/TReCS on 2004 March 15, overlaid with a dashed arc showing the location of a prominent dust arc plume. Bottom: 13.04 µm image of WR48a taken by VLT/VISIR on 2016 June 15 overlaid with the 2004 dust arc position (dashed) as above and its new position (dotted–dashed) at the date of the VISIR observation.</figcaption></figure>

clusters Danks 1 and 2 (~3.3–4.6 kpc; Danks et al. 1984). However, van der Hucht (2001) determine a distance of 1200 pc based on its WC8 spectral type. We therefore adopt the Gaia distance for our dust models, because it falls between the two distance estimates. We also adopt a line-of-sight interstellar extinction toward WR48a of $A_v = 8.28$ (Rate & Crowther 2020).

In Figure 4, we present mid-IR imaging observations taken by Gemini-South/TReCS (Marchenko & Moffat 2007; PID—GS-2004A-Q-63, PI—A. Moffat) and VLT/VISIR at similar N-band wavelengths at 2004 March 15 and 2016 June 15, respectively. The overlaid arc segment in Figure 4 shows the ~1.1' proper motion of eastern dust arc over the 12.3 years separating the two observations. Given the low signal-to-noise ratio of this dust arc in the VISIR images, as well as its extended nature, we assume position uncertainty of 0.023, which is $0.5 \times$ the PSF FWHM (0.45). Adopting the Gaia distance of 2270 pc, the dust expansion velocity of WR48a is $v_{exp} = 1000_{-240}^{+390} \pm 200$ km s⁻¹, where the first and second uncertainties are from the distance and proper motion, respectively. The derived expansion velocity is consistent with the $1200 \pm 170$ km s⁻¹ wind velocity measured for WR48a from its 10830 Å emission line width (Williams et al. 2012).

We apply the two-component dust SED model to fit the ISO/SWS spectrum of WR48a between 2 and 27.5 µm

<figure><img src="image_8.png" /><figcaption>**Figure 5.** Best-fit two-component DustEM model of WR48a using the small grain size distribution fit to ISO/SWS spectroscopy taken in 1996 February.</figcaption></figure>

(Figure 5). From the two-component model, we derive separation distances between the central system and dust components of $r_1 = 350_{-130}^{+100}$ au (110$_{-40}^{+70}$ au) and $r_2 = 2370_{-900}^{+640}$ au (720$_{-150}^{+440}$ au) for the small (large) dust distribution models. Historic mid-IR light curves show an IR brightening consistent with enhanced dust formation around 1994.5 (Williams et al. 2012). Using the estimated dust expansion velocity of 1000 km s⁻¹, this dust formed in 1994.5 would be located at ~340 au from the central binary at time of the ISO/SWS observation (1996 February). This distance estimate is consistent with the small grain distribution, which we therefore infer for the dust composition.

Based on our small grain distribution SED model (Figure 5), we derive dust masses for components 1 and 2 of $M_{d1} = (1.42_{-0.99}^{+1.12}) \times 10^{-7} M_\odot$ and $M_{d2} = (9.82_{-5.93}^{+5.84}) \times 10^{-6} M_\odot$. Using the component 1 distance and mass and a dust expansion velocity of $v_{exp} = 1000$ km s⁻¹, we estimate a dust production rate of $\dot{M}_d = (8.46_{-4.38}^{+3.48}) \times 10^{-8} M_\odot$ yr⁻¹. The carbon dust condensation fraction is therefore $\chi_c \approx 0.1\%$, assuming the radio-derived mass-loss rate of $1.7 \times 10^{-4} M_\odot$ yr⁻¹ (Zhekov et al. 2014) and a 40% carbon mass composition.

### 3.4.4. IR-luminous Dustars

WR98a. WR98a is a WC8 or WC9+OB dustar exhibiting continuous dust formation and appears like a “rotating pinwheel” that “revolves” with a 1.54 period (Monnier et al. 1999). The Gaia distance toward WR98a derived by Rate & Crowther (2020) is $d = 1260_{-410}^{+870}$ pc. Despite its close distance, WR98a is buried behind over 10 mag of visual extinction from the ISM along the line of sight based on a deep silicate absorption feature (Chiar & Tielens 2006). Rate & Crowther (2020) derive a much lower visual extinction of $A_v = 1.25$, but assign it a “b” flag that denotes it as highly implausible. Additionally, they note that the Gaia parallax measurement of WR98a exhibits large astrometric noise (>1 mas).

Due to the uncertainties in the Rate & Crowther (2020) interstellar extinction and distance measurements, we adopt $A_v = 13.79$ from van der Hucht (2001) and a distance of $d = 1900$ pc derived by Monnier et al. (1999) that is consistent with the upper distance uncertainties from Gaia. Monnier et al. (1999) derive this distance to WR98a based on its dust proper motion and adopting an expansion velocity of 900 km s⁻¹
<figure><img src="image_9.png" /><figcaption>**Figure 6.** Best-fit DustEM models of WR98a, WR118, WR19, WR70, WR137, and WR125. A large grain size distribution two-component dust model was used for the WR98a and WR118, and a small grain size distribution one-component dust model was used for the other four dustars.</figcaption></figure>

(Williams et al. 1995). The stellar luminosity of WR98a is not well-characterized given the high extinction. We therefore adopt a system luminosity of $L_* = 1.5 \times 10^5 L_\odot$ that is consistent with the value adopted in the WR98a models by Hendrix et al. (2016). Given its lower wind velocity, WC9 stellar properties are assumed for WR98a.

We fit double-component dust SED models to the ISO/SWS spectrum of WR98a between 2 and 27.5 µm (Figure 6, top left) and derive dust component distances of $r_1 = 200_{-70}^{+50}$ au ($60_{-10}^{+10}$ au) and $r_2 = 860_{-400}^{+1120}$ au ($370_{-130}^{+190}$ au) for the small (large) grain size distribution models. Spatially resolved images of the circumstellar dust around WR98a demonstrate that the majority of the mid-IR dust emission is concentrated in the central ~70 mas, or ~100 au, of the central binary (Monnier et al. 1999, 2007), with no prominent extended emission beyond ~150 mas, or ~200 au. As with WR104, we prefer the large grain

size distribution for WR98a. Notably, based on the dust component distances derived from the large grain SED model and adopting a dust expansion velocity of 900 km s⁻¹ (Williams et al. 1995), the expansion timescale between the two components is ~1.6 yr, which is consistent with the 1.54 yr orbital period. The separation between the two dust components is therefore consistent with the distance between the dust spirals in WR98a.

From the large grain dust SED model, we derive dust masses for components 1 and 2 of $M_{d1} = (1.82_{-0.68}^{+1.08}) \times 10^{-7} M_\odot$ and $M_{d2} = (9.25_{-5.29}^{+11.73}) \times 10^{-6} M_\odot$. Using the distance and mass of component 1, the dust production rate is $\dot{M}_d = (6.10_{-1.38}^{+1.77}) \times 10^{-7} M_\odot$ yr⁻¹. Adopting the mean WC9 star mass-loss rate of $\dot{M} = 10^{-4.66} M_\odot$ yr⁻¹ and a carbon mass fraction of 40% (Sander et al. 2019) for WR98a, we estimate a carbon dust condensation fraction of $\chi_C \approx 7\%$.
We note that there is a discrepancy between our derived dust production rate and the value determined from the hydrodynamic simulations by Hendrix et al. (2016), who obtain a dust production rate of ~$10^{-8} M_{\odot} yr^{-1}$. This discrepancy is primarily due to their assumption of a fixed value of 0.2% for the dust condensation fraction of the total WC mass-loss rate.

WR118. WR118 is a continuous WC9 dust-maker with a possible pinwheel morphology (Millour et al. 2009). The Gaia distance toward WR118 determined by Rate & Crowther (2020) is $d = 2490_{-680}^{+780}$ pc, which exhibited astrometric noise in excess >1 mas and large parallax uncertainty. However, this distance is consistent with previous estimates based on photometry and spectral type ($d=3130$ pc; van der Hucht 2001). We therefore adopt the Gaia distance for the dust SED models. WR118 is also heavily reddened by interstellar extinction, where $A_{\nu} = 13.68$ (Rate & Crowther 2020). We adopt a combined stellar luminosity of $L_{*} = 2.5 \times 10^{5} L_{\odot}$, consistent with the mean WC9 luminosity derived by Sander et al. (2019).

We fit double-component dust SED models to the ISO/SWS spectrum of WR118 between 2 and 27.5 µm (Figure 6, top right) and derive separation distances for components 1 and 2 of $r_1 = 120_{-50}^{+30}$ au ($30_{-10}^{+10}$ au) and $r_2 = 680_{-300}^{+320}$ au ($260_{-80}^{+120}$ au) for the small (large) grain size distributions. High spatial resolution interferometric IR observations of WR118 indicate that dust is located with ~15 mas (~40 au) of the central binary (Yudin et al. 2001; Monnier et al. 2007; Millour et al. 2009), which is comparable to the size of component 1 in the large grain model. We therefore favor the large grain model for WR118, which is also consistent with the large grain composition trend for the other IR-luminous WC9 dustars in our sample.

From the large grain SED models of WR118, we determine dust masses of $M_{d1} = (5.58_{-2.95}^{+6.24}) \times 10^{-8} M_{\odot}$ and $M_{d2} = (6.61_{-3.42}^{+6.99}) \times 10^{-6} M_{\odot}$ for components 1 and 2. Adopting an expansion velocity consistent with WC9 stars of $v_{\text{exp}} = 1390$ km s$^{-1}$, a WC9 mass-loss rate of $\dot{M} = 10^{-4.66} M_{\odot} yr^{-1}$ (Sander et al. 2019), and a 40% carbon fraction by mass, we find that the dust production rate and carbon dust condensation fraction are $\dot{M}_d = (6.27_{-1.94}^{+2.78}) \times 10^{-7} M_{\odot} yr^{-1}$ and $\chi_C \approx 7\%$, respectively.

### 3.4.5. Selected Single-component Dustars

WR19. WR19 is a periodic WC5+O9 dustar with an orbital period of 10.1 yr and an eccentricity of $e = 0.8$ (Williams et al. 2009b). This system hosts one of the earliest-type WC stars among the known WC dustars. The Gaia distance to WR19 is $d = 4330_{-580}^{+780}$ pc (Rate & Crowther 2020), which is slightly larger than previous distance estimates based on optical photometry and spectral type ($d \sim 3300$ pc; van der Hucht 2001). However, because Gaia observations of WR19 were taken before the onset of dust formation in 2017, it exhibited a mostly dust-free environment that likely contributed to the well-constrained distance. We therefore adopt the Gaia-derived distance for our dust SED modeling. WR19 is also affected by line-of-sight interstellar extinction, where $A_{\nu} = 5.19$ (Rate & Crowther 2020). We adopt a combined stellar luminosity of $L_{*} = 4 \times 10^{5} L_{\odot}$, the mean luminosity for WC5 stars as derived by Sander et al. (2019).

We fit single-component dust SED models to the stellar wind-subtracted 5–37 µm Spitzer/IRS spectrum of WR19 (Figure 6, center left) and derive a dust separation distance from the central binary of $r = 6870_{-3630}^{+3130}$ au ($2680_{-1530}^{+2030}$ au) for the small (large) grain size models. Based on the orbital phase of WR19 when Spitzer observed it (2005 June, $\phi \approx 0.84$), we can estimate the extent of dust formed at periastron passage in early 1997 (Williams et al. 2009b). Adopting the mean WC5 wind velocity from Sander et al. (2019) of $v_{\text{exp}} = 2780$ km s$^{-1}$ for the dust expansion velocity, the dust formed in the 1997 periastron passage should be located ~5000 au away in 2005 June. We therefore favor the small grain composition given the consistency with the dust expansion distance estimate.

For the small grain SED models of WR19, we derive a dust mass of $M_d = (3.98_{-2.84}^{+3.76}) \times 10^{-8} M_{\odot}$, which implies a dust production rate of $\dot{M}_d = (3.94_{-2.81}^{+3.72}) \times 10^{-9} M_{\odot} yr^{-1}$ (Equation (7)) given its orbital period of 10.1 yr. Adopting the mean WC5 mass-loss rate of $\dot{M} = 10^{-4.39} M_{\odot} yr^{-1}$ (Sander et al. 2019) and a 40% carbon mass fraction, the carbon dust condensation fraction is $\chi_C \approx 0.02\%$.

WR70. WR70 is a continuous but variable WC9+B0I dustar with a proposed period of 2.8 yr derived from its IR light curve (Niemela 1995; Williams et al. 2013b). The Gaia-derived distance to WR70 is $d = 3010_{-340}^{+440}$ pc, which is consistent with the recently revised extinction-based distance estimate of $d \approx 3500$ pc (Williams et al. 2013b). We therefore adopt the Gaia distance for our dust SED models. The adopted line-of-sight extinction and combined stellar luminosity of the WR70 system are $A_{\nu} = 5.20$ and $7 \times 10^5 L_{\odot}$, respectively (Williams et al. 2013b; Rate & Crowther 2020).

We fit single-component dust SED models to the stellar wind-subtracted 2–27.5 µm ISO/SWS spectrum of WR70 (Figure 6, center right) and determine a separation distance between dust and the central binary of $r = 550_{-70}^{+80}$ au ($150_{-50}^{+80}$ au) for the small (large) grain size models. Extended emission around WR70 has not yet been resolved, so it is difficult to distinguish between a small or large grain size distribution. However, we choose to adopt the small grain size distribution in order to be consistent with the WC dust formation analysis by Zubko (1998).

The small grain SED models of WR70 provide a dust mass of $M_d = (6.60_{-1.78}^{+2.40}) \times 10^{-8} M_{\odot}$. Adopting the mean WC9 mass-loss rate and velocity of $\dot{M} = 10^{-4.66} M_{\odot} yr^{-1}$ and $v_{\text{exp}} = 1390$ km s$^{-1}$ (Sander et al. 2019), we derive a dust production rate and carbon dust condensation fraction of $\dot{M}_d = (3.51_{-0.56}^{+0.65}) \times 10^{-8} M_{\odot} yr^{-1}$ and $\chi_C \approx 0.4\%$, where we have assumed the winds are composed of 40% carbon by mass.

WR137. WR137 is a periodic WC7+O9 dustar with an orbital period of 13.05 yr derived from both radial velocity measurements and IR light curve variations (Williams et al. 2001; Lefèvre et al. 2005). The Gaia-derived distance of $d = 2100_{-160}^{+180}$ pc is consistent with previous distance estimates toward the system (e.g., $d \approx 1820$ pc; Nugis & Lamers 2000). We therefore adopt the Gaia distances for our dust SED models of WR137. We also adopt an interstellar extinction and a heating source luminosity of $A_{\nu} = 1.70$ (Rate & Crowther 2020) and $L_{*} = 5.4 \times 10^{5} L_{\odot}$, respectively, where the stellar luminosity was derived by a pseudo-spectral fit from the PoWR models by Sander et al. (2012).
<figure><img src="image_10.png" /><figcaption>**Figure 7.** Best-fit DustEM models of WR95 and WR106 with the large grain size distribution.</figcaption></figure>

We fit single-component dust SED models to the stellar wind-subtracted contemporaneous 1.2–20 µm photometry taken in 1985 May by Williams et al. (2001) (Figure 6, bottom left) and derive a separation distance between dust and the central binary of $r = 660_{-220}^{+240}$ au (210$_{-70}^{+110}$ au) for small (large) grain size models. The mid-IR emission from dust formation in WR137 peaked in mid-1984 (Williams et al. 2001); therefore, assuming the expansion velocity of 2000 km s$^{-1}$ (Sander et al. 2012), we expect the dust to be located ~400 au at the time of the WR137 observations. Since dust formation from WR137 started several years before the mid-1984 peak and may therefore have propagated further, we suggest that the dust is most likely composed of small grains. We therefore favor the small grain model. This is also supported by WR137’s similarity in both orbital properties and spectral subtype to WR140, which exhibits dust consistent with small grains.

Based on the small grain SED models and the 13.05 yr dust formation/orbital period of WR137, we derive a dust mass of $M_d = (1.20_{-0.57}^{+0.64}) \times 10^{-8} M_\odot$ and a dust production rate of $\dot{M}_d = (9.21_{-4.36}^{+4.90}) \times 10^{-10} M_\odot \text{ yr}^{-1}$ (Equation (7)). Adopting the mean WC7 mass loss of $\dot{M} = 2.7 \times 10^{-5} M_\odot \text{ yr}^{-1}$ and a 40% carbon composition by mass, the carbon dust condensation fraction is $\chi_C = 0.009\%$.

WR125. WR125 is a WC7+O9III system (Williams et al. 1994; Midooka et al. 2019) that recently exhibited a second observed IR rebrightening, implying a 28.3 yr period (Williams 2019). The Gaia-derived distance toward WR125 is $3360_{-650}^{+990}$ pc (Rate & Crowther 2020), which is closer than the previous distance estimate of 4700 pc (Williams et al. 1992) based on a near-IR flux comparison to the similar system WR140 and its well-defined distance. Given the distance uncertainties arising from a relative comparison to WR140, we use the Gaia distance for our dust SED models of WR125. The adopted line-of-sight extinction and combined stellar luminosity of WR125 are $A_v = 6.48$ (Rate & Crowther 2020) and $L_* = 1.6 \times 10^5 L_\odot$, respectively, where the stellar luminosity was derived by a pseudo-spectral fit from the PoWR models by Sander et al. (2012).

We fit single-component dust SED models to the stellar wind-subtracted 1.2–20 µm photometry of WR125 taken in 1993.47–1993.5 (Figure 6, bottom right), several years after the onset of observed dust formation between 1990 and 91 (Williams et al. 1994), and determine a dust separation distance between the central binary of $r = 490_{-100}^{+170}$ au (170$_{-40}^{+70}$ au) assuming small (large) grains. Assuming a dust expansion

<figure><img src="image_11.png" /></figure>

velocity of 2000 km s$^{-1}$, consistent with the mean WC7 wind velocity derived by Sander et al. (2019), the dynamical timescale of the modeled dust is 1.2 yr (0.4 yr) for the small (large) dust grains distributions. Because the small grain model timescale is consistent with formation during the observed peak IR emission in mid-1992 (Williams et al. 1994), we favor the small grain dust model.

For the small grain SED models, we derive a dust mass of $M_d = (1.00_{-0.35}^{+0.61}) \times 10^{-7} M_\odot$ and a dust production rate of $\dot{M}_d = (3.54_{-1.22}^{+2.15}) \times 10^{-9} M_\odot \text{ yr}^{-1}$ (Equation (7)), assuming an orbital period of 28.3 yr. Our derived dust mass is roughly consistent with the $M_d = 1.64 \times 10^{-7} M_\odot$ derived from the isothermal dust model by Williams et al. (1994) after correcting for different distance assumptions. We adopt a mass-loss rate of $1.6 \times 10^5 M_\odot \text{ yr}^{-1}$ from the spectral model fit by Sander et al. (2012), and a 40% carbon mass fraction. We then determine a carbon dust condensation fraction of $\chi_C \approx 0.03\%$.

### 3.4.6. Additional WC Dustars with Angular Size Constraints

WR95 and WR106 are continuous WC9 dustars (van der Hucht 2001) that have constraints on their extended dust emission from IR high spatial resolution observations (Monnier et al. 2007; Rajagopal et al. 2007). The Gaia-derived distances toward WR95 and WR106 are $2070_{-310}^{+430}$ pc and $3070_{-430}^{+560}$ pc, respectively (Rate & Crowther 2020). The dust separation distances determined from our SED models are $r = 130_{-40}^{+110}$ au ($40_{-20}^{+70}$ au) for WR95 and $r = 170_{-60}^{+220}$ au ($50_{-20}^{+120}$ au) for WR106, assuming small (large) grains.

Rajagopal et al. (2007) measure a 28.4 and a 45.0 mas Gaussian FWHM at 10.5 µm for WR95 and WR106, respectively, consistent with the FWHM measured by Monnier et al. (2007) for these systems at shorter IR wavelengths (2.2 and 3.08 µm) near the emission peaks in the SEDs (Figure 7). These angular sizes correspond to a linear size of 60 au for WR95 and 100 au for WR106. Assuming a circular morphology, the extent of the dust shells should be approximately half of these linear sizes: ~30 au for WR95 and ~50 au for WR106. This is consistent with our large grain size models, and agrees with the interpretation from Rajagopal et al. (2007). Those authors also favor dust models with larger (~1 µm) grains for WR95 and WR106, and determine dust shell radii of a few tens of astronomical unit from the stars.

Assuming the large grain SED models, the dust production rates are $\dot{M}_d = (9.36_{-3.36}^{+2.00}) \times 10^{-8} M_\odot \text{ yr}^{-1}$ and $\dot{M}_d = (2.97_{-1.01}^{+1.14}) \times 10^{-7} M_\odot \text{ yr}^{-1}$ for WR95 and WR106, respectively. Given the
<figure><img src="image_12.png" /><figcaption>**Figure 8.** Best-fit DustEM models of WR48b, WR53, WR59, WR96, WR103, and WR119.</figcaption></figure>

adopted mass-loss rates (Table 3) from PoWR models by Sander et al. (2012) and assuming a 40% carbon composition in the winds by mass, the carbon dust condensation fraction is ~1.2% for WR95 and ~4.6% for WR106.

For the dustars without dust radius constraints (e.g., WR70), we assume a small grain size distribution. In addition to WR70, the eight other dustars that fall into this category are WR48b, WR53, WR59, WR96, WR103, WR119, WR121, and WR124-22. DustEM SED models of these eight dustars are shown in Figures 8 and 9.

### 3.4.7. Dust Properties and Comparison to Previous Studies

The results of the DustEM SED models to the 19 WC dustars are presented in Table 4. In Figure 10, the dust condensation fraction $\chi_C$ and the DPR are plotted against the total IR luminosity of each dustar. Both plots show that $\chi_C$ and the DPR

increase with increasing with IR luminosity. The periodic or highly variable dust producers (WR19, WR48a, WR125, WR137, and WR 140) show slight discrepancies from this trend because their IR luminosity varies depending on the orbital phase. These results demonstrate that WC dustars exhibit a range of DPRs and $\chi_C$ values ranging from $1 \times 10^{-10} - 4 \times 10^{-6} M_{\odot} \text{yr}^{-1}$ and ~0.002%–40%, respectively.

The most recent and relevant study to compare our results to is that of Zubko (1998), who conducted a homogeneous modeling analysis of the SEDs and dust shells of 17 Galactic WC dustars incorporating grain physics and dynamics. A majority of their sample overlaps with ours. Zubko (1998), however, modeled the dust formation physics and dynamics under the assumption of spherically symmetric winds, because the complex, colliding-wind morphology of WC dustars was not known at that time. Despite the difference in analytical
<figure><img src="image_13.png" /></figure>

<figure><img src="image_14.png" /><figcaption>Figure 9. Best-fit DustEM models of WR121 and WR124-22.</figcaption></figure>

<figure><img src="image_15.png" /></figure>

<figure><img src="image_16.png" /><figcaption>Figure 10. Left: Distribution of derived dustar values for $\chi_C$ vs. $L_{IR}$. Markers with thick borders represent dustars where the large grain model was favored, markers with thin borders represent dustars where the small grain model was favored, and markers without borders show dustars where the small grain model was adopted by default. Open, unfilled markers show the highly variable or periodic dustars WR19, WR48a, WR125, WR137, and WR140. Right: Distribution of DPR vs. $L_{IR}$ overlaid with the best-fit power-law relation DPR $\propto L_{IR}^{1.21\pm0.06}$. Highly variable or periodic dustars (WR19, WR48a, WR125, WR137, and WR140) were excluded from the fit.</figcaption></figure>

techniques, Zubko (1998) derive a range of mass-loss rates and carbon dust fractions consistent with our results. They determine condensed carbon fractions of 0.002%–20% and dust production rates of $5 \times 10^{-10}-6 \times 10^{-6} M_{\odot} yr^{-1}$.

### 3.4.8. IR Luminosity, Dust Production Rate, and Grain Size

The distribution of $L_{IR}$ and the DPR of WC dustars provides a valuable empirical relation for estimating DPRs from IR luminosities that can be directly measured from IR observations. Excluding the periodic and highly variable dustars, the DPR versus $L_{IR}$ (Figure 10, right) relation can be characterized by the following power-law fit:

$$ \text{Log}\left(\frac{\dot{M}_d}{M_{\odot} \text{yr}^{-1}}\right) = 1.21_{-0.06}^{+0.06} \text{Log}\left(\frac{L_{IR}}{L_{\odot}}\right) - 11.55_{-0.18}^{+0.32} \quad (8) $$

In Figure 10 (left), it is apparent that the efficient dust- makers, where $\chi_C \gtrsim 1\%$, are consistent with large grain dust models. This suggests that there is no homogeneous grain size distribution exhibited by all WC dustars, but rather that larger grain sizes ($\gtrsim 0.1 \mu\text{m}$) are produced by dustars with higher dust condensation fractions. However, if the WC subtype influences the grain size, it is possible that WC9 dustars preferentially

form large grains. Regardless, the results suggest that efficient WC9 dustars form larger dust grains that will be more robust against destruction and sputtering via SN shocks and survive longer in the ISM.

### 3.5. Investigating the Colliding-wind Dust Formation Model

Usov (1991) investigated a theoretical model of dust formation in the wind collision between WR and OB binaries. They derived an analytical relation for the fraction of the WR wind that is strongly cooled and compressed in the shock layer between the two stars:

$$ \alpha \propto \eta^2 (1 + \eta^{1/2})^2 \dot{M}_{\text{WR}}^2 v_{\text{WR}}^{-6} D^{-2}, \quad (9) $$

where $\eta$ is the wind momentum ratio between the OB and WR stars ($\eta = \dot{M}_{\text{OB}}/v_{\text{OB}}$), $\dot{M}_{\text{WR/OB}}$ is the mass-loss rate of the WR/ OB star, $v_{\text{WR/OB}}$ is the terminal wind velocity of the WR/OB star, and $D$ is the separation distance between the WR and OB star. Usov (1991) assume that most of the carbon from the WC wind in this cold, compressed shock layer condenses into dust,
<figure><img src="image_17.png" /><figcaption>**Figure 11.** Graph of $\chi_C \dot{M}_{WR}^{-2} v_{WR}^6$ normalized to the $\chi_C$, velocity, mass-loss rate of WR104 plotted against $P_{orb}$ for the eight dustars in our sample with orbital periods. Dashed line shows the best-fit power-law relation $\chi_C \dot{M}_{WR}^{-2} v_{WR}^6 \propto P_{orb}^{-1.2 \pm 0.4}$, where the 32.5 yr orbital period dustar WR48a is the outlier.</figcaption></figure>

which implies that

$$ \alpha \propto \chi_C. \qquad (10) $$

We can therefore test the theoretical model from Usov (1991) by investigating the relation between $\chi_C$ and the WC binary mass-loss rates, wind velocities, and orbital separation.

We focus on the eight WC dustars in our sample with known orbital periods ($P_{orb}$): WR19, WR48a, WR70, WR98a, WR104, WR125, WR137, and WR140. Using Kepler's Third Law and assuming a similar total mass for each dustar system, it follows that the orbital separation $D$ is proportional to $P_{orb}^{2/3}$. Eccentricity is also an important parameter that regulates the dust formation and orbital separation for the periodic, noncontinuous dustars (e.g., WR19, WR125, WR137, and WR140); however, their dust production rates and thus their carbon condensation efficiencies were calculated as an average over the orbital period (Equation (7)). Since the time-averaged separation distance is $a(1 + e^2/2)$ (see Williams 2003), where $a$ is the semimajor axis and $e$ is the eccentricity of the system, the eccentricity should only contribute variations from $a$ by a factor of 1.5 at most. Given the minor effect of eccentricity in this relation and the eccentricity uncertainty in most dustars, we assume that $D \propto P_{orb}^{2/3}$ for all eight systems.

Due to uncertainties in the mass-loss rate and wind velocity of the binary OB companions of the WC stars, we make the assumption that $\eta$ is similar for the dustars. If the observed dustar properties are consistent with the Usov (1991) colliding-wind model (Equation (9)), they should exhibit the following relation with our assumptions:

$$ \chi_C \dot{M}_{WR}^{-2} v_{WR}^6 \propto P_{orb}^{-4/3}. \qquad (11) $$

In Figure 11, we plot $\chi_C \dot{M}_{WR}^{-2} v_{WR}^6$ versus $P_{orb}$ for the eight dustars normalized to the carbon dust condensation efficiency, wind velocity, and mass-loss rate of WR104. Values for $\dot{M}_{WR}$, $v_{WR}$, and $P_{orb}$ are provided from Table 3, and $\chi_C$ is provided from Table 4. We note that the power-law relation is insensitive to the dustar system used to provide the normalization factors. We determine a $P_{orb}^{-1.2 \pm 0.4}$ power-law fit, where WR48a is excluded as an outlier since it exhibits a higher mass-loss rate than the other systems (see Table 3). This derived power-law fit is consistent with the theoretical power-law relation from the

Usov (1991) model, with our assumption of constant $\eta$ across the dustars (Equation (11)). The consistency of the observational results and theoretical model supports the colliding-wind interpretation of WC dustars and demonstrates the effect of the binary orbital parameters on dust formation efficiency

# 4. BPASS Dust Production Model Results

## 4.1. Binary Population and Spectra Synthesis Models

The Binary Population and Spectral Synthesis (BPASS; Eldridge et al. 2017) suite of binary stellar evolution and stellar population models provides a valuable tool for investigating dust production from stellar populations. The presence of a binary companion is a key factor that influences both the formation of WC stars and their dust production via colliding winds, which underscores the utility of the binary evolution tracks utilized in BPASS for modeling the dust contribution from WC binaries in stellar populations. Notably, BPASS evolutionary models may also be used to study the formation history of individual WC binary system. This work presents the first use of BPASS to model dust production from stellar populations.

We use the BPASS v2.2.1 binary models as described in Stanway & Eldridge (2018), with the default “135_300” BPASS initial mass function (IMF). This IMF is based on Kroupa et al. (1993), with a power-law slope from 0.1 to 1.0 $M_{\odot}$ of $-1.30$ that steepens to $-2.35$ above this to an upper stellar mass of 300 $M_{\odot}$. Stellar mass loss is determined by the difference between the stellar mass in the current and previous time steps.

In addition to the WC binaries, three more sources of dust are considered in our BPASS models, based on the current understanding of leading dust producers: dust condensation in SN II ejecta and in stellar winds of AGB and RSG stars (Matsuura et al. 2009; Dwek & Cherchneff 2011; Boyer et al. 2012; Temim et al. 2015; Srinivasan et al. 2016). We do not include the dust input from H-poor type Ia or Ib/c SNe since type II SNe are considered to be the most efficient and prolific source of SN dust (e.g., Sarangi et al. 2018 and references therein).

Luminous Blue Variables (LBVs), massive evolved stars that exhibit giant outbursts of extreme mass loss (Humphreys & Davidson 1994; Smith 2014), may also be important sources of dust. Based on observations of their circumstellar dusty nebulae, LBVs are capable of forming between ~0.001 and 0.1 $M_{\odot}$ of dust (Kochanek 2011 and references therein). This corresponds to a DPR of ~$10^{-7}-10^{-5} M_{\odot}$ yr$^{-1}$ for an LBV, assuming a lifetime of ~$10^4$ yr (Smith 2014), which is comparable to the upper DPR range of WC dustars. However, given their rarity, the evolutionary formation pathways for LBVs and the physics of their eruptive mass loss are not well-understood. Due to these uncertainties, it is difficult to identify LBVs and quantify their dust contribution via eruptive mass loss in current BPASS models (Eldridge et al. 2017). We therefore we omit them from the BPASS models in this study, but acknowledge their potential as important dust producers that warrant further investigation.

The dust formation prescriptions and the assumed parameters for identifying the four dust inputs sources (type II SNe, AGBs, RSGs, and WC binaries) in the BPASS models are described as follows:
1. Type II SNe. For stars with a final mass greater than 2 $M_\odot$, CO core mass greater than 1.38 $M_\odot$, and H surface mass fraction greater than 1%, the terminal explosion is considered a Type II SN. A 10% dust condensation efficiency by mass is adopted for the dust-forming elements in the SN ejecta (Dwek et al. 2014; Sarangi et al. 2018). Dust-forming elements include carbon, oxygen, magnesium, silicon, and iron. The hydrogen, helium, nitrogen, and neon synthesized in the SN ejecta are not considered as dust-forming elements.

2. Type II SNe Dust Destruction. When an SN explodes, its shock interacts with the surrounding ISM and destroys the ISM dust grains. Given a total ISM mass that the SN clears of dust, $m_g$, the destroyed quantity of dust per SN explosion can be estimated as

$$m_d = m_g \chi_{d/g}, \qquad (12)$$

where $\chi_{d/g}$ is the gas-to-dust mass ratio of the ISM (Dwek et al. 2014; Temim et al. 2015). Note that this does not include destruction of SN ejecta dust by the reverse shock, which is not considered in our work. The following fitting formula derived by Yamasawa et al. (2011) is adopted for determining $m_g$:

$$m_g = 1535 n_{\text{ISM}}^{-0.202} ((Z/Z_{\odot}) + 0.039)^{-0.298} M_{\odot}, \qquad (13)$$

where $n_{\text{ISM}}$ is the density of the surrounding ISM and $Z_{\odot}$ is the solar metallicity ($Z = 0.02$). Given the weak dependence on $n_{\text{ISM}}$, an ISM density of $n_{\text{ISM}} = 1 \text{ g cm}^{-3}$ is assumed for all BPASS models. Total SN dust input rates are determined with and without incorporating this dust destruction prescription.

3. AGBs. An AGB star is identified as a star with an effective temperature less than $10^{3.66}$ K, a luminosity less than $10^{4.4} L_\odot$, a CO core mass greater than 0.5 $M_\odot$, and a difference between the He and CO core masses less than 0.1 $M_\odot$. This luminosity cut off is consistent with the most luminous AGB in the Large Magellanic Cloud (Riebel et al. 2012). The AGB dust production is calculated assuming a 0.5% condensation fraction of the total mass-loss and is scaled by the metallicity ratio of the model and the solar metallicity, $Z/Z_\odot$ (van Loon 2000).

4. RSGs. A star is identified as an RSG if its effective temperature is less than $10^{3.66}$ K and its luminosity is greater than $10^{4.4} L_\odot$. The empirically derived $M_{\text{bol}}$-DPR relation from Massey et al. (2005) is adopted and scaled by the metallicity ratio of the model and the solar metallicity, $Z/Z_\odot$, to quantify the dust input from the RSGs. An upper limit was set on the RSG dust production requiring the dust input of each time step to be less than the total metals in the mass loss.

5. WC Binaries. WC+OB binaries are identified as systems where the WC star has no hydrogen, a combined C and O mass fraction higher than N, a He mass fraction higher than 40%, a C mass fraction higher than 25%, and a companion with an effective temperature greater than $10^{4.38}$ K. The metallicity dependence for the total mass-loss rates from WR stars in BPASS are described as a power-law relation $\dot{M} \propto Z^\alpha$, where $\alpha = 0.5$ (Eldridge et al. 2017). The dust production rate from WC binaries is $\dot{M}_d = \dot{M} X_C \chi_C$, where $X_C$ is the surface mass fraction of carbon that is output from BPASS. The $\chi_C$ relation

(Equation (11)) that we verified in Section 3.5 is adopted and scaled to the properties of WR104:

$$\chi_C \approx 40\% \left( \frac{\dot{M}}{\dot{M}_{\text{WR104}}} \right)^2 \left( \frac{v_\infty}{v_{\text{WR104}}} \right)^{-6} \left( \frac{P_{\text{orb}}}{P_{\text{WR104}}} \right)^{-4/3} . \qquad (14)$$

The dust condensation fraction is therefore treated as a function of the mass-loss rate and terminal wind velocity of the primary star (i.e., the WC star) and the orbital period of the system. The terminal wind velocity ($v_\infty$) is defined by the stellar mass-loss rate ($\dot{M}$), luminosity ($L_*$), and the momentum transfer efficiency ($\eta_{\text{mom}}$) via the relation $\eta_{\text{mom}} = \dot{M}v_\infty/(L_*/c)$, where c is the speed of light. Here, $\dot{M}$ and $L_*$ are output parameters from BPASS, and we adopt a fixed value of $\eta_{\text{mom}} = 5.3$, the median value of the 11 Galactic dusty WC stars from the optical spectroscopic analysis by Sander et al. (2019). In the BPASS models, we only consider the dust input from WC binaries with orbital periods consistent with the range of orbital periods in our sample: $P_{\text{orb}} = 0.66-32.5$ yr. Finally, since not all WC binaries are dust producers, we conservatively assume that 28% of the BPASS WC systems that meet these criteria are dustars, which is based on the observed WC-dustar/WC ratio in the Galaxy (Rosslowe & Crowther 2015). This is percentage is probably an underestimate, given the recent mid-IR variability analysis by Williams (2019).

Since the shocks from end-of-life SN explosions may destroy the newly formed circumstellar dust from RSGs and WC binaries, their dust input is not included if they end their lives as SNe. It is assumed that RSGs and WC stars end their lives as explosive SNe if their remnant core mass is less than 3 $M_\odot$, whereas a core with a mass greater than 3 $M_\odot$ leads to a direct collapse into a black hole with minimal ejecta energy and circumstellar dust survival.

## 4.2. BPASS Dust Models at $Z = 0.001, 0.008, \text{ and } 0.020$

BPASS dust models were run at low ($Z = 0.001$), LMC-like ($Z = 0.008$), and solar ($Z = 0.020$) metallicities in order to explore the relative dust contribution of WC binaries, AGBs, RSGs, and SNe at different epochs in cosmic time as characterized by different metallicities. The star formation history will impact the stellar populations and their dust input. For example, in a scenario with constant star formation (SF), the dust sources associated with massive stars have short lifetimes but will be continuously replenished by ongoing star formation. Two star formation scenarios were therefore considered for the BPASS models at each metallicity: a coeval $10^6 M_\odot$ stellar population and constant SF. These two different star formation scenarios represent contrasting star-forming histories, i.e., a single starburst event versus continuous star formation.

The dust-to-gas mass ratios, which are relevant for the SN destruction rates, are assumed to be $\chi_{d/g} = 0.0008, 0.0018, \text{ and } 0.007$ (Zubko et al. 2004) for the $Z = 0.001, 0.008, \text{ and } 0.020$ metallicity models, respectively. The low-Z and LMC-like dust-to-gas mass ratios are adopted from the values derived for the SMC and LMC, respectively, by Temim et al. (2015). Based on Equation (13), the total ISM masses swept up by SNe in the low-Z, LMC-like, and solar metallicity models are $m_g = 3088, 1947, \text{ and } 1518 M_\odot$, respectively.
<figure><img src="image_18.png" /><figcaption>**Figure 12.** Left column: Constant SF BPASS models at low (Z = 0.001), LMC-like (Z = 0.008), and solar (Z = 0.020) metallicities showing the cumulative dust mass input as a function of time. SNe models with dust destruction do not appear in any of the plots, because SNe are net destroyers of dust at all time steps. Right column: Coeval $10^6 M_{\odot}$ stellar population BPASS models at low, LMC-like, and solar metallicities showing the total cumulative dust mass input (solid lines) as a function of time. Dotted lines show the dust input in a single time bin for each dust source. SNe + dust destruction model is omitted for the same reason as described above.</figcaption></figure>

The cumulative dust production rates for the constant SF and coeval $10^6 M_{\odot}$ stellar population models are shown in Figure 12. Results from the constant SF and $10^6 M_{\odot}$ stellar population models at a time $t = 1.0$ Gyr are provided in Tables 5 and 6, respectively.

### 4.2.1. Constant SF Models

The results from the constant SF models are shown in Figure 12 (left column) and Table 5. In order to check the DPR prescriptions of the four dust sources, we compare the model results against observationally measured values and previous studies. The average DPRs for AGB and RSG stars in the LMC-like metallicity model are $\langle \dot{M}_{d,\text{AGB}} \rangle = 1.1 \times 10^{-9} M_{\odot} \text{ yr}^{-1}$

and $\langle \dot{M}_{d,\text{RSG}} \rangle = 1.7 \times 10^{-9} M_{\odot} \text{ yr}^{-1}$, which are consistent within factors of a few to the average DPRs measured from mid-IR observations of the LMC (Riebel et al. 2012; Srinivasan et al. 2016), where $\langle \dot{M}_{d,\text{AGB}} \rangle = 6.9 \times 10^{-10} M_{\odot} \text{ yr}^{-1}$ and $\langle \dot{M}_{d,\text{RSG}} \rangle = 4.0 \times 10^{-10} M_{\odot} \text{ yr}^{-1}$. Temim et al. (2015) calculate a theoretical IMF-averaged SN dust yield of 0.65 $M_{\odot}$ for a 100% dust condensation efficiency in SNe. A 10% condensation efficiency, which we adopt for the BPASS models, applied to the Temim et al. (2015) SN dust yields is therefore consistent within a factor of two to the BPASS SN dust yields ($\sim 0.1 M_{\odot}$). The average dust mass destroyed per SN in the LMC measured by Temim et al. (2015) is $6.1 \pm 2.6 M_{\odot}$ for silicates and $1.6 \pm 0.7 M_{\odot}$ for carbon grains, which is consistent within factors of a few to the average dust mass destroyed per SN in the
<table><caption>**Table 5**
BPASS Constant SF Dust Models (at $t = 1.0$ Gyr)</caption><thead><tr><th></th><th>Low-Z<br>Z = 0.001</th><th>LMC-like<br>Z = 0.008</th><th>Solar<br>Z = 0.020</th></tr></thead><tbody><tr><td>χ<sub>d/g</sub></td><td>0.0008</td><td>0.0018</td><td>0.007</td></tr><tr><td>m<sub>g</sub> (M<sub>☉</sub>)</td><td>3088</td><td>1947</td><td>1518</td></tr><tr><td>N<sub>AGB</sub></td><td>1739 (100%)</td><td>5465 (100%)</td><td>6613 (100%)</td></tr><tr><td>N<sub>RSG</sub></td><td>2332 (65%)</td><td>6070 (30%)</td><td>6900 (28%)</td></tr><tr><td>N<sub>WC</sub></td><td>3 (100%)</td><td>27 (100%)</td><td>47 (56%)</td></tr><tr><td>R<sub>SN</sub> (yr<sup>-1</sup>)</td><td>2.6 × 10<sup>-2</sup></td><td>1.7 × 10<sup>-2</sup></td><td>1.2 × 10<sup>-2</sup></td></tr><tr><td>M<sub>d,AGB</sub> (M<sub>☉</sub>)</td><td>4.4 × 10<sup>-7</sup></td><td>5.7 × 10<sup>-6</sup></td><td>1.9 × 10<sup>-5</sup></td></tr><tr><td>yr<sup>-1</sup>)</td><td>[2.5 × 10<sup>-10</sup>]</td><td>[1.1 × 10<sup>-9</sup>]</td><td>[2.9 × 10<sup>-9</sup>]</td></tr><tr><td>M<sub>d,RSG</sub> (M<sub>☉</sub>)</td><td>3.1 × 10<sup>-7</sup></td><td>3.0 × 10<sup>-6</sup></td><td>7.4 × 10<sup>-6</sup></td></tr><tr><td>yr<sup>-1</sup>)</td><td>[2.0 × 10<sup>-10</sup>]</td><td>[1.7 × 10<sup>-9</sup>]</td><td>[3.8 × 10<sup>-9</sup>]</td></tr><tr><td>M<sub>d,WC</sub> (M<sub>☉</sub>)</td><td>1.7 × 10<sup>-10</sup></td><td>1.7 × 10<sup>-6</sup></td><td>2.6 × 10<sup>-5</sup></td></tr><tr><td>yr<sup>-1</sup>)</td><td>[6.5 × 10<sup>-11</sup>]</td><td>[6.2 × 10<sup>-8</sup>]</td><td>[9.9 × 10<sup>-7</sup>]</td></tr><tr><td>M<sub>d,SN</sub> (M<sub>☉</sub>)</td><td>3.1 × 10<sup>-3</sup> [0.12</td><td>1.8 × 10<sup>-3</sup> [0.10</td><td>1.1 × 10<sup>-3</sup> [0.10</td></tr><tr><td>yr<sup>-1</sup>)</td><td>M<sub>☉</sub>]</td><td>M<sub>☉</sub>]</td><td>M<sub>☉</sub>]</td></tr><tr><td>M<sub>SN,dest</sub> (M<sub>☉</sub>)</td><td>-6.5 × 10<sup>-2</sup></td><td>-5.9 × 10<sup>-2</sup></td><td>-8.2 × 10<sup>-2</sup></td></tr><tr><td>yr<sup>-1</sup>)</td><td>[-2.5 M<sub>☉</sub>]</td><td>[-3.5 M<sub>☉</sub>]</td><td>[-10.6 M<sub>☉</sub>]</td></tr><tr><td>Log (t<sub>0,AGB</sub>)</td><td>7.9</td><td>7.9</td><td>7.8</td></tr><tr><td>Log (t<sub>0,RSG</sub>)</td><td>6.3</td><td>6.4</td><td>6.5</td></tr><tr><td>Log (t<sub>0,WC</sub>)</td><td>6.3</td><td>6.3</td><td>6.4</td></tr><tr><td>Log (t<sub>0,SN</sub>)</td><td>6.5</td><td>6.6</td><td>6.7</td></tr></tbody></table>

**Note.** Dust source populations from the constant SFR BPASS dust model for AGBs, RSGs, and WC binaries, as well as their SN rates ($R_{SN}$), and the total and average (in brackets) DPR for each source at time $t = 1.0$ Gyr. Here, $\chi_{d/g}$ is the adopted dust-to-gas mass ratio, and $m_g$ is the total ISM mass swept up by each SN as determined by Equation (13) at each metallicity. The value of $N$ shows the total numbers of sources present at this time, and the percentage indicates the fraction that contribute to the dust input and do not end their lives as SN. The AGB, RSG, and WC binary DPRs are determined from the population that do not end their lives as SN. The SN DPR value in brackets is the average mass of dust formed/destroyed per SN. The time, $t_0$ (yr), corresponds to the onset time of each dust source.

<table><caption>**Table 6**
BPASS $10^6 M_{\odot}$ Stellar Population Dust Models (at $t = 1.0$ Gyr)</caption><thead><tr><th></th><th>Low-Z<br>Z = 0.001</th><th>LMC-like<br>Z = 0.008</th><th>Solar<br>Z = 0.020</th></tr></thead><tbody><tr><td>M<sub>d,AGB</sub> (M<sub>☉</sub>)</td><td>1.0 × 10<sup>1</sup></td><td>1.0 × 10<sup>2</sup></td><td>2.7 × 10<sup>2</sup></td></tr><tr><td>M<sub>d,RSG</sub> (M<sub>☉</sub>)</td><td>2.9 × 10<sup>-1</sup></td><td>2.1 × 10<sup>0</sup></td><td>4.3 × 10<sup>0</sup></td></tr><tr><td>M<sub>d,WC</sub> (M<sub>☉</sub>)</td><td>1.9 × 10<sup>-5</sup></td><td>2.3 × 10<sup>-1</sup></td><td>3.7 × 10<sup>0</sup></td></tr><tr><td>M<sub>d,SN</sub> (M<sub>☉</sub>)</td><td>1.4 × 10<sup>3</sup></td><td>9.1 × 10<sup>2</sup></td><td>5.4 × 10<sup>2</sup></td></tr><tr><td>Δt<sub>AGB</sub> (yr)</td><td>1.5 × 10<sup>9</sup></td><td>1.9 × 10<sup>9</sup></td><td>2.4 × 10<sup>9</sup></td></tr><tr><td>Δt<sub>RSG</sub> (yr)</td><td>2.0 × 10<sup>8</sup></td><td>1.6 × 10<sup>8</sup></td><td>2.5 × 10<sup>8</sup></td></tr><tr><td>Δt<sub>WC</sub> (yr)</td><td>5.2 × 10<sup>5</sup></td><td>4.3 × 10<sup>6</sup></td><td>2.5 × 10<sup>6</sup></td></tr><tr><td>Δt<sub>SN</sub> (yr)</td><td>9.7 × 10<sup>7</sup></td><td>7.5 × 10<sup>7</sup></td><td>5.8 × 10<sup>7</sup></td></tr></tbody></table>

**Note.** Model results show the cumulative dust mass, $M_d$, produced by each source and the timescale, $\Delta t$, on which each source is actively forming dust. These results show the population of RSG and WC stars that do not end their lives as SNe

LMC-like BPASS models ($3.5 M_{\odot}$). We discuss the comparisons to the LMC in further detail in Section 4.3.

The average dust production rate from WC binaries in the solar metallicity models is $\langle \dot{M}_{d,WC} \rangle = 9.9 \times 10^{-7} M_{\odot} \text{yr}^{-1}$, which is slightly greater than the average DPR from our sample of 19 Galactic dustars $\langle \dot{M}_{d,WC} \rangle = 3.3 \times 10^{-7} M_{\odot} \text{yr}^{-1}$. However, we find that the WC binary dust production from the

BPASS models is largely dominated by systems with the shortest orbital period: ~50% of the total WC binary dust input in the solar metallicity model originates from three systems with $P_{\text{orb}}$ between 0.66 and 1.0 yr. These three systems have an average DPR of $3.7 \times 10^{-6} M_{\odot} \text{yr}^{-1}$, which is consistent with our measured DPR of the $P_{\text{orb}} = 0.66$ yr system WR104, $\dot{M}_d = 4.4 \times 10^{-6} M_{\odot} \text{yr}^{-1}$. We therefore conclude that our BPASS dust prescriptions provide reasonable estimates for the dust input from the four sources.

The constant SF BPASS models show that dust production from SNe without dust destruction dominates AGBs, RSGs, and WC binaries by over an order of magnitude for all three metallicities. However, SNe are consistently net dust destroyers at all time steps when incorporating SN dust destruction, and the SN dust destruction rate exceeds the production rate by over an order of magnitude. The average dust produced per SN is ~0.1 $M_{\odot}$ for all three metallicities. We note that dust formation in SN ejecta is still uncertain and measurements range from ~$10^{-6}$ to $1 M_{\odot}$ of dust per SN (Temim et al. 2015; Gall & Hjorth 2018; Sarangi et al. 2018), but SNe would have to form dust beyond the highest mass range of the observed values in order to compensate for the dust destruction rate.

In the solar metallicity models, WC binaries are the dominant source of dust for ~60 Myr until the onset of the AGB star population that forms dust at a similar rate to the WC binaries. We note that the WC binaries are significant dust producers in the solar metallicity model despite almost half of the WC stars ending their lives as SNe with their dust input removed (Table 5).

The LMC-like metallicity model shows slightly reduced dust input from AGB and RSG stars compared to the solar metallicity models. The WC binaries show a much greater reduction in dust production due to having a smaller population (by a factor of two) and a lower average dust production rate (by over an order of magnitude) versus those of the solar metallicity model. AGB stars lead the dust production in the LMC-like metallicity models, but only by factors of 2–3 over RSGs and WC binaries. WC binaries notably form within the first 2 Myr after the onset of star formation and are the first sources of dust in the LMC-like (and solar) models.

In the low-Z models, RSGs, AGBs, and WC binaries show smaller populations and form dust at a much lower rate than the LMC-like and solar models. WC binaries exhibit the most substantial decrease in dust production and population. This is because only the most massive and luminous stars evolve to the WC phase at lower metallicities. Such luminous stars drive faster winds and therefore form dust much less efficiently, given the sensitivity of $\chi_C$ to wind velocity (Equation (9)). At low Z, if SNe are indeed net dust destroyers, AGBs and RSGs are likely the dominant dust sources, with dust input rates exceeding WC binaries by over three orders of magnitude.

### 4.2.2. Coeval $10^6 M_{\odot}$ Stellar Population Models

Figure 12 (right column) shows the dust mass input in each time bin and the cumulative dust mass input as a function of time for WC binaries, AGBs, RSGs, and SNe. The most notable difference in the dust production between the coeval population and the constant SF models is the relatively higher dust contribution by AGBs. For example, the total dust mass from AGB stars at a time $t = 1.0$ Gyr is over an order of magnitude greater than the total dust mass input from RSGs
and WC binaries in all three metallicity models (Table 6). Like the constant SF models, SNe are consistently net dust destroyers at all three metallicities. The RSG dust input shows two peaks, where the secondary peak around ~100 Myr is consistent with the onset of AGBs. This secondary, late-time peak arises due to the difficulties in distinguishing the population of luminous AGB stars from faint RSGs (Eldridge et al. 2017). SN models with dust destruction show that the AGBs dominate the dust input at all three metallicities.

The dust production timescales, which we define as the length of time during which the sources are forming dust, are provided in Table 6 and show that WC binaries only produce dust for $\Delta t_{WC} = 0.5-4.3$ Myr. Note that the RSG and WC binary timescales in Table 6 only include sources that do not end their lives as SNe. When including sources that end their lives as SNe, the first RSG peak is broader and bridges the second peak in all three metallicity models, and the WC binary dust formation timescale in the solar metallicity model is broader by factor of ~2 ($\Delta t_{WC} = 5.4$ Myr) than the 2.5 Myr timescale reported in Table 6.

As expected, the WC dust production timescale is substantially shorter than the dust production timescales of AGB stars, where $\Delta t_{AGB} = 1-2.4$ Gyr. The short dust production timescale of WC binaries explains the difference in their cumulative dust input rate relative to the longer-lived AGB stars between the constant SF and coeval population models. Given their short timescales, the presence of WC binaries should also reflect the intensity of recent star formation.

The results from the coeval $10^6 M_\odot$ stellar population BPASS models indicate that, for a “bursty” SFH, AGB stars will dominate the dust budget at late times ($t \gtrsim 70$ Myr) over WC binaries and RSGs with SNe as net dust destroyers. However, the LMC-like and solar models show that WC binaries are the first and the dominant source of dust for ~1 Myr before the RSG dust production starts to increase. Notably, in the solar metallicity mode, the WC binaries are the leading dust source until the onset of AGB stars.

### 4.3. LMC BPASS Dust Model and Massive Star Population Discrepancy

The well-studied dust and supernova remnant properties of the nearby LMC ($d \sim 50$ kpc) highlight its utility as a laboratory to test dust input rates in its ISM (e.g., Matsuura et al. 2009; Temim et al. 2015). Importantly, the stellar populations and dust input from AGBs, RSGs, and SNe have been previously measured or estimated observationally, and thus provide a comparison for the BPASS model output (Section 4.2.1). In this section, we apply the $Z = 0.008$ constant star formation BPASS models to the LMC and demonstrate the discrepancy in the massive star populations due to its complex and variable star formation history (SFH; Harris & Zaritsky 2009). The LMC dust-to-gas mass ratio, which is relevant for the SN dust destruction rates, is assumed to be $\chi_{d/g} = 0.0018$ (Temim et al. 2015). The total ISM mass swept up per SN adopted in the BPASS models based on Equation (13) is $m_g = 1947 M_\odot$, consistent with the LMC value measured by Temim et al. (2015) of $m_g \sim 2000 M_\odot$.

In order to model the dust input the LMC, the dust production and stellar population output from the constant SF BPASS models at a metallicity of $Z = 0.008$ were scaled by a constant factor such that the total number of AGB stars at a time $t = 12.6$ Gyr matches the current observed number of

<table><caption>**Table 7**
BPASS LMC Dust Model Table</caption><thead><tr><th></th><th>BPASS (t = 12.6 Gyr)</th><th>From Obs.</th></tr></thead><tbody><tr><td>N<sub>AGB</sub></td><td>20000 (100%)</td><td>20040</td></tr><tr><td>N<sub>RSG</sub></td><td>11887 (30%)</td><td>3589</td></tr><tr><td>N<sub>WC</sub></td><td>52 (100%)</td><td>2</td></tr><tr><td>R<sub>SN</sub> (yr<sup>-1</sup>)</td><td>3.3 × 10<sup>-2</sup></td><td>2.6 × 10<sup>-3</sup></td></tr><tr><td>M&#775;<sub>d,AGB</sub> (M<sub>☉</sub> yr<sup>-1</sup>)</td><td>1.4 × 10<sup>-5</sup> [7.1 × 10<sup>-10</sup>]</td><td>1.4 × 10<sup>-5</sup> [6.9 × 10<sup>-10</sup>]</td></tr><tr><td>M&#775;<sub>d,RSG</sub> (M<sub>☉</sub> yr<sup>-1</sup>)</td><td>5.9 × 10<sup>-6</sup> [1.6 × 10<sup>-9</sup>]</td><td>1.4 × 10<sup>-6</sup> [4.0 × 10<sup>-10</sup>]</td></tr><tr><td>M&#775;<sub>d,WC</sub> (M<sub>☉</sub> yr<sup>-1</sup>)</td><td>3.3 × 10<sup>-6</sup> [6.2 × 10<sup>-8</sup>]</td><td>≥3 × 10<sup>-8</sup></td></tr><tr><td>M&#775;<sub>d,SN</sub> (M<sub>☉</sub> yr<sup>-1</sup>)</td><td>3.4 × 10<sup>-3</sup> [0.10 M<sub>☉</sub>]</td><td>1.7 × 10<sup>-3</sup> [0.65 M<sub>☉</sub>]<sup>a</sup></td></tr><tr><td>M&#775;<sub>SN,dest</sub> (M<sub>☉</sub> yr<sup>-1</sup>)</td><td>-1.1 × 10<sup>-1</sup> [-3.5 M<sub>☉</sub>]</td><td>~3 × 10<sup>-2</sup> [-8 M<sub>☉</sub>]</td></tr><tr><td>Log (t<sub>0,AGB</sub>)</td><td>7.9</td><td>...</td></tr><tr><td>Log (t<sub>0,RSG</sub>)</td><td>6.4</td><td>...</td></tr><tr><td>Log (t<sub>0,WC</sub>)</td><td>6.3</td><td>...</td></tr><tr><td>Log (t<sub>0,SN</sub>)</td><td>6.6</td><td>...</td></tr></tbody></table>

**Notes.** Results from the LMC BPASS dust models at time $t = 12.6$ Gyr compared to observationally derived values (Temim et al. 2015; Srinivasan et al. 2016; Williams 2019). Values in this table are presented in the same format as in Table 5. Note that the observed value of $N_{WC}$ corresponds to the number of known dust-forming WC systems in the LMC. For the SN dust destruction in the BPASS model, it is assumed that the dust-to-gas mass ratio in the LMC is $\chi_{d/g} = 0.0018$ and that the SN shock clears dust from a total surrounding ISM mass of $m_g = 1947 M_\odot$. The dust destruction rates in the right column are from Temim et al. (2015), who derive $m_g \sim 2000 M_\odot$.

<sup>a</sup> Temim et al. (2015) estimate the SN dust production rate in the LMC assuming 100% dust condensation efficiency, and present these as absolute upper limits on the injected dust mass.

AGBs in the LMC, $N_{AGB} \approx 20000$ (Riebel et al. 2012; Zhukovska & Henning 2013; Srinivasan et al. 2016). The current age of the LMC is assumed to be $t = 12.6$ Gyr.

In Table 7, we compare the dust input, populations, and SN rates from the LMC BPASS models with values measured from mid-IR observations of the LMC by Riebel et al. (2012), Srinivasan et al. (2016), and the Magellanic Cloud SN remnant analysis by Temim et al. (2015). The total dust input rate from the AGB stars at time $t = 12.6$ Gyr is closely consistent with the observed rate. However, the SN rate and the total WC binary and RSG dust input rates from the models are larger than the observed values.

Both the BPASS models and results from previous studies show an agreement that the dust production from SNe without dust destruction dominates the other three sources by over an order of magnitude. However, when incorporating the effects of ISM dust destruction from SN shocks, SNe are net destroyers of dust (Temim et al. 2015).

It is apparent that the population of sources associated with massive stars is overpredicted by the constant SF BPASS models: the total number of RSGs at $t = 12.6$ Gyr ($N_{RSG} = 11887$), which includes RSGs that eventually will explode as SNe and not input dust, is a factor of ~3 greater than the observed number of RSGs, and both the SN rate and the number of dust-producing WC binaries exceed the observed value by over an order of magnitude.

For the WC binaries, there are only two currently known dust producers in the LMC, as opposed to the 52 predicted by the models. These two systems, HD 36402 and HD 38030, both host WC4 stars and exhibit dust formation periods of 5.1 and >20 yr, respectively (Williams 2019). The estimated lower limit on the current observed WC binary dust input rate in the
LMC, $3 \times 10^{-8} M_{\odot} \text{yr}^{-1}$ (Table 7) is derived by taking the measured amorphous carbon dust mass of $M_d \sim 1.5 \times 10^{-7} M_{\odot}$ around HD 36402 by Williams (2011) divided by its 5.1 yr orbital period. Interestingly, the DPR of HD 36402 is consistent within a factor of two with the mean dust input rate for the WC binaries in the BPASS model where $\langle \dot{M}_{d,\text{WC}} \rangle = 6.2 \times 10^{-8} M_{\odot} \text{yr}^{-1}$. The dust input rate from the longer-period HD 38030 WC system is not yet known, because it was only recently identified as a periodic/episodic dust former. However, given their longer orbital period ($>5$ yr) and lower IR luminosities (Williams et al. 2013a), these two LMC WC systems are unlikely to be significant dust producers like WR104 where $\dot{M}_d \sim 10^{-6} M_{\odot} \text{yr}^{-1}$ (Table 4).

We attribute the number discrepancy of the massive star population to the assumption of a constant SF in the BPASS models, since the SFH of the LMC is complex and variable (Harris & Zaritsky 2009). Analysis of the SFH in the LMC by Harris & Zaritsky (2009) indicates that the current total stellar mass is significantly composed of stars formed in an initial epoch of star formation that ceased 10 Gyr ago; this was followed by a quiescent period until 5 Gyr ago with recent star formation peaks occurring 12, 100, 500 Myr, and 2 Gyr ago. Normalizing the number of sources from the constant SF BPASS models to the observed AGB stars in the LMC would therefore overpredict the number of RSGs and WC binaries and SN rates, all of which are products of shorter-lived massive stars. Notably, the lifetime of the WC binary dust producers in the BPASS models is $<5 \times 10^6$ Myr, which is shorter than one of the most recent star formation peaks in the LMC that occurred 12 Myr ago (Harris & Zaritsky 2009).

The results from this LMC analysis highlight the importance of considering the SFH in modeling the dust input rates, especially from massive stars. Although this is beyond the scope of the work presented here, future work with BPASS dust models incorporating variable SFH is crucial for obtaining a more complete understanding of galactic dust production rates.

# 5. Discussion and Conclusions

## 5.1. Astrophysical Implications

The WC dustar SED analysis and the predicted dust input from BPASS models indicate that WC binaries are leading sources of dust, comparable to AGB stars, in environments with constant star formation at solar metallicities if SNe are indeed net dust destroyers. For a stellar population formed in an instantaneous starburst at solar metallicities, WC binaries will dominate the dust contribution until the onset of AGBs. WC binaries are also among the earliest sources of dust in LMC and solar metallicities, which implies they are an early reservoir of carbon-rich dust.

Due to the net dust destruction rates from SNe, additional dust input is still necessary to account for the observed dust in galactic ISM (Dwek et al. 2014). For example, in the LMC, an additional dust source with an input rate more than an order of magnitude greater than the combined input from AGB stars and SNe is needed to balance out the SN dust destruction (Temim et al. 2015). The net SN dust destruction rate for all three metallicity models also shows that the dust destruction exceeds the input from AGBs, RSGs, and WC binaries by several orders of magnitude.

Although WC dustars may not provide enough dust to account for this deficit, the large 0.1–1.0 µm sized grains produced by the heavy WR104-like dust-makers would be more robust against destruction from SN shocks. Such large grains would exhibit grain lifetimes a factor of ~3 longer relative to small ~100 Å sized grains (Jones et al. 1996). Therefore, the destruction timescales of the ISM grains would be lengthened if a significant fraction of the ISM was composed of dust from WC binaries.

In the context of cosmochemistry, the formation of large dust grains supports the hypothesis of WR stars as sources of short-lived radionuclides (SLRs) in the solar system. The longer destruction timescales of larger grains suggest WR winds were able to transport and seed SLRs like $Al^{26}$ in the presolar nebula (Dwarkadas et al. 2017).

Given the impact of WC dustars on the ISM, it is important to understand their dust chemistry and the grain formation physics in their hostile wind collision environment. The C-rich dust that forms in the H-deficient and fast winds ($\gtrsim 1000$ km s$^{-1}$) of WC binaries likely follows a chemical formation pathway different from that of C- and H-rich AGBs (Le Teuff & Cherchneff 1998; Cherchneff et al. 2000). Distinguishing emission features of the molecular precursors of C-rich dust condensed in WC binaries from other sources could provide a means of observationally tracing their contribution to galactic ISM. Interestingly, the ISO/SWS spectroscopy of WC dustars reanalyzed in this work also show evidence of hydrocarbon features (Marchenko & Moffat 2017).

Our analysis of the carbon dust condensation fraction (Section 3.5) strengthens our understanding of dust formation in colliding-wind binaries; however, it is important to consider that not all WC+OB binaries form dust. For example, the well-studied colliding-wind system $\gamma^2$ Velorum (a.k.a WR11) hosts a WC8+O7.5III binary with a 78.5 day orbital period, but does not exhibit any dust formation (Schmutz et al. 1997; Lamberts et al. 2017). Systems with shorter orbital periods—and thus shorter orbital separations—may be affected by “sudden radiative braking” of the WR winds due to deceleration as winds approach the luminous OB companion (Gayley et al. 1997). In their analysis, Gayley et al. (1997) suggest that the wind–wind collision $\gamma^2$ Velorum is indeed significantly affected by radiative braking. Tuthill et al. (2008) also explore the impact of radiative braking in WR104. They find that it may play an important role in the wind–wind interaction and shock cone geometry, which are both related to the dust formation. Continued high spatial-resolution observations combined with hydrodynamic simulations of the wind–wind collisions are therefore crucial for understanding the conditions leading to dust formation in colliding-wind systems.

From a dust chemistry perspective, it is interesting to consider why dust formation via colliding winds has only been observed for carbon-rich WRs but not for nitrogen-rich WR (WN) stars. As pointed out by Mauerhan et al. (2015), a C-rich chemistry is not a requirement for dust formation in a colliding-wind binary. For example, the N-rich LBV system $\eta$ Car exhibits episodic dust formation via colliding winds, similarly to WR140, but dust around $\eta$ Car is likely composed of metal-rich silicates, corundum ($Al_2O_3$), and/or iron as opposed to amorphous carbon (Smith 2010; Morris et al. 2017). However, the high mass-loss rate of $\sim 10^{-3} M_{\odot} \text{yr}^{-1}$ from the primary star in $\eta$ Car (Hillier et al. 2001) exceeds the mass-loss rates of WN and WC stars by an order of magnitude. We therefore
speculate that the density of refractory elements such as Al and O may be too low in colliding-wind WN binaries to facilitate dust formation, whereas C in WC winds is readily abundant (40% by mass; Sander et al. 2019).

## 5.2. Conclusions

We have conducted an SED analysis of a sample of 19 Galactic WC dustars with archival ground- and space-based IR photometry and spectroscopy obtained over the past several decades and revised Gaia distance estimates by Rate & Crowther (2020). The SEDs were modeled using DustEM with single or double dust component models.

Results from the SED modeling show that WC dustars can produce small (0.01–0.1 µm) or large (0.1–1.0 µm) dust grains, where large grains are favored for the efficient ($\chi_C \gtrsim 1\%$) dust-makers WR 95, WR98a, WR104, WR106, and WR118. We find that WC dustars exhibit a range of dust production rates and carbon dust condensation fractions, which increase with increasing IR luminosity (Section 3.4.8). We also fit the empirical relation between the DPR and IR luminosity that applies to the nonperiodic dustars (Equation (8)). For the WC dustars with known periods and under the assumption of consistent wind momentum ratios ($\eta$), we successfully demonstrate the predicted power-law relation between $\chi_C$, $M_{WR}$, $v_{WR}$, and $P_{orb}$ (Figure 11) for theoretical models of dust-forming colliding-wind WR systems by Usov (1991).

In order to study the dust contribution of WC dustars in comparison to other leading stellar dust producers (AGBs, RSGs, SNe), we performed binary population synthesis models using BPASS incorporating dust production prescriptions for the four sources in addition to dust destruction from SNe. For the WC binaries, the dust production rate in the BPASS models was defined via the verified dust condensation relation in Equation (11) from the SED analysis results and scaled to the properties of WR104.

We performed two sets of BPASS dust models at three different metallicities ($Z = 0.001, 0.008$, and $0.020$) representing different phases of cosmic time, assuming constant star formation and a coeval $10^6 M_\odot$ stellar population. Both the constant SF and coeval models show that SNe produce the most dust out of all the sources, but are net dust destroyers at all metallicities when considering the ISM dust destruction from SN shocks. At low ($Z = 0.001$) metallicities, AGB stars input a comparable amount of dust to RSGs and outproduce WC binaries by over three orders of magnitude. At LMC-like ($Z = 0.008$) metallicities, the AGB stars slightly outproduce WC binaries and RSGs by factors of 2–3, whereas at solar metallicities, WC binaries are the dominant source of dust until the onset of AGBs, which form dust at a comparable rate (Figure 12, Table 5).

Results from the coeval population models indicate that, for a “bursty” star formation history, AGB stars become significant sources of dust at late times after their onset ($t > 100$ Myr) but still do not produce as much dust as SNe without considering dust destruction. However, at LMC-like and solar metallicities, the AGBs dominate the dust production over WC binaries and RSGs, as well as SNe, if they are indeed net dust destroyers at these metallicities (Figure 12, Table 6).

Constant star formation BPASS models of the LMC ($Z = 0.008$) and comparisons to the observationally derived populations of AGBs, RSGs, WC binaries, and SNe demonstrated the inconsistencies between the observed and model

massive star populations. The SN rates and the population of WC binaries and RSGs from the constant SF BPASS model overpredicted the observed values, which we attribute to the complex and variable star-forming history of the LMC (Harris & Zaritsky 2009).

Finally, our conclusions regarding WC binaries as important but overlooked dust sources strongly motivates further pursuits of BPASS dust models with variable SFH and future studies addressing the open questions pertaining to WC dust formation physics and dust chemistry. With the combination of high sensitivity and spatial resolution in the mid-IR, observations with the upcoming James Webb Space Telescope will transform our understanding of WC dustars and dust formation in hostile environments.

We thank T. Onaka, I. Endo, M. Buragohain, E. Lagadec, T. Nozawa, A. Soulain, T. Bakx, T. Takeuchi, and M. Corcoran for enlightening discussions on dust, SNe, AGBs, RSGs, and WR binaries. We also thank Greg Sloan for guidance in assessing the quality of ISO/SWS data. We would also like to thank the anonymous referee for their careful review and insightful comments that have improved the focus and clarity of our work. P.M.W. thanks the Institute for Astronomy for continued hospitality and access to the facilities of the Royal Observatory. R.M.L. acknowledges the Japan Aerospace Exploration Agency’s International Top Young Fellowship. This work is based in part on observations made with the Spitzer Space Telescope, obtained from the NASA/IPAC Infrared Science Archive, both of which are operated by the Jet Propulsion Laboratory, California Institute of Technology under a contract with the National Aeronautics and Space Administration. Based in part on observations collected at the European Organisation for Astronomical Research in the Southern Hemisphere under ESO program 097.D-0707(A). Based in part on observations with ISO, an ESA project with instruments funded by ESA Member States (especially the PI countries: France, Germany, the Netherlands, and the United Kingdom) and with the participation of ISAS and NASA. This publication makes use of data products from the Wide-field Infrared Survey Explorer, which is a joint project of the University of California, Los Angeles, and the Jet Propulsion Laboratory/California Institute of Technology, funded by the National Aeronautics and Space Administration.

**Facilities:** Spitzer(IRS, IRAC, and MIPS), WISE, ISO (SWS), VLT(VISIR).

**Software:** BPASS, DustEM, MIRA.

## ORCID iDs

J. J. Eldridge https://orcid.org/0000-0002-1722-6343

Matthew J. Hankins https://orcid.org/0000-0001-9315-8437

Itsuki Sakon https://orcid.org/0000-0001-7641-5497

Peredur M. Williams https://orcid.org/0000-0002-8092-980X

## References

Ardila, D. R., Van Dyk, S. D., Makowiecki, W., et al. 2010, ApJS, 191, 301
Bailer-Jones, C. A. L., Rybizki, J., Fouesneau, M., Mantelet, G., & Andrae, R. 2018, AJ, 156, 58

Boyer, M. L., Srinivasan, S., Riebel, D., et al. 2012, ApJ, 748, 40

Carey, S. J., Noriega-Crespo, A., Mizuno, D. R., et al. 2009, PASP, 121, 76

Cherchneff, I. 2015, in Wolf-Rayet Stars, ed. W.-R. Hamann, A. Sander, & H. Todt (Potsdam: Universitätsverlag Potsdam), 269

Cherchneff, I., Le Teuff, Y. H., Williams, P. M., et al. 2000, A&A, 357, 572
Chiar, J. E., & Tielens, A. G. G. M. 2001, ApJL, 550, L207

Chiar, J. E., & Tielens, A. G. G. M. 2006, ApJ, 637, 774

Churchwell, E., Babler, B. L., Meade, M. R., et al. 2009, PASP, 121, 213

Cohen, M., Barlow, M. J., & Kuhi, L. V. 1975, A&A, 40, 291

Compiègne, M., Verstraete, L., Jones, A., et al. 2011, A&A, 525, A103

Conti, P. S. 1975, MSRSL, 9, 193

Conti, P. S., & Vacca, W. D. 1990, AJ, 100, 431

Crowther, P. A. 1997, MNRAS, 290, L59

Crowther, P. A. 2007, ARA&A, 45, 177

Crowther, P. A., De Marco, O., & Barlow, M. J. 1998, MNRAS, 296, 367

Crowther, P. A., Morris, P. W., & Smith, J. D. 2006, ApJ, 636, 1033

Cutri, R. M. 2013, yCat, 2328, 0

Cutri, R. M., Skrutskie, M. F., van Dyk, S., et al. 2003, yCat, 2246, 0

Danks, A. C., Wamsteker, W., Shaver, P. A., & Retallack, D. S. 1984, A&A, 132, 301

de Graauw, T., Haser, L. N., Beintema, D. A., et al. 1996, A&A, 315, L49

De Marco, O., & Izzard, R. G. 2017, PASA, 34, e001

Dougherty, S. M., Beasley, A. J., Claussen, M. J., et al. 2005, ApJ, 623, 447

Dwarkadas, V. V., Dauphas, N., Meyer, B., et al. 2017, ApJ, 851, 147

Dwek, E., & Cherchneff, I. 2011, ApJ, 727, 63

Dwek, E., Staguhn, J., Arendt, R. G., et al. 2014, ApJL, 788, L30

Eldridge, J. J., Stanway, E. R., Xiao, L., et al. 2017, PASA, 34, e058

Fazio, G. G., Hora, J. L., Allen, L. E., et al. 2004, ApJS, 154, 10

Gall, C., & Hjorth, J. 2018, ApJ, 868, 62

Gayley, K. G., Owocki, S. P., & Cranmer, S. R. 1997, ApJ, 475, 786

Gehrz, R. D., & Hackwell, J. A. 1974, ApJ, 194, 619

Gräfener, G., Koesterke, L., & Hamann, W.-R. 2002, A&A, 387, 244

Gutermuth, R. A., & Heyer, M. 2015, AJ, 149, 64

Hankins, M. J., Lau, R. M., Morris, M. R., et al. 2016, ApJ, 827, 136

Harries, T. J., Monnier, J. D., Symington, N. H., & Kurosawa, R. 2004, MNRAS, 350, 565

Harris, J., & Zaritsky, D. 2009, AJ, 138, 1243

Hendrix, T., Keppens, R., van Marle, A. J., et al. 2016, MNRAS, 460, 3975

Herter, T. L., Vacca, W. D., Adams, J. D., et al. 2013, PASP, 125, 1393

Hillier, D. J., Davidson, K., Ishibashi, K., et al. 2001, ApJ, 553, 837

Houck, J. R., Roellig, T. L., van Cleve, J., et al. 2004, ApJS, 154, 18

Howarth, I. D., & Schmutz, W. 1992, A&A, 261, 503

Humphreys, R. M., & Davidson, K. 1994, PASP, 106, 1025

Jones, A. P., Tielens, A. G. G. M., & Hollenbach, D. J. 1996, ApJ, 469, 740

Kanarek, G., Shara, M., Faherty, J., et al. 2015, MNRAS, 452, 2858

Kessler, M. F., Steinz, J. A., Anderegg, M. E., et al. 1996, A&A, 315, L27

Kochanek, C. S. 2011, ApJ, 743, 73

Kroupa, P., Tout, C. A., & Gilmore, G. 1993, MNRAS, 262, 545

Lagage, P. O., Pel, J. W., Authier, M., et al. 2004, Msngr, 117, 12

Lamberts, A., Millour, F., Liermann, A., et al. 2017, MNRAS, 468, 2655

Le Teuff, Y. H., & Cherchneff, I. 1998, in ASP Conf. Proc. 147, Abundance Profiles: Diagnostic Tools for Galaxy History, ed. D. Friedli et al. (San Francisco, CA: ASP), 188

Lebouteiller, V., Barry, D. J., Goes, C., et al. 2015, ApJS, 218, 21

Lefèvre, L., Marchenko, S. V., Lépine, S., et al. 2005, MNRAS, 360, 141

Marchenko, S. V., & Moffat, A. F. J. 2007, in ASP Conf. Ser. 367, Massive Stars in Interactive Binaries, ed. N. St.-Louis & A. F. J. Moffat (San Francisco, CA: ASP), 213

Marchenko, S. V., & Moffat, A. F. J. 2017, MNRAS, 468, 2416

Marchenko, S. V., Moffat, A. F. J., Vacca, W. D., et al. 2002, ApJL, 565, L59

Massey, P., Plez, B., Levesque, E. M., et al. 2005, ApJ, 634, 1286

Matsuura, M., Barlow, M. J., Zijlstra, A. A., et al. 2009, MNRAS, 396, 918

Mauerhan, J., Smith, N., Van Dyk, S. D., et al. 2015, MNRAS, 450, 2551

Midooka, T., Sugawara, Y., & Ebisawa, K. 2019, MNRAS, 484, 2229

Millour, F., Driebe, T., Chesneau, O., et al. 2009, A&A, 506, L49

Monnier, J. D., Tuthill, P. G., & Danchi, W. C. 1999, ApJL, 525, L97

Monnier, J. D., Tuthill, P. G., Danchi, W. C., et al. 2007, ApJ, 655, 1033

Monnier, J. D., Zhao, M., Pedretti, E., et al. 2011, ApJL, 742, L1

Morris, P. W., Brownsberger, K. R., Conti, P. S., Massey, P., & Vacca, W. D. 1993, ApJ, 412, 324

Morris, P. W., Gull, T. R., Hillier, D. J., et al. 2017, ApJ, 842, 79

Niemela, V. S. 1995, in Proc. IAU Symp. 163, Wolf-Rayet Stars: Binaries; Colliding Winds; Evolution, Vol. 163, ed. K. A. van der Huch & P. M. Williams (Dordrecht: Kluwer), 223

Nugis, T., & Lamers, H. J. G. L. M. 2000, A&A, 360, 227

Rajagopal, J., Menut, J.-L., Wallace, D., et al. 2007, ApJ, 671, 2017

Rate, G., & Crowther, P. A. 2020, MNRAS, 493, 1512

Riebel, D., Srinivasan, S., Sargent, B., et al. 2012, ApJ, 753, 71

Rieke, G. H., Young, E. T., Engelbracht, C. W., et al. 2004, ApJS, 154, 25

Roche, P. F., & Aitken, D. K. 1984, MNRAS, 208, 481

Roslowe, C. K., & Crowther, P. A. 2015, MNRAS, 447, 2322

Roslowe, C. K., & Crowther, P. A. 2018, MNRAS, 473, 2853

Sana, H., de Mink, S. E., de Koter, A., et al. 2012, Sci, 337, 444

Sander, A., Hamann, W.-R., & Todt, H. 2012, A&A, 540, A144

Sander, A. A. C., Hamann, W.-R., Todt, H., et al. 2019, A&A, 621, A92

Sarangi, A., Matsuura, M., & Micelotta, E. R. 2018, SSRv, 214, 63

Schmutz, W., Schweickhardt, J., Stahl, O., et al. 1997, A&A, 328, 219

Shara, M. M., Moffat, A. F. J., Smith, L. F., et al. 1999, AJ, 118, 390

Skrutskie, M. F., Cutri, R. M., Stiening, R., et al. 2006, AJ, 131, 1163

Sloan, G. C., Kraemer, K. E., Price, S. D., & Shipman, R. F. 2003, ApJS, 147, 379

Smith, L. F. 1968, MNRAS, 140, 409

Smith, L. F., Shara, M. M., & Moffat, A. F. J. 1990, ApJ, 358, 229

Smith, N. 2010, MNRAS, 402, 145

Smith, N. 2014, ARA&A, 52, 487

Soulain, A., Millour, F., Lopez, B., et al. 2018, A&A, 618, A108

Srinivasan, S., Boyer, M. L., Kemper, F., et al. 2016, MNRAS, 457, 2814

Stanway, E. R., & Eldridge, J. J. 2018, MNRAS, 479, 75

Sugawara, Y., Maeda, Y., Tsuboi, Y., et al. 2015, PASJ, 67, 121

Temim, T., Dwek, E., Tchernyshyov, K., et al. 2015, ApJ, 799, 158

Tuthill, P. G., Monnier, J. D., & Danchi, W. C. 1999, Natur, 398, 487

Tuthill, P. G., Monnier, J. D., Lawrance, N., et al. 2008, ApJ, 675, 698

Usov, V. V. 1991, MNRAS, 252, 49

van der Hucht, K. A. 2001, NewAR, 45, 135

van der Hucht, K. A., Hidayat, B., Admiranto, A. G., et al. 1988, A&A, 199, 217

van der Hucht, K. A., Morris, P. W., Williams, P. M., et al. 1996, A&A, 315, L193

van Loon, J. T. 2000, A&A, 354, 125

Werner, M. W., Roellig, T. L., Low, F. J., et al. 2004, ApJS, 154, 1

Williams, D. M. 2003, AmJPh, 71, 1198

Williams, P. M. 2011, BSRSL, 80, 195

Williams, P. M. 2019, MNRAS, 488, 1282

Williams, P. M., Chu, Y.-H., Gruendl, R. A., et al. 2013a, MNRAS, 431, 1160

Williams, P. M., Cohen, M., van der Hucht, K. A., Bouchet, P., & Vacca, W. D. 1995, MNRAS, 275, 889

Williams, P. M., Kidger, M. R., van der Hucht, K. A., et al. 2001, MNRAS, 324, 156

Williams, P. M., Marchenko, S. V., Marston, A. P., et al. 2009a, MNRAS, 395, 1749

Williams, P. M., Rauw, G., & van der Hucht, K. A. 2009b, MNRAS, 395, 2221

Williams, P. M., & van der Hucht, K. A. 2000, MNRAS, 314, 23

Williams, P. M., van der Hucht, K. A., Bouchet, P., et al. 1992, MNRAS, 258, 461

Williams, P. M., van der Hucht, K. A., Kidger, M. R., et al. 1994, MNRAS, 266, 247

Williams, P. M., van der Hucht, K. A., Pollock, A. M. T., et al. 1990, MNRAS, 243, 662

Williams, P. M., van der Hucht, K. A., & The, P. S. 1987, A&A, 182, 91

Williams, P. M., van der Hucht, K. A., van Wyk, F., et al. 2012, MNRAS, 420, 2526

Williams, P. M., van der Hucht, K. A., van Wyk, F., et al. 2013b, MNRAS, 429, 494

Wright, E. L., Eisenhardt, P. R. M., Mainzer, A. K., et al. 2010, AJ, 140, 1868

Yamasawa, D., Habe, A., Kozasa, T., et al. 2011, ApJ, 735, 44

Yudin, B., Balega, Y., Blöcker, T., et al. 2001, A&A, 379, 229

Zhekov, S. A., Tomov, T., Gawronski, M. P., et al. 2014, MNRAS, 445, 1663

Zhukovska, S., & Henning, T. 2013, A&A, 555, A99

Zubko, V., Dwek, E., & Arendt, R. G. 2004, ApJS, 152, 211

Zubko, V. G. 1998, MNRAS, 295, 109

Zubko, V. G., Mennella, V., Colangeli, L., et al. 1996, MNRAS, 282, 1321