The long-period spectroscopic orbit and dust creation in the Wolf-Rayet binary system WR 125

Several long-period binaries with a carbon-rich Wolf-Rayet star and an O star produce dust in their wind collisions. In eccentric binaries, this is seen most strongly near periastron passage. The exact conditions leading to dust creation require orbital properties to be determined, which is difficult owing to their long periods. Recently, the binary system WR 125 (WC7+O9III) began a dust creation episode seen through an infrared outburst first detected by NEOWISE-R, which was the first outburst detected since 1991. We present new near- and mid-infrared photometry, which we use to show consistency between the two outbursts and derive an orbital period of 28.12$^{+0.10}_{-0.05}$ yr. We use a long time-series of optical spectra to place the first constraints on its orbital elements, on the assumption that this system will produce dust near periastron. The orbit has a mild eccentricity of 0.29$\pm$0.12 and is only derived for the Wolf-Rayet component, as the O star's radial velocities have noise that is likely larger than the expected semi-amplitude of the orbit. We also present SOFIA/FORCAST grism spectroscopy to examine the infrared spectral energy distribution (SED) of the dust during this outburst, comparing its properties to other WCd binaries, deriving a dust temperature of 580 K in 2021. This collection of observations will allow us to plan future observations of this system and place the system in the context of dust-creating Wolf-Rayet binaries.


INTRODUCTION
A classical Wolf-Rayet star is a massive, evolved star that has lost its outer envelope to reveal a compact, hydrogen-depleted, hot star with a strong stellar wind.These stars are usually either nitrogen-rich (spectral type WN) or carbon-rich (WC), and could form through either strong stellar winds (Conti 1975) which would likely also include episodic mass loss (Smith & Owocki 2006), or through binary interactions (e.g., Vanbeveren et al. 1998;Eldridge 2009).The WC stars are sometimes observed to form dust, which is inferred through excess emission at infrared wavelengths first discovered by Allen et al. (1972).The dust production is often seen as persistent or periodic, but episodic bursts are also seen in some systems.Periodic and persistent dust makers are often seen in binary systems as first discovered with WR 140 (WC7+O5.5I;Williams et al. 1978) and WR 48a (WC9; Danks et al. 1983).Both of these systems were observed to show rapid increases in the infrared attributed to dust formation.WR 140 was subsequently shown to be a high-eccentricity binary with dust formation near periastron passage (Williams et al. 1990).It is predicted that near the time of a periastron passage, the density of wind collisions between the WC star and a companion OB star trigger dust formation.In contrast to the long-period episodically-producing dustar WR 140, persistent dust makers tend to be in low eccentricity orbits.Curiously, some such systems like γ 2 Velorum have never been observed to form dust, although very similar to ones that do.
The dust production in WC binaries could be of cosmological importance.Both Marchenko & Moffat (2007) and Lau et al. (2021) suggest that the formation of WC binaries could be the first source of dust in a low metallicity environment.They can form dust earlier than supernovae because the higher mass and hence more rapidly evolving WC star can form through binary interactions before the supernova takes place, where the primary star loses its envelope through Roche lobe overflow onto its OB companion.Then, as the Roche lobe overflow ends, the strong stellar winds of the two stars will form a shock interface conducive to dust formation down stream from the heated apex.
Of the WC binary "dustars", the prototype is often considered to be WR 140.This system has a very well-established orbit, both spectroscopically (Fahed et al. 2011;Thomas et al. 2021) and with a visual orbit established through interferometry (Monnier et al. 2011;Thomas et al. 2021).It has a long 7.93 yr period with high eccentricity (e = 0.8993 ± 0.0013) with the dust formation triggered by a changing gas density in the shock front near periastron.At this binary phase, the increase in the density of the shocked gas prompts increased X-ray production, but this seems to switch to cooling via optical emission lines at phases very close to periastron, perhaps allowing the conditions for dust production to occur (Pollock et al. 2021).Williams et al. (2009) found that the dust around this prototype "dustar" (WR 140) survived at least two cycles from imaging fossil dust around the system at midinfrared wavelengths.The overarching question of dust survivability in these hostile environments was still open however in part due to limited spatial resolution and sensitivity with ground-based imaging, although Marchenko et al. (2003) showed from ground-based mid-IR imagery that the dust in WR140 very likely reached the ISM.With a short exposure using JWST +MIRI at mid-infrared wavelengths, Lau et al. (2022) showed that the dust survives at distances out to 70,000 AU from WR 140, likely implying that the dust should survive and be included in the dust budget of galaxies.WR 125 is a near spectroscopic twin to WR 140 (Abbott et al. 1986), and is the topic of this paper, first noted as IC 14-36 by Iriarte & Chavira (1956).Williams et al. (1987) observed this star with infrared photometry but found it not to show variation in its infrared flux.However, its near twin status prompted continued infrared monitoring until Williams et al. (1992) discovered an infrared outburst in the years 1990-1991.Williams et al. (1994) reported on additional infrared photometry of WR 125 during this outburst and presented some of the first optical spectroscopy of the system.They confirmed that the emission line dilution is caused by an O9III companion star based on the optical absorption lines present.Furthermore, Williams et al. (1994) presented infrared spectroscopy of the system and the dust cloud, finding no signatures of the WR wind lines in the 10 µm region, nor the 11.52 µm graphite feature, suggesting that the dust had to be amorphous carbon.
Very little progress was made on WR 125 for many years following this initial discovery of dust creation and that the star had remarkable similarities to WR 140.While the star had been occasionally included in infrared sky surveys, infrared monitoring was largely absent.The launch of the Wide-field Infrared Survey Explorer (WISE) in 2009, which was reactivated (Mainzer et al. 2014) as the Near-Earth Object Wide-field Infrared Survey Explorer (NEOWISE-R) in 2013 December following a hibernation of the satellite that began in 2011 February.While NEOWISE-R was unable to observe at the longer wavelengths from the WISE mission, its continued sky monitoring at 3.4 and 4.6 µm showed that the WR 125 system began a second dust creation outburst in 2018 as reported by Williams (2019), indicating that the binary was once again approaching the phases amenable to dust production.
The increase in infrared flux observed with the NEOWISE-R mission (Williams 2019) has prompted a few new studies into WR 125.Midooka et al. (2019) reported on a few public X-ray observations of WR 125 taken with the Swift and XMM-Newton X-ray observatories.The observations, largely taken prior to the start of the current dusty outburst, were fairly constant in their flux and shape.The flux was the same as observed in 1981 with the Einstein satellite, but the flux was lower during the previous dust outburst in 1991 as observed with ROSAT.This is likely similar to that of the prototype of the dust-making WC binaries, WR 140, which also shows an X-ray dip at some phases close to periastron as described by Pollock et al. (2021).
In addition to the X-ray observations, Endo et al. (2022) presented a mid-infrared spectrum of the system taken with Subaru and the COMICS instrument.They found that the system shows a broad 8 µm feature that was also seen in several WCd stars observed with ISO and the SWS instrument.This feature is a typical "unidentified infrared" (UIR) band, which often correlates to other infrared features.These UIR features are also seen in spatially-resolved spectroscopy of the dust surrounding WR 140 and imaged with JWST (Lau et al. 2022).Endo et al. (2022) presented MIR spectroscopic observations of WR 125 obtained in 2019 that were characterized with a blackbody temperature of nearly 800 K after removing the underlying stellar flux.This was hotter dust than seen in the other WR binaries they used as a comparison from ISO observations, likely as it was more recently formed.They also derived a period of ∼28.1 yr by a comparison of a single infrared flux point with the light curve presented by Williams et al. (1992Williams et al. ( , 1994)).
Many of these results show that WR 125 is a prime target for additional observations and understanding the dust production around these systems.To date, no spectroscopic orbit has been measured, and a time-series analysis of the infrared photometry is strongly needed.To that end, we present new observations in Section 2 and analyze the time-series photometry in Section 3 to obtain a period with higher confidence than elsewhere.In Section 4, we present the first spectroscopic orbit of the system.We present an infrared spectrum taken with SOFIA and FORCAST in Section 5. We discuss these findings in Section 6, and then conclude our study.

Infrared Photometry
Near-and mid-infrared photometry of WR 125 has been collected with the Sternberg Astronomical Institute's (SAI) Crimean Laboratory of Moscow State University using a JHKLM photometer with an InSb photovoltaic detector cooled with liquid nitrogen (Shenavrin et al. 2011).Observations have been taken regularly since the beginning of the current dust creation episode that was observed with NEOWISE-R (Williams 2019) and are shown in Fig. 1; subsequent observations of WR 125 with NEOWISE-R are severely saturated.The new data allow comparison with the previously observed outburst.One of these observations was also used by Endo et al. (2022).Each data point consists of multiple sub-exposures of 30-60 s with a total integration time of 5-10 minutes in each of the JHKL filters and 20-25 minutes in the M filter.The observations are compared to the standard star BS 7488 (Johnson et al. 1966) observations taken before or after each observation of the target.The light curve is tabulated as an online data file.In addition to these new data, along with the archival measurements from Williams et al. (1992Williams et al. ( , 1994)), we include the two measurements from 2007 and 2008 taken at UKIRT and reported by Endo et al. (2022) in this analysis.For our analysis, we have combined L and L ′ data into one light curve, as this is also how archival data (e.g., Williams et al. 1994) have been presented.

Optical Spectroscopy
We have collected spectroscopy of WR 125 with three instruments.These spectra typically have a signal-to-noise of at least 100 at wavelengths longer than ∼ 5000 Å and were all reduced using standard pipelines and techniques for the instruments used.
The first set of spectra were obtained with the Keck I telescope and the Low Resolution Imaging Spectrometer (LRIS; Oke et al. 1995;Rockosi et al. 2010) on seven independent nights, which provides spectra with wavelength coverage from ∼ 4900 Å to ∼ 6200 Å.The spectrograph had a dispersion of 0.64 Å pixel −1 and we typically obtained a signal-to-noise of 100 per sub-exposure with a 5-minute exposure.
We also collected spectra with the Keck II telescope and the Echellette Spectrograph and Imager (ESI; Sheinis et al. 2002).These data have a resolving power of 13,000 and cover several orders of the spectrum between 3900 Å and 1.1 µm.Near the C III λ5696 line, our data had a typical signal-to-noise of 100 in an order of the spectrograph covering 5100 to 5970 Å with a dispersion of 0.21 Å pixel −1 .A typical spectrum was obtained with a 5-minute exposure.
Lastly, we obtained a few spectra with Gemini-North and the Gemini Multi-Object Spectrograph (GMOS; Hook et al. 2004).We used the B600 grating with the G5307 blocking filter to obtain spectra between 4115 Å and 7410 Å, with a dispersion of 1.0 Å pixel −1 .In the vicinity of the C III λ5696 line, our data had a typical signal-to-noise of 100 for an exposure of 100 s.

SOFIA Infrared Spectroscopy
We obtained mid-IR grism spectroscopy using the FORCAST instrument (Herter et al. 2013(Herter et al. , 2018) ) onboard SOFIA (Temi et al. 2018) during the peak of the recent dust creation outburst.These observations, taken on 2021 April 8, used both the G063 and G111 grating setups.The G063 setup covers the wavelength range of 4.9-8.0µm, but we found the spectra near the edge of the chip to be too noisy for analysis, effectively reducing the useful range to 5.1-7.9µm.The G111 setup was used for a spectrum in the range of 4.9-8.0µm.Both of these spectra were taken with the 2.4 ′′ slit, yielding a spectral resolving power of ∼ 180 and ∼ 260 respectively.There was an additional spectrum taken with the G227 setup that covers the 17.6--27.7 µm range.This setup provided very low signal-to-noise, so we use this primarily to provide a flux point at 23 µm.All data were processed through the typical SOFIA pipelines.The FORCAST grism data are known to have variable slit-losses, adding uncertainty to the flux calibration, but these data should be accurate to a few percent (e.g., Gehrz et al. 2021).Williams et al. (1992Williams et al. ( , 1994) ) as well as new measurements of WR 125.NEOWISE-R W 1 and W 2 after 2018 are too heavily saturated to use.We have also incorporated measurements from 2MASS.For each bandpass, we have shifted the points to the same quiescent level to highlight the outbursts and different amplitudes as a function of wavelength.The new observations are tabulated in a machine-readable format as data behind the figure.The vertical lines indicate the time of the infrared spectrum reported by Endo et al. (2022) in 2019, the 2020 spectrum reported by Endo (2022, PhD diss) and our SOFIA spectrum obtained in 2021.

THE ORBITAL PERIOD FROM THE INFRARED PHOTOMETRY
Fig. 1 shows the infrared light curve of WR 125.In this figure, we see that the near-infrared J-band flux does not show much excess during the two outbursts, while H-band shows a slight excess.The K-band flux reaches nearly a magnitude of excess during an outburst, while L/L ′ -band and M-band reach about 2.5 and 3 magnitudes of excess compared to the quiescent flux.We note that the amplitude of the L/L ′ and W 1 outbursts are heterogeneous, but the amplitude of variation at these wavelengths is so large that we have not attempted to put the data on a common scale.The addition of the two points in 2007-2008 from UKIRT and knowing the K-band magnitude reported in 2MASS (taken between 1997.5 and 2001.1 Skrutskie et al. 2006) was 8.214±0.017,provides enough information to know that these two recorded outbursts were consecutive outbursts for the system despite poor coverage in the time-series as they rule out the presence of another outbust mid-way between the two observed outbursts.
Given the nature of WR 125, we can use the infrared light curve to derive a binary period.We began our time-series analysis with Fourier techniques using Period04 (Lenz & Breger 2005).Fourier techniques assume a sinusoidal function, which is not the observed shape of the light curve we see in Fig. 1, so the results tended to lie on harmonics of the fundamental (orbital) period that is seen to be ∼ 28.1 yr (Endo et al. 2022).Determining the true period based on harmonics introduces extra errors because of the long time-scales involved.As we wish to use this period for a spectroscopic orbit, we needed a different approach.
We then used the phase dispersion minimization routines described by Stellingwerf (1978) to determine the period independent of the light curve shape.This method calculates a dispersion statistic at each potential period, and a minimum in this quantity represents a candidate period.We show this in Fig. 2, where we see a minimum at presumed period at 28.12 +0.10 −0.05 years, as well as the harmonic at 2/3 the period (3/2 the frequency) of ∼ 18 years.The error on this period comes from the statistical uncertainty in the pdm statistic.We calculated this statistic for the Kand L-band light curves.The shorter-wavelength bandpasses had a smaller variability amplitude, casting doubt on the minimum found, while the longer-wavelength Mband has a sparser light curve for the calculation.
From the K-and L-band light curve analyses, we derive a period of 28.12 +0.10 −0.05 years by fitting the minima with Gaussian curves, overlaid as red and blue curves in Fig. 2.This period is consistent with other recent papers (Arora et al. 2021;Endo et al. 2022), and we show the phased light curves in Fig. 3.In the phased light curves, we assume that the outburst begins as the WISE flux began going into an excess, and that the plateau in the K-band peak represents the end of new dust production.We then take the mid-point time between these two epochs to be the time of a periastron passage, assuming that the dust-creation event should be centered on the periastron passage.We checked this period with the Lafler-Kinman approach for sparse timeseries (e.g., Saha & Vivas 2017) and found the same result.In order to measure the radial velocities of WR 125, we employed similar techniques as that done in the WR 140 system by Thomas et al. (2021).We examined the fairly isolated C III λ5696 emission line to measure radial velocities of the WR star.We examined the normalization of each spectrum in the region surrounding the emission line and applied local corrections as necessary.Unlike WR 140, we did not see any strong evidence of excess emission in the C III line although Williams et al. (1992) observed some excess red emission in the He I λ10830 line.We measured the line using bisectors at five points above the continuum level.A radial velocity measurement was taken to be the average of these five points, with a standard deviation of the individual points used as an error measurement, and are shown in Table 1 and illustrated in Fig. 4. For our data taken with LRIS, we often had several sub-exposures which we measured independently to verify consistency in our measurements.
We attempted to measure the O star's absorption line velocities using the He I λ5876 line, but found that this did not produce realistic measurements that would show an orbit with time.This line is seen in emission in late sub-type WC stars and it is possible that this occurs in the spectrum of WR 125, making interpretation difficult.Our spectra did not extend to the higher Balmer absorption lines attributed to the companion (Williams et al. 1994).
Given the sparse nature of the spectroscopic measurements, we found that fitting an orbit was not a straight-forward task.We utilized the orbit-fitting procedures used for the γ 2 Velorum system (Richardson et al. 2017).The RV data span about 20 years so we did not search for a period in the RV data.Spectroscopy over at least another decade will be required to cover the dust-formation cycle and longer to establish a more definitive period.Meanwhile, therefore, by analogy with WR 140 and WR 137 whose dust-formation and spectroscopic periods are the same, we fixed the period of WR 125 to 28.12 years based on the infrared light curves (see previous section).
Similarly, to estimate the date of periastron passage, we look to the infrared light curve.It is evident that the dust condensation by WR 125 occurred over a relatively longer interval than that of WR 140 but in both cases the only way the wind-collision region could "know" the orbital phase is through the variation of the ambient preshock wind density and stellar radiation at the collision, which vary as the inverse square of the stellar separation.The conditions necessary for dust condensation must occur at critical stellar separations, and likely symmetrical in phase about periastron passage, as demonstrated in the case of WR 140 by Han et al. (2022).For the present we therefore assume that periastron passage occurs between the beginning and end of dust condensation.To that end, we used the first observation with NEOWISE-R as the system began to brighten to estimate the onset of dust production (Williams

2019)
. We estimate that the dust production ends when the infrared light curve reaches a maximum in the near-infrared K-band, after which dust being carried away by the stellar winds is no longer being replenished by the condensation of new dust.
We then estimate that the midpoint of these two times was the periastron passage to fit an orbit.The result of this fit is shown in Fig. 5, with the orbital elements  given in Table 2.The resulting orbit has a moderate eccentricity, but will need better spectroscopic coverage in the coming decades to be refined.

THE MID-INFRARED SPECTRUM AS OBSERVED BY SOFIA
The spectrum of WR 125 observed by SOFIA is shown in Fig. 6.The three grism spectra are of moderate signal-to-noise.The raw flux shows a peak near 6-7 µm, and the data taken from 17.6-27.7 µm with the G227 grating are too noisy for a thorough analysis given the lower flux and sensitivity of the instrument at these wavelengths.Following the procedures used by Endo et al. (2022), we aimed to correct the fluxes for interstellar extinction using the extinction laws from Weingartner & Draine (2001) (WD01) and Gordon et al. (2021) (G21).Endo et al. (2022) noted that the value of A V was calculated to be 5.89 ± 0.75 from A v = 6.48 ± by Rate & Crowther (2020) from the relation A v = 1.1AV (Turner 1982).We show the extinction correction for A V = 5.89 as well as the values at the extrema of the value reported by Rate & Crowther (2020) in Fig. 6.
We also used the photometry obtained at SAI converted to flux1 in our analysis of the SOFIA data.We began our analysis by subtracting the stellar contribution that was calculated by Endo et al. (2022).The stellar contribution was calculated as a power-law with two forms based on the extinction law adopted.The power-law fit held the form of a × λ −b , where the coefficients were derived from dust-free epochs of infrared photometry and correspond to (a, b) = (2.08,1.58) for the WD01 extinction law and (a, b) = (1.79,1.40) for the G21 extinction law.In both cases, the units of this relationship provide the flux in Jy and use a wavelength given in µm and was fit with the infrared flux observed in the 1980s when no dust was seen in the SED.The J-band light curve shows little or no variability, so this relationship will provide the normal J-band flux.In both extinction situations, a blackbody fit provides a temperature for the dust of 400 K (Fig. 7) if we neglect the emissivity of the dust.This is cooler than the temperature derived by Endo et al. (2022), but also was from data taken two years later when the dust cloud should have expanded and cooled and is past the time of active dust formation as the KLM light curves were dimming at the time of the SOFIA observations (see Fig. 3).
The SOFIA spectra were relatively short exposures and thus have limited signal-tonoise.However, there is a small rise at the long-wavelength end of the G063 spectrum and a general downward slope at the short-wavelength end of the G111 spectrum.We think that these spectral features are part of the 8 µm feature reported by Endo et al. (2022).This feature is often associated with a feature near 6-6.5 µm.We see a noisy feature that is reminiscent of the feature present in the dust of WR 137 that is discussed by Peatt et al. (2023).Unlike the SOFIA data for WR 137, we see a shorter wavelength of the feature indicating a higher hydrogen-content than that of WR 137.The WR wind lines observed in WR 137 early in its dust creation episode  2022) and a follow-up Subaru/COMIC spectrum presented by Endo (2022, PhD dissertation).Each grism setting from FORCAST is shown in a different color, and we also show the averaged flux of the G227 grating as a dark point with the standard deviation of all flux-calibrated points for this observation overplotted.In the second set of plots, we show the extinction-corrected flux of the SOFIA spectra for the values of A V discussed by Endo et al. (2022), with the range corresponding to the error in A V .(Peatt et al. 2023) are not easily seen in the spectrum of WR 125 likely due to the larger dust contribution for WR 125.After correcting the infrared SOFIA data for interstellar extinction (WD01 correction, A V = 5.89) and subtracting the stellar contribution as approximated by a powerlaw, we fit a blackbody distribution to the SOFIA data and SAI photometry including the emissivity of the dust and found the temperature of the dust to be ∼ 500 K as shown in the red dashed line.A warmer or cooler temperature can change the fit to better compare to the short-or long-wavelength ends of the distribution, but this is similar to the other fits we examined.
WR 125 is a member of the unusual class of dust-forming carbon-rich binaries.It has spectral properties similar to the nearby and well-studied prototype of these systems, WR 140.The O star in the WR 125 system is a later type and less luminous than in WR 140, but the spectrum does not show a disk-like geometry around the O star like in the WR 137 system (St-Louis et al. 2020).Given that the dense disk around the O star in WR 137 may help promote the dust production, WR 125 represents a longer period system with stellar parameters closer to that of WR 140.
Our derived orbit of WR 125 shows a moderate eccentricity compared to the extreme value for WR 140.WR 140 has an eccentricity of 0.8993±0.0013(Thomas et al. 2021).In contrast, the unusual WR 137 system, while also possessing a WC7 primary star has a very low eccentricity of 0.178 according to the spectroscopic orbit of Lefèvre et al. (2005), although the decretion disk around the Oe star likely adds to the dust production (Peatt et al. 2023) and the narrow geometry of the dust plume (Lau et al. 2023).We also note that a new visual orbit of WR 137 (Richardson et al. in prep) reports a larger eccentricity of 0.314±0.001and casts additional doubts on the O star radial velocities of Lefèvre et al. (2005).
Recently, WCd binary systems have received much attention due to their photogenic appearance in imaging from the recent JWST observations first presented by Lau et al. (2022).The imaging-based geometries can be reconciled with models for dust production based on orbital elements and the balance of momenta from the stellar winds.With the orbital kinematics for the WR star derived for WR 125, future observations of the dust from JWST will be able to create models for the dust production with fewer free parameters.
With the long period for WR 125, we can also consider the types of observations that could better pin down the orbit and evolution of the system.Despite the ∼ 8 yr period of WR 140, Thomas et al. (2021) presented a binary evolution model for the system that showed that the stars likely interacted to produce the modern-day WR star and stellar masses we observe.WR 125 has a much longer period and lower eccentricity, thus begging the question of its evolutionary status.If the stellar masses were similar to those of WR 140, then even with the larger distance to WR 125, the system would be resolvable with the CHARA Array, as Richardson et al. (2021) resolved the WN+O binary WR 133 with an apparent semi-major axis of less than one milliarcsecond.With a new higher-sensitivity beam combiner coming online (Lanthermann et al. 2022), the visual orbit of the fainter WR 125 may be possible to track in the coming decade.
Another context in which WR 125 resembled WR 140 was its non-thermal radio emission pointing to the presence of shock-accelerated electrons in its wind (Abbott et al. 1986).The fading of this (Williams et al. 1992), like that of WR 140 suggesting the burial of the non-thermal source in the wind (Williams et al. 1990) intensified the infrared monitoring of WR 125 leading to the discovery of the first dust-formation episode.While Midooka et al. (2019) presented a very limited X-ray data set, the similarities of WR 125 to WR 140 demands further analysis with the addition of more data.In particular, the X-ray light curve of WR 140 is explained through the radiative transfer based on hydrodynamical simulations based on the orbital geometry (Pollock et al. 2021).With additional data from X-ray satellites and the orbit we presented in Section 4, similar computations may allow for a complete understanding of the gas and its cooling and dust formation for WR 125.Every astrophysical laboratory in the population of dust-forming WC binaries allows us additional information into the physics of dust formation in the extreme environments surrounding the WCd binaries.
The infrared light curve shows a remarkable repeatability in its shape between the two observed outbursts, as seen in the phased light curves in Fig. 3. Our observations with SOFIA show a cooler dust cloud than seen earlier in the dust outburst with Subaru/COMICS (Endo et al. 2022).As our observations were taken after the nearinfrared peak in the light curve and the observations of Endo et al. (2022) were taken prior to the peak, this is not surprising and supports a dust cloud that forms and then expands and cools as it moves into the interstellar medium.
The contrast of the dusty spectral features such as the 8 µm and 6.2 µm features compared to the continuum and even in comparison to the WR emission lines appears to be less than the earlier observations taken in 2019.This shows that the continuum emission processes are dominant by the time the dust begins to cool in these systems.Likely, the changes in the 6.2 µm emission line seen in the early phases of the dust outburst of WR 137 (Peatt et al. 2023) when combined with the low contrast seen in these later phases of dust production for WR 125 give us a means to fully understand the chemistry of dust production in these binaries.Namely, we must be able to observe the spectral features at high signal-to-noise at several epochs early in the dust formation episodes in order to model the dust production and understand the chemistry of these systems.Given the potential importance of WC dust in our understanding of the modern-day and early-Universe dust budget (e.g.Marchenko & Moffat 2017;Lau et al. 2020), such observations and the associated laboratory measurements are crucial to understand the relative importance of WC dust.
Lastly, we can use the SOFIA observations and contemporaneous infrared photometry to estimate the mass of the dust produced in this binary system.From the G21 analysis shown in Fig. 7, we first note a few parts of this analysis.The fluxes at the lowest and highest wavelengths are seen to lie above the blackbody fit to the spectrum shown that has a temperature of ∼400 K.We divided the residual spectra shown in Fig. 6 by the emissivity for amorphous carbon grains (Zubko et al. 1996, ACAR sample).Using all the IR data and the dust emissivity, we can best fit the distribution with a temperature of 500 K, which is consistent with the temperature fitted to the multi-wavelength IR photometry of Williams et al. (1994) in 1993, at a similar orbital phase to the SOFIA observations.The dust temperature is a reasonable fit over 2-20 µm, but deviates towards the short and long wavelength regions of the range, indicating that the data are not well represented with a single dust temperature.As the dust has been moving away from the star for roughly three years since the beginning of the present dust-formation episode at this time and the stellar radiation heating it diluted, there will be a range of dust temperatures, with the oldest dust being coolest.Fitting multiple components to the flux distribution at this time is beyond the scope of this paper.
If we adopt an isothermal model for the dust with T = 500 K, a distance to the system of 5.88 kpc (Bailer-Jones et al. 2021), we get a dust mass of 1.5 × 10 −6 M ⊙ created, a bit higher than that inferred from the data in 1993 reported by Williams et al. (1994) but with a similar temperature.If the temperature is increased from 500 K to 550 K, then the long wavelength fluxes no longer fit the blackbody distribution, whereas if we decrease the temperature from 500 K to 450 K, the short wavelength fluxes are no longer fit, so we estimate an error of the temperature for an isothermal model to be ∼ 25 K.With the ∼3-year dust formation prior to our observation, we can then estimate the amount of carbon available for dust formation.For this, we can consider the mass-loss parameters of Sander et al. (2012), where we find that Ṁ = 3.05 × 10 −5 M ⊙ yr −1 and a carbon fraction of 0.4.We estimate that ∼10% of the WC wind enters the shock region based on an opening angle of 35-40 • dependent on the wind parameters of the two stars.From these considerations, we would see that the available carbon that had gone into the shock by the time of our SOFIA observation is about 3.5 × 10 −6 M ⊙ , with about a quarter of this being condensed into dust.As a comparison, if the dust temperature was 400 K, the dust mass becomes 1.1 × 10 −5 M ⊙ , which is higher than the amount of available carbon, further pointing to a need for multiple temperatures in the dust formation region years after periastron.While this seems high at first glance as WR 140 produces about 2 − 6 × 10 −8 M ⊙ each periastron (Williams et al. 2009), WR 140's dust production happens over a very small time window in comparison to that of WR 125.Furthermore, we note that planned observations in the coming years will solve this much better with high spatial resolution than with our spatially unresolved spectra presented here.

CONCLUSIONS
Our analysis of WR 125 has provided several key insights to the system.These include: • The infrared light curve has been fully consistent across two separate outbursts near the periastron passage of the system.The repeatability is reminiscent of the prototype of the dusty WC binaries: WR 140.A phase-dispersion minimization routine allowed us to measure the period of the binary to be 28.12 +0.10 −0.05 years.
• From archival spectroscopy spanning ∼20 years, we have measured a first orbit of WR 125.Only the WR star had measurable motion with our data.The eccentricity of the binary is 0.29±0.12based on an analysis with the periastron passage held constant between the rise of the NEOWISE-R mid-infrared flux and the peak of the K−band near-infrared flux.The O star kinematics was not able to be measured from our spectra, which did not cover the range of the higher Balmer absorption lines.
• The infrared spectrum from SOFIA shows a weak signature of a 6.2 µm UIR band, similar to that of WR 137.Unlike the recent analysis of WR 137, and similar to the results of Williams et al. (1994), no Wolf-Rayet emission lines are seen in the mid-infrared when the dust emission is near maximum for this system.Our analysis of the dust emission shows that it cooled from 800 K in 2019 to 580 K in 2021, likely consistent with the dust cloud expanding and cooling as it is radiatively driven away from the binary.Future studies detailing the expansion of the dust will allow this to be fully tested.
• Our data support the idea of a non-isothermal dust around WR 125, although the best fit temperature seems to be near 580 K, resulting in a dust mass of 9.5 × 10 −7 M ⊙ from the recent periastron.
The current outburst of WR 125 highlights the importance of long-term spectroscopic monitoring of such systems.The next few years will be especially crucial for WR 125, including the radial velocity maximum for the WR star.Likely, the orbit will be crucial to modeling efforts for any dust images that are made with infrared imaging with either ground-based telescopes or JWST.Long-term monitoring of WCd stars with ground-based photometry is instrumental in determining orbital periods in studies like these, and with a large sample of such systems, we could create an accurate assessment of the amount of dust created in such systems, which could be of cosmological importance.
Based in part on observations obtained at the international Gemini Observatory with programs GN-2019B-Q-410 and GN-2020A-Q-309 and then processed using the Gemini IRAF package, a program of NSF's NOIRLab, which is managed by the Association of Universities for Research in Astronomy (AURA) under a cooperative agreement with the National Science Foundation on behalf of the Gemini Observatory partnership: the National Science Foundation (United States), National Research Council (Canada), Agencia Nacional de Investigación y Desarrollo (Chile), Ministerio de Ciencia, Tecnología e Innovación (Argentina), Ministério da Ciência, Tecnologia, Inovações e Comunicações (Brazil), and Korea Astronomy and Space Science Institute (Republic of Korea).This work was enabled by observations made from the Gemini North telescope, located within the Maunakea Science Reserve and adjacent to the summit of Maunakea.We are grateful for the privilege of observing the Universe from a place that is unique in both its astronomical quality and its cultural significance.Some of the data presented herein were obtained at the W. M. Keck Observatory, which is operated as a scientific partnership among the California Institute of Technology, the University of California and the National Aeronautics and Space Administration.The Observatory was made possible by the generous financial support of the W. M. Keck Foundation.The authors wish to recognize and acknowledge the very significant cultural role and reverence that the summit of Maunakea has always had within the indigenous Hawaiian community.We are most fortunate to have the opportunity to conduct observations from this mountain.
Based in part on observations made with the NASA/DLR Stratospheric Observatory for Infrared Astronomy (SOFIA).

Figure 1 .
Figure1.Infrared photometry fromWilliams et al. (1992Williams et al. ( , 1994) ) as well as new measurements of WR 125.NEOWISE-R W 1 and W 2 after 2018 are too heavily saturated to use.We have also incorporated measurements from 2MASS.For each bandpass, we have shifted the points to the same quiescent level to highlight the outbursts and different amplitudes as a function of wavelength.The new observations are tabulated in a machine-readable format as data behind the figure.The vertical lines indicate the time of the infrared spectrum reported byEndo et al. (2022) in 2019, the 2020 spectrum reported by Endo (2022, PhD diss) and our SOFIA spectrum obtained in 2021.

Figure 2 .
Figure2.Results of the phase dispersion minimization for determination of the period based on the K-and L/L ′ -band photometry.The resulting period is 28.12 years.

Figure 3 .
Figure3.Phased infrared light curves based on the phase-dispersion minimization routines used in Section 3. We show photometry from both outbursts here, illustrating the repeatability of the photometry from one cycle to the next.As in Fig.1, the vertical lines represent the timing of the three mid-infrared spectra: fromEndo et al. (2022) in 2019,  Endo (PhD dissertaion, 2022)  in 2020, and then our SOFIA spectrum obtained in 2021.

Figure 4 .
Figure 4.An example spectrum of WR 125 in the region around C III λ5696 (reference for the velocity scale) from Gemini/GMOS on UT 2020 July 7.We show the measured blue and red measured velocities on the line wings with cross hairs on the central line indicating the measured velocity.The small absorption at velocity +9500 km s −1 is the He I λ5876 from the O star's photosphere, which we point out with a vertical arrow.

Figure 5 .
Figure 5.The single-lined spectroscopic orbit representing the motion of the WR star in the WR 125 binary.The data are represented with black circles for Keck/ESI velocities, blue diamonds for Gemini-N/GMOS, and red + signs for Keck/LRIS.

Figure 6 .
Figure 6.The flux-calibrated SOFIA data taken in 2021 April, compared to the spectrum presented by Endo et al. (2022) and a follow-up Subaru/COMIC spectrum presented byEndo (2022, PhD dissertation).Each grism setting from FORCAST is shown in a different color, and we also show the averaged flux of the G227 grating as a dark point with the standard deviation of all flux-calibrated points for this observation overplotted.In the second set of plots, we show the extinction-corrected flux of the SOFIA spectra for the values of A V discussed byEndo et al. (2022), with the range corresponding to the error in A V .
Figure7.After correcting the infrared SOFIA data for interstellar extinction (WD01 correction, A V = 5.89) and subtracting the stellar contribution as approximated by a powerlaw, we fit a blackbody distribution to the SOFIA data and SAI photometry including the emissivity of the dust and found the temperature of the dust to be ∼ 500 K as shown in the red dashed line.A warmer or cooler temperature can change the fit to better compare to the short-or long-wavelength ends of the distribution, but this is similar to the other fits we examined.
SOFIA was jointly operated by the Universities Space Research Association, Inc. (USRA), under NASA contract NNA17BF53C, and the Deutsches SOFIA Institut (DSI) under DLR contract 50 OK 2002 to the University of Stuttgart.Financial support for this work was provided in part by NASA through award #08-0220 issued by USRA.ARD is grateful to Embry-Riddle Aeronautical University's Undergraduate Research Institute for financial support through their IGNITE program and the NASA Space Grant program.NDR is grateful for support from the Cottrell Scholar Award #CS-CSA-2023-143 sponsored by the Research Corporation for Science Advancement as well as support from NASA through award #09-0163 issued by USRA.