© The Authors 2026

1. Introduction

Massive stars profoundly influence their surrounding interstellar medium (ISM) through mechanical, radiative, and chemical feedback, particularly in their late evolutionary stages. In these phases, their intense radiation field drives powerful, supersonic stellar winds. As these winds propagate, they can either interact with the ambient ISM or, in the case of a binary system, collide with the wind from a companion star. This leads to two possible interaction scenarios: wind–wind and wind–ISM interactions, which generate high Mach number shocks — suitable sites for relativistic particle acceleration through diffusive shock acceleration (DSA; e.g., Drury 1983). Relativistic particles in the shocked region emit nonthermal (NT) radiation. In particular, relativistic electrons produce synchrotron emission when they interact with the magnetic fields present in the shocks (e.g., White 1985). The detection of this synchrotron emission in the radio domain constitutes indirect but unambiguous evidence of particle acceleration.

The energy budget for particle acceleration depends on the stellar wind kinetic power (P_kin¹). Consequently, in recent years, studies have focused on O-type and Wolf-Rayet (WR) systems, which typically have P_kin > 10³⁶ erg s⁻¹. These investigations detected synchrotron emission in the colliding wind region of binary systems, revealing a strong correlation between particle acceleration and binarity (Dougherty & Williams 2000; De Becker & Raucq 2013).

To date, approximately 50 NT radio-emitting binary systems have been identified, which are classified as particle-accelerating colliding-wind binaries (PACWBs) (De Becker & Raucq 2013; De Becker et al. 2017). On the other hand, for single stars, stellar bow shocks produced by runaway massive stars are also a conducive environment, although clear detections of synchrotron emission have been reported in only three such systems to date (Benaglia et al. 2010; Moutzouri et al. 2022; van den Eijnden et al. 2024). For non-runaway stars, NT emission has so far only been detected in one bubble, G2.4+1.4 around WR 102 (Prajapati et al. 2019). This observational breakthrough opens a new window for studying particle acceleration and the physical processes associated with shocks in stellar bubbles. These observational studies, along with several theoretical works (e.g., Seo et al. 2018; Meyer 2024; De Becker et al. 2024, and reference therein), highlight the ability of massive stars to accelerate particles, positioning them as potential sources of Galactic cosmic rays (GCRs). However, the potential of single stars as particle accelerators and their significance as relevant sources of GCRs remains largely unexplored.

In this context, studying the bubbles surrounding WR stars is especially important for understanding the interaction between the stellar wind and the ambient ISM, given that WR winds typically exhibit P_kin > 10³⁷ erg s⁻¹. The termination shocks in these wind-blown bubbles (WBBs) can create the required conditions to drive efficient particle acceleration. Studies by Toalá & Arthur (2011, 2016) on the evolution of circumstellar material around WR stars provide a detailed overview of WBB formation due to wind–ISM interaction. Recent gamma-ray detections highlight stellar wind termination shocks in single and isolated WR environments as potential particle acceleration sites, where emission likely arises from GCR interactions with nearby matter or radiation fields (Inventar et al. 2025). In the early evolutionary phase, thermal overpressure and stellar winds drive the expansion of the H II region. As the stellar wind propagates through the surrounding ionized gas, it generates shocks that sweep up material, sculpting a bubble. By the end of the main sequence phase, the massive star is typically encompassed by a low-density, hot bubble bounded by a dense shell of neutral material, which expands at a few kilometers per second (Arthur 2007). After the main sequence phase, the mass-loss rate increases and the star evolves into a red supergiant, yellow supergiant, or luminous blue variable stage, depending on its initial mass (Meynet & Maeder 2003). These stages are characterized by very high mass-loss rates (10⁻⁴–10⁻³ M_⊙ yr⁻¹; Nugis & Lamers 2000) and slow winds (v_∞ = 10–100 km s⁻¹). As a result, the star can lose a significant fraction of its mass (up to 50%) through episodic, nonspherical outbursts of material (Humphreys 2010). In the final stage of stellar evolution, as the star enters the WR phase, it generates a powerful wind with terminal velocities between 700 and 5000 km s⁻¹ and a mass-loss rate of ∼(1–5)×10⁻⁵ M_⊙ yr⁻¹ (Nugis & Lamers 2000; Sander et al. 2012, 2019). This fast wind sweeps up and compresses the material expelled during earlier evolutionary stages, thereby forming a WR nebula.

We investigated the role of the local ISM and its interaction with stellar winds as a necessary ingredient for accelerating particles in single massive stars and assessed their potential as relevant sources of GCRs. For this, we carried out low-frequency (250–500 MHz and 550–950 MHz) radio observations of the WR bubble NGC 2359 encompassing WR 7.

This paper is organized as follows. Section 2 gives a brief overview of the WR bubble NGC 2359. Section 3 describes the upgraded Giant Meterwave Radio Telescope (uGMRT) observations, the complementary archival radio datasets used in this study, along with the data reduction process. The main results are presented in Sect. 4, and a detailed discussion of these results and their implications follows in Sect. 5. Finally, Sect. 6 summarizes our key conclusions.

2. WR bubble NGC 2359

With an angular size of 7.4′, this bubble presents a complex morphology in the IR and optical (shown in Fig. 1), indicative of wind–ISM interaction regions. As discussed in Schneps et al. (1981), NGC 2359 presents a filamentary, nearly spherical shell morphology with multiple streamers and a prominent bright region in its southern part. WR 7 is classified as a WN4b type (Smith et al. 1996) and exhibits a terminal velocity of 1600 km s⁻¹ with a mass-loss rate of 4 × 10⁻⁵ M_⊙ yr⁻¹ (Hamann & Koesterke 1998), resulting in a wind kinetic power of 3.3 × 10³⁷ erg s⁻¹. In the literature, different kinematic distances have been adopted for NGC 2359, ranging from 4.1 to 6.3 kpc. Furthermore, the photometric distance to WR 7 varies between 3.5 and 6.9 kpc (for a more in-depth discussion, we refer to Goudis et al. 1994). In addition, the parallax measured by Gaia (Gaia Collaboration 2023) places the WR star at a distance of 4.8 ± 0.6 kpc. Considering these estimates, we adopted a distance of 4.8 kpc for our study. Several parameters of WR 7 and the bubble are listed in Table 1.

Fig. 1. Maps of NGC 2359 in the optical and IR. (a) Optical morphology of NGC 2359 from the Digitized Sky Survey 2 (DSS2). (b) Spitzer color-composite image of the 3.6 μm (blue), 4.5 μm (green), and 8.0 μm (red) bands. The arrow indicates the direction of motion of WR 7 with respect to the local medium (see Sect. 4.1 for details).

Several studies have examined the distribution and characteristics of the molecular material associated with WR 7 and NGC 2359 and probed the interaction between ionized and neutral gas using radio and spectral line observations. Cappa et al. (1999) reported a strong similarity between the continuum and the line emissions of the nebula using Very Large Array (VLA) radio continuum (1465 MHz) and HI 21 cm observations of NGC 2359, with angular resolutions of approximately 30″ and 45″, respectively. Assuming that the radio emission is free–free, the ionized mass of the filamentary shell region was estimated to be 70 M_⊙, and for the entire complex, including the filaments, the estimated mass was in the range of 900–1100 M_⊙, adopting a distance of 5 kpc. This large ionized mass suggests that the shell predominantly consists of interstellar matter. The gas kinematics of NGC 2359 have also been studied using CO molecular transitions at different resolutions (Schneps et al. 1981; Cappa et al. 2001; Rizzo et al. 2001, 2003). The CO (1–0) observations from these studies revealed similar results of enhanced molecular line emission at velocities of 37, 54, and 67 km s⁻¹. Rizzo et al. (2001) reported that the surrounding molecular cloud at a velocity of 67 km s⁻¹ was accelerated (from 67 to 53 km s⁻¹) by at least 14 km s⁻¹, likely due to interaction with stellar winds. In a follow-up work, using higher-resolution CO and ¹³CO data, Rizzo et al. (2003) identified three distinct velocity components. They analyzed their spatial distribution and properties and suggested that multiple layers of shocked molecular gas in the WBB are associated with episodic interactions of the central star with the circumstellar medium, likely created during the earlier luminous blue variable stage and/or the current WR phase of WR 7.

The complex filamentary bubble morphology and multiple velocity components observed in CO spectra are likely related to episodic mass loss during previous evolutionary phases and shock-driven interactions. Additionally, X-ray observations of NGC 2359 reveal high plasma temperatures (∼10⁶–10⁷ K; Toalá et al. 2015). This is consistent with the predictions of models for hot shocked stellar wind bubbles, where the stellar wind impacts wind material from a previous evolutionary phase or the ISM and forms an adiabatically shocked region of gas with temperatures ranging from 10⁷ K to 10⁸ K. Together with the high wind kinetic power of WR 7, these factors suggest that shocks in NGC 2359 can provide the necessary conditions for accelerating particles.

3. Radio datasets

We probed the nature of the radio emission and the complex morphologies associated with NGC 2359 using our uGMRT observations at Band 3 (250–500 MHz) and Band 4 (550–950 MHz), along with complementary archival radio datasets at different frequencies. The properties – the central frequency, resolution, rms, and the largest angular structure (LAS) – of the radio maps are compiled in Table 2. To address the “missing flux” problem in interferometric observations, we estimated the LAS, which is the largest coherent structure that can be imaged by an interferometer and is determined by the length of the shortest baseline. In addition, for all radio observations, we adopted a conservative systematic flux uncertainty due to absolute flux calibration of 10%.

3.1. uGMRT observations

We measured the continuum emission associated with NGC 2359 using the uGMRT located in Pune, India. The uGMRT employs a hybrid Y-shaped array configuration with 30 antennas of 45-m diameter; 12 antennas are randomly placed within a compact central region spanning 1 km², while the remaining 18 antennas are uniformly distributed along the three arms. This hybrid configuration, with baselines ranging from 100 m to 25 km, enables the study of ionized emission across both large and small spatial scales. An overview of the GMRT system is available in Swarup et al. (1991), with details of the upgrades presented in Gupta et al. (2017). Dedicated observations were conducted in Bands 3 and 4 using the GMRT Wideband Backend (GWB) correlator, configured to process 2048 and 4096 channels, respectively. The primary calibrator (3C147) was observed at the beginning and end of the observation sessions for flux and bandpass calibration. The phase calibrator (0744–064) was observed after each 30-minute scan of the target to correct for phase and amplitude variations throughout the observing period.

We processed the full-band data in Bands 3 and 4 using the CAPTURE continuum imaging pipeline for uGMRT (Kale & Ishwara-Chandra 2021). It uses tasks from Common Astronomy Software Applications (CASA; McMullin et al. 2007) for the flagging, calibration, imaging, and self-calibration processes. The Perley & Butler (2017) scale was used to set the flux density calibration. After the initial rounds of editing and calibration, we used the multi-term multifrequency synthesis (see Rau & Cornwell 2011) algorithm in the tclean task to account for possible deconvolution errors in wide-band imaging. Five rounds of phase-only and one round of amplitude-phase self-calibration were conducted to obtain the final images of the full-band data. We note that the residual map for Band 3 (not shown here) indicates that we could not fully recover the flux from the larger scale structures (due to large phase variation and broad radio frequency interference), suggesting that the flux estimate in this band is a lower limit.

To investigate the in-band spectral index in Band 4, we employed the standard procedures of data reduction using the Astronomical Image Processing System (AIPS) for generating sub-band images. A recent paper by Rashid et al. (2024) elaborates on the reliability and accuracy of spectral index estimation using sub-band splitting method. The Band 4 uv data were initially divided into six sub-bands, each with a bandwidth of approximately 35 MHz, centered at 622, 657, 692, 728, 761, and 801 MHz. Each sub-band was then carefully examined for corrupted data, including nonfunctional antennas, bad baselines, and radio frequency interference (RFI). Calibration and imaging were performed separately for each sub-band. Channel averaging was optimized to have negligible bandwidth smearing. Further, wide-field imaging techniques were used to account for w-term effects (Cornwell & Perley 1992). Several iterations of self-calibration were performed to correct for phase errors and improve the image quality. We did not carry out sub-band imaging of Band 3 data as it (i) suffered from broad RFI, (ii) showed large phase variations, and (iii) had a smaller bandwidth of ∼200 MHz compared to Band 4. All images used in our analysis were primary beam-corrected.

3.2. Archival radio data

To complement our targeted uGMRT observations, we utilized several archival radio datasets at different frequencies. These include radio maps at (i) 150 MHz from the TIFR GMRT Sky Survey (TGSS; Intema et al. 2017), (ii) 887 and 943 MHz from the Rapid Australian Square Kilometre Array Pathfinder (ASKAP) Continuum Survey (RACS; McConnell et al. 2020), and (iii) 1425, 4860, and 8689 MHz from the National Radio Astronomy Observatory VLA Archive Survey (NVAS). For these datasets, we extracted pipeline-processed images from the respective survey portals². The parameters of the archival radio datasets are listed in Table 2.

Regarding the RACS data, it is worth noting that we have not included available data at frequencies 1367 and 1655 MHz in our analysis. As discussed in Duchesne et al. (2023, 2025), (i) both frequencies are significantly affected by broadband RFI and (ii) there is a decrease in flux density with increasing source size due to incomplete uv coverage. The authors mention that baselines below 75 m (for 1367 MHz) and 100 m (for 1655 MHz) were excluded while creating the maps. Additionally, Duchesne et al. (2023) highlight the challenges in imaging fields containing sources with extended emission in the Galactic Plane (as is the case with NGC 2359). These include an increased number of artifacts localized to extended sources, higher rms, and a significant portion of the flux density not recovered from them.

4. Results

In this study, we primarily focused on determining the nature of the radio emission in NGC 2359 by estimating the spectral index and searching for evidence of NT emission. The spectral index helps us understand the radiation mechanisms and decipher the underlying physical processes.

4.1. Morphology of the source

Radio continuum maps from GMRT, ASKAP and VLA are shown in Figs. 2 and 3. All maps are convolved with a common circular beam size of 22″ to compare the emission at different frequencies, except for the map at 150 MHz that has a poorer resolution. Almost all radio maps trace a similar complex morphology as seen in optical and IR (Fig. 1). The southern part of the complex shows strong radio emission. In the maps from the six sub-bands, the diffuse emission near WR 7 is observed to be weak, especially at 761 and 801 MHz, compared to the full-band 735 MHz image, even though the resolution and uv coverage are similar. This is possibly because the sub-band images have much smaller bandwidths than the full-band CAPTURE maps, which reduces their sensitivity to fainter, extended emission. The radio morphology reported in Cappa et al. (1999), where this source was mapped using VLA data at 1465 MHz with a lower resolution of ∼39″ × 25″, is similar to the uGMRT maps. Due to the lower resolution, more extended diffuse emission is traced in their map (see Fig. 2 of Cappa et al. 1999). The maps reveal an asymmetrical bubble structure around WR 7, which is offset from the center. A similar picture is also seen in the bubble, G2.4+1.4, where the star WR 102 is observed to be offset from the bubble center (Prajapati et al. 2019). In a previous study of G2.4+1.4, Dopita & Lozinskaya (1990) attributed this asymmetry to Rayleigh-Taylor instabilities, proposing that the location of WR 102 near a molecular cloud restricts expansion to that side while allowing expansion on the other side into a low-density ISM, resulting in a “scalloping” morphology. For NGC 2359, CO observations (Schneps et al. 1981; Rizzo et al. 2001, 2003) have revealed bright extended emission in the southern part, supporting the above explanation.

Fig. 2. Radio maps of NGC 2359 obtained using our uGMRT data. (a) Maps obtained using the full-band data. (b) Sub-band images for Band-4 data. All maps are convolved to a circular beam of 22″. For each map, the central frequency and rms noise (σ) in units of mJy beam−1 are mentioned in each panel. The “X” marks the location of the star in all the panels. Identified apertures (E1–E4) and the central region used to obtain the SEDs are shown in the 735 MHz map. Contour levels are as follows: (i) For maps with central frequencies of 402 and 735 MHz, contour levels start from 3σ and increase in steps of 4σ. (ii) The maps with central frequencies of 622, 657, 692, and 728 MHz have contour levels of [3, 8, 13, 28, 43, 58, 73, 88, 103, 118] × σ. (iii) For maps with central frequencies of 761 and 801 MHz, the contour levels are [3, 8, 13, 28, 43, 58, 73] × σ. In all panels, the contours are smoothed over 3 pixels using a Gaussian kernel.

Fig. 3. Radio maps of NGC 2359 obtained using archival data. For each map, the central frequency and rms noise (σ) in units of mJy beam−1 are mentioned in each panel. Barring the image at 150 MHz, all maps are convolved to a circular beam of 22″. The 150 MHz map has a beam size of 29.1″ × 24.9″ and is smoothed across 3 pixels using a Gaussian kernel. The “X” marks the location of the star in all the panels. Identified apertures (E1–E4) and the central region used to obtain the SEDs are shown in the 1425, 4860, and 8689 MHz maps. Contour levels are as follows: (i) For maps with the central frequency of 150 MHz, contour levels start from 2σ and increase in steps of 1σ. (ii) For maps with central frequencies of 887 and 943 MHz, contour levels start from 3σ and increase in steps of 5σ. (iii) The 1425 MHz map has contour levels of [3, 7, 11, 15, 23, 31, 39, 47, 54, 62, 70, 78] × σ. (iv) For maps with central frequencies of 4860 and 8689 MHz, contour levels start from 3σ and increase in steps of 4σ. In all panels, the contours are smoothed over 3 pixels using a Gaussian kernel.

Additionally, several simulation studies (e.g., Brighenti & D’Ercole 1995b; Brighenti & D’Ercole 1995a; Reyes-Iturbide et al. 2019) have shown that a high proper motion of the star leads to an asymmetry in the bubble structure. To investigate this possibility, we estimated the stellar velocity with respect to its surrounding medium by correcting the proper motions measured by Gaia for Galactic rotation as in Martínez et al. (2023)³. The radial (V_r) and tangential (V_t) velocities, and the direction of motion (θ), are provided in Table 1. This implies that WR 7 is moving with a supersonic velocity in the direction shown in Fig. 1, which is the projection of its motion onto the plane of the sky and explains the off-centered location of WR 7. Moreover, such a high peculiar velocity should shape the WBB into a more paraboloidal structure known as a bow shock (Weaver et al. 1977). The apex of the bow shock is located in the direction of motion of the star, at a stagnation radius given by the balance of the momentum flux of the wind with the ram pressure of the surrounding ISM. The fact that the radio maps do not show any discernible bow shock morphology could be due to the radial velocity being higher than the tangential velocity, though the uncertainty in the velocity estimates is large. This results in a projected structure that is more circular, which would be challenging to identify from our radio maps. Furthermore, the medium through which the WR wind propagates might deviate significantly from the standard homogeneous medium assumed in bow shock models. Moreover, no bow shocks have been identified around any WR stars, although some of them are classified as runaway stars (Meyer et al. 2020). This non-detection of bow shocks is consistent with numerical simulations by Meyer et al. (2020) that study the circumstellar material surrounding a 60 M_⊙ runaway WR star. The simulations suggest that the interaction of the high-velocity WR phase wind with the slower wind from earlier evolutionary stages leads to the formation of spherical shells. These shells seem to co-move with the star regardless of its proper motion or the properties of the surrounding ISM, and consequently, no bow-shaped region is identified. Furthermore, Meyer et al. (2020) proposed that a low-density medium around these stars can explain the absence of bow shocks in runaway stars like WR 71, WR 124, and WR 148. For NGC 2359, without knowledge of the density, it is difficult to draw a conclusion on the reason for the non-detection of a bow shock. Additional studies to estimate the density of the ISM are required for further investigation.

4.2. Nature of the radio emission: Spectral index

The spectral index (α) is defined as S_ν ∝ ν^α, where S_ν is the flux density at frequency ν. For thermal free–free emission from ionized gas, α ranges from +2 (optically thick) to −0.1 (optically thin). The free–free opacity depends on the frequency, which causes sources to shift from being optically thick at lower frequencies to optically thin at higher frequencies. This transition point (turnover frequency) depends on the plasma density. Moreover, the slope of the spectrum can vary depending on the density distribution of the plasma. For instance, in the case of stellar winds of massive stars (with a density distribution that scales with the inverse of the distance to the star), an intermediate value of α ∼ 0.6 is expected (Wright & Barlow 1975). On the other hand, optically thin synchrotron radiation generally has a steep spectral index of α ≈ −0.5, assuming high Mach number adiabatic shocks and monoatomic gas. However, absorption or suppression processes can modify the observed spectrum, making the spectral index less negative or even positive (e.g., Melrose 1980). A spectral index of α < −0.1 strongly indicates NT emission, while positive values typically point to a thermal origin. In many astrophysical environments, thermal plasma can coexist with a population of relativistic electrons, leading to spectra that include both thermal and NT components. Despite this complexity, it is possible to differentiate between these components by analyzing the spectral energy distribution (SED) across a wide range of frequencies. NT emission, with a negative spectral index, usually dominates at lower frequencies, while thermal radiation contributes mainly at higher frequencies (e.g., De Becker 2018). Therefore, a comprehensive sampling of the SED at low frequencies is necessary to investigate the interplay between absorption and emission processes.

To estimate the global spectral index, we calculated the integrated flux density within the 5σ contour of the 735 MHz map. This aperture encompasses the central region of the bubble, as shown in Figs. 2 and 3, and matches the extent of the field of view (FOV) of the 4860 MHz map. However, as the FOV at 8689 MHz is much smaller than this region, flux was not obtained at that frequency. The resultant SED shows an interesting trend, where the slope is seen to vary: (i) between 150 MHz and 402 MHz, the slope is positive (α = 0.27 ± 0.15); (ii) between 402 MHz and 1425 MHz, the slope becomes approximately flat within the error bars (α = −0.01 ± 0.17); and (iii) from ∼1.4 GHz to higher frequencies, the slope turns negative (α = −0.93 ± 0.12). To examine the variation of spectral index across the bubble, we identified apertures with bright emissions at all frequencies from 150 MHz to 8.7 GHz, labeled E1–E4 (see Figs. 2 and 3), and estimated the spectral index from the integrated flux density within each aperture. The SEDs shown in Fig. 4 follow the same trend as the global SED of the central region: (i) between 150 MHz and 402 MHz the spectral index is α = 0.4 to 0.8; (ii) between 402 MHz and 1425 MHz, the slope is almost flat; and (iii) from ∼1.4 GHz to higher frequencies, the slope ranges from α = −0.2 to −0.6. The spectral indices for the above frequency ranges are given in Table 3. Consistent with this, the radio SED presented in Fig. 4b, derived from the Band 4 uGMRT data for the six sub-bands (622–801 MHz) shows a nearly flat trend. As discussed previously, the sub-band images do not fully recover the diffuse emission, and the average flux densities of the apertures for the sub-band images are slightly lower than the values at 735 MHz, especially for E2. This can be attributed to the smaller bandwidth of the sub-band maps, as discussed previously in Sect. 4.1. Hence, in our spectral analysis, we used the flux densities derived from the full-band CAPTURE maps and archival radio maps.

Fig. 4. Radio SED for the identified apertures (E1–E4) shown in Figs. 2 and 3. (a) SED determined using the maps at 150, 402, 735, 887, 943, 1425, 4860, and 8689 MHz. The data point at 402 MHz (Band 3) is a lower limit (see Sect. 3.1 for details). (b) SED determined using the maps of the six sub-bands (622, 657, 692, 728, 761, and 801 MHz) of the Band 4 GWB uGMRT data. The uncertainty in the flux density measurements was estimated using (Lal 2021) ΔS = [(Sν × f)2 + rms2 × Nbeams]0.5, where Sν is the flux density, f is an absolute flux density calibration uncertainty (taken as 10%), and Nbeams is the number of synthesized beams in the aperture.

The accuracy of the spectral index estimates strongly depends on the spatial scales probed in the frequency bands. Given the large spatial extent of the bubble, we explored whether the flux densities at higher frequencies are affected by missing flux because of filtering out of structures larger than the LAS. For this, we compared the FOV, aperture sizes and LAS of the observations at 4860 and 8689 MHz. For VLA 4860 MHz (observed in VLA-C configuration), the native resolution, FOV, and LAS are 6.5″ × 3.9″, 9.2′, and 6′, respectively. For VLA 8689 MHz (observed in VLA-D configuration), they are 11.9″ × 7.6″, 5.2′, and 3′, respectively. On the other hand, the diameter of the entire bubble is 7.4′, while the central region and the largest aperture (E1), selected for deriving the SEDs have sizes of 5.8′ and 1.3′, respectively. In the VLA 4860 MHz map, the size of the entire bubble is slightly larger than the LAS; the size of the central region is comparable to the LAS, and aperture E1 is around five times smaller than the LAS. Given the small FOV of the VLA 8689 MHz map, we could only obtain the flux for aperture E1, which has a size much smaller than the LAS. Hence, it is unlikely that there is any significant loss of flux in these higher frequencies for the selected apertures.

5. Discussion

The observed radio SEDs of the bubble, as well as regions within, present distinct slopes in different frequency domains. Evidence of NT emission is seen above 1 GHz, where the spectral index values are negative, even steeper than the canonical value of −0.5 for optically thin synchrotron radiation from relativistic electrons accelerated by DSA in high Mach number adiabatic shocks. As we move below 1 GHz to ∼400 MHz, the flattening is indicative of a turnover due to absorption/suppression processes at work (e.g., Melrose 1980). To elucidate the turnover mechanism and derive physical information about the system, we need to introduce a radiative model.

5.1. SED modeling

We introduced an SED model that is physically motivated. The emission region is most likely composed of a combination of ionized gas and relativistic particles. Thus, we included both free–free and synchrotron emission components, although the negative spectral indices observed already point to synchrotron being the dominant component. The observed low-frequency synchrotron SED can be affected by: (1) synchrotron self-absorption (SSA), (2) the Razin-Tsytovich (RT) effect, and (3) free–free absorption (FFA) by thermal plasma mixed with the relativistic particles (hereafter, internal FFA). Each of these processes leads to specific spectral features. Under the assumption of a homogeneous region, we discuss the emission and absorption processes in more detail below.

5.1.1. Emission processes

5.1.1.1. Synchrotron radiation:

The optically thin synchrotron emission arising from a relativistic electron population in the presence of a magnetic field can be written as

$$\begin{array}{c}\hfill S_{\mathrm{sy}}=A_{\mathrm{sy}}{\left(\frac{\nu }{{\nu }_0}\right)}^{{\alpha }_{\mathrm{sy}}},\end{array}$$(1)

where S_sy is the synchrotron flux density, A_sy is a normalization factor, α_sy is the synchrotron spectral index and ν₀ is a reference frequency. This parameterization encodes information on the relativistic particle population and the magnetic field intensity. Namely, for a power-law electron energy distribution, N(E) = N₀ E^−p, A_sy is a function of N₀, p, and B, and α_sy = −(p − 1)/2 (e.g., Blumenthal & Gould 1970). We can rewrite these quantities in terms of energy density in NT electrons, U_e, and in magnetic fields, U_B, as

$$\begin{array}{c}\hfill U_{\mathrm{e}}\approx \frac{E_{\min }^{2-p}}{p-2}\frac{N_0}{V},U_B=\frac{B^2}{8\pi },\end{array}$$(2)

where V is the volume of the emitting region and E_min ∼ 1 MeV is the minimum energy of the relativistic electrons.

5.1.1.2. Thermal free–free emission:

The thermal free–free flux density can be expressed as (Mezger & Henderson 1967)

$$\begin{array}{c}\hfill S_{\mathrm{ff}}=3.07\times {10}^4{T}_e{\nu }^2\mathrm{\Omega }(1-e^{-\tau (\nu )}),\end{array}$$(3)

where the optical depth (τ) is given by

$$\begin{array}{c}\hfill \tau (\nu )=8.21a\times {10}^{-2}{\nu }^{-2.1}EM{T}_{\mathrm{e}}^{-1.35},\end{array}$$(4)

and the emission measure (EM) and solid angle (Ω) are

$$\begin{array}{c}\hfill EM=2n_{\mathrm{e}}^{2}R,\mathrm{\Omega }=\pi {\left(\frac{R}{d}\right)}^2.\end{array}$$(5)

In the above equations, S_ff is the integrated flux density in Jy; T_e is the electron temperature in K; ν is in GHz; EM is in cm⁻⁶ pc; n_e is the electron number density in cm⁻³; Ω is the solid angle subtended by the source in steradian; R is the radius of the apertures in parsecs, and d is the distance to the bubble in parsecs. Also, based on Table 6 of Mezger & Henderson (1967), we used a value of 0.99 for the correction factor (a).

5.1.2. Turnover processes

5.1.2.1. SSA:

This is the process by which synchrotron photons are re-absorbed by the same relativistic electron population that emits them. It leads to an optically thick spectrum with a positive slope (α = 2.5) below a specific turnover frequency, causing a notable reduction in synchrotron flux density at lower frequencies. To note, SSA is only effective in very compact and luminous radio sources in which the number density of relativistic electrons is high. In contrast, NGC 2359 is an extended structure that spans several parsecs and exhibits relatively low luminosities. Hence, SSA is unlikely to be a significant absorption mechanism in this particular environment and is not considered in our model.

5.1.2.2. RT effect:

The synchrotron radiation produced by relativistic electrons embedded in a thermal plasma is suppressed when the refractive index of the medium is less than unity (Melrose 1972). This is known as the RT effect, and it depends on the plasma frequency, which in turn depends on n_e and B. Observationally, it causes a turnover below the Razin frequency (ν_Razin) given by (Pacholczyk 1970)

$$\begin{array}{c}\hfill {\nu }_{\mathrm{Razin}}=20\frac{n_{\mathrm{e}}}{B}(\mathrm{MHz}),\end{array}$$(6)

where n_e is the electron number density in cm⁻³ and B is the magnetic field strength in μG. The effect results in an exponential cutoff in the synchrotron radio spectrum at lower frequencies (e.g., Bloot et al. 2022):

$$\begin{array}{c}\hfill S_{\mathrm{sy}}=A_{\mathrm{sy}}{\left(\frac{\nu }{{\nu }_0}\right)}^{{\alpha }_{\mathrm{sy}}}{e}^{-\frac{{\nu }_{\mathrm{Razin}}}{\nu }},\end{array}$$(7)

where the parameters have the same meaning as in Eq. 1.

5.1.2.3. Internal FFA:

Where a thermal plasma and relativistic electrons coexist, the ionized gas can significantly absorb the synchrotron photons as they propagate within the emitting region, causing a turnover in the synchrotron spectrum referred to as internal FFA. In this case, the drop in radio flux can be expressed as (Tingay & de Kool 2003)

$$\begin{array}{c}\hfill S_{\mathrm{sy}}=A_{\mathrm{sy}}{\left(\frac{\nu }{{\nu }_0}\right)}^{{\alpha }_{\mathrm{sy}}}\frac{1-e^{-\tau (\nu )}}{\tau (\nu )},\end{array}$$(8)

where the parameters have the same meaning as in Eqs. 1 and 4. The occurrence of internal FFA in astrophysical synchrotron sources has been addressed by de Bruyn (1976).

Although that study was motivated by the need to provide an adequate interpretation framework for extragalactic sources, the basic principle deserves to be considered for other synchrotron sources. Additionally, internal FFA has more recently been proposed to be responsible for the morphology of the radio SED of the PACWB WR 147 (Tasseroul et al. 2025).

5.1.3. Composite SED

Considering the above processes, we fit the observed SEDs with the following composite model:

$$\begin{array}{cc}\hfill S_{\mathrm{total}}& =A_{\mathrm{sy}}{\left(\frac{\nu }{{\nu }_0}\right)}^{{\alpha }_{\mathrm{sy}}}{e}^{-\frac{{\nu }_{\mathrm{Razin}}}{\nu }}\left(\frac{1-e^{-\tau (\nu )}}{\tau (\nu )}\right)\hfill \\ \hfill & +3.07\times {10}^4{T}_e{\nu }^2\mathrm{\Omega }(1-e^{-\tau (\nu )}).\hfill \end{array}$$(9)

We fixed the reference frequency at ν₀ = 1 GHz. To reduce the number of free parameters, the electron temperature was fixed at T_e = 9000 K, a typical value for ionized regions. We note that varying T_e within plausible bounds has a very small impact on the SED.

To ensure that all possible solutions are explored and obtain reliable fits, we utilized a Bayesian inference approach that employs a Markov chain Monte Carlo (MCMC) algorithm. In particular, we used the EMCEE sampler (Foreman-Mackey et al. 2013) implemented in the Bilby Python package (Ashton et al. 2019). We adopted uniform priors within a broad range of physically reasonable values for all fitted parameters. In particular, for α_sy we used a prior range of [ − 1.0, −0.5]. However, in regions for which its posterior distribution peaks at the hard limit of −0.5, we fixed α_sy = −0.5. This value corresponds to the theoretical prediction of DSA in high Mach number adiabatic shocks (Drury 1983), as expected in NGC 2359. Values harder than this (α_sy > −0.5) would imply that particles follow a power-law energy distribution with index p < 2, which is not favored by particle-acceleration theory. The list of parameters and the corresponding explored ranges are given in Table 4.

The fitted SEDs are shown in Fig. 5, and the corner plots presenting the complete details of the posterior distributions are provided in Appendix A. The best-fit parameters along with the chi-square (χ²) values and degree of freedom (d.o.f) are given in Table 5.

Fig. 5. SED fitting for apertures E1–E4 and the central region using the composite model described in Eq. (9). The shaded regions show the 1σ confidence interval. The bottom panels of each image present the residuals of the fit. The corresponding corner plots showing the posterior distribution of each parameter are presented in Appendix A.

5.2. Parameter estimates and caveats

The composite phenomenological model can closely reproduce the observed SEDs and the behavior of the spectral index (see Table 3). To explain the observed turnover, ν_Razin has to be ∼400 MHz. Given its dependence on electron density and magnetic field (see Eq. 6), this constrains the allowed values of n_e and B (namely the ratio between these two quantities, as can be seen in the corner plots in Figs. A.1 and A.2). Additionally, the electron density cannot be too high, as this would increase the free–free contribution and result in an incorrect slope at higher frequencies, inconsistent with the observed SED (which can also be seen in the corner plots in Figs. A.1 and A.2). Therefore, the turnover can only be reproduced if the magnetic field is sufficiently weak, given the range of electron densities allowed by the free–free emission constraints.

The values of n_e that we derive are a factor of four less than the density estimated using the ionized mass and size of the whole complex reported by Cappa et al. (1999); if we consider the entire bubble, it is an order of magnitude lower. However, these authors assumed that the emission arises entirely from thermal free–free radiation (which our analysis shows is not the case); hence, their ionized mass and number densities are likely overestimated.

Of interest is the low magnetic field strengths of ∼0.5–4 μG derived by the model. In comparison, Prajapati et al. (2019) estimated the magnetic field strength in the synchrotron-emitting region of the bubble G2.4+1.4 to lie in the range of 120–770 μG. This range corresponds to the assumed energy density ratio between NT particles and the magnetic field, ranging between 10 and 0.01. The synchrotron luminosity is proportional to the energy density in the magnetic fields and to the energy density in relativistic electrons. In this case, the weak magnetic field implies that an extremely large amount of energy must be injected into the relativistic particles to reproduce the observed synchrotron emission. Specifically, for a magnetic field of only a few microgauss, using Eq. 2 we infer that the energy density in relativistic particles must be around 10⁶ times greater than the magnetic energy density. Such an extreme energy distribution may not be physically plausible, as it is unlikely to be sustained by the available stellar wind kinetic power.

The main caveat in our analysis is the assumption of a homogeneous emitting region, which is clearly an oversimplification for NGC 2359. Assuming a different density distribution would change the exact shape of the absorbed spectrum. For instance, in the case of FFA, for a homogeneous cloud, the spectral index is 2, while for a massive stellar wind where the density decreases with the square of the distance from the star, the index is close to 0.6 (Wright & Barlow 1975). This exemplifies that internal FFA can lead to a range of spectral indices in the observed spectrum.

A physically plausible scenario for NGC 2359 could involve an inhomogeneous medium composed of dense, irregular substructures embedded within a lower-density environment. In such a configuration, FFA is not uniform across the region, as it depends on both the opacity and spatial distribution of the clumps. Radiation passing through denser regions can undergo significant absorption, while emission that travels through more diffuse inter-clump zones may remain largely unaffected. This selective filtering of emission, governed by the spatial variation in density and opacity, provides a more flexible framework for interpreting the observed SED features. A similar scenario of clumpy ionized regions is considered in different astrophysical environments (e.g., Lacki 2013; Conway et al. 2018; Mutie et al. 2025). However, a quantitative assessment of this possibility would require an appropriate description of such complex structures, which can only be obtained through detailed magnetohydrodynamic simulations that are beyond the scope of this work.

6. Conclusions

Nonthermal emission in massive stars is commonly detected in the colliding-wind region of binary systems under the wind–wind interaction scenario (De Becker & Raucq 2013). In contrast, only a few studies have investigated isolated, non-runaway WR stars with high wind kinetic power using low-frequency radio data. Among these, WR 114 and WR 142 (Saha et al. 2023) are not associated with extended structures and no radio emission is detected from them, whereas WR 102 (Prajapati et al. 2019) and WR 7 (this work) are surrounded by bubble-like structures and exhibit clear evidence of synchrotron emission arising from wind–ISM interaction. These results suggest that, for single stars, high wind kinetic power is a necessary but not sufficient condition for particle acceleration.

The interaction between the stellar wind and the ambient ISM is crucial. Termination shocks in the stellar winds can provide the physical conditions required for efficient particle acceleration, leading to detectable synchrotron radiation. This depends on the properties of the local ISM, such as density and temperature, together with stellar characteristics, including how long the star remains in the WR phase, the stellar proper motion, the mass of the swept-up material, and how efficiently wind energy is dissipated. However, it is unclear how these factors regulate the formation of strong shocks on parsec scales, underscoring the need for further observational and theoretical studies.

We investigated the WR bubble NGC 2359 surrounding WR 7 using uGMRT observations at low frequencies (250–500 MHz and 550–950 MHz), along with complementary archival datasets from GMRT, VLA, and ASKAP. The observed SEDs exhibit a steep negative slope above 1 GHz, a turnover below 1 GHz with a relatively flat slope between 400 MHz and 1 GHz, and a steep positive slope at frequencies below 400 MHz. To extract physical information and understand the processes responsible for the spectral turnover, we modeled the observed SEDs using a composite model. This model includes contributions from both synchrotron and free–free emission, with two turnover mechanisms: internal FFA and the RT effect. We fit the SEDs using the MCMC algorithm and estimated the electron number density and magnetic field strength. Our analysis revealed synchrotron radiation from the bubble characterized by spectral indices that are typically steeper than the canonical value of −0.5. The observed turnover is mainly due to the RT effect, with a small contribution from internal FFA. This is only the second detection of NT radiation from a stellar bubble associated with a single massive star, following the first discovery in G2.4+1.4. These findings confirm that wind–ISM interactions in such bubbles can create the necessary conditions for accelerating particles to relativistic speeds, establishing that isolated massive stars can be considered potential sources of GCRs. Our analysis enabled us to place tighter constraints on the electron number density, which is significantly lower than the previous estimates (Cappa et al. 1999). The model also infers a very low magnetic field strength of a few microgauss. A potential caveat of our analysis is the assumption of a homogeneous emitting region. Based on the observed SEDs and derived model parameters, a more physically realistic scenario for NGC 2359 involves an inhomogeneous medium consisting of dense, irregular substructures embedded within a lower-density environment. Overall, it is crucial to investigate particle acceleration in a larger sample of WR bubbles using data across a broad radio frequency range to gain insights into the necessary conditions that drive particle acceleration, the underlying physical processes, and various aspects of shock physics. With the upcoming Square Kilometre Array facility, significant breakthroughs are expected due to its very high sensitivity across different scales and broad frequency coverage below 1 GHz.

Acknowledgments

We thank the referee for valuable suggestions that helped us to improve the manuscript. A.S acknowledges support from the National SKA Program of China (2025SKA0140100) and the National Natural Science Foundation of China (No. 12573025). A.S. and A.T. thank Dr. Dharam V. Lal from the National Centre for Radio Astrophysics for fruitful discussions. S.d.P. acknowledges support from ERC Advanced Grant 789410. R.S. and I.C.CH. acknowledge the support of the Department of Atomic Energy, Government of India, under project no. 12-R&D-TFR-5.02-0700. This research is part of the PANTERA-Stars collaboration (https://www.astro.uliege.be/~debecker/pantera), an initiative aimed at fostering research activities on the topic of particle acceleration associated with stellar sources. We thank the staff of the GMRT that made these observations possible. GMRT is run by the National Centre for Radio Astrophysics of the Tata Institute of Fundamental Research. This research has made use of NASA’s Astrophysics Data System Bibliographic Services.

References

Appendix A: Posterior distribution: Corner plots

In Figs. A.1 and A.2 we present the corner plots showing the posterior distributions of the model parameters for all SED fits. The diagonal panels display the 1D distributions of each parameter, with the orange line indicating the best-fit value (median), and the dashed lines marking the 16th and 84th percentiles. The off-diagonal panels show the 2D joint distributions of parameter pairs, with contours representing the 1, 2, and 3σ confidence intervals.

Fig. A.1. Corner plot for the posteriors of the MCMC fitting for apertures E1 (top left), E2 (top right), E3 (lower left), and E4 (lower right). The normalization of the synchrotron component (Asy), electron number density (ne), and magnetic field strength (B) are in units of mJy, cm−3, and μG, respectively.

Fig. A.2. Corner plot for the posteriors of the MCMC fitting for the central region of the bubble. The normalization of the synchrotron component (Asy), electron number density (ne), and magnetic field strength (B) are in units of mJy, cm−3, and μG, respectively.

All Tables

All Figures

	Fig. 1. Maps of NGC 2359 in the optical and IR. (a) Optical morphology of NGC 2359 from the Digitized Sky Survey 2 (DSS2). (b) Spitzer color-composite image of the 3.6 μm (blue), 4.5 μm (green), and 8.0 μm (red) bands. The arrow indicates the direction of motion of WR 7 with respect to the local medium (see Sect. 4.1 for details).
In the text

Fig. 2. Radio maps of NGC 2359 obtained using our uGMRT data. (a) Maps obtained using the full-band data. (b) Sub-band images for Band-4 data. All maps are convolved to a circular beam of 22″. For each map, the central frequency and rms noise (σ) in units of mJy beam−1 are mentioned in each panel. The “X” marks the location of the star in all the panels. Identified apertures (E1–E4) and the central region used to obtain the SEDs are shown in the 735 MHz map. Contour levels are as follows: (i) For maps with central frequencies of 402 and 735 MHz, contour levels start from 3σ and increase in steps of 4σ. (ii) The maps with central frequencies of 622, 657, 692, and 728 MHz have contour levels of [3, 8, 13, 28, 43, 58, 73, 88, 103, 118] × σ. (iii) For maps with central frequencies of 761 and 801 MHz, the contour levels are [3, 8, 13, 28, 43, 58, 73] × σ. In all panels, the contours are smoothed over 3 pixels using a Gaussian kernel.

In the text

Fig. 3. Radio maps of NGC 2359 obtained using archival data. For each map, the central frequency and rms noise (σ) in units of mJy beam−1 are mentioned in each panel. Barring the image at 150 MHz, all maps are convolved to a circular beam of 22″. The 150 MHz map has a beam size of 29.1″ × 24.9″ and is smoothed across 3 pixels using a Gaussian kernel. The “X” marks the location of the star in all the panels. Identified apertures (E1–E4) and the central region used to obtain the SEDs are shown in the 1425, 4860, and 8689 MHz maps. Contour levels are as follows: (i) For maps with the central frequency of 150 MHz, contour levels start from 2σ and increase in steps of 1σ. (ii) For maps with central frequencies of 887 and 943 MHz, contour levels start from 3σ and increase in steps of 5σ. (iii) The 1425 MHz map has contour levels of [3, 7, 11, 15, 23, 31, 39, 47, 54, 62, 70, 78] × σ. (iv) For maps with central frequencies of 4860 and 8689 MHz, contour levels start from 3σ and increase in steps of 4σ. In all panels, the contours are smoothed over 3 pixels using a Gaussian kernel.

In the text

Fig. 4. Radio SED for the identified apertures (E1–E4) shown in Figs. 2 and 3. (a) SED determined using the maps at 150, 402, 735, 887, 943, 1425, 4860, and 8689 MHz. The data point at 402 MHz (Band 3) is a lower limit (see Sect. 3.1 for details). (b) SED determined using the maps of the six sub-bands (622, 657, 692, 728, 761, and 801 MHz) of the Band 4 GWB uGMRT data. The uncertainty in the flux density measurements was estimated using (Lal 2021) ΔS = [(Sν × f)2 + rms2 × Nbeams]0.5, where Sν is the flux density, f is an absolute flux density calibration uncertainty (taken as 10%), and Nbeams is the number of synthesized beams in the aperture.

In the text

	Fig. 5. SED fitting for apertures E1–E4 and the central region using the composite model described in Eq. (9). The shaded regions show the 1σ confidence interval. The bottom panels of each image present the residuals of the fit. The corresponding corner plots showing the posterior distribution of each parameter are presented in Appendix A.
In the text

	Fig. A.1. Corner plot for the posteriors of the MCMC fitting for apertures E1 (top left), E2 (top right), E3 (lower left), and E4 (lower right). The normalization of the synchrotron component (Asy), electron number density (ne), and magnetic field strength (B) are in units of mJy, cm−3, and μG, respectively.
In the text

	Fig. A.2. Corner plot for the posteriors of the MCMC fitting for the central region of the bubble. The normalization of the synchrotron component (Asy), electron number density (ne), and magnetic field strength (B) are in units of mJy, cm−3, and μG, respectively.
In the text