A global view on star formation: The GLOSTAR Galactic plane survey

Context. While over 1000 supernova remnants (SNRs) are estimated to exist in the Milky Way, only less than 400 have been found to date. In the context of this apparent deficiency, more than 150 SNR candidates were recently identified in the D-configuration Very Large Array (VLA-D) continuum images of the 4–8 GHz global view on star formation (GLOSTAR) survey, in the Galactic longitude range − 2 ◦ < l < 60 ◦ . Aims. We attempt to find evidence of nonthermal synchrotron emission from 35 SNR candidates in the region of Galactic longitude range 28 ◦ < l < 36 ◦ , and also to study the radio continuum emission from the previously confirmed SNRs in this region. Methods. Using the short-spacing corrected GLOSTAR VLA-D+Effelsberg images, we measure the ∼ 6 GHz total and linearly polarized flux densities of the SNR candidates and the SNRs that were previously confirmed. We also attempt to determine the spectral indices by measuring flux densities from complementary Galactic plane surveys and from the temperature-temperature plots of the GLOSTAR-Effelsberg images. Results. We provide evidence of nonthermal emission from four candidates that have spectral indices and polarization consistent with a SNR origin, and, considering their morphology, we are confident that three of these (G28.36+0.21, G28.78-0.44, and G29.38+0.10) are indeed SNRs. However, about 25% of the candidates (8 out of 35) have spectral index measurements that indicate thermal emission, and the rest of them are too faint to have a good constraint on the spectral index yet. Conclusions. Additional observations at longer wavelengths and higher sensitivities will shed more light on the nature of these candidates. A simple Monte Carlo simulation reiterates the view that future studies must persist with the current strategy of searching for SNRs with small angular size to solve the problem of the Milky Way’s missing SNRs.


Introduction
The structure formed from the expelled material and the shockwave of a supernova explosion interacting with the surrounding ⋆ The movie associated to Fig. 11 is available at https://www. aanda.org ⋆⋆ Member of the International Max Planck Research School (IMPRS) for Astronomy and Astrophysics at the Universities of Bonn and Cologne.
⋆⋆⋆ Jansky Fellow of the National Radio Astronomy Observatory. interstellar medium (ISM) is known as a supernova remnant (SNR). The interactions of expanding SNRs and the ISM are important feedback mechanisms that may trigger star formation or, on the contrary, disperse gas and thus suppress the star formation rate in a galaxy. Gas can be blown out of the Galactic plane, and turbulent pressure is produced and maintained on both small (molecular cloud) and large (galaxy-wide) scales (e.g. Efstathiou 2000; Ostriker & Shetty 2011;Dubner & Giacani 2015;Bacchini et al. 2020). In order to fully understand and quantify the impact SNRs have on the dynamics of star formation in the Milky Way from an observational point of view, having a complete catalog of Galactic SNRs is highly desirable. The most recent Galactic SNR catalogs (Ferrand & Safi-Harb 2012;Green 2019) contain fewer than 400 objects. This number is, however, significantly smaller than the expected ∼1000 discussed by Li et al. (1991), who arrived at this estimate by using statistical arguments primarily based on our knowledge of the distances to SNRs in the Milky Way. Ranasinghe & Leahy (2022), using a similar statistical analysis but with improved distances to the currently known SNRs, estimate that the number quoted by Li et al. (1991) must be a lower limit, and that there could be over 3000 SNRs in the Galaxy. It is believed that this apparent discrepancy is only due to SNRs that are fainter and smaller than the currently known sample of SNRs not being detected, rather than insufficient knowledge of the local universe (Brogan et al. 2006;Anderson et al. 2017: hereafter A17).
In an attempt to improve this situation, recent radio Galactic plane surveys have been carried out with good sensitivity to both compact and extended emission leading to the identification of well over one hundred SNR candidates (A17; Hurley-Walker et al. 2019;Dokara et al. 2021: D21 from here on). These studies used the radio and mid-infrared (MIR) anti-correlation property of SNRs (Fürst et al. 1987b;Haslam & Osborne 1987). While H II regions emit brightly at both radio and MIR wavelengths, SNRs are typically bright radio emitters on the one hand but weak MIR emitters on the other hand. Fürst et al. (1987b) found that the ratio (R) of 60 µm MIR to 11 cm radio flux density is much higher for H II regions than SNRs (R HII ∼ 1000 and R SNR ∼ 15). Subsequently multiple other studies confirmed this anti-correlation property (Broadbent et al. 1989;Whiteoak & Green 1996;Pinheiro Gonçalves et al. 2011).
Most of these SNR candidates are yet to be confirmed as genuine SNRs with clear nonthermal radio emission. In addition, some objects in the Galactic SNR catalogs either do not have good radio measurements (such as G32.1-0.9 and G32.4+0.1), or, worse, the evidence that they emit nonthermal synchrotron radiation is rather weak (e.g., G31.5-0.6; Mavromatakis et al. 2001). It is not uncommon for H II regions, which emit thermally, to be confused as SNRs due to their similar radio morphology (e.g., A17 and D21). The presence of nonthermal synchrotron radio emission is thus vital to determine whether an object is truly a SNR. Synchrotron radiation is linearly polarized and typically has a negative spectral index at frequencies where synchrotron losses do not occur (typically over 1 GHz; Wilson et al. 2013), where the spectral index α is determined via a power law fit to the flux density spectrum, as S ν ∝ ν α , S ν being the flux density and ν is the frequency. In this work, we focus on confirming the status of the SNR candidates and the sample of objects that were catalogued as SNRs (hereafter called 'known SNRs') in the region of the Galactic longitude range 28 • < l < 36 • and |b| < 1 • (hereafter called 'the pilot region') by measuring linearly polarized flux densities (LPFD) and spectral indices.
The rest of the paper is organized as follows: Sect. 2 contains the descriptions of the data and the methods used for this study. The results for known and candidate SNRs are presented in Sects. 3 and 4 respectively. Their implications are discussed in Sect. 5, and we provide a summary of this work in Sect. 6.

The GLOSTAR survey
The global view on star formation in the Milky Way (GLOSTAR) survey 1 is a 4-8 GHz sensitive, unbiased, large scale continuum and spectral line survey of the first quadrant of the Galaxy, covering the region bounded by the Galactic longitudes −2 • < l < 60 • and Galactic latitudes |b| < 1 • , in addition to the Cygnus X star forming complex (76 • < l < 83 • and −1 • < b < 2 • ). The observations were done using the Karl Jansky Very Large Array (VLA) in B-and D-configurations, as well as the 100 meter Effelsberg telescope. Full details of the observations and data reduction are presented in Medina et al. (2019) and Brunthaler et al. (2021). The catalogs of continuum sources in the GLOSTAR pilot region (28 • < l < 36 • and |b| < 1 • ) from the VLA images, which contain 1575 sources in the D-configuration images including several supernova remnants and 1457 in the B-configuration images, are discussed in Medina et al. (2019) and Dzib et al. (2023), respectively. An overview of the survey and initial results are described in Brunthaler et al. (2021). Further results are presented in D21 (SNRs), Ortiz-León et al.
(2021, Cygnus X methanol masers), Nguyen et al. (2021, Galactic center continuum sources), and Nguyen et al. (2022, methanol maser catalog). Here, we give a brief overview of the data that we use for this study.
We focus on the GLOSTAR pilot region, which contains numerous extended and compact sources overlapping with a strong Galactic background (see Fig. 1 and also Medina et al. 2019). The W43 mini-starburst complex located at l ∼ 30 • (where the bar of the Milky Way meets the Scutum-Centaurus Arm, e.g., Zhang et al. 2014) and the W44 supernova remnant at l ∼ 35 • are among the brightest objects observed in this region. In our previous work (D21), we detected over 150 SNR candidates in the D-configuration VLA images of the full survey, with the pilot region containing 35 of them. We use only the D-configuration VLA and the Effelsberg continuum images in this work.
Since the D-configuration images do not completely recover emission on scales larger than a few arcminutes, they are not suitable to accurately measure the total flux densities and spectral indices of extended objects such as many SNRs. For this purpose, we use images from the single-dish 100 m Effelsberg telescope and their combination with the VLA D-configuration images. The Effelsberg images of this survey do not contain information on the very large scales (>1 • ) due to the baseline subtraction and limited coverage in Galactic latitude (|b| < 1 • ). The large scale information had been 'restored' using the Urumqi 6 cm survey images (Sun et al. 2007(Sun et al. , 2011b to produce the GLOSTAR-Effelsberg maps with the correct intensities (see Brunthaler et al. 2021, for details). However, all the objects we study are smaller than half a degree, and since we need to filter out the large scale background in any case, we use the original maps directly (i.e., without restoration using the Urumqi maps) to avoid a source of uncertainty.
The calibration and imaging of the VLA and the Effelsberg data, along with their feathering for the total power Stokes I images, are described by Brunthaler et al. (2021). Feathering is a method to combine two images with emission from different angular scales, where the two images are co-added in the Fourier domain (uv-space) (Anderson et al. 2014) and the GLOSTAR VLA D-configuration catalog (Medina et al. 2019), are marked using grey circles.
response (e.g., Vogel et al. 1984;Cotton 2017). In this work, we exclusively use the combination of VLA-D and Effelsberg data, and hereafter these images are called 'the combination images'. Since the frequency coverage is not exactly the same on both the VLA and the Effelsberg telescopes, producing the combination images is not straightforward. The final VLA continuum data from the GLOSTAR survey are binned into nine subbands centered on frequencies from ∼4.2 to ∼7.5 GHz, whereas the Effelsberg continuum maps have two subbands centered at frequencies f E,lo ∼ 4.9 GHz and f E,hi ∼ 6.8 GHz (see Brunthaler et al. 2021, for more details). The procedure that we followed to combine the D-configuration VLA and the Effelsberg data for the different Stokes parameters is described below.

Image combination: Stokes I
We average the VLA images from the first five subbands at lower frequencies and the next three subbands at higher frequencies separately to form two VLA images, I V,f V,lo and I V,f V,hi respectively. We discard the ninth subband since it is mostly corrupted by radio frequency interference. The first five subbands have an average frequency f V,lo ∼ 4.7 GHz and the next three subbands have f V,hi ∼ 6.9 GHz, which are already close to the central frequencies of the Effelsberg continuum data, f E,lo ∼ 4.9 GHz and f E,hi ∼ 6.8 GHz respectively, but they are not exactly equal. To bring the two VLA images (I V,f V,lo and I V,f V,hi ) to the frequencies of the Effelsberg images, we use a pixel-by-pixel VLA spectral A145, page 3 of 21 A&A 671, A145 (2022) index α pix to scale the intensities of each pixel: where I pix V,f E,lo and I pix V,f E,hi are the VLA pixel values estimated at the Effelsberg central frequencies. For pixels with intensities above a signal-to-noise threshold of three, α pix is measured from the two Stokes I images I V,f V,lo and I V,f V,hi . For pixels below the threshold, we take a spectral index of zero, i.e., we use the average intensity. After bringing the VLA images to the central frequencies of the Effelsberg images, we feather the VLA and the Effelsberg maps I V,f E,lo + I E,f E,lo to produce the low frequency combination image, and I V,f E,hi + I E,f E,hi to produce the high frequency combination image. Finally, the low and high frequency images are averaged to form the 5.85 GHz GLOSTAR combination image. These combination images will be made available on the GLOSTAR image server 2 before the publication of this work.

Image combination: Stokes Q and U
Similar to the Stokes I procedure, we make the low-and highfrequency VLA images by averaging the first five and next three subbands. We then directly feather each of these averaged images with their respective Effelsberg maps: P V,f V,lo + P E,f E,lo and P V,f V,hi + P E,f E,hi , where P represents Stokes Q or U. We do this without any intensity scaling applied to bring them to the exact same frequency as we did for Stokes I. This is because Stokes Q and U have both positive and negative features, and a direct spectral index calculation is not possible. The Stokes Q and U images at each of the two frequencies are then combined to form the linearly polarized intensity maps Q 2 + U 2 . The low and high frequency maps are then averaged to form the 5.85 GHz GLOSTAR combination image of linearly polarized intensity.
We note that this method may introduce a bias in the measured polarized intensities and the polarization vectors due to the different central frequencies. However, we find that this bias is negligible since the frequencies are quite close ( f V,lo ≈ f E,lo and f V,hi ≈ f E,hi ). Assuming a spectral index of −0.7 for synchrotron emission, the different central frequencies of the feathered VLA and Effelsberg images of linearly polarized emission introduce an error of approximately 4%, which is close to the calibration uncertainty. For the polarization vector to change by just five degrees from f V,lo to f E,lo , the rotation measure must be greater than about 2500 rad m −2 , which is unlikely to be seen in any typical Galactic source. Nonetheless, to introduce this bias in the uncertainty measurement of flux densities and also the instrumental polarization (≲2% in both VLA and Effelsberg data), we adopt a conservative 10% error that will be added in quadrature to the usual uncertainty we obtain from the measurement of flux density of an extended source. In addition, we observe that the LPFD measured in the combination images may be lower than the values measured in the VLA D-configuration only images. This can happen due to the depolarization that occurs when the polarization vectors in the small scale structure detected by the VLA are misaligned with the polarization vectors measured from the Effelsberg data. It is worth noting that, in this study, the exact degree of polarization is not exceptionally important except to the degree it establishes whether the source is or is not polarized, i.e., we only use it as a tool to identify nonthermal emission.

Supernova remnant catalogs
In D21, we presented the list of the SNR candidates that are detected in the D-configuration VLA images of the GLOSTAR survey. It contains 77 objects that were noted as potential SNRs by earlier studies (their Table 3), and 80 new identifications as well (their Table 4). These candidates were identified using the MIR-radio anti-correlation property of SNRs as discussed earlier.
In the GLOSTAR pilot region, there are 21 candidates discovered in the 1-2 GHz HI/OH/Recombination line survey (THOR; A17) and 14 from the GLOSTAR survey (D21). These 35 candidates, in addition to the 12 confirmed SNRs in the Galactic SNR catalogs by Ferrand & Safi-Harb (2012) and Green (2019), are the targets of this study.

Ancillary data
In addition to the GLOSTAR survey continuum images previously described, we also use other complementary radio surveys that are able to recover emissions at the scale of several arcminutes: the 1-2 GHz HI/OH/Recombination line survey (Beuther et al. 2016;Wang et al. 2020) combined with the VLA Galactic Plane Survey (VGPS; Stil et al. 2006), which is called the THOR+VGPS 3 , the 80-230 MHz GaLactic and Extragalactic All-sky Murchison Widefield Array survey (GLEAM; Hurley-Walker et al. 2019) 4 , the Effelsberg 11 cm (∼2.7 GHz) survey of the Galactic plane by Reich et al. (1984) 5 , and the 3 cm (10 GHz) survey of the Galactic plane with the Nobeyama telescope by Handa et al. (1987) 6 .

Flux density and spectral index measurements
We use the GLOSTAR combination images to measure the flux densities at 5.85 GHz, in addition to other surveys as mentioned earlier (Sect. 2.3). We note that we do not measure the flux densities from the two sub-bands to derive an 'in-band' GLOSTAR spectral index, since each of those images depends uponthough only partly -the pixel-by-pixel spectral index from the VLA data, which suffer from the problem of the undetected large scale flux density (see Sect. 2.1.1).
The presence of background emission may bias the value of the measured spectral index. This is particularly true for extended objects in the Galactic plane since the nonthermal Galactic background is strong and ubiquitous at low radio frequencies. In addition, the intensity of this background is dependent on frequency and position (e.g., Paladini et al. 2005). The method of 'unsharp masking' (Sofue & Reich 1979) is generally used to filter out the large scale Galactic emission, but it is not appropriate for smaller scale background emission across an object with the size of a few arcminutes. In this work, we fit a 'twisted plane' that removes the background contribution up to a first order variation. Points are chosen around an object such that they represent the background emission in that area, and a two-dimensional least-squares linear fit is performed to the pixel intensities to measure the background variation. The uncertainty from this background subtraction operation is determined by choosing multiple sets of vertices. We subtract the local background in both the total intensity and the polarized intensity images, and we mask pixels typically below a 3σ-level, where the noise is determined locally by a sigma-clipping algorithm.
While several objects we discuss in this paper already have their low frequency flux densities derived in multiple previous studies, for the sake of consistency with regards to spectral index, we make our own measurements of the flux densities of these objects using the images directly from their survey data, performing background subtraction in the same manner as we do for the GLOSTAR images. We also mask the point sources that are clearly unrelated (e.g., Tian & Leahy 2005) to keep the measurement as accurate as possible. In addition, since radio interferometric artefacts such as radial 'spokes' are common near bright sources, we do not measure the total or linearly polarized flux densities if we are unable to disentangle such effects from the emission of an object. Due to such artefacts, polarization measurements are not possible for about a third of the objects studied in this work.
The spectral index of an object is usually measured by fitting a straight line to the relation between flux densities and frequencies in logarithmic space: However, the values determined in this manner are sensitive to the presence of background emission. Turtle et al. (1962) introduced the concept of temperature-temperature (TT) plots, in which a spectral index is extracted from the slope of a straight line fit to the pixel intensities at one frequency against the pixel intensities at another frequency. In essence, we integrate over the whole area to measure the flux density spectral index (α FD ), whereas the TT-plot spectral index (α TT ) is calculated by measuring the variation of each pixel at different frequencies.
The intensities on TT-plots can be represented by brightness temperatures in Kelvin, or pixel intensities in Jy beam −1 . In this work, we exclusively use pixel intensities, and the spectral index is calculated using: where m S is the slope of the line that is fit to pixel intensities. This is a more reliable measurement of spectral index of an extended object because the flux density bias introduced by a constant large scale background emission moves all the points equally, and hence does not affect the slope of the fit. Since the combination images are produced using the spectral index derived from the D-configuration GLOSTAR-VLA images, they are not suitable to measure the TT-plot spectral index (α TT ). We only use the GLOSTAR-Effelsberg images for this purpose. We also measure the flux density spectral index (α FD ); this serves as a useful consistency check since we subtract the background regardless, as described above. We note that, at low radio frequencies such as the regime of the GLEAM survey (≲200 MHz), absorption effects become important, either via synchrotron self-absorption or free-free absorption (e.g., Wilson et al. 2013;Arias et al. 2019). This lowers the emitted flux at low frequencies and increases the power-law spectral index compared to values determined at higher frequencies. Such a 'spectral break' effect had been noted in several SNRs before (e.g., Sun et al. 2011a). Spectral breaks are also observed in pulsar wind nebulae due to the central pulsar's time-dependent energy injection (Reynolds & Chevalier 1984;Woltjer et al. 1997). When we calculate the flux density spectral index in this work, if we clearly see a break, we split the spectrum into two and calculate two spectral indices; if the break is not obvious, we calculate only a single spectral index.
The nonthermal emission from the Galactic disk is polarized, and it may have structure on small scales that is not filtered out by an interferometer. While this is more significant at longer wavelengths, it might affect the GLOSTAR images too (see D21). We verify that in the objects we study in this work, there exist no features with no Stokes I counterparts when measuring LPFD. In addition, a Ricean polarization bias might introduce a positive offset. This occurs because LPFD is the square root of the sum of squares ( Q 2 + U 2 ), and any positive or negative noise in Q and U will always add up and result in a non-zero LPFD. We find that this effect is an order of magnitude smaller than the flux densities we report, and in fact there is no need to correct for this bias due to the background subtraction procedure and the 3σ-level mask we use (see Wardle & Kronberg 1974). Nonetheless, the twisted-plane background subtraction procedure is applied to the linearly polarized intensity images as well, which accounts for the Galactic plane polarized background and also any possible Ricean polarization bias.

Known SNRs
The Galactic SNR catalogs of Green (2019) and Ferrand & Safi-Harb (2012) list 12 confirmed SNRs in the region we study.
In Table 1, we present the GLOSTAR 5.8 GHz integrated flux densities of these SNRs along with their spectral indices. If possible, overlapping H II regions and clearly unrelated point sources are masked while measuring the flux densities, taking special care in crowded regions. If it is unclear whether a particular region of emission belongs to the SNR, we do not remove that region. We find that the flux densities and spectral indices are generally consistent with previous studies. We present the GLOSTAR combination images of some interesting known SNRs and discuss them below. The total intensity images and the linearly polarized intensity images of all the known SNRs studied in this work are shown in Appendix A.
3.1. G29.6+0.1 While we had already detected linear polarization in the SNR G29.6+0.1 using the VLA D-configuration images (in D21), the emission in the combination images seems to be depolarized due to the addition of large scale information from the GLOSTAR-Effelsberg data. We do not measure any polarized emission over a 1σ upper limit of ∼0.1 Jy. The flux densities we measure (see Table 1) appear to be lower than what is expected from the lower limits reported by Gaensler et al. (1999): ∼0.41 and ∼0.26 Jy at 5 and 8 GHz, respectively. The reason for this inconsistency is unclear. Nonetheless, the broadband spectral index we derive from our measurements (approximately −0.5) is in line with the TT-plot spectral indices derived by Gaensler et al. (1999). We show the GLOSTAR combination image of the SNR G29.6+0.1 in Fig. 2. The spectrum of this SNR shows that it might be falling more rapidly from 1.4-5.8 GHz than from 0.2-1.4 GHz, suggesting the presence of a spectral break around 1 GHz. But given the uncertainties, we reserve judgment on the changing spectral index.
this is a SNR-H II region complex. The Stokes I flux densities we measure are consistent with those given by Fürst et al. (1987a) within uncertainties, and we also find a morphology similar to their image. However, the spectral index we derive from 200 MHz to 10 GHz is essentially zero, which is consistent with our TT-plot result (Fig. 2), but in slight tension with the value of approximately −0.2 given by Fürst et al. (1987a). Even after separating from the region the thermal emission that they reported, we find no evidence for synchrotron emission. In the 24 µm images of MIPSGAL (Multiband Infrared Photometer for Spitzer Galactic plane survey; Carey et al. 2009), we find weak emission following the radio morphology, hinting that the emission may be thermal. Based on sulfur and Hα optical lines, Mavromatakis et al. (2001) also suggest that this may be an H II region instead of a SNR. High resolution deeper observations at lower frequencies will shed more light on the nature of the emission from this object, but the evidence so far suggests that G31.5-0.6 is not a SNR.
3.3. G32.1-0.9 Folgheraiter et al. (1997) discovered the SNR G32.1-0.9 in the X-ray regime, and they found a possible faint radio counterpart in the 11 cm Effelsberg images. A17 reported a possible detection in the THOR+VGPS data too, but no radio spectral index was ever determined. While we cannot confidently identify any counterpart in the GLOSTAR data, the 200 MHz GLEAM image shows a shell that corresponds to the 11 cm Effelsberg and THOR+VGPS detections (Fig. 3). Using these three images, we derive a radio spectral index for this unusually faint SNR for the first time, α FD = −0.68 ± 0.12. Its average 1 GHz surface brightness is approximately 3 × 10 −22 W m −2 Hz −1 sr −1 , which makes it one of the faintest radio SNRs currently known: it is only three times brighter than the faintest SNR known in the Milky Way (G181.1+9.5; Kothes et al. 2017).
3.4. G32.4+0.1 G32.4+0.1 was discovered in the X-ray regime by Yamaguchi et al. (2004), who also noted a possible counterpart in the images of the 1.4 GHz NRAO VLA Sky Survey (Condon et al. 1998).
The radio emission from this SNR is faint but clearly visible in the GLEAM, the THOR+VGPS and the GLOSTAR combination images, allowing us to measure, for the first time for this SNR, a spectral index of −0.21 ± 0.07 (from brightness values) to −0.39 ± 0.10 (from a TT-plot). The GLOSTAR combination image and the plots for spectral index determination are shown in Fig. 2. As noted in Sect. 2.4, the low frequency emission detected in GLEAM may be self-absorbed which brings the spectral index close to zero; hence we favor the TT-plot spectral index (approximately −0.4) for higher frequencies where the effects of synchrotron self-absorption are not important. Linear polarization is undetected, with an upper limit on the linearly polarized flux density of ∼0.3 Jy.
3.5. G32.8-0.1 Green (2019) lists the SNR G32.8-0.1 with an uncertain spectral index of −0.2 based on the work of Caswell et al. (1975), who report a flux density of 12.8 Jy at 408 MHz. Unfortunately, no uncertainties were quoted, but they reported that their error might be large. Later, Kassim (1992) observed this SNR with the VLA at a similar frequency of 330 MHz, and their results are in dispute with the result from Caswell et al. (1975). They measured a significantly higher flux density of ∼32 Jy and consequently a more negative spectral index of approximately −0.5, but no uncertainties were quoted once again. This SNR is clearly visible in the GLOSTAR survey, in addition to the GLEAM and the THOR surveys (Fig. 2), which helps us resolve the tension. Our measurements of flux density (16.3 ± 1.7 Jy at 200 MHz) and spectral index (α TT = −0.27 ± 0.04) are consistent with the values given by Caswell et al. (1975), which is confirmed by our TT-plot spectral index (α TT = −0.32 ± 0.05).
3.6. G35.6-0.4 The nature of emission from G35.6-0.4 had been a subject of discussion for a long time. It was included in the early SNR catalogues (e.g., Downes 1971;Milne 1979 radio recombination line by Lockman (1989) among other studies, had cast doubts that the emission is nonthermal (see Green 2009, for an overview). Finally, using higher quality radio continuum data, Green (2009) "re-identified" this as a SNR with a spectral index of approximately −0.5. This object is clearly visible in the GLOSTAR survey, where we also unambiguously detect linearly polarized emission (Fig. 4). Its spectrum appears to be broken; from 200 MHz to 1.4 GHz the flux density has no significant change (α ∼ 0), and from 1.4 to 10 GHz it falls with a spectral index of α = −0.31 ± 0.07. This is confirmed with the GLOSTAR Effelsberg TT-plot spectral index as well (Fig. 4). This result ( index derived by Rennie et al. (2022): α = −0.34 ± 0.08 from 2.7-30 GHz. Green (2009) derives a slightly steeper spectral index of −0.47 ± 0.07. This is probably because of the different choice of polygons used for measuring the flux density and also the subtraction of background emission in this complex region, but we note that the values are consistent within 2σ.
Given the presence of radio recombination lines that indicate thermal emission and a spectral index approximately −0.3, this region appears to be a complex of thermal and nonthermal emissions. Paredes et al. (2014) suggest that there may be two circularly shaped extended objects present in this complex (marked by two red dotted circles in Fig. 4), and with one of them with thermal and the other one with nonthermal emission. We find that MIR emission is detected from the southern part in the GLIMPSE and MIPSGAL images (Churchwell et al. 2009;Carey et al. 2009), providing further evidence of thermal emission from this region. The linearly polarized emission detected in the GLOSTAR combination data (see Fig. 4) also hints at the presence of two shells, one centered at G35.60-0.40 and the other at G35.55-0.55, similar to those reported by Paredes et al. (2014). However, since we find polarization from both these regions, it is likely that emissions from these regions have both thermal and nonthermal components.

Candidate SNRs
In the pilot region, we had discovered 14 new candidate SNRs from the GLOSTAR survey in our previous work (D21), in addition to the 21 candidates discovered by A17 using THOR+VGPS images. The continuum flux densities of these candidates are presented in Table 2 (from THOR+VGPS) and Table 3 (from GLOSTAR). We derived flux density spectral indices whenever possible, and these are plotted in Fig. 5. We discuss five objects for which there is good evidence of nonthermal emission in detail in the following sections. We also find that 14 other candidates possibly have a negative spectral index. But since they are quite faint and the morphology of these candidates is not clear (see Figs. B.1 and B.2), we do not discuss them further.      We detect this object in the Stokes I images of our GLOSTAR survey (Fig. 6), with the same morphology as observed in the THOR+VGPS images. Its fractional polarization is about 2%, which is not unusual in SNRs (e.g., Sun et al. 2011a). The linearly polarized intensity map from GLOSTAR shows a faint structure, close to the noise level in this region, that resembles the total intensity of the shell of this object. From the Effelsberg images A145, page 10 of 21  of our survey, we made a TT-plot and obtained a spectral index of −0.33 ± 0.14. By measuring the background-subtracted flux densities in the images of GLOSTAR combination, THOR+VGPS and GLEAM, we obtain a brightness spectral index of −0.28 ± 0.11. These measurements and the morphology we observe in the total and linearly polarized intensity images provide ample evidence of nonthermal emission from this object, and hence we conclude that G28.36+0.21 is indeed a SNR.

G28.78-0.44
The candidate SNR G28.78-0.44 (Fig. 7) had previously been identified in the MAGPIS and the THOR+VGPS surveys (Helfand et al. 2006;A17). Hurley-Walker et al. (2019) derive a spectral index of approximately −0.7 in their GLEAM survey (70-230 MHz), consistent with the spectral index from the TIFR-GMRT Sky Survey and the NRAO VLA Sky Survey (de Gasperin et al. 2018;Dokara et al. 2018). While the polarization from this object was already clearly visible in the VLA images of the GLOSTAR survey (D21), the addition of the Effelsberg data allows us to measure its flux densities at 5.8 GHz. The fractional polarization we measure in the combination images is about 4%. We also detect this object in the Effelsberg 11 cm survey (Reich et al. 1984) and the Nobeyama 10 GHz survey (Handa et al. 1987). These give us a broadband flux density spectral index of −0.42 ± 0.04, which is consistent with the TTplot spectral index from the Effelsberg images of the GLOSTAR survey alone (−0.52 ± 0.12, see Fig. 7). Thus we find strong evidence that this filled-shell object is a SNR.

G29.38+0.10
This source appears to have a complex structure with a bright pulsar wind nebula (PWN) and a faint SNR shell in the GLOSTAR combination image (Fig. 8). The central structure of this complex is bright and highly polarized in the combination images, with the degree of linear polarization reaching as high as 30% in some pixels. For the whole complex, this value is 5.5 ± 0.8%. We had detected strong linear polarization from this object in our previous work as well (D21), which was based only on the D-configuration VLA images. Its low frequency spectral index measured using the GLEAM images by Hurley-Walker et al. (2019) for the whole complex, and for the central source by Dokara et al. (2018) using the TGSS-NVSS spectral index map (de Gasperin et al. 2018) is approximately zero, which is typical of PWNe. We calculate a similar spectral index using the THOR+VGPS and GLEAM images as well. However, between the THOR+VGPS and the GLOSTAR combination images, the flux density falls with a spectral index of α FD ∼ − 0.34. Constructing a TT-plot using images from the two bands of the GLOSTAR-Effelsberg data, we measure a value α TT ∼ − 0.46. This implies that there is a break in the spectrum of this source near 2 GHz, or a gradual turnover. Such a varying spectral index at these frequencies is once again typical of PWNe (see Pacini & Salvati 1973;Reynolds & Chevalier 1984;Sun et al. 2011a). These facts provide further evidence that G29.38+0.10 is a PWN+SNR shell complex. 4.4. G031.256-0.041 A17 cataloged G31.22-0.02 as a shell-shaped SNR candidate based on the THOR+VGPS images. It lies in a crowded field with a strong background, due to which the determination of the TT-plot spectral index from the Effelsberg images (α TT ) was not possible. This region is better resolved in the GLOSTAR combination images, in which we identify the brightest part of the supposed shell of G31.22-0.02 (at l ∼ 31.26 • , b ∼ −0.02 • ) to be inside another shell (Fig. 9). We believe that this is a PWN+shell complex, and named it as a GLOSTAR SNR candidate G031.256-0.041 in our previous work (D21). The flux densities we measured in the THOR+VGPS and the GLOSTAR combination images are similar within uncertainties (S ∼ 0.35 Jy), giving a spectral index close to zero between 1.4 and 5.8 GHz. In the 200 MHz GLEAM images, what we believe is the center of the PWN (at l ∼ 31.26 • , b ∼ − 0.02 • ) is barely resolved, with a peak brightness of nearly 0.8 Jy beam −1 . The background level in this region is about 0.4 Jy, implying that the flux density of the peak is ∼0.4 Jy, similar to the flux densities from the GLOSTAR combination and the THOR+VGPS images. Unfortunately, the linearly polarized intensity images from GLOSTAR in this region are contaminated with sidelobe artefacts of nearby bright sources, prohibiting us from measuring its degree of polarization. The morphology and the estimated spectral index are, however, consistent with our PWN+SNR shell interpretation.

G034.524-0.761
We discovered the SNR candidate G034.524-0.761 in our previous GLOSTAR work, where we had identified clear linear polarization from the VLA data (see Fig. 11 of D21). With the addition of the Effelsberg data to the VLA images, we now obtain a degree of polarization ∼10% from this candidate. In addition, we obtain a TT-plot spectral index of approximately −0.6 using the Effelsberg images, although with a large A145, page 12 of 21 R. Dokara et al.: A&A proofs,   uncertainty of ∼0.5. We measured flux densities in the 200 MHz GLEAM and the 1.4 GHz THOR+VGPS images, which give us a spectral index of approximately −0.9. While all these facts point to a nonthermal origin of the emission from this region, the morphology of this candidate (Fig. 10) indicates that this might be a filament. For this reason, we cannot conclude that this object is a SNR.

Discussion
It is evident from Fig. 5 that the spectral indices of several SNR candidates are not well constrained yet. Most of them have a small angular size and a low surface brightness, and they lie in crowded regions with a strong background; these conditions result in large uncertainties in the measurement of their spectral indices. Moreover, the polarization signals from several SNRs may remain undetected because of limited sensitivity (the linearly polarized flux density is typically only a few percent of the total flux density, e.g., Sun et al. 2011a). Deeper observations of these candidates across the radio band are necessary to constrain their spectral indices and linear polarization better. However the current results do not look very promising since the rate of confirmation appears to be quite low, and we are forced to ponder over the strategy to identify new SNRs.
Since most of the bright SNRs are likely to have been discovered already, it might progressively get more difficult to find the remaining fainter ones. H II regions are more numerous in the Galaxy, and there is a chance that the fainter H II regions contaminate the sample of the faint SNRs. However, the SNR candidates identified by A17 and D21 do not have any significant coincident MIR emissions detected in the Spitzer MIR surveys, which can detect H II regions anywhere in the Galaxy (Anderson et al. 2014). Hence, we believe that, if the SNR candidates do not turn out to be SNRs, the confusion must be due to radio emitters other than H II regions, although it is unclear what kind of objects they might be. A17 and D21 suggest that the remaining undetected SNRs must be faint and also have a small angular size. We turn our attention toward these properties of the sample of the SNR candidates.

Angular radius
One question that needs to be answered before starting the search for the remaining SNRs is whether most of them are indeed small, since that would determine what resolution is necessary to detect the 'missing' SNRs. To estimate their apparent angular extents, we ran a simple Monte Carlo simulation of evolution of SNRs in the Milky Way. SNRs are evolved in a locally uniform ISM using the expressions from Draine (2011), which are based on the four classical stages as proposed by Woltjer (1972): 1. The earliest part of the evolution is known as the freeexpansion or the ejecta-dominated phase. We assume that the mass of the swept up ISM (m sw ) is negligible compared to the mass of the SN ejecta (m ej ) in this stage. 2. Sedov-Taylor phase begins when the shocked and swept up mass is comparable to the ejecta mass m sw ∼ m ej , during which the explosion can be approximated as a point source injecting only energy. 3. Snowplow phase begins when the radiative cooling losses become important and the matter behind the SNR shock cools rapidly to form a cold and dense shell. In the hot and tenuous medium that is interior to the shock, however, the energy losses do not yet play a role, and the pressure from this hot central volume drives the momentum of the dense outer shell. 4. The final phase is 'dispersion' as the SNR merges into the surrounding ISM and fades away when the shock speed drops to the ambient velocity dispersion levels. We derive the radius of each SNR based on the time since explosion and the position in the Galaxy. Following are the main parameters and inputs of the simulation: -Galactic supernova rate of one per 40 yr, with the corecollapse and thermo-nuclear types being 85 and 15% respectively (Tammann et al. 1994;Reed 2005). -Three dimensional gas density model of the Milky Way from Misiriotis et al. (2006). -A random Monte Carlo model of the two-dimensional distribution of supernova events in a disk with a central hole and a two-arm spiral following Li et al. (1991). The central hole is to account for the dearth of massive star formation, and by extension SNRs, near the Galactic center (see Nguyen et al. 2021;Ranasinghe & Leahy 2022, for example). -Core-collapse SN events, which trace massive star formation, are chosen to have a scale height of 80 pc, the same as the scale height of the molecular gas (from Misiriotis et al. 2006). -Type Ia SNe arise due to mass accretion onto old degenerate stars; accordingly we use the thick disk scale height of 0.7 kpc from Kordopatis et al. (2011). -The maximum lifetime of SNRs is fixed at 80 000 yr (Frail et al. 1994). -For a Type Ia supernova, the kinetic energy of the ejecta is fixed at 10 51 erg and the ejecta mass is normally distributed from 0.8 M ⊙ -1.8 M ⊙ (following Scalzo et al. 2014). -For the more numerous core-collapse supernova events, the ejecta mass (8 M ⊙ -11 M ⊙ ) and the kinetic energy (0.2-1.3 times 10 51 erg) are randomly drawn from distributions adapted from the results of Martinez et al. (2022). There are, however, some caveats to consider: -Realistically, the properties of the ISM are not smoothly varying functions of position as the model given by Misiriotis et al. (2006). The ISM number density can drastically change depending on the environment, especially in the case of previous mass-loss events such as stellar winds. These affect the evolution of SNRs in a crucial and nontrivial manner (e.g., Yasuda et al. 2021). -The distribution of supernova events follows the model of Li et al. (1991), which is quite simplistic. But similar to their findings, we also observe that the results are insensitive to parameters of the disk and the spiral arms. The inverse dependence of angular radius with distance makes our result even more robust than that of Li et al. (1991). -The distributions of ejecta mass we used (from Scalzo et al. 2014;Martinez et al. 2022) may not hold for the Milky Way accurately, since those results are from the nearby local universe with supernovae from several galaxies. However, we find that even if the ejecta mass for core-collapse supernovae was only 1 M ⊙ instead of 8 M ⊙ -11 M ⊙ , the results are mostly the same. -There is evidence that the explosion energies of supernovae can have a range wider than that we have taken, for both Type Ia and core-collapse, from ∼10 49 to ∼10 52 erg (e.g., Benetti et al. 2005;Fisher & Jumper 2015;Pejcha & Thompson 2015;Murphy et al. 2019;Leahy et al. 2020). Even with a wider range, we find that the resultant radius distribution does not significantly change. -We do not take into account the effects of clustering. This is the main drawback of this simulation. A significant fraction of massive star formation -and the number of SN events by extension -happens in clusters (e.g., Krumholz 2014). Ferrière (2001) estimates that ∼60% of O stars probably remain in their natal group, while the rest of them end up in the 'field'. If multiple supernovae occur in succession in such clusters, this might result in the formation of a super-bubble (e.g., Ehlerová & Palouš 2013). We ran the simulation for two million years, which is several generations of SNRs. A snapshot at a time of 1.8 Myr is presented in Fig. 11, and a movie of the whole two million years is available online. Given that the lifetime of a SNR and the SN rate are fixed at 80 000 yr and one for every 40 yr respectively, about 2000 SNRs exist at the end of the simulation. It is clear that most of the SNRs are quite small with angular radii of only a few arcminutes, similar to the THOR and GLOSTAR SNR candidates. Even if the lifetime of a typical SNR is longer than 80 000 yr as we had used, the resultant distribution does not shift to higher angular scales significantly. This is due to the fact that the expansion is considerably slower in the later stages of SNR evolution. While this simulation only serves as a first A145, page 15 of 21 A&A 671, A145 (2022) approximation since we do not consider several effects such as those mentioned above, it is nevertheless useful to give us an idea of what to expect. And the result reiterates the views of A17 and D21 that SNR searches must focus on small angular sized objects to make the most gains.

Radio surface brightness
In the simulation described above, we also measured the area of overlap of SNRs. We find it to be typically less than 10% of the total sky area covered by SNRs, suggesting that the confusion due to SNRs overlapping themselves may not be important. However, the SNRs originating from core-collapse events are located near massive star forming complexes, which also contain other extended structures emitting at radio wavelengths. H II regions are the most likely sources of positional overlapping confusion: they are probably over 8000 in number (Anderson et al. 2014), and the range of the values their radio surface brightness is similar to that of SNRs.
Currently, the faintest SNR known has a brightness temperature of about 0.33 K at 1 GHz (Kothes et al. 2017), and, by extrapolating to 1 GHz assuming a nonthermal spectral index, we find that the SNR candidates from A17 and D21 are at a similar or lower surface brightness. On the other hand, the background emission from the diffuse gas in the Milky Way is at a level of a few Kelvin in the inner Galactic plane at 1 GHz (e.g., Reich et al. 1990), and it is even higher in regions such as the mini-starburst W43 where one expects many SNRs due to recent massive star formation activity. This implies that the diffuse background emission is a critical source of confusion, and finding new SNRs will probably be more difficult from now on. Interferometric surveys at lower frequencies, such as MeerKAT, appear promising in the search for new SNRs (e.g., Heywood et al. 2022), but the nonthermal Galactic background emission is also stronger at lower frequencies and may contribute to the confusion.

Summary and conclusions
We derived spectral indices of previously confirmed SNRs in the Galactic longitude range 28 • < l < 36 • , using the VLA-D+Effelsberg combination images of the 4-8 GHz GLOSTAR survey in addition to other complementary and archival survey data. These include the first radio spectral index determinations for SNRs G32.1-0.9 and G32.4+0.1, along with the first reported spectral break for SNR G35.6-0.4. We showed that G31.5-0.6 may not be a SNR, and we provided further evidence of nonthermal emission from the SNR candidates G28.36+0.21, G28.78-0.44, G29.38+0.10, and G034.524-0.761. We find that G28.36+0.21 and G28.78-0.44 are typical SNR shells, and G29.38+0.10 is a PWN+shell complex. Based on a simple Monte Carlo simulation of SN events in the Milky Way, we find that most of the SNRs yet to be discovered must have angular sizes smaller than half a degree. Hence, despite the low rate of confirmation, we believe that future studies must focus on small angular sized objects such as the THOR and GLOSTAR SNRs. The forthcoming Effelsberg images from the GLOSTAR survey for the rest of the coverage will be analyzed in the coming months, which will undoubtedly help us study more SNRs and candidates in the near future.