© The Authors 2025

1 Introduction

Young massive stellar clusters (YMCs) provide ideal environments to study the evolution of massive stars, planet and star formation processes, and stellar feedback to the Galaxy (Portegies Zwart et al. 2010). In recent times, they have also been proposed as sources of cosmic rays due to detections of gamma rays (e.g. H.E.S.S. Collaboration 2022). Westerlund 1 is thought to be the most massive YMC in the Galaxy, with mass estimates of the order of 5 × 10⁴ to 10⁵ M_⊙ (Clark et al. 2005; Andersen et al. 2017), although the inferred mass is dependent on the derived luminosity (and hence the age) of the massive stars in the cluster, which is quite uncertain because of the large extinction to the cluster (Beasor et al. 2021). Westerlund 1 hosts a large and diverse population of massive stars at various evolutionary stages, with 24 Wolf–Rayet (WR) stars (Crowther et al. 2006) as well as OB supergiants, yellow super-giants, red supergiants (RSGs), luminous blue variables, and an sgB[e] star (Clark et al. 2020). The WR stars drive a cluster wind of several thousand km s⁻¹ (Härer et al. 2023) that sculpts a superbubble cavity around the cluster. The region around Westerlund 1 is a bright gamma-ray emitter to very high energies (H.E.S.S. Collaboration 2022), likely due to the combined action of stellar winds and supernovae in accelerating cosmic rays (e.g. Vieu & Reville 2023; Härer et al. 2023).

The concept of mass loading of a fast wind by hydrodynamic (HD) mixing with cold and slow material was introduced by Hartquist et al. (1986), who applied it to the expansion of WR nebulae and the consequent X-ray and nebular emission that should arise according to their model. Mass loading in stellar clusters was considered theoretically by Stevens & Hartwell (2003) in the context of protostars and protoplanetary disks and their effects on diffuse X-ray emission. The effects of stellar winds and supernovae on molecular material remaining from the formation of a YMC were modelled by Rogers & Pittard (2013), who found the cluster wind to be mass loaded by a factor of up to several hundred. Arthur (2012) modelled the X-ray emission from the mass loading of protoplanetary disks in Orion. Mass loading of the cluster wind by the dense and slow winds of cool supergiants reduces the thermal pressure driving the cluster wind (Silich et al. 2011), potentially limiting its effectiveness in creating a superbubble and accelerating cosmic rays.

Observations of diffuse thermal X-ray emission in Westerland 1 show a hard component (Kavanagh et al. 2011; Haubner et al. 2025). While pre-main-sequence stars are thought to contribute to this (Clark et al. 2008), the majority is attributed to thermalisation of the winds from WR stars in the cluster core. These winds must mix with cooler material to explain the gas temperature derived from X-ray observations, but the origin of this material is not clear. Dust and ejecta from supernovae are thought to not contribute significantly to this process (Kavanagh et al. 2011; Haubner et al. 2025) since the supernova rate of 1 per ~10 kyr would imply the cluster is free of supernova remnant material for ~90% of the time (Muno et al. 2006); however, we note that this estimated supernova rate is very uncertain. RSG winds inject cool material that can be ionised by an external radiation field, such as that found in a YMC (Mackey et al. 2014). This has already been demonstrated for W26, one of the four RSGs in Westerlund 1 (Wright et al. 2014; Mackey et al. 2015). Furthermore, the potential for the enriched RSG wind material to contribute to future star formation within the cluster depends on whether it is heated and ejected rapidly, or can remain confined, cool, and dense (Mackey et al. 2015).

Recent JWST images of Westerlund 1 from the Extended Westerlund 1 and 2 Open Clusters Survey (EWOCS) show resolved shells and outflows from the four RSGs and other evolved cool massive stars in the cluster (Guarcello et al. 2025). The emission mechanism for gas and dust in Westerlund 1 is not clear, but its proximity to the core of such a massive cluster implies that the material was recently ejected from the cluster member stars. If the mid-IR emission is produced by dust, then the most likely source is winds and mass ejections from the cool supergiants and hypergiants in the cluster core.

In this paper, we explore, in a sense, the opposite case to that of Rogers & Pittard (2013), that of a cool RSG wind embedded in a hot, fast collective cluster wind from a YMC. We present a 3D HD simulation of our scenario to explore the efficiency of mass loading from evolved stellar winds, and use it to create synthetic dust-emission maps. There is some similarity between our work and previous simulations of the G2 cloud orbiting the Galactic Centre: De Colle et al. (2014) and Ballone et al. (2018) modelled a cool star emitting a slow wind and moving rapidly through the hot coronal gas around the supermassive black hole Sgr A^*. To our knowledge, no multi-dimensional simulations have been made of the interaction of stellar wind from a RSG embedded in a hot star-cluster wind, though such a situation must be common in YMCs.

Our paper is organised as follows: in Sect. 2 we discuss our simulation setup and the values adopted for the RSG and cluster winds. In Sect. 3, we present the results of our simulation, including slices through the XY domain (Sect. 3.1), phase diagrams for the RSG wind and cluster wind material (Sect. 3.2), mass loading of the cluster wind over time (Sect. 3.3), synthetic dust continuum-emission maps for 11 µm and 24 µm derived from our simulation (Sect. 3.4), and how those maps correspond to observations from JWST (Sect. 3.5). In Sect. 4 we discuss the uncertainties in the winds of RSGs and clusters (Sect. 4.1) and in our modelling (Sect. 4.2), as well as the implications for this work in understanding the role of mass loading (Sect. 4.3), its connection to super star clusters (Sect. 4.4), and the need for follow-up observations (Sect. 4.5). We present our conclusions in Sect. 5, and a resolution study in Appendix A.

2 Methods

2.1 HD simulations

We performed our 3D Cartesian simulations using the (magneto-) HD code PION (Mackey et al. 2021). We employed the static mesh refinement option in PION to achieve high resolution only where required. For our radiative heating and cooling prescription, we used model 8 in PION as described in Green et al. (2019), which assumes a sufficiently intense extreme-ultraviolet (EUV) radiation field to keep H and He photoionised everywhere on the domain, and which includes three cooling sources. The first source is the maximum of the Solar metallicity collisional ionisation equilibrium curve described by Wiersma et al. (2009), and the cooling function for forbidden lines described by Henney et al. (2009). This ensures we account for the significant rates of cooling that occur in photoionisation equilibrium at T ~ 10⁴ K. The second source included is ionised hydrogen (Hummer 1994) and helium (Rybicki & Lightman 1979) Bremsstrahlung. The third source is the cooling rate of hydrogen recombination from Hummer (1994), where we assumed hydrogen is fully ionised. Heating is calculated by assuming each H recombination (excluding recombinations directly to the ground state, i.e. the Case B approximation) rapidly results in a re-ionisation by the intense Lyman-continuum radiation field found in YMCs, with associated photo-heating of 5 eV per photoionisation. We obtain an equilibrium gas temperature T ≈ 8300 K for this prescription.

In Westerlund 1, the EUV radiation field is provided by the OB and WR stars in the cluster, which are observed to ionise the circumstellar nebulae of the cool stars to a high degree in the vicinity of the cluster (Andrews et al. 2018). This validates our above assumption that the medium is photoionised. Consistent with this observation, postprocessing of simulation snapshots in Sect. 3.4 with the TORUS radiation transfer code calculates the ionisation state of the gas as well as dust emission. From this we find that H is more than 99% ionised in the full domain except for within a radius of ≈ 3 × 10¹⁶ cm around the RSG where the wind is dense enough to self-shield from ionising radiation.

The simplified microphysics treatment is quite standard in HD modelling of nebulae (e.g. Garcia-Segura et al. 1996; Meyer et al. 2014; Mackey et al. 2025), where calculation of radiative transfer and non-equilibrium ionisation would make the simulations prohibitively expensive (Mathew et al. 2025). Haworth et al. (2015) investigated the differences between a simplified microphysics treatment and more complicated calculations with full radiative transfer and ionisation of H, He, and metals, applied to expanding H II regions. They found that differences were predominantly driven by the equilibrium temperature of the photoionised gas, determined by the gas composition and spectrum of the radiation field.

We performed our simulation using a 3D Cartesian domain of (x, y, z) ϵ [−10.0 × 10¹⁸,10.0 × 10¹⁸] cm, using five static mesh refinement levels. At the coarsest level, n = 0, the entire domain D is covered, and more refined levels with n > 0 cover a sub-domain of diameter D/2ⁿ, i.e. each finer level has a domain two times smaller than the level preceding it. We included 256³ grid cells per level of refinement, such that the maximum grid resolution is ∆x = 4.88 × 10¹⁵ cm.

To simulate the stellar wind, we placed the RSG star at the origin, and the nested grids were centred on the star. Our aim is to study the ablation of the RSG wind, and thus a constant mass loss source term is applied, with a surface temperature T_eff = 3500 K, a mass-loss rate Ṁ = 2 × 10⁻⁵ M_⊙ yr⁻¹ (Mackey et al. 2015) and wind terminal velocity υ_∞ = 25 km s⁻¹ (Mauron & Josselin 2011). The wind injection region is 5 × 10¹⁶ cm, or ten cells. This is much larger than the stellar radius, and thus we injected the wind at the terminal velocity.

To simulate the cluster wind, we imposed a constant x-velocity V_x = −1000 km s⁻¹, density ρ = 1 × 10⁻²⁴ g cm⁻³, and pressure 5 × 10⁻⁹ dyn cm⁻². While direct observational constraints are challenging, the cluster wind velocities cannot significantly exceed the adiabatic sound speed. This is of order 500–1000 km s⁻¹, using the diffuse X-ray gas temperatures measured by Kavanagh et al. (2011). Our velocity is also consistent with that of Silich et al. (2011), but lower than theoretical predictions (e.g. Stevens & Hartwell 2003; Härer et al. 2023), which are of order 2000–3000 km s⁻¹. Our density value is based on the density of the X-ray emitting gas from Kavanagh et al. (2011). The upstream boundary condition was set to the fixed inflow conditions of the cluster wind, and the other external boundaries are set to zero-gradient (outflow).

The simulation was run for 200 000 years with an explicit, finite-volume integration scheme that is 2nd-order accurate in time and space. The radiative heating and cooling source terms are included by the operator-splitting method. The timescale for the cluster wind to advect from the RSG to the downstream boundary is 10¹¹ s (about 3000 yr), much shorter than the simulation duration. The advection timescale for the RSG wind to reach the downstream boundary is 4 × 10¹² s (about 125 000 yr), and so we can expect the simulation to have reached a time-stationary state by 200 000 yr.

3 Results

3.1 XY slices

In Fig. 1, we show x − y slices of gas density, temperature and velocity magnitude at 50 kyr intervals. The cluster wind ablates the slowly expanding RSG wind. The contact discontinuity of the RSG wind and cluster wind becomes dynamically unstable in the upstream direction, with a combination of Rayleigh-Taylor (RT) and Kelvin-Helmholtz (KH) instability shearing clumps of shocked RSG wind into the wake downstream from the star. At early times, there is a dense and contiguous tail of ablated RSG material, probably influenced by initial conditions. This tail breaks up into clumps, and in the last two columns we have a clumpy and turbulent wake behind the star. These clumps are themselves being slowly ablated and their gas absorbed into the hot phase.

The ablation flow bears some similarity to the 2D von Karman vortex street, generated by subsonic flow past a solid object. Direct estimates of Reynolds numbers (Re) in simulations of this type are difficult due to the requirement for artificial viscosity. Since we are in between a laminar flow (Re ~ 10) and fully turbulent (Re ~ 1000 s), we can assume our Re is of order ~100 s. With this assumption, for an idealised von Karman vortex street (where the characteristic length scale, D, is the diameter of the cylinder), the shedding frequency, f, is related to the Strouhal number (St) by the ratio of D (in our case, twice the stand-off distance to the apex of the bow shock) and flow velocity υ. For flow past a cylinder υ is well defined and, importantly, sub-sonic, but in our case we have the complication of two different flow speeds (the cluster wind and the RSG wind) with widely differing sound speeds. The difficulty in assigning a characteristic velocity in such a multi-scale shear flow is addressed by Tan et al. (2021). If we nevertheless assume υ corresponds to the RSG wind velocity, we can calculate a characteristic timescale. As St only varies slightly for Re = 10²−10⁵, we assumed St = 0.2 and find a typical shedding timescale of ~3 kyr. However, due to the instability of the bow shock, D is neither fixed nor constant over the bow shock, and so the vortices are only quasi-periodic. This is similar to the behaviour found by Wareing et al. (2007) in simulations of asymptotic giant branch stars, albeit with much lower relative velocities.

Fig. 1 x − y slices through our simulation of the plane z = 0 at 50, 100, 150, and 200 kyr, showing distributions of density (upper), temperature (middle), and velocity magnitude (lower). The cluster wind is flowing from top to bottom, and the RSG is positioned at the origin. The colour scale is chosen to highlight dense, cool, and slow-moving material from the RSG.

3.2 Phase diagrams

In Fig. 2, we show the density-temperature phase diagrams for only the RSG wind, only the cluster wind, and the sum of these components at the end of our simulation. We considered the RSG wind to be ‘cool’, material above the photoionisation equilibrium temperature ~8300 K to be ‘mixed’, and material significantly hotter than this to be ‘hot’. In particular, RSG wind material that is over 10⁴ K must have been heated by interactions, shocks and other processes.

Focusing on the RSG wind first, there is a horizontal feature at low temperature that corresponds to material immediately surrounding the star, mainly unshocked stellar wind. At the bow shock, the RSG wind is ablated by the cluster wind and clumps break off. As seen in Fig. 1, these clumps can remain dense and cool. This explains the ‘hotspot’ around ~10⁻²⁰ g cm⁻³ and the equilibrium temperature of ~8300 K. As these clumps continue to be mixed, some of the material is heated and produces the long tail in that panel.

For the cluster wind in the centre panel of Fig. 2, there is a hotspot at the ambient cluster wind values of ρ = 1 × 10⁻²⁴ g cm⁻³ and T ~ 4 × 10⁷ K. As the cluster wind mixes with the RSG wind, some of this material is cooled. Much of the mixed material cools until it reaches the equilibrium temperature (~8300 K), and this is seen as a horizontal feature in the phase diagram. The presence of a feature at the equilibrium temperature in all three panels reflects the dynamical mixing of the two gas phases, which both heat the RSG wind and cool the cluster wind.

In Fig. 3, we show phase diagrams of V_x versus temperature for the RSG wind material at 50 kyr intervals. There is a hotspot at the wind injection temperature, visible as the isolated point in the upper left in all panels. As the RSG wind material expands spherically, we see a vertical feature where the material reaches the equilibrium temperature with a dispersion in velocity at the bow shock, which widens over time as the simulation approaches a time-stationary state.

The RSG wind material is ablated by the cluster wind, and rapidly accelerated as it is heated and then exits the domain, and the effect increases with time. The similarity of the last two panels implies that the simulation is approaching a time-stationary state.

Fig. 2 Temperature-density phase diagrams for the 3 pc sphere centred at the origin of our model for the RSG wind (left), cluster wind (middle), and all material (right) at the end of our simulation.

Fig. 3 Vx–temperature phase diagrams for the 3 pc sphere centred at the origin of our model for only the RSG wind material for t = 50, 100, 150, and 200 kyr.

3.3 Mass loading

To get an idea of the efficiency of the mass loading over time, we compared the time evolution of ‘hot’ M_h (T > 10⁴ K) and ‘cool’ M_c (T < 10⁴ K) RSG wind material (see Fig. 4). As expected, the cool RSG wind mass increases rapidly at early times. It stays relatively constant over time for small radii, but appears to peak and then decrease for larger radii (Fig. 4a). As the contiguous tail is broken up and ablated, the cool material is removed from the domain more rapidly over time. The mass of hot RSG wind material changes rapidly at first (Fig. 4b), likely due to initial conditions. At later times it tends to increase to a quasi-steady state for all radii. As this material is accelerated and exits the domain, we only find a lower limit for the amount of RSG material that is heated. This behaviour is consistent across different heating temperatures at small and large radii (Fig. 4c). The small fluctuations arise from the turbulent nature of the instabilities themselves.

The ratio of hot to cool RSG wind material tends to increase over time, again with turbulent fluctuations. Given that some of the heated wind is escaping from the domain, this represents a lower limit on the mass loading. Taking M_h/M_c ≈ 0.5, this represents a mass-loading efficiency of 33% at late times, albeit a more detailed numerical study with a larger simulation domain is required to measure this number accurately. We study the numerical convergence of these results in Appendix A.

Second generation star formation from rapidly cooling shocked stellar winds has been studied in detail by Wünsch et al. (2017). However, that requires material to be shielded (i.e. not be photoionised), which is never the case in the central part of Westerlund 1 that we are considering due to the magnitude of the incident ionising radiation. We therefore do not expect there to be second generation star formation resulting from the ablation flows in this case.

3.4 Dust maps

We produced synthetic dust maps using the TORUS Monte Carlo radiative transfer code (Harries et al. 2019) with the procedure outlined in Mackey et al. (2021) to post-process our simulation. To account for the external radiation field that contributes to dust heating, we included the RSG itself as one source, and placed two hot-star sources at positions [1, 1, 1] × 10¹⁸ cm and [1, −1, −1] × 10¹⁸ cm.

The spectral energy distribution (SED) of the two hot-star sources are assumed to be blackbodies with an effective temperature T_eff = 40 kK, radius 20 R_⊙, and a luminosity L_★ = 9.1 × 10⁵ L_⊙. While this is more luminous than estimated values for individual stars in Westerlund 1 (the supergiants are at least a factor of three less luminous; Clark et al. 2005; Beasor et al. 2021), the motivation for these parameters is to generate a radiation field comparable to that found in the inner regions of the Westerlund 1 cluster. In this way we avoid the complications of including tens of radiation sources with different positions and photospheric properties, making the simplest possible assumptions for this first study of the subject. At the location of the RSG, the resulting radiation energy flux from the two sources is 92.4 erg cm⁻² s⁻¹ and radiation energy density is 1.9 × 10³ eV cm⁻³. Further down-stream the radiation flux decreases with the inverse square law, so the radiative heating is reduced.

For the RSG we took a Kurucz (1992) stellar atmosphere model with T_eff = 3750 K and stellar radius 1000 R_⊙ with a stellar mass of 25 M_⊙. This corresponds to a luminosity of 1.78 × 10⁵ L_⊙. The RSG radiation field heats the dust in the immediate vicinity of the star, but in the downstream region the heating is dominated by the two hot radiation sources, particularly because the dust absorption cross-section peaks in the UV. We assumed silicate grains (Draine 2003) with a dust-to-gas ratio of 0.01, and a Mathis et al. (1977) grain-size distribution with power-law index of −3.3 from 0.005 to 0.2 µm. We further assumed that grains do not exist in gas with temperatures T > 10⁶ K, i.e. instantaneous grain destruction. This is a gross simplification of the rich diversity of dust observed in winds of RSGs (e.g. Decin 2021), but serves to obtain a reasonable estimate of the dust temperature and emissivity in the RSG wind. It would be relatively easy to obtain intensity variations of a factor of two or more by changing the dust composition, size distribution and dust-to-gas ratio within reasonable values, and so the predicted intensity maps have relatively large modelling uncertainties.

We present synthetic dust continuum-emission maps for 11 µm in Fig. 5 and 24 µm in Fig. 6, as these are of interest to compare to current and possible future observations of RSG outflows in Westerlund 1 (Guarcello et al. 2025). The upper panel in each figure shows a projection with the line of sight being the z-axis, and the lower panel shows the same but along the y-axis. The logarithmic intensity scale brings out some of the faint features further in the downstream wake of the flow. The faintest features are starting to show some of the photon-sampling noise that is inherent in Monte Carlo radiative transfer modelling.

The brightest emission at both wavelengths is concentrated around the RSG, where the radiation field is most intense and the dust has the highest temperature. Because the peak of the dust SED is at λ > 11 µm, and the 11 µm band is in the Wien tail of the SED, the emission in Fig. 5 is about 10× fainter than at 24 µm (Fig. 6). Moving to the downstream wake (more negative x coordinates) the dust temperature decreases further with distance from the radiation sources, and so there is very little 11 µm emission for x < −0.75 pc. The 24 µm emission is less sensitive to distance but eventually the dust temperature decreases such that this band is also in the Wien tail of the SED and so the brightness decreases significantly.

We see a range of clumps, globules with and without tails, and some filamentary emission connecting clumps that are being ablated from the RSG wind. These features are somewhat reminiscent of models of bow shocks from winds of cool runaway stars (Wareing et al. 2007; Mohamed et al. 2012), although here the much larger flow velocity of the cluster wind (and the higher resolution of our simulations) leads to a more highly clumped wake behind the star. The half opening-angle of the IR-emitting clumps in the downstream is about 30° – 40°, with the majority of the dusty wind material enclosed within a cylinder of radius ≈ 0.5 pc. It should be noted, however, that we imposed a planar inflow rather than a spherically expanding cluster wind (which is relevant on scales of a few parsecs), and so we probably underestimate the dispersion of the clumps in directions perpendicular to the cluster-wind flow ($-\widehat{x}$). The optical depth of the clumps to EUV/far-UV and optical stellar radiation is τ ≈ 0.075 – 0.3 (depending on clump size and photon energy), calculated using the opacity tables from TORUS, and so they are not self-shielding and are photoionised throughout to a high degree.

Projections of the 11 µm and 24 µm emission from different viewing perspectives are shown in Fig. 7, progressing from face-on (θ = 0°, where θ is the angle between the line of sight and the cluster-wind flow velocity) to edge-on (θ = 90°) from left to right. The face-on images have a circular symmetry with a ‘fried-egg’ morphology that results from a peak of emission in the termination-shock of the RSG wind. As θ increases, the nebula acquires the head-tail morphology of an ablation flow, with dynamical instabilities and a turbulent wake. At 11 µm, Fig. 7 shows an emission feature that resembles a bow shock from the edge-on perspective, although the cluster wind is only trans-sonic and no emission from it is visible in dust maps. This feature is the shocked RSG wind, which accumulates in a dense isothermal layer behind the RSG wind termination-shock and then advects downstream around the unshocked RSG wind region.

Fig. 4 Time evolution of stellar wind masses of different temperatures enclosed within spheres of selected radii. Panel a: cool (T < 104 K) material. Panel b: Hot (T > 104 K) material. Panel c: Material with T > 104 K (dotted), T > 3 × 104 K (dashed), and T > 105 K (dot-dashed) enclosed within spheres of 0.5 pc (blue) and 3.0 pc (purple). Panel d: ratio of cool (T < 104 K) to hot (T > 104 K) material.

Fig. 5 Thermal dust emission maps from the final snapshot of our simulation, showing surface brightness at 11 µm on a logarithmic scale. The colour bar shows log10 I11 µm/(MJy ster−1). The radiation sources heating the dust are the RSG at the origin (orange star) and two sources indicated by the blue stars at x = 0.32 pc and z = ±0.32 pc.

Fig. 6 Same as Fig. 5 but showing surface brightness at 24 µm on a logarithmic scale. The colour bar shows log10 I24 µm/(MJy ster−1).

3.5 Comparison to observations

We compared our synthetic observations with the JWST/MIRI F1130W observations towards Westerlund 1 by Guarcello et al. (2025), taken as part of GO-1905 (PI Guarcello). Since our objective is a phenomenological comparison rather than bespoke modelling of a specific system, we used the level 3 data products from the MAST¹ science archive without further processing or reduction.

We focused on three M-type supergiants: W20, W26, and W237, which are plotted in Fig. 8. In each case there is a substantial saturated inner zone, simply due to the high sensitivity of JWST and bright nature of the sources. However, exterior to that is a plethora of unsaturated features.

Qualitatively and quantitatively, the 11 µm synthetic image from Fig. 7 at θ = 90° is very similar to the JWST F1130W image of the RSG W237 (Guarcello et al. 2025; see also our Fig. 8). The synthetic images at θ = 60° and 90° are also remarkably similar to the F1130W image of the RSG W20, although there the JWST image seems to show brighter emission that our prediction. For both W20 and W26, a detailed comparison with simulations is difficult because large parts of the nebula are saturated in the F1130W image.

The synthetic 24 µm emission maps in Fig. 7 show emission that is more than ten times brighter than at 11 µm, as discussed already for Fig. 6. Plotting with a linear rather than logarithmic intensity scale brings out the chaotic and filamentary morphology of the wake behind the star, with clumps connected by ridges of dusty gas as the ablation process rips them away from the bulk of the RSG wind material.

The fluxes in the unsaturated parts of the extended features associated with the M-type supergiants are typically in the range 3000–8000 MJy ster⁻¹, comparable to (i.e. same order of magnitude as) the values in the simulated maps of Figs. 5, 6, and 7. The morphology of the flow is also similar, with a primary flow, multiple cometary secondary objects and even extended features that connect to the primary flow that likely result from some combination of RT and/or KH instabilities, as illustrated in Fig. 9. Such features have previously been discussed in the context of runaway asymptotic giant branch stars such as Mira (Wareing et al. 2007), although here the numerical resolution is significantly higher and a more complex flow structure is apparent.

It is notable that the model predicts similar fluxes to the observations without any contribution from polycyclic aromatic hydrocarbon (PAH) emission, despite there being a significant PAH C-H bending feature within this band (e.g. Tappe et al. 2006; Candian & Sarre 2015; Wesson et al. 2024). This suggests that the PAH abundance in the ablation flow is depleted, most likely due to destruction by shocks or radiation (e.g. Micelotta et al. 2010; Siebenmorgen & Krügel 2010). We note that PAHs are not universally present in RSG winds (Verhoelst et al. 2009) and so it is possible that they are not formed in the first place. Future observations may be able to resolve this question.

This appears to be supported by ALMA continuum observations of free-free emission from the nebulae around the cool supergiants in Westerlund 1 (Fenech et al. 2018), ATCA radio observations that also detect extended nebulae around the cool supergiants (Andrews et al. 2018), and Hα emission from the nebula around W26 (Wright et al. 2014), all implying that the winds of the cool supergiants are externally photoionised. However, confirming a deficit of PAH abundance will require a multi-wavelength study with dedicated models that treat PAHs, which is beyond the scope of this work.

Fig. 7 Synthetic thermal dust emission maps at 11 µm (top row) and 24 µm (bottom row) viewed from different angles, with emission shown on a linear scale in units of MJy ster−1. The location of the RSG is shown by the orange star, and the other two synthetic radiation sources by the blue stars projected at z = ±0.32 pc. The angle θ is measured between the observer’s line of sight and the cluster-wind flow velocity.

4 Discussion

4.1 RSG and cluster wind uncertainties

The mass-loss rates of RSGs are poorly constrained (Beasor et al. 2020) and can vary rapidly over time, including strongly clumped mass-loss episodes (Montargès et al. 2021). In this work we chose a constant ‘canonical’ value, but a time-varying mass-loss rate would increase the variation of the mass-loading rate and shedding timescale, as well as the morphology of dust emission. The cluster wind properties are also poorly constrained from observations. To fully account for this, a model of the cluster wind’s velocity, density, and temperature as a function of radius would be required. We also assumed a planar inflow here, whereas in reality it is spherically expanding, with potentially strong local perturbations from nearby individual stars. The cluster wind expansion radius is of the order of the cluster radius itself, and therefore could impact the cluster wind at the positions of the RSGs. We see the effects of our initial conditions for the first 25–50 kyr of our simulation, as we ‘switch on’ our star. In reality, the RSG would have evolved from a main sequence star with a hot, fast wind. The RSG would expand into the bubble excavated during the main sequence as opposed to the uniform ISM we assumed here.

4.2 Modelling uncertainties

Many previous works have focused on the physics of radiative turbulent mixing layers, which arise in many different astrophysical contexts (e.g. Fielding et al. 2020; Tan et al. 2021; Mackey et al. 2025). Fielding et al. (2020) show that while the phase structure is resolution dependent, the total cooling is well converged even for moderate resolutions. We did not consider the effects of magnetic fields in this work. We would expect the morphology of the dust filaments to be elongated along the cluster magnetic field lines (Gronke & Oh 2020) but the contribution to the turbulent mixing is likely not significant (Li et al. 2020). We did not consider thermal conduction in this work either. The importance of thermal conduction is debated in the literature, but recent work suggests that its effects are only significant if larger than the turbulent diffusivity (Ji et al. 2019; Tan et al. 2021).

Our dust modelling has several uncertainties. We assumed a generic hot radiation field, but the dust emission is sensitive to the specifics of this radiation field. For example, a hot star in the negative-x direction would heat up the dust and change the observed dust map considerably. We also made assumptions about the properties of the dust grains, and the ratio of gas to dust in the ablation flow. Given the number of uncertainties and the fact that we did not tune any parameters, the qualitative and quantitative agreement between our synthetic 11 µm emission maps and the observed emission around the RSGs in Westerlund 1 is encouraging.

4.3 Mass loading

We obtain a lower limit for mass-loading efficiency due to RSG wind material being accelerated and exiting the domain before it can be heated, as well as due to resolution as we discuss in Appendix A. Using the cluster wind velocity formalism given in Eq. (1) of Härer et al. (2023), and assuming a total current RSG mass-loss rate in Westerlund 1 of 1 × 10⁻⁴ M_⊙ yr⁻¹, we would expect the cluster wind to be slowed by a factor of ~10% due to RSG winds. Westerlund 1 is a young cluster with a small number of RSGs at present, so here this contribution is probably not important. At later times however, when a cluster wind is weaker due to fewer WR stars and mass loading is stronger due to more RSG stars, this effect could become more important. Furthermore, the JWST observations (Guarcello et al. 2025) show that there is significant mid-IR nebular emission in all directions around the star cluster, some of which cannot have come from the RSGs unless the cluster-wind geometry is very unusual. This material also contributes to the mass loading of the cluster wind, and so the combined effects of all of the cool, dusty gas is larger than the 10% effect we estimate for just the RSG winds.

In this work, we only considered a single set of steady-state stellar wind parameters, and future work is required to examine the effects of initial stellar mass and stellar evolution on mass loading in clusters. We also neglected other possibilities such as bursts and asymmetries in mass-loss, which would affect the morphology of the ablation structures and dust emission.

Fig. 8 JWST/MIRI F1130W observations towards three M-type supergiants in Westerlund 1. These are W26, W2O, and W237 from top to bottom and in order of decreasing declination. The colour bar units are MJy ster−1. In each case the central white pixels are saturated. The stars denote the positions of the three target supergiants. The colour scale is the same as the synthetic observations in Fig. 7.

4.4 Relevance to super star clusters

Super star-clusters are extremely massive YMCs that are sufficiently massive and tightly bound that they could be progenitors of globular clusters, which are often observed to have multiple stellar populations (Gratton et al. 2012; Bastian & Lardo 2018). The enrichment patterns of the second-generation stars imply that they formed from gas that was polluted by nuclear-processed gas ejected from the first generation of (massive) stars. Feedback from (hot) massive stars is one explanation for these multiple populations (e.g. Wünsch et al. 2008; de Mink et al. 2009; Krause et al. 2013), and cool supergiants have also been proposed (e.g. Szécsi et al. 2018). Rather specific conditions are required in order for winds of hot stars to cool efficiently and remain bound to the star cluster to form a second generation (Wünsch et al. 2017). RSG winds are already cool, dense and dusty, and thus readily suited to forming a second generation of stars if the wind material can be confined within the cluster. While Westerlund 1 may not be massive enough to retain its gas and form a 2nd generation in situ, it is very interesting to study (presumably) cool and dusty gas right in the core of the cluster. Measuring the thermal and chemical properties of this gas, and estimating its lifetime in the cluster, may give insights into how the multiple populations in globular clusters could have formed. This is especially interesting because some targeted searches for intra-cluster gas and dust failed to detect either in significant quantities in other star clusters (Bastian & Strader 2014).

4.5 Follow-up observations

Given the high velocities associated with the ablation flow, of order 1000 km s⁻¹, they should be spectrally resolvable even at medium resolution. VLT/ERIS-SPIFFER (Davies et al. 2023) is a near-IR integral field spectrograph facility with a high resolution K-middle filter including Br γ (e.g. Reiter et al. 2024). This should therefore be well suited to confirming those high velocities without the same susceptibility to saturation that JWST suffers. That would immediately rule out a primordial origin for the nebular emission seen towards Westerlund 1 and confirm that it instead results from stellar ejecta. It would also provide the ability to make empirical mass-loss rate estimates, which is a key factor in interpreting these ablation flows.

Westmoquette et al. (2010, 2013) obtained spectra of dynamical mixing layers at the boundary between hot cluster winds and pillars of dense gas and dust. They were able to measure the enhanced velocity dispersion in the mixing layer of about 100 km s⁻¹. Similar observations of these dusty clumps in Westerlund 1 should give observational constraints on the degree of mixing between hot and cool phases, and hence on the mass loading of the cluster wind. Comparison of such observations with synthetic spectra from HD simulations will be valuable for constraining physical processes operating in the mixing layer (e.g. HD mixing vs thermal conduction).

Fig. 9 Left: JWST/MIRI F1130W observations towards W20. Right: simulation density map. In addition to the primary flow and cometary knots, there are extended features, indicated by the arrows, that connect to the primary flow and arise due to RT and/or KH instabilities.

5 Conclusion

In this work, we computed a 3D HD model of a RSG embedded in a hot cluster wind with values typical for those in Westerlund 1. We summarise our findings below:

• A clumpy tail of ablated RSG wind material forms through dynamical instabilities at the contact discontinuity, with the cluster wind shearing and entraining cold gas into the cluster wind;

• The RSG wind material is heated as the clumps are further ablated in this tail, while some of the cluster wind is cooled due to dynamical mixing;

• The RSG wind material is accelerated as it is heated, and this acceleration tends towards a steady state over time;

• After any effects of the initial conditions have left the simulation domain, the fraction of cool RSG wind mass decreases over time while the fraction of hot RSG wind approaches a steady state;

• We determine a mass-loading efficiency of about 33% for the RSG wind, although this should be confirmed with higher-resolution simulations and for different RSG-wind and cluster-wind parameters;

• We produced synthetic dust continuum-emission maps for 11 µm and 24 µm for several viewing angles for comparison with current and future observations of RSG outflows in Westerlund 1 and other clusters;

• Our dust maps closely match observations of W20, W26, and W237 without any fine-tuning, in terms of both observed fluxes and morphology. Our omission of PAH emission in our dust map modelling suggests that PAHs are depleted in the ablation flow due to the effects of winds and radiation and/or a lack of PAH formation in the RSG winds in the first place;

• We encourage further observational studies of the fate of RSG winds and dusty clumps in clusters. This will help ascertain the role RSG winds play in mass loading, as well as their potential to seed future populations of stars, as inferred from observations of globular clusters;

• Follow-up observations with IR integral field spectroscopy such as with VLT/ERIS-SPIFFIER would enable detailed kinematic studies of the RSG-cluster wind interaction and resulting ablation flow. This would provide insights into the origin of Westerlund 1’s nebular emission, a direct probe of the so-far elusive RSG mass-loss rates, and help disentangle the various physical processes in the subsequent mixing.

Acknowledgements

We thank the anonymous referee for their constructive feedback which has improved the article. The authors thank Brian Reville for useful discussions on hydrodynamics. The simulations presented here were performed on the HPC system Raven at the Max Planck Computing and Data Facility. CJKL gratefully acknowledges support from the International Max Planck Research School for Astronomy and Cosmic Physics at the University of Heidelberg in the form of an IMPRS PhD fellowship. This publication results from research conducted with the financial support of Taighde Éireann – Research Ireland under Grant number 20/RS-URF-R/3712 (JM). JM is grateful to B. Reville and the Max Planck Institute for Nuclear Physics (MPIK) for welcoming him as a regular visitor during the preparation of this paper. AACS is supported by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) in the form of an Emmy Noether Research Group - Project-ID 445674056 (SA4064/1-1, PI Sander). TJH acknowledges funding from a Royal Society Dorothy Hodgkin Fellowship and UKRI guaranteed funding for a Horizon Europe ERC consolidator grant (EP/Y024710/1). This research has made use of the Astrophysics Data System, funded by NASA under Cooperative Agreement 80NSSC21M00561. This study made use of the following software packages: Astropy (Astropy Collaboration 2018), Numpy (Harris et al. 2020), matplotlib (Hunter 2007), yt (Turk et al. 2011), PION (Mackey et al. 2021), PYPION (Green & Mackey 2021), TORUS (Harries et al. 2019). Contributions. C.J.K.L. and J.M. conceived the project. C.J.K.L. performed and analysed the simulations and led the writing of the article. J.M. computed the TORUS dust maps, assisted with performing the simulations and interpreting the results. T. J.H. created the JWST figures, assisted with TORUS and interpreting the results. A.A.C.S. provided feedback on the manuscript and theoretical support for the article.

Appendix A Resolution test

To check the effects of resolution on our results, we ran an otherwise identical model using 192³ grid cells per level of refinement. PION has a minimum injection region radius of 7 cells, making this the lowest achievable resolution without having an artificially large injection region. We compare the time evolution of the RSG wind mass for the two models in Fig. A.1.

For a radius of 0.5 pc, both models have converged to similar values in all panels, and the differences likely originate from instabilities in the ablation flow. At 3.0 pc, the instabilities are more poorly resolved, reducing the efficiency of gas mixing. This suggests we do not have sufficient resolution at this radius, making our mass-loading estimate a lower limit. We present additional figures of our low resolution model as Figs. A.2, A.3, A.4, and A.5.

Figure A.2 shows that the downstream wake contains larger and more contiguous clumps that are more closely confined to the x-axis, when comparing this lower-resolution calculation with the higher-resolution simulation in Fig. 1. Despite this qualitative difference between the two resolutions, the phase diagrams (Figs. A.3 and A.4) nevertheless appear quite similar between the two resolutions. The mass ratio of hot-to-cold gas in Fig. A.5(d) is significantly lower than in Fig. 4. This shows that dynamical mixing is not sufficiently well resolved at this lower numerical resolution, and that higher resolution is required to demonstrate whether or not our results for mass loading have converged.

Fig. A.1 Time evolution of stellar wind masses of different temperatures for the low (dashed) and high (dotted) resolution models: Panel a: Mass of cool (T < 104 K) material vs time enclosed within spheres of 0.5 pc (blue) and 3.0 pc (purple). Panel b: Idem for hot (T > 104 K) material. Panel c: Plot of material with T > 104 K (blue) and T > 105 K (purple) enclosed within spheres of 0.5 pc (blue) and 3.0 pc (purple). Panel d: Ratio of cool (T < 104 K) to hot (T > 104 K) material.

Fig. A.2 Same as Fig. 1 but for the lower-resolution model.

Fig. A.3 Same as Fig. 2 but for the lower-resolution model.

Fig. A.4 Same as Fig. 3 but for the lower-resolution model.

Fig. A.5 Same as Fig. 4 but for the lower-resolution model.

References

All Figures

	Fig. 1 x − y slices through our simulation of the plane z = 0 at 50, 100, 150, and 200 kyr, showing distributions of density (upper), temperature (middle), and velocity magnitude (lower). The cluster wind is flowing from top to bottom, and the RSG is positioned at the origin. The colour scale is chosen to highlight dense, cool, and slow-moving material from the RSG.
In the text

	Fig. 2 Temperature-density phase diagrams for the 3 pc sphere centred at the origin of our model for the RSG wind (left), cluster wind (middle), and all material (right) at the end of our simulation.
In the text

	Fig. 3 Vx–temperature phase diagrams for the 3 pc sphere centred at the origin of our model for only the RSG wind material for t = 50, 100, 150, and 200 kyr.
In the text

Fig. 4 Time evolution of stellar wind masses of different temperatures enclosed within spheres of selected radii. Panel a: cool (T < 104 K) material. Panel b: Hot (T > 104 K) material. Panel c: Material with T > 104 K (dotted), T > 3 × 104 K (dashed), and T > 105 K (dot-dashed) enclosed within spheres of 0.5 pc (blue) and 3.0 pc (purple). Panel d: ratio of cool (T < 104 K) to hot (T > 104 K) material.

In the text

	Fig. 5 Thermal dust emission maps from the final snapshot of our simulation, showing surface brightness at 11 µm on a logarithmic scale. The colour bar shows log10 I11 µm/(MJy ster−1). The radiation sources heating the dust are the RSG at the origin (orange star) and two sources indicated by the blue stars at x = 0.32 pc and z = ±0.32 pc.
In the text

	Fig. 6 Same as Fig. 5 but showing surface brightness at 24 µm on a logarithmic scale. The colour bar shows log10 I24 µm/(MJy ster−1).
In the text

Fig. 7 Synthetic thermal dust emission maps at 11 µm (top row) and 24 µm (bottom row) viewed from different angles, with emission shown on a linear scale in units of MJy ster−1. The location of the RSG is shown by the orange star, and the other two synthetic radiation sources by the blue stars projected at z = ±0.32 pc. The angle θ is measured between the observer’s line of sight and the cluster-wind flow velocity.

In the text

	Fig. 8 JWST/MIRI F1130W observations towards three M-type supergiants in Westerlund 1. These are W26, W2O, and W237 from top to bottom and in order of decreasing declination. The colour bar units are MJy ster−1. In each case the central white pixels are saturated. The stars denote the positions of the three target supergiants. The colour scale is the same as the synthetic observations in Fig. 7.
In the text

	Fig. 9 Left: JWST/MIRI F1130W observations towards W20. Right: simulation density map. In addition to the primary flow and cometary knots, there are extended features, indicated by the arrows, that connect to the primary flow and arise due to RT and/or KH instabilities.
In the text

Fig. A.1 Time evolution of stellar wind masses of different temperatures for the low (dashed) and high (dotted) resolution models: Panel a: Mass of cool (T < 104 K) material vs time enclosed within spheres of 0.5 pc (blue) and 3.0 pc (purple). Panel b: Idem for hot (T > 104 K) material. Panel c: Plot of material with T > 104 K (blue) and T > 105 K (purple) enclosed within spheres of 0.5 pc (blue) and 3.0 pc (purple). Panel d: Ratio of cool (T < 104 K) to hot (T > 104 K) material.

In the text

	Fig. A.2 Same as Fig. 1 but for the lower-resolution model.
In the text

	Fig. A.3 Same as Fig. 2 but for the lower-resolution model.
In the text

	Fig. A.4 Same as Fig. 3 but for the lower-resolution model.
In the text

	Fig. A.5 Same as Fig. 4 but for the lower-resolution model.
In the text