# Radiation-driven acceleration in the expanding WR140 dust shell

Yinuo Han, Peter G. Tuthill, Ryan M. Lau, Anthony Soulain

► To cite this version:

Yinuo Han, Peter G. Tuthill, Ryan M. Lau, Anthony Soulain. Radiation-driven acceleration in the expanding WR140 dust shell. Nature, 2022, 610, pp.269-272. 10.1038/s41586-022-05155-5. insu-03860271

HAL Id: insu-03860271

https://insu.hal.science/insu-03860271v1

Submitted on 3 Oct 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

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

Distributed under a Creative Commons Attribution 4.0 International License
# Radiation-driven acceleration in the expanding WR140 dust shell

https://doi.org/10.1038/s41586-022-05155-5

Received: 15 December 2021

Accepted: 27 July 2022

Published online: 12 October 2022

Open access

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

Check for updates

Yinuo Han<fnref n="1" />,²✉, Peter G. Tuthill², Ryan M. Lau³,⁴ & Anthony Soulain²,⁵

The Wolf–Rayet (WR) binary system WR140 is a close (0.9–16.7 mas; ref. <fnref n="1" />) binary star consisting of an O5 primary and WC7 companion² and is known as the archetype of episodic dust-producing WRs. Dust in WR binaries is known to form in a confined stream originating from the collision of the two stellar winds, with orbital motion of the binary sculpting the large-scale dust structure into arcs as dust is swept radially outwards. It is understood that sensitive conditions required for dust production in WR140 are only met around periastron when the two stars are sufficiently close²⁻⁴. Here we present multiepoch imagery of the circumstellar dust shell of WR140. We constructed geometric models that closely trace the expansion of the intricately structured dust plume, showing that complex effects induced by orbital modulation may result in a ‘Goldilocks zone’ for dust production. We find that the expansion of the dust plume cannot be reproduced under the assumption of a simple uniform-speed outflow, finding instead the dust to be accelerating. This constitutes a direct kinematic record of dust motion under acceleration by radiation pressure and further highlights the complexity of the physical conditions in colliding-wind binaries.

The WR140 was observed on six occasions between 2001 and 2017 with the near-infrared cameras NIRC and NIRC2 at the Keck Observatory. The near-infrared imagery, displayed in Fig. 1, spans orbital phases between $\phi = 0.043$ and $0.592$, clearly showing an expanding dust shell with an evolving apparent morphology. Although the images span two orbital cycles, features appear to be consistent, implying a high degree of cycle-to-cycle replication of the same underlying morphology. Prominent structures at earlier orbital phases ($\le 0.111$) include an eastern and western dust arc. As the dust shell expands, these structures give rise to an ‘eastern arm’ and a prominent ‘southern bar’ at later orbital phases, following the nomenclature of structures identified in previous imaging⁴.

The well-resolved detailed form of the plume motivates the construction of a geometric model to explain the structural variation and expansion over time. We modelled the geometry of the dust plume as a linearly expanding spiral based on the ‘pinwheel mechanism’⁵,⁶. We assumed that dust production turns on and off episodically as the binary approaches, and departs the region near periastron. The wind-collision region can be approximated as the surface of a cone at large distances from the stars where the velocity of the compressed wind has reached its asymptotic value⁷. Dust assumed to form on this conical interface between the WR star and the O-star primary subsequently expands radially outwards as the binary continues in its orbit. Originally modelled on WR104 (ref. ⁵), such pinwheel models have been shown to accurately reproduce dust structures in several Wolf–Rayet (WR) binaries including Apep⁸ and WR112 (ref. ⁹).

The orbital parameters of WR140 are well constrained<fnref n="1" />,¹⁰ and form the basis for generation of the geometric dust plume model. By fitting to the location and geometry of dust structures across the multiple

epochs of observations, our model suggests that the half-opening angle of the conical shock front is $\theta_w = 40 \pm 5°$. This value is in close agreement with that estimated by Fahed et al.¹¹ (42 ± 3°), which is derived by fitting a cone model¹² to the 5696 A C III emission line, and Williams et al.¹³ (34 ± 1°) by modelling the 10830 He I subpeak. The $\theta_w$ value estimated in this study implies a momentum ratio between the O and WR star of $\eta = 0.043^{+0.021}_{-0.015}$ assuming radiative postshock conditions⁷,¹⁴, which is a larger value than that expected from the mass-loss rate and wind speed of the stars (Extended Data Table 2). Upon fitting to the turn-on and turn-off values independently, we find that dust production occurs over a period of $0.7^{+0.3}_{-0.1}$ yr centred at the periastron passage, with key parameters of the dust production thresholds displayed in Extended Data Table 3.

An image simulated under this geometric model at $\phi = 0.592$ is displayed in Fig. 2a. When comparing with the corresponding observations (Fig. 1), this model successfully reproduces a number of prominent features, with very accurate registration between the structural edges in the model and the data. The eastern arm represents a segment of the ellipse that corresponds to the earliest dust produced in the dust production episode (when dust production turned on), whereas the southern bar corresponds to the most recently produced dust in the episode (just before dust production turned off).

However, not all predicted features are apparent in the data. Although the geometric model is primarily designed to reproduce structural features, its physical interpretation implies a pinwheel system producing isothermal, optically thin dust at a constant rate. The clear absence of structures predicted by the geometric model, most noticeably dust features to the north and west (an outer-western arc) and below the southern bar, cannot be explained by simple density variations due

<footnote n="1">¹Institute of Astronomy, University of Cambridge, Cambridge, UK. ²Sydney Institute for Astronomy, School of Physics, The University of Sydney, Sydney, New South Wales, Australia. ³NSF’s NOIR Lab, Tucson, AZ, USA. ⁴Institute of Space & Astronautical Science, Japan Aerospace Exploration Agency, Sagamihara, Japan. ⁵Université Grenoble Alpes, CNRS, IPAG, Grenoble, France.</footnote>

<footnote>✉e-mail: yinuo.han@ast.cam.ac.uk</footnote>
<figure><img src="image_3.png" /><figcaption>**Fig. 1** | Near-infrared imagery of WR140's expanding circumstellar dust structure. The observing bands (top left), instrument (top right) and orbital phase $\phi$ (bottom right) are labelled in each panel. Coronagraphic (cor.) and aperture masking (mask) images are labelled in parentheses after the observing bands. All other images were observed with the full aperture. The scale bar in each panel indicates 0.3 arcsec (corresponding to 501 AU at a distance of 1.67 kpc). The images are stretched linearly by the following</figcaption></figure>

amounts to accentuate the dust structures that are otherwise faint relative to the bright stellar core: $\phi = 0.077$ Kp-band image stretched between 0 and 0.8 relative to the brightest pixel; $\phi = 0.111$ Kp-band image between 0 and 0.5; $\phi = 0.111$ H-band image between 0 and 0.8; $\phi = 0.592$ Lp-band image between 0 and $1.5 \times 10^{-3}$; $\phi = 0.592$ Ms-band image between $1 \times 10^{-4}$ and $5 \times 10^{-3}$. All other images are not stretched. An observing log is presented in Extended Data Table 1.

to orbital velocity, and deeper astrophysical insight into the plume generative model is required.

We found that two modifications to the geometric model enable it to account for the lack of the north-western outer dust arcs, both of which have physical motivation in the mechanics of the system. First, as dust nucleation and condensation sensitively depend on the physical conditions mediated by stellar winds³, a change in the shock structure will occur over the orbit, resulting in a dust production rate that varies continuously within the dust-producing segment of the orbit. We found that when the dust production rate is smoothly reduced to a local minimum at periastron in the model, we obtain a significantly improved fit that removes the western dust arc from the image, in accordance with observation.

Second, asymmetries introduced at the wind base can be amplified into asymmetric large-scale structures. Williams et al.⁴ found that an increased dust density at the trailing edge of the shock cone relative to the orbital direction produced a better fit than a uniform model. This asymmetry could be introduced by the fast orbital motion near periastron, which results in an intrinsic 'headwind' against the instantaneous

direction of orbit. By allowing dust to preferentially form on the trailing edge of the conical shock front and smoothly decrease in density towards the leading edge, we again obtained an improved fit with the relative brightness of the northern dust structures reduced. A schematic diagram illustrating the two effects is displayed in Fig. 3a.

Figure 2b,c displays model images under the two modifications to the original model (see also Extended Data Fig. 1). We find that a very close fit to the data is achieved when both the orbital and azimuthal modulation of dust production are combined (Fig. 2d-f). Overlaying the outline of the model on the corresponding $\phi = 0.592$ Lp-band observation shows that the geometry predicted by the model excellently reproduces the dust structures and their location (Fig. 2e).

Previous interpretations of dusty WR binaries have focused on identifying the upper limit on the binary separation that enables sufficiently high densities in the colliding winds to form dust. As gas compression is derived from wind-wind collision, it is natural to expect an upper distance threshold above which a sufficiently high density cannot be reached to facilitate dust nucleation. However, our model appears to
<figure><img src="image_4.png" /><figcaption>**Fig. 2 | Geometric model of WR140's circumstellar dust structure.**
**a–d.** Model images at $\phi = 0.592$ simulated under four different variants of the geometric model. The four models assume that the dust production rate is uniform (a), orbitally modulated (b), azimuthally asymmetric (c) and both orbitally modulated and azimuthally asymmetric (d). An animation showing the evolution of the model in time is available in Supplementary Video 1.</figcaption></figure>

**e–h.** The outline of the 'orbital + azimuthal' model overlaid on the observations at three different epochs, in which **e–g** show models at orbital phases of 0.592, 0.111 and 0.043 adjusted for dust acceleration, whereas **h** shows a model at an orbital phase of 0.043 that is not adjusted for acceleration, resulting in a clear misfit to the data. The scale bar in each panel indicates 0.3 arcsec.

suggest that there may also exist a lower limit, implying that dust is only formed within a Goldilocks zone where just the right conditions of density and temperature are met.

Observationally, several other WR binaries with infrared-imaged dust plumes, such as WR112 (ref. ⁹), WR98a¹⁵,¹⁶ and WR104 (refs. ⁵,⁶), appear to be continuous dust producers that do not show apparent signs of such a threshold effect. The Apep system appears to oscillate in and out of a single threshold, with dust production occurring near periastron⁸. If dust production in WR140 indeed reaches two maxima over a single orbit as suggested by our model, it would be the first WR binary known to exhibit the full range of this Goldilocks effect.

Previous work by Usov³ has shown that the degree of WR wind compression and cooling is sensitively dependent on the wind velocity, and so it is plausible that, at very short distances, the WR wind is significantly slowed by its binary companion via radiative braking, thereby reducing the rate of dust production. Lau et al.¹⁷ proposed that this mechanism may be responsible for impeding dust formation in the Gamma Velorum system. Alternatively, near periastron, the wind of the O star may not have accelerated to its terminal velocity before reaching the shock front.

With all other parameters fixed, the final parameter of interest in the model is the dust expansion speed. Dust grains are expected to form in a population of postshock gas where the two stellar winds of the binary collide. This dust may then be accelerated by stellar radiation pressure until reaching the terminal dust drift velocity. The assumption of a constant expansion speed similar to the stellar wind speed has been shown to accurately reproduce the expanding dust shell in similar systems such as WR112 (ref. ⁹).

However, the multiepoch observations of WR140 suggest that a uniformly expanding dust shell cannot simultaneously reproduce the time-evolving spatial extent of the dust shell observed in all epochs. Fitting to the $\phi = 0.592$ epoch, which shows the best-resolved dust structures, the model suggests a dust expansion speed of 2,400 ± 100 km s⁻¹,

which is broadly consistent with a streaming velocity of 2,170 ± 100 km s⁻¹ along the shock cone determined by Fahed et al.¹¹. However, extending this model to earlier epochs shows a clear misfit to the location of the dust shell. Figure 2h shows the outline of such a uniform expansion model overlaid on the earliest epoch of imagery at $\phi = 0.043$. The dust shell predicted by uniform expansion fails to fit, being significantly larger than that observed, implying that the average expansion speed up to $\phi = 0.043$ must be lower than in subsequent epochs.

Dropping the assumption of uniform expansion, the locations of the dust structures were fitted independently to each epoch. The resulting expansion speeds yielded a clear accelerating trend, with most of the impulse imparted onto the dust by the first two epochs. Given the uncertainties associated with fitting the location of dust features, direct derivation of the magnitude of acceleration as a function of distance from the star was not possible. To constrain the basic physics, we therefore posited a simple model based on expectations for radiatively accelerated dust as illustrated in Fig. 3b. Dust is not expected to form very close to the star owing to high temperatures and harsh ultraviolet radiation, and so the postshock gas is assumed to originate at a constant drift velocity, $v_0$. At a distance of $r_{nuc}$, dust nucleates and condenses to form an optically thick sheet, which experiences maximal radiation pressure with a constant acceleration, $a_{max}$. The dust continues to expand and eventually becomes optically thin at $r_t$, at which stage the acceleration decreases as $1/r^2$.

Under these working assumptions, we find that dust grains form at a speed of $v_0 = 1,810^{+140}_{-170}$ km s⁻¹, which is then accelerated in the optically thick regime by $a_{max} = 900^{+700}_{-400}$ km s⁻¹ yr⁻¹ up to $r_t = 220^{+150}_{-80}$ AU before becoming optically thin. The dust nucleation radius is fitted to be $r_{nuc} = 50 \pm 30$ AU, although this value is not well constrained.

Figure 3c shows the acceleration and velocity as a function of distance resulting from the best-fit model. Radiation pressure has long been suggested to play a major role in accelerating material near WR stars¹⁸. We find that the best-fit acceleration in the optically thin regime may be
<figure><img src="image_5.png" /><figcaption>**Fig. 3 |** *A model for dust production and acceleration in WR140. a, The 'orbital' and 'azimuthal' directions referred to in this paper relative to the orbit of a binary consisting of a WR star and an O or B star. A dust production rate that varies in the both the orbital and azimuthal directions is found to reproduce the observed structures. b, The acceleration regimes near the binary. No dust is produced interior to the dust nucleation radius, where material flows at a constant rate. Upon nucleation, dust forms an optically thick shell accelerated*</figcaption></figure>

reproduced by radiation-driven acceleration expected from the grain and stellar properties of the system (Methods). It is therefore plausible, if not expected, for radiation pressure to drive dust expansion at the acceleration derived from the observations.

The acceleration of matter by radiation pressure is ubiquitous throughout astrophysics, and is central to a range of astrophysical theatres ranging from the clearing of primordial material in proto-planetary disks around young stars¹⁹ to the formation of planetary nebulae towards the end of stellar evolution²⁰. These multiepoch observations of WR140 constitute the first direct detection of dust structures undergoing acceleration in the circumstellar radiation field: more commonly proper motions recorded are under the influence of gravity. Acceleration due to interactions other than gravity, however, are rarely witnessed because typically radiative acceleration zones are very short and terminal velocity rapidly attained. However, as we have discovered, WR binaries produce postshock material that reaches a significant distance from the photosphere before dust condensation, providing a rare laboratory in which the acceleration zone may be probed. This motivates further high angular resolution observations of such dust structures close to their nucleation point to test theories of radiatively driven acceleration physics in the circumstellar environment.

## Online content

Any methods, additional references, Nature Research reporting summaries, source data, extended data, supplementary information, acknowledgements, peer review information; details of author contributions and competing interests; and statements of data and code availability are available at https://doi.org/10.1038/s41586-022-05155-5.

at a uniform rate before the acceleration decreases as the dust becomes optically thin. c, The best-fit acceleration (dashed) and resulting dust expansion speed (solid) and the corresponding 1σ interval of models drawn from the posterior distribution fitted to the data. The theoretically expected acceleration via radiation pressure in the optically thin regime assuming a dust-to-gas mass ratio of 0.02 is overplotted (dotted).

1. Monnier, J. D. et al. First visual orbit for the prototypical colliding-wind binary WR 140. Astrophys. J. **742**, L1 (2011).

2. Williams, P. M. et al. Multi-frequency variations of the Wolf-Rayet system HD 193793 - I. Infrared, X-ray and radio observations. Mon. Not. R. Astron. Soc. **243**, 662–684 (1990).

3. Usov, V. V. Stellar wind collision and dust formation in long-period, heavily interacting Wolf-Rayet binaries. Mon. Not. R. Astron. Soc. **252**, 49 (1991).

4. Williams, P. M. et al. Orbitally modulated dust formation by the WC7+O5 colliding-wind binary WR140. Mon. Not. R. Astron. Soc. **395**, 1749–1767 (2009).

5. Tuthill, P. G., Monnier, J. D. & Danchi, W. C. A dusty pinwheel nebula around the massive star WR104. Nature **398**, 487–489 (1999).

6. Tuthill, P. G. et al. The prototype colliding-wind pinwheel WR 104. Astrophys. J. **675**, 698–710 (2008).

7. Cantó, J., Raga, A. C. & Wilkin, F. P. Exact, algebraic solutions of the thin-shell two-wind interaction problem. Astrophys. J. **469**, 729 (1996).

8. Han, Y. et al. The extreme colliding-wind system Apep: resolved imagery of the central binary and dust plume in the infrared. Mon. Not. R. Astron. Soc. **498**, 5604–5619 (2020).

9. Lau, R. M. et al. Resolving decades of periodic spirals from the Wolf-Rayet dust factory WR 112. Astrophys. J. **900**, 190 (2020).

10. Thomas, J. D. et al. The orbit and stellar masses of the archetype colliding-wind binary WR 140. Mon. Not. R. Astron. Soc. **504**, 5221–5230 (2021).

11. Fahed, R. et al. Spectroscopy of the archetype colliding-wind binary WR 140 during the 2009 January periastron passage. Mon. Not. R. Astron. Soc. **418**, 2–13 (2011).

12. Lührs, S. A colliding-wind model for the Wolf-Rayet system HD 152270. Publ. Astron. Soc. Pacific **109**, 504–513 (1997).

13. Williams, P. M. et al. Conditions in the WR 140 wind-collision region revealed by the 1.083-µm H I line profile. Mon. Not. R. Astron. Soc. **503**, 643–659 (2021).

14. Gayley, K. G. Asymptotic opening angles for colliding-wind bow shocks: the characteristic-angle approximation. Astrophys. J. **703**, 89–95 (2009).

15. Monnier, J. D., Tuthill, P. G. & Danchi, W. C. Pinwheel nebula around WR 98A. Astrophys. J. **525**, L97–L100 (1999).

16. Monnier, J. D., Tuthill, P. G., Danchi, W. C., Murphy, N. & Harries, T. J. The Keck aperture-masking experiment: near-infrared sizes of dusty Wolf-Rayet Stars. Astrophys. J. **655**, 1033–1045 (2007).

17. Lau, R. M. et al. Revisiting the impact of dust production from carbon-rich Wolf-Rayet binaries. Astrophys. J. **898**, 74 (2020).

18. Beals, C. S. On the nature of Wolf-Rayet emission. Mon. Not. R. Astron. Soc. **90**, 202–212 (1929).

19. Hughes, A. M., Duchêne, G. & Matthews, B. C. Debris disks: structure, composition and variability. Annu. Rev. Astron. Astrophys. **56**, 541–591 (2018).

20. Peimbert, M. Planetary nebulae. Rep. Prog. Phys. **53**, 1559–1619 (1990).

**Publisher's note** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

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

**Open Access** This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/.

© The Author(s) 2022
# Methods

## Image reduction

Imagery from Keck's NIRC2²¹ were taken using sparse-aperture interferometry techniques as described in Tuthill et al.²², with an observing log describing specific configurations given in Extended Data Table 1. After data reduction of sequences of short-exposure data cubes, interleaved with observations of a point-source reference star, our codes deliver calibrated standard interferometric observables (visibilities and closure phases). Image recovery was performed with the BSMEM algorithm²³ at a pixel scale of 10 mas pixel⁻¹. The BSMEM algorithm reconstructs a model-independent image iteratively, taking into account both the difference between the model and data (χ²) to assess the goodness of fit, and the entropy function (H) to minimize information and avoid overfitting as the normalized χ² converges towards unity. The balance between the two terms during the iterative minimization is set by the regularization weight (a), which is adaptively determined by the algorithm. The value of the regularization weight is chosen such that it maximizes the Bayesian evidence (that is, the probability of obtaining the data given the regularization weight) given a Jeffreys prior²⁴.

Data analysis for the NIRC2 instrument, also operational at the Keck observatory, proceeded by centring and stacking data frames common to the same epoch, observing band and aperture configuration (full pupil or coronagraph) and interpolating outlying pixels in the stacked image. For the $\phi = 0.592$ images that showed the faintest dust structure, we subtracted the sky background determined using the dithered data frames.

## Geometric model

**Original model.** Our model assumes that dust is produced on a thin conical shock front, after which it expands radially away from the binary, generating a spiral form as a consequence of rotation of the nose of the shock with the orientation of the orbiting binary. Noting WR140's well-constrained orbital parameters, we adopted values from Monnier et al.¹ presented in Extended Data Table 2, which were derived using a combination of knowledge from radial velocity¹¹ and optical interferometry. The model is available at github.com/yinuohan/WR140 (ref. ⁸).

Parameters additional to the orbit that are required to define the plume model are the half-opening angle of the conical wind-shock interface ($\theta_w$), the terminal dust velocity and the range of orbital phases over which dust production occurs. We constrained these parameters by fitting the geometry generated by the model to the observed dust structures across the multiepoch images. The images generated by the geometric model assume isothermal dust. This is a reasonable assumption given that the temperature gradient in the orbital direction is expected to be small; however, a more exact model of the emission properties would require further radiative transfer simulations. As the primary objective is to reproduce the geometric structure of the dust plume, the features of interest are the shape and location of dust structures. These properties are difficult to programmatically extract from the observations, thereby hindering the construction of a metric and subsequently using an automated sampling algorithm to explore the parameter space, as found by Han et al.⁸ in the case of the Apep system. We therefore fitted the structures manually with an understanding of the physical effects of parameters and visual inspection of the model overlaid on the observations.

**Orbital modulation.** Even though the abrupt turn-on and turn-off are able to closely reproduce corresponding features such as the southern bar and the southern segment of the eastern arm, the orbital modulation of dust production rate once dust production begins appears to require a gradual variation, as suggested by the gradual dimming of the eastern arm towards the north. We modelled the variation of the dust production rate over the orbit as the superposition of two distinct peaks, each of which is modelled as a half-Gaussian function varying

over the true anomaly. The turning on and off of dust production corresponds to the peaks of the two half-Gaussians, with dust production gradually suppressed near periastron as shown in Extended Data Fig. 2. We find that a reasonable fit is obtained when the half-Gaussians peak at the same turn-on and turn-off orbital phase as previously fitted, each with a standard deviation of $40_{-10}^{+30}°$ in true anomaly. This results in an interval between half powers of $0.15_{-0.09}^{+0.05}$ yr along the orbit.

**Azimuthal asymmetry.** We assumed that the density of dust being produced at the conical shock front varies over the azimuthal angle about the WR–OB axis as a Gaussian function (with periodic boundary condition). Fixing the peak dust density at the trailing side of the shock cone (where the shock cone intersects the orbital plane), the model was best fitted when the dust density distribution exhibited a standard deviation of $80 \pm 20°$ in the azimuthal angle measured from the trailing edge. This results in a density contrast of $13_{-8}^{+77}$ between the trailing and leading sides of the dust plume.

## Expansion speed fitting

The image of the spiral dust plume is a projection of its true three-dimensional structure, and so measuring the true displacement in three-dimensional space relies on the geometric model.

Different dust features in the same observation are produced at different points in time. The effect of this time offset is particularly important towards earlier orbital phases. Although the construction of the geometric model assumes uniform expansion across the entire dust plume, it is still possible to measure the displacement of specific features within an observation by fitting only to the location of that feature. We therefore measured the displacement of the eastern arm and southern bar independently by fitting to their location separately for each epoch of observation, thereby dropping the assumption of uniform expansion speed. Extended Data Fig. 3 displays a series of model outlines overplotted on the corresponding data for models with and without accounting for acceleration. A misfit between the model-predicted location of features and their observed location in the absence of acceleration, which is particularly prominent at earlier epochs, implies that the inclusion of acceleration is required to explain their motion.

The geometric model predicts the range of orbital phases over which dust is produced, and so it is possible to predict the time when the two dust features were produced for each epoch. Here we assumed that the eastern arm and southern bar are produced at the very beginning and end of a dust production episode, respectively. Subtracting this time point from the time of observation gives the total time since a given feature was produced, giving two sets of displacements as a function of time, one for the eastern arm and one for the southern bar. We then combined the two sets of measurements so that they share the same time axis and used the combined data and a point at the origin for subsequent fitting.

We fitted the model for acceleration driven by radiation pressure described in the main text using a Markov chain Monte Carlo, solving for the differential equations

$$a(r) = \frac{d^2 r}{dt^2} = \begin{cases} 0 & \text{if } 0 \le r < r_{\text{nuc}} \\ a_{\text{max}} & \text{if } r_{\text{nuc}} \le r < r_{\text{t}} \\ \frac{a_{\text{max}} r_{\text{t}}^2}{r^2} & \text{if } r_{\text{t}} \le r \end{cases} \quad (1)$$

at each iteration to obtain $r(t)$, which is then compared with the data. The reported fitted results are derived from the median of the marginalized posterior distributions, with uncertainties estimated using the 16th and 84th percentiles. Extended Data Fig. 4 compares the measured location of dust with the best-fit model to these measurements with and without accounting for acceleration. The systematic offset
in residuals for the model, assuming that the radial distance of dust is proportional to time since production, suggests that uniform expansion is not a suitable model.

We note that acceleration fitting relies on knowledge of the time of production of the dust features being fitted to, which in turn is derived from the orbital phase of the earliest and most recently produced dust, which carries uncertainties. Even though the displacement-time measurements at most orbital phases appear to be approximately linear, and even if we were to drop the assumption that displacement is strictly proportional to time in a model without acceleration, a best-fit straight line fitted through these points does not intercept the origin, implying the presence of acceleration at least at very early stages of the expansion. Such an offset from the origin arises from the combined measurements of both the eastern arm and southern bar, and any systematic errors resulting in such an offset across both sets of measurements is unlikely.

## Acceleration from radiation pressure

To estimate the expected effect of radiation pressure, we assumed stellar and grain properties adopted by the dust grain model by Lau et al. (manuscript in preparation), which includes a combined stellar luminosity of $L_{\text{bin}} = 1.43 \times 10^6 L_{\odot}$, a grain size of $a_g = 0.04 \mu\text{m}$, a grain density of $\rho_g = 1.6 \text{ g cm}^{-3}$ (ref. 25) and a radiation-pressure efficiency of $Q_{\text{pr}} = 1$ as appropriate for hot stellar sources given the relatively large grain size$^{26}$. We calculated the radiation force experienced by an individual dust grain as

$$F_{\text{rad}} = Q_{\text{pr}} \sigma_g \frac{L_{\text{bin}}}{4\pi r^2 c} \quad (2)$$

where $\sigma_g$ is the cross-sectional area of a dust grain, $r$ is the distance of the dust grain from the star and $c$ is the speed of light. Assuming that any gas is comoving with the dust as it accelerates, we found that a dust-to-gas ratio of 0.019 reproduces the best-fit acceleration model derived from the observations.

Lau et al.$^{17}$ estimated a mean dust-to-gas ratio of approximately $4 \times 10^{-5}$ in the system by comparing the total dust mass derived from observations with the stellar mass-loss rate; however, the larger dust-to-gas ratio found in this study is not surprising given that dust production exhibits significant spatial and temporal variation. As the geometric model suggests, dust is only produced over a fraction of the orbit, and further orbital modulation implies that the dust features (the eastern arm; the southern bar) probed in this study correspond to only a very small fraction of the orbit. Furthermore, dust production occurs only in a confined stream. The two effects together imply that dust in the system is expected to be highly concentrated in the two structures probed in this study, which may account for the more than two orders of magnitude increase in dust density in these structures relative to the mean density, assuming a uniform dust distribution across the system.

Additional contributions to the acceleration could theoretically include further wind collision downstream, which is not modelled here. However, the two stellar winds rapidly become parallel farther away from the apex of the shock front and so this effect is not expected to be a major contribution given the geometry.

The interval of constant acceleration ($a_{\text{max}}$) models a scenario in which the dust is initially optically thick, harvesting all photon momentum available over the area that it spans. The extent to which this assumption may hold true motivates an important point of observational follow-up, such as by monitoring variations in the spectral energy distribution of newly formed dust over time.

## Instrumental and optical depth effects

Critical constraints on the accelerating motion of the dust plume are provided by data at early orbital phases, which were observed with aperture masking interferometry. It is therefore important to ensure that the plate scale of images recovered with this technique is not adversely

affected by systematic biases, for example due to different observing bandpasses. Aperture masking observations of reference binary stars of known configuration spanning a diverse range of filters and bandwidths at each epoch (Alp UMi with CH$_4$ and Beta Del with KCont in June 2001; HD 2150123 with H, K in July 2001; Zeta Aqr with H, K, Beta Del with KCont and HD 2150123 with H, K in July 2002) all recovered separations within errors of anticipated locations, which suggest that bandwidth effects do not adversely affect proper motion tracking as performed in this study.

Furthermore, the same apparatus and pipelines have also been applied as part of a relatively well-tested and mature experiment that has been benchmarked, delivering numerous multiepoch astrophysical studies of similar resolved structures in other contexts. Notable science targets whose structures were registered and tracked over these observing epochs include IRC+10216 (ref. 27) and LkHα 101 (refs. 28,29). In particular, IRC+10216 had outflow motions tracked by tagging structures in recovered images over many epochs with both broad and narrow filters, with the resulting motions consistent with a physically plausible, monotonically expanding flow, and no biases introduced with bandwidth or other settings$^{27}$. In the WR140 dataset presented in this study, the close agreement between features taken with different filters at the same epoch also argue against plate scale biases being a significant effect.

It is also possible that optical depth effects can alter the apparent spatial extent of structures, as inner dust could potentially shadow outer dust at earlier orbital phases if the dust is optically thick, making the structure appear smaller. However, this is unlikely to alter the observed spatial extent by an amount that could explain the apparent acceleration. If shadowing were to be a prominent effect, the thickness of the dust plume required would result in a significantly different structure at earlier orbital phases that resemble a much smaller half-opening angle (approximately 10° smaller than the best-fit value), and at later orbital phases this would produce thick and fuzzy structures that deviate from the fine structures observed. These observations therefore imply that flux gradients and/or shadowing within a thick dust plume are not large enough to bias our kinematic analysis of WR140.

## Data availability

The NIRC and NIRC2 data underlying this article are publicly available on the Keck Observatory Archive (koa.ipac.caltech.edu/cgi-bin/KOA/nph-KOalogin) under programme IDs U59N, U2N, U45N, U14N2, U044 and Z273. Source data are provided with this paper.

## Code availability

The dust plume model used in this study is available at github.com/yinuohan/WR140.

21. Matthews, K. & Soifer B. T. The near infrared camera on the W.M. Keck telescope. *Exp. Astron.* **3**, 77–84 (1994).

22. Tuthill, P. G., Monnier, J. D., Danchi, W. C., Wishnow, E. H. & Haniff, C. A. Michelson interferometry with the Keck I telescope. *Publ. Astron. Soc. Pacific* **112**, 555–565 (2000).

23. Buscher, D. F. in *Very High Angular Resolution Imaging* (eds Robertson, J. G. & Tango, W. J.) 91–93 (Kluwer, 1994).

24. Baron, F. & Young, J. S. in *Optical and Infrared Interferometry Vol. 7013* of Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series (eds Schöller, M. et al.) 70133X (SPIE, 2008).

25. Jones, A. P. et al. The evolution of amorphous hydrocarbons in the ISM: dust modelling from a new vantage point. *Astron. Astrophys.* **558**, A62 (2013).

26. Hoang, T. Effect of dust rotational disruption by radiative torques on radiation pressure feedback from massive protostars. *Astrophys. J.* **921**, 21 (2021).

27. Tuthill, P. G., Monnier, J. D. & Danchi, W. C. Multiwavelength diffraction-limited imaging of the evolved carbon star IRC+10216. II. *Astrophys. J.* **624**, 352–358 (2005).

28. Tuthill, P. G., Monnier, J. D. & Danchi, W. C. A dusty torus around the luminous young star LkHα 101. *Nature* **409**, 1012–1014 (2001).

29. Tuthill, P. G., Monnier, J. D., Danchi, W. C., Hale, D. D. S. & Townes, C. H. Imaging the disk around the luminous young star LkHα 101 with infrared interferometry. *Astrophys. J.* **577**, 826–838 (2002).

30. Foreman-Mackey, D., Hogg, D. W., Lang, D. & Goodman, J. emcee: the Markov chain Monte Carlo hammer. *Publ. Astron. Soc. Pacific* **125**, 306 (2013).
31. Van Der Walt, S., Colbert, S. C. & Varoquaux, G. The numpy array: a structure for efficient numerical computation. *Comput. Sci. Eng.* **13**, 22–30 (2011).

32. Hunter, J. D. Matplotlib: A 2D graphics environment. *Comput. Sci. Eng.* **9**, 90–95 (2007).

33. Astropy Collaboration et al. Astropy: a community Python package for astronomy. *Astron. Astrophys.* **558**, A33 (2013).

34. Pollock, A. M. T. et al. Competitive X-Ray and optical cooling in the collisionless shocks of WR 140. *Astrophys. J.* **923**, 191 (2021).

35. Setia Gunawan, D. Y. A. et al. Multi-frequency variations of the Wolf–Rayet system HD193793 (WC7pd+O4-5) III. IUE observations. *Astron. Astrophys.* **376**, 460–475 (2001).

**Acknowledgements** Y.H., P.G.T. and A.S. acknowledge the traditional owners of the land, the Gadigal people of the Eora Nation, on which the University of Sydney is built and this work was carried out. Y.H. acknowledges funding from a Gates Cambridge Scholarship enabled by Grant No. OPP1144 from the Bill and Melinda Gates Foundation. P.G.T. acknowledges the longstanding support for innovative imaging at the W.M. Keck observatory, and development of the techniques particularly from J. Monnier. 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 financial support of the W. M. Keck Foundation. The authors recognize and acknowledge the significant cultural role and reverence that the summit of Maunakea has always had within the indigenous Hawaiian

community. This research made use of NASA's Astrophysics Data System; the emcee package<sup>30</sup>; NUMPY<sup>31</sup>; MATPLOTLIB<sup>32</sup>; and Astropy, a community-developed core Python package for Astronomy<sup>33</sup>.

**Author contributions** P.G.T. prepared and performed several sets of observations and reduced NIRC data. Y.H. and P.G.T. devised and fitted the geometric model. Y.H. reduced NIRC2 data, fitted the plume motion and, with input from P.G.T., R.M.L. and A.S., wrote the manuscript. R.M.L. analysed dust properties and interpreted the interaction between dust and stellar winds. A.S. analysed possible optical thickness variations for different grain sizes. All authors contributed to analysing the data and synthesizing interpretations of the physics.

**Competing interests** The authors declare no competing interests.

**Additional information**

**Supplementary information** The online version contains supplementary material available at https://doi.org/10.1038/s41586-022-05155-5.

**Correspondence and requests for materials** should be addressed to Yinuo Han.

**Peer review information** Nature thanks Peredur Williams and the other, anonymous, reviewer(s) for their contribution to the peer review of this work. Peer reviewer reports are available.

**Reprints and permissions information** is available at http://www.nature.com/reprints.
<figure><img src="image_7.png" /><figcaption>**Extended Data Fig. 1 | Face-on view of the geometric model for WR140.**</figcaption></figure>

The four panels correspond to model images at $\phi = 0.592$ simulated under the assumption that the dust production rate is uniform (original), orbitally

modulated (orbital), azimuthally asymmetric (azimuthal) and both orbitally modulated and azimuthally asymmetric (orbital + azimuthal), respectively.
a

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

**Extended Data Fig. 2** | Orbital and azimuthal variation of dust production rate in the geometric model. (a) Dust production rate as a function of orbital phase near periastron. Dust is not produced at all other orbital phases. (b) Dust

b

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

production rate as a function of azimuthal angle (about the WR–O star axis) relative to the trailing edge. The $1\sigma$ uncertainty regions in both plots are shaded.
<figure><img src="image_10.png" /><figcaption>Extended Data Fig. 3 | Comparison between the model-predicted and observed location of the Eastern Arm for models with and without acceleration. All observations in which the Eastern Arm is visible are reproduced twice, one of which is overplotted with the outline of the model</figcaption></figure>

with acceleration and fitted to the location of the Eastern Arm (left panel of each pair) and the other without acceleration (right panel of each pair). The model without acceleration is fitted to the final epoch of observations. Lack of acceleration results in a misfit at earlier epochs.
<figure><img src="image_11.png" /></figure>

**Extended Data Fig. 4 | Comparison between observed and model-predicted location of dust features.** (a) The data points correspond to the location of the Eastern Arm and Southern Bar obtained by fitting the geometric model to these features at epochs where the features are visible. The error bars correspond to 1σ uncertainties, which were estimated in an Markov chain Monte Carlo fit under the model with acceleration assuming independently

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

and identically distributed Gaussian errors. The solid lines show 100 models randomly drawn from the posterior distribution. The dotted line shows the best-fit model assuming instead that the radial distance of dust features is proportional to time (that is, no acceleration). (b) residuals of the models with and without acceleration. The model without acceleration resulted in a systematic bias in the residuals.
Extended Data Table 1 | A log of NIRC and NIRC2 observations of WR140 presented in this study

<table><thead><tr><th>Instrument</th><th>Filter</th><th>Obs Mode</th><th>Obs date</th><th>Orbital phase</th></tr></thead><tbody><tr><td rowspan="4">NIRC</td><td>CH<sub>4</sub>, K</td><td>Masking</td><td>12 Jun 2001</td><td>0.043</td></tr><tr><td>K</td><td>Masking</td><td>30 Jul 2001</td><td>0.059</td></tr><tr><td>PAHCS</td><td>Full pupil</td><td>30 Jul 2001</td><td>0.059</td></tr><tr><td>PAHCS</td><td>Full pupil</td><td>24 Jul 2002</td><td>0.183</td></tr><tr><td rowspan="4">NIRC2</td><td>Lp, Ms</td><td>Full pupil</td><td>22 Oct 2005</td><td>0.592</td></tr><tr><td>Kp, Lp, Ms</td><td>Full pupil</td><td>31 Jul 2017</td><td>0.077</td></tr><tr><td>H, Kp, Lp</td><td>Full pupil</td><td>5 Nov 2017</td><td>0.111</td></tr><tr><td>Kp</td><td>Coronagraph200</td><td>5 Nov 2017</td><td>0.111</td></tr></tbody></table>

Filters [name / central wavelength / fractional bandpass] for the NIRC camera: [CH<sub>4</sub> / 2.269 µm / 6.8%]; [K / 2.2135 µm / 19%]; [PAHCS / 3.0825 µm / 3.3%]; and for the NIRC2 camera: [H / 1.633 µm / 18%]; [Kp / 2.124 µm / 16.5%]; [Lp / 3.776 µm / 18.5%]; [Ms / 4.67 µm / 5.2%]. The "Masking" observing mode corresponds to a partially redundant aperture mask as described in Tuthill et al.<sup>22</sup>. "Full pupil" indicates normal imagery, while one epoch also employed an occulting "Coronagraph200" spot. Orbital phases were calculated using orbital parameters derived by Monnier et al.<sup>1</sup>.
Extended Data Table 2 | Properties of the WR140 system

<table><thead><tr><th></th><th>Value</th><th>Reference</th></tr></thead><tbody><tr><td>Orbital period (P)</td><td>2896.35 ± 0.20 days</td><td>(1)</td></tr><tr><td>Time of periastron (T<sub>0</sub>)</td><td>46154.8 ± 0.8 MJD</td><td>(1)</td></tr><tr><td>Eccentricity (e)</td><td>0.8964<sup>+0.0004</sup><sub>-0.0007</sub></td><td>(1)</td></tr><tr><td>Inclination (i)</td><td>119.6 ± 0.5°</td><td>(1)</td></tr><tr><td>Longitude of ascending node (Ω)</td><td>53.6 ± 0.4°</td><td>(1)</td></tr><tr><td>Argument of periastron (ω)</td><td>46.8 ± 0.4°</td><td>(1)</td></tr><tr><td>Semi-major axis (a)</td><td>8.82 ± 0.05 mas</td><td>(1)</td></tr><tr><td>Distance (d)</td><td>1.67 ± 0.03 kpc</td><td>(1)</td></tr><tr><td>WR star mass loss rate (M&#775;<sub>WR</sub>)</td><td>1.7 × 10<sup>-5</sup> M<sub>☉</sub>/yr</td><td>(2)</td></tr><tr><td>O star mass loss rate (M&#775;<sub>O</sub>)</td><td>3.7 × 10<sup>-7</sup> M<sub>☉</sub>/yr</td><td>(2)</td></tr><tr><td>WR star terminal wind speed (v<sub>∞,WR</sub>)</td><td>2800 km/s</td><td>(3)</td></tr><tr><td>O star terminal wind speed (v<sub>∞,O</sub>)</td><td>3100 km/s</td><td>(3)</td></tr><tr><td>O/WR momentum ratio</td><td>0.020</td><td>M&#775;<sub>O</sub>v<sub>∞,O</sub>/(M&#775;<sub>WR</sub>v<sub>∞,WR</sub>)</td></tr></tbody></table>

References: (1) Monnier et al.¹, (2) Pollock et al.³⁴ and (3) Setia Gunawan et al.³⁵.
<table><caption>Extended Data Table 3 | Dust production thresholds in the original model</caption><thead><tr><th></th><th>Turn-on</th><th>Turn-off</th></tr></thead><tbody><tr><td>True anomaly (°)</td><td>-135<sup>+5</sup><sub>-5</sub></td><td>135<sup>+10</sup><sub>-5</sub></td></tr><tr><td>Separation (AU)</td><td>7.9<sup>+1.3</sup><sub>-1.1</sub></td><td>7.9<sup>+3.0</sup><sub>-1.1</sub></td></tr><tr><td>Orbital phase</td><td>-0.041<sup>+0.008</sup><sub>-0.011</sub></td><td>0.041<sup>+0.025</sup><sub>-0.008</sub></td></tr><tr><td>Duration of dust production (yr)</td><td>0.7<sup>+0.3</sup><sub>-0.1</sub></td><td></td></tr><tr><td>Separation at closest approach (AU)</td><td>1.53<sup>+0.04</sup><sub>-0.05</sub></td><td></td></tr></tbody></table>