2.1. Askap Emu

The object WR16 has been seen in three Australian Square Kilometre Array Pathfinder (ASKAP) Evolutionary Map of the Universe (EMU) observations. SB46953 observed the tile EMU $\_1017$ -60 on December 12th 2022, SB51428 observed the tile EMU $\_0954$ -55 on July 13th 2023, and SB53568 observed the tile EMU $\_0936$ -60 on October 7th 2023. SB46953 is excluded from analysis due to its proximity to the tile’s edge, which can cause a drop in sensitivity.

The data were reduced using the standard Australian Square Kilometre Array Pathfinder (ASKAP) pipeline, ASKAPSoft, using multi-frequency synthesis imaging, multi-scale cleaning, self-calibration and convolution to a common beam size (Guzman et al. Reference Guzman2019). As the source was observed in two usable Evolutionary Map of the Universe (EMU) fields, the images were combined with the radio imaging software Miriad (Sault, Teuben, & Wright Reference Sault, Teuben, Wright, Shaw, Payne and Hayes1995) using the task imcomb. The final image was created by combining the observations with equal weighting. We measure a number of background regions around WR16 to find the background difference in the images. We find a median rms value of 37.5 $\unicode{x03BC}$ Jy beam $^{-1}$ and a mean value of 37 $\unicode{x03BC}$ Jy beam $^{-1}$ , compared with the previous values of median = 36 $\unicode{x03BC}$ Jy beam $^{-1}$ and mean = 43.5 $\unicode{x03BC}$ Jy beam $^{-1}$ (SB51428), and median = 46.5 $\unicode{x03BC}$ Jy beam $^{-1}$ and mean = 50.5 $\unicode{x03BC}$ Jy beam $^{-1}$ (SB53568). This background noise level is slightly higher than typical for Evolutionary Map of the Universe (EMU) ( $\sigma=20-30\,\unicode{x03BC}$ Jy beam $^{-1}$ ; Hopkins et al. Reference Hopkins2025), which can be attributed to WR16 being in an area of higher background Galactic emission ( $b =$ $-$ 2.55). The final image is shown in Figure 1 with a resulting sensitivity of $\sigma$ =37 $\unicode{x03BC}$ Jy beam $^{-1}$ and a synthesised beam of $15^{\prime\prime}\times15^{\prime\prime}$ .

2.2. Gaia

We use measurements made by the Gaia astrometry satellite (Gaia Collaboration et al. Reference Collaboration2016) for WR16. We obtain proper motion (PM) and parallax values from Data Release 3 (DR3) (Gaia Collaboration et al. Reference Collaboration2023) and use them to determine the star’s movement and distance. Discussion surrounding the use of this specific data release is found in Section 3.1.

2.3. Nanten CO

In order to explore the distribution of molecular clouds surrounding WR16, we analysed the archival $^{12}$ CO(J = 1–0) data taken by the NANTEN 4-m radio telescope (Mizuno & Fukui Reference Mizuno, Fukui, Clemens, Shah and Brainerd2004). The angular resolution of the data cube was $2.^{\prime}6$ . The typical noise fluctuations were $\sim$ 0.3 K at a velocity resolution of 0.65 km s $^{-1}$ .

2.4. Other data

We include observations from the SuperCOSMOS (Hambly et al. Reference Hambly2001) and Wide-Field Infrared Survey Explorer (WISE) (Wright et al. Reference Wright2010) sky surveys to compare the nebulosity surrounding WR16 at 656.281 nm and 22 $\unicode{x03BC}$ m (Figure 2).

Figure 2. Detection of WR16 and its shell at 22 $\unicode{x03BC}\!$ m from and Wide-Field Infrared Survey Explorer (WISE) infrared (Left), and SuperCOSMOS H $\alpha$ (Right). Both sub-images are linearly scaled.

We also use previous flux density measurements using Australia Telescope Compact Array (ATCA) from Leitherer et al. (Reference Leitherer, Chapman and Koribalski1995), Leitherer et al. (Reference Leitherer, Chapman and Koribalski1997), and Chapman et al. (Reference Chapman, Leitherer, Koribalski, Bouter and Storey1999). These measurements are included in Table 1, and discussed in Section 3.1.

Table 1. Measured flux densities of the Wolf-Rayet star WR16 and its nebulous shell.

^aMeasurements made as part of this paper.

^bMeasurements using Australia Telescope Compact Array (ATCA) observations from Chapman et al. (Reference Chapman, Leitherer, Koribalski, Bouter and Storey1999).

^cMeasurements using Australia Telescope Compact Array (ATCA) observations from Leitherer et al. (Reference Leitherer, Chapman and Koribalski1995).

We have searched for high-energy associations to WR16 and its shell. Fermi-LAT Data Release 4 (Abdollahi et al. Reference Abdollahi2022; Ballet et al. Reference Ballet, Bruel, Burnett and Lott2023) shows no corresponding Gamma-Ray emission. The SRG/eROSITA all-sky survey (eRASS) Data Release 1 (Merloni et al. Reference Merloni2024), and observations taken with XMM-Newton (Strüder et al. Reference Strüder2001; Turner et al. Reference Turner2001) show no corresponding X-Ray emission.

3.1. Measurements

Using a Gaia parallax of 0.4380 $\pm$ 0.0168 mas (Gaia Collaboration Reference Collaboration2020), we estimate the distance of WR16 to be 2.28 $\pm$ 0.09 kpc. This is in accordance with the photometric distance calculated in van der Hucht (Reference van der Hucht2001) of 2.37 kpc, and with a distance of 2.23 $\pm$ 0.39 kpc derived using Ca ii lines (Megier et al. Reference Megier, Strobel, Galazutdinov and Krełowski2009). Previously, Gaia Data Release 2 (DR2) (Gaia Collaboration Reference Collaboration2018) provided a distance of 2.66 $\pm$ 0.23 kpc (Bailer-Jones et al. Reference Bailer-Jones, Rybizki, Fouesneau, Mantelet and Andrae2018). Considering that Gaia DR3 provides a distance that is concordant with other measurements, we adopt 2.28 $\pm$ 0.09 kpc to be the most accurate distance value.

We obtain proper motion (PM) from Gaia DR3 (Gaia Collaboration Reference Collaboration2020), for WR16, which are $\mu_{\alpha}$ : $-9.458\pm 0.021$ mas yr $^{-1}$ and $\mu_{\delta}$ : $5.054\pm 0.018$ mas yr $^{-1}$ . Using the Gaia parallax and proper motion values, we calculate the peculiar velocity with respect to the local interstellar medium following Comerón & Pasquali (Reference Comerón and Pasquali2007) and Cichowolski et al. (Reference Cichowolski2020). We determine the stellar peculiar velocity to be $V_{\alpha}(pec) =$ –45.3 $\pm$ 5.4 km s $^{-1}$ and $V_{\delta}(pec) =$ 22.8 $\pm$ 4.7 km s $^{-1}$ , indicating that the star is moving in a north-west direction, which is discussed further in Section 3.3. This matches the brighter north-western component of the circular nebulosity observed in radio, infrared, and H $\alpha$ (Figures 1 and 2). We then calculate the peculiar tangential velocity to be 50.7 $\pm$ 6.9 km s $^{-1}$ .

We measure the inner WR16 nebula to have a flux density of 72.2 $\pm$ 7.2 mJy at 943.5 MHz, taking an error of 10% (as discussed in Smeaton et al. Reference Smeaton2024b; Filipović et al. Reference Filipović2024). We also measure the flux density of the central star at 943.5 MHz, and find it to be 0.35 $\pm$ 0.04 mJy. We combine this with the Australia Telescope Compact Array (ATCA) flux density measurements listed in Table 1 and Section 2.4 to compute a spectral index for the star (See Figure 3). The radio spectral index is defined as $S \propto \nu^{\alpha}$ , where S is flux density, $\nu$ is the frequency and $\alpha$ is the spectral index. We calculate the spectral index using the linregress Footnote ^a function in the Python scipy (Virtanen et al. Reference Virtanen2020) library. The function uses the linear least-squares regression method to find a best fit line of $\alpha\,=\,+0.74\pm0.02$ for WR16. This is very close to the canonical spectral index of an isothermal, spherical stellar wind ( $\alpha=+0.6$ , Panagia & Felli Reference Panagia and Felli1975; Wright & Barlow Reference Wright and Barlow1975), and is similar to that of Dougherty & Williams (Reference Dougherty and Williams2000) for similar Wolf-Rayet (WR) stars that also exhibit thermal characteristics.

Figure 3. Spectral index plot of WR16 star, using flux density values from Table 1.

Figure 4. (a) – Distribution of the NANTEN $^{12}$ CO(J = 1–0) line emission toward WR16 (Mizuno & Fukui Reference Mizuno, Fukui, Clemens, Shah and Brainerd2004). The integration velocity range is from $-7$ to 0 km s $^{-1}$ . The incomplete circular shell indicates a wind-blown bubble detected in CO (see Section 3.2). The dashed circle shows the position of the inner circular shell of WR16. (b,c) – Position–velocity diagrams of CO. The integration range of Right Ascension is from $149{{{.\!\!^\circ}}}33 to 149{{{.\!\!^\circ}}}92$ for (b) and from $148{{{.\!\!^\circ}}}42$ to $149{{{.\!\!^\circ}}}01$ for (c). The dashed semicircles represent the expanding motion of CO due to the stellar wind from WR16 (see Section 3.2).

We fit a circular region to the outer edge of the shell and measure an angular diameter of $8.^{\prime}42$ , which is similar to the angular size estimated in Cichowolski et al. (Reference Cichowolski2020). We use the derived distance and the angular size measured on the Evolutionary Map of the Universe (EMU) image (Figure 1) to determine the shell’s true linear size and find that the shell has a diameter of 5.57 $\pm$ 0.22 pc. Due to the shell’s presence in near–infrared and H $\alpha$ (Figure 2) we infer that the shell also has thermal origin. Because the shell has not yet been observed at other frequencies, we cannot confirm thermal origin from the spectral index.

3.2. Large-scale wind-blown bubble of CO

Figure 4(a) shows the velocity-integrated intensity map of CO toward WR16 with an integration range of from $-7$ to 0 km s $^{-1}$ . We found a new candidate for a large-scale wind-blown bubble whose radius is approximately $0\rlap{.}^{\circ}75$ or $\sim$ 30 pc at the distance of WR16. Since the inner circular shell of WR16 is approximately placed at the geometric center of the CO bubble, they appear to be physically associated to each other.

Figure 4(b) and (c) show the position–velocity diagrams of CO. We found clear evidence of an expanding gaseous bubble whose expanding velocity is $\sim$ 7 km s $^{-1}$ . The fact that the spatial extent in the Declination direction on the position–velocity diagram varies depending on the integrated Right Ascension range is also roughly consistent with the three-dimensional expansion motion of the molecular clouds.

We argue that the dynamical timescale of this expanding gaseous bubble can be naturally understood if it is assumed to have been driven by WR16. The dynamical time scale of the expanding bubble is estimated to be (bubble radius)/(expanding velocity) $\sim$ 30 pc/7 km s $^{-1}$ $\sim$ 4 Myr. This time scale is roughly consistent with the typical lifetime of a massive star, suggesting that the CO expanding bubble was likely formed by stellar feedback effects such as strong stellar winds from the progenitor of WR16.

3.3. Shell expansion

Figure 5 shows several on-image measurements of the shell surrounding WR16. We find that the tangential peculiar velocity of the star passes through the geometric center of the $8{{{.\mkern-4mu^\prime}}}42$ circle (RA (J2000) = 09:54:57.8, Dec (J2000) = –57:43:58.9), as indicated by the dashed black line in the figure. The distance of the star to the center of the nebulosity is $44\rlap{.}^{\prime\prime}81$ . We convert this distance from angular size to linear size (0.49 pc, or $1.53\times 10^{13}$ km), and multiply the calculated velocity in Section 3.1 by seconds in a year ( $3.154\times 10^{7}$ s) to obtain a yearly velocity of $1.58\times 10^{9}$ km yr $^{-1}$ , $5.12\times 10^{-5}$ pc yr $^{-1}$ . We can determine the amount of time it has taken for WR16 to travel to its current point by dividing the linear size by the yearly velocity, and return a value of $\sim 9\,500\pm 1\,300$ yr.

Assuming that the shell originated at the point that WR16 passed through its geometric centre coordinates, we can measure its average expansion velocity based on the distance travelled by WR16. The radius of the shell is $4.^{\prime}21$ , which converts to a linear size of $8.62\times 10^{13}$ km (2.79 pc). Dividing this by the years travelled ( $\sim 9\,500\pm 1\,300$ yr) and dividing by seconds in a year, we determine the average expansion velocity of the circular shell to be $280\pm 40$ km s $^{-1}$ . It is important to note that this only assumes a 2–dimensional plane, and is therefore a lower limit on age values and an upper limit on expansion values.

We take our calculated expansion velocity to be the average rate of expansion because it is presumed that the shell is slowing down the more it expands and interacts with the surrounding interstellar medium. WR16 also has a wind terminal velocity ( $v_{\infty}$ ) of 630 km s $^{-1}$ (Toalá et al. Reference Toalá, Guerrero, Ramos-Larios and Guzmán2015), so it is likely that the wind-driven shell started with velocities similar to this, and has slowed over the $\sim 9\,500\pm 1\,300$ yr timescale. This velocity sits conveniently between the typical expansion speed of Luminous Blue Variable (LBV) shells of $\sim$ 50 km s $^{-1}$ (Weis Reference Weis, Neiner, Wade, Meynet and Peters2011), like that of AG Car (Smith Reference Smith, van der Hucht and Hidayat1991), and the faster ejecta of Eta Car (600 km s $^{-1}$ , Steffen et al. Reference Steffen2014). This supports the idea that the shell originated in a previous Luminous Blue Variable (LBV) phase.

3.4. Mass-loss rate

Based on the observed parameters of WR16 at 943.5 MHz, we are able to determine a mass-loss rate ( $\dot{M}$ ) at this frequency using the equations outlined in (Wright & Barlow Reference Wright and Barlow1975, Equation 20) and (Leitherer et al. Reference Leitherer, Chapman and Koribalski1995, Equation 3). We calculate a Gaunt factor, $g_v$ , at 943.5 MHz to be 5.99, assuming a gas temperature of 10 $^4$ K and a Z of 1. Using the same constants for our calculation as Leitherer et al. (Reference Leitherer, Chapman and Koribalski1995), and substituting our flux, distance, and frequency, we derive a mass-loss rate of $1.753\times 10^{-5}\,{\rm M}_\odot\,$ yr $^{-1}$ .

This is lower than the calculated mass-loss rate from Leitherer et al. (Reference Leitherer, Chapman and Koribalski1995) of $3.981\times 10^{-5}\,{\rm M}_\odot\,$ yr $^{-1}$ . However, recalculating their value using our Gaia distance measurement of 2.28 kpc, we get a result of $2.188\times 10^{-5}\,{\rm M}_\odot\,$ yr $^{-1}$ . These measurements are quite similar, and may indicate that the mass-loss rate has decreased slightly over the 30 yr between observations. It is also important to note that the canonical value for mass-loss rates in WN8 stars like WR16 (Section 1) is $1.995\times 10^{-5}\,{\rm M}_\odot\,$ yr $^{-1}$ (Crowther Reference Crowther2007).

Figure 5. Measurements of the circular nebulosity surrounding WR16 and peculiar velocity mapped over the EMU image, similar to Figure 1. The $8.^{\prime}42$ diameter shell is shown as a black circle. The magenta circle indicates the position of WR16, and the red arrow shows its distance from the shell’s geometric center, as well as its direction. The black dashed line represents the projected path of WR16’s peculiar velocity, which points toward the north-west (top-right).

Using the calculated mass-loss rate, we are also able to determine a lower limit of ionising photons that WR16 are released per second as in (Wright & Barlow Reference Wright and Barlow1975, Equation 23). We use a radius of 19 R $_\odot$ (Schmutz, Hamann, & Wessolowski Reference Schmutz, Hamann and Wessolowski1989) to calculate a lower limit value of $N_{UV} \gt 1.406\times 10^{47}\,\textit{s}^{-1}$ , which is in concordance with the canonical value for WN8 stars of $1.259\times 10^{49}\,\textit{s}^{-1}$ (Crowther Reference Crowther2007).

3.5. Similarities to other objects

WR16 presents an almost symmetrical circular shell similar to the Galactic Supernova Remnants (SNRs) Teleios (Greek $\tau\epsilon\lambda\epsilon\iota o\varsigma$ – meaning perfect, Filipović et al. submitted). Although they share similarities, they are unique objects for their given type. This shows that different objects can share the same morphologies, and can possibly be equated to similar formation processes.

Another object similar to the WR16 shell is the radio ring Kýklos, discussed in Bordiu et al. (Reference Bordiu2024).They discuss the likelihood of the ring-like object being an outburst of a Wolf-Rayet (WR) star. This object is also situated in our Galaxy, but is significantly smaller than WR16 (WR16: $8.^{\prime}42$ , Kýklos: 80 $^{\prime\prime}$ ). Both systems show evidence of thermal emission, however, Kýklos presents a ring-like structure rather than a diffuse circle like WR16. The ring-like structure and smaller size may indicate that Kýklos is at an earlier developmental stage compared to WR16.

We see significant brightening toward the north–west part of the WR16 circular shell. This is likely due to the closer proximity of WR16, which also explains the much fainter south–east portion. This edge brightening is similar to the Lagotis H ii region discussed in Bradley et al. (Reference Bradley2025). In the case of Lagotis, edge brightening is caused by a cluster moving into a molecular cloud, whereas for WR16, it can be seen as the star catching up with its shell edge.