© ESO 2018

1. Introduction

The radio emission from normal star-forming galaxies traces the underlying distributions of thermal and relativistic plasmas, cosmic ray (CR) electrons, and magnetic fields, thus providing vital information about the physical processes at work in galaxies. Studying the radio emission at different frequencies via the radio continuum spectra of galaxies allows us to understand radio emission processes and the structure and local properties of the galaxy’s interstellar medium (ISM).

The shape of radio continuum spectra can be characterised to first order by their power-law spectral index α (S_ν ≈ ν^α), the value of which can be related to the various radiation processes responsible for the emission. An optically thin plasma yields α = −0.1 for thermal bremsstrahlung radiation, while synchrotron radiation gives α ≈ −0.5 for freshly accelerated CR electrons injected from supernova remnants into the star-forming disk. These CR electrons can sustain considerable synchrotron and inverse Compton radiation losses, giving rise to steeper spectra, with α also dependent on the structure and strength of the magnetic field and the confinement of CRs (Beck & Wielebinski 2013; Han 2017).

The transport of CR electrons away from supernova remnants can take the form of either diffusion, which depends on the magnetic field structure, or advection in a galactic wind (Pohl et al. 1991; Heesen et al. 2016, 2018). Therefore, galaxy spectra, particularly the integrated (global) ones, depend on a complex interplay between thermal and nonthermal components, CR electron energy losses, and propagation effects (Lisenfeld & Völk 2000).

At low radio frequencies, the spectra of galaxies are expected to be modified by additional mechanisms. For instance, H II regions become optically thick at low frequencies, which affects not only the propagation of free-free emitted photons but also the transmission of photons generated from synchrotron emission. At low frequencies, relativistic bremsstrahlung and especially ionisation losses can also be much more important than at higher frequencies, in particular in starburst galaxies (Murphy 2009). Finally synchrotron self-absorption and Razin effects can further suppress the radio emission from regions of dense ISM and produce breaks in the galaxy spectra below 10 MHz (Lacki 2013).

Our current understanding of low-frequency spectra of galaxies and the processes that shape them is limited, mostly due to the lack of high-quality radio observations at these frequencies. In the well-known work of Israel & Mahoney (1990), flux densities of 68 galaxies were observed with the Clarke Lake Telescope at 57 MHz and compared with flux densities extrapolated from high-frequency measurements. The authors found that flux densities at 57 MHz were systematically lower than expected, with the difference largest for highly inclined objects, which they interpreted as evidence that the low-frequency flattening was due to free-free absorption from thermal gas concentrated in the plane of galaxy disks. A new ISM phase was proposed, containing an ionised but relatively cool (T <  1000 K) gas, to account for this absorption. However, other authors have subsequently come to different conclusions: for example Hummel (1991) after reanalysing the same data, confirmed the reduction of radio emission but questioned whether this flattening is correlated with galaxy inclination. The authors instead proposed that the observed spectral breaks were due to the steepening of high-frequency spectra caused by energy losses of the CR electrons propagating within the galaxy. Due to the difficulties of obtaining high-quality low-frequency measurements, the discussion on the shape of galaxy spectra and the role of thermal absorption has not yet been resolved (Marvil et al. 2015; Basu et al. 2015; Mulcahy et al. 2018).

It is clear that thermal absorption is required to account for the observed spectra of local regions in the centres of galaxies, including the Milky Way (e.g. Roy & Pramesh Rao 2006). Such an effect can also clearly be seen in the centre of M 82 at 408 MHz (Wills et al. 1997) and at 150 MHz (Varenius et al. 2015) and likely also affects the integrated spectrum of M 82, being visible as a weak flattening of the spectrum at frequencies below 300 MHz (Condon 1992; Yoast-Hull et al. 2013; Lacki 2013; Adebahr et al. 2013). The Milky Way also shows a spectral turnover of the global spectrum, but at a much lower frequencies of about 3 MHz (Brown 1973). The frequency of such spectral breaks is related to the amount of ionised gas, specifically of the warm ionised medium, and therefore also to the age of star-forming regions, with recent star formation on a 10 Myr timescale being most important. Complex free-free absorption features can serve as an indicator of an early evolutionary state of a starburst, as shown by Clemens et al. (2010).

Free-free absorption can also strongly affect the relation between radio and far-infrared (FIR) emission (Schleicher & Beck 2013). For the same reason it can modify the radio emission of high-redshift galaxies, thus influencing the source counts. Therefore, a better understanding of the role played by free-free absorption in galaxies is essential for studies of cosmological galaxy evolution.

A new observational facility, the Low Frequency Array (LOFAR; van Haarlem et al. 2013), opens up the possibility for systematic studies of nearby galaxies at low frequencies and allows us to reinvestigate the problems related to their low-frequency spectra. The Multifrequency Snapshot Sky Survey (MSSS; Heald et al. 2015) covers the entire northern sky, enabling the detection of many catalogued nearby galaxies, which span a large range of star formation rates (SFRs), sizes, and morphological types.

In this paper, we determine flux densities at 150 MHz from both the MSSS source catalogue and image mosaics for a large number of galaxies (Sect. 2). We then obtain global spectra of galaxies, combining our new MSSS flux densities with flux densities from the literature (Sect. 3). A three-dimensional (3D) model of galaxy radio emission is introduced (Sect. 4) for the purpose of analysing and interpreting these spectra. As cases in point, we analyse the spectra of M 51 and M 82 in detail, using spatially resolved information on the observed radio emission at different frequencies to constrain our model parameters. Synthetic radio intensity maps are then produced for various galaxy inclinations at different frequencies, which enable us to predict global spectra of galaxies and assess the role that free-free absorption plays in shaping them. In Sect. 5 we provide a summary and conclusions.

2. Data and galaxy selection

2.1. Selection of the galaxy sample

As the parent sample for our study we selected the compilation of Yun et al. (2001), which contains radio counterparts to the IRAS Redshift Survey galaxies detected in the NRAO VLA Sky Survey (NVSS; Condon et al. 1998). The catalogue lists over 1800 IRAS flux-density-limited (S ^60 μm ≤ 2 Jy) objects with known radio properties at 1.4 GHz, and constitutes the largest sample of this type within the local Universe. It has been used to investigate the radio-luminosity function of galaxies, the radio-FIR correlation, and the extinction-free star-formation density for the local volume (Yun et al. 2001; Condon et al. 2002).

The sample is not complete at low Galactic latitudes of |b_Gal| < 10°. Therefore, we included a number of well-known galaxies from similar studies by Condon (1987) and Condon et al. (1990), namely: IC 10, NGC 628, UGC 12914, NGC 3646, NGC 4217, NGC 4449, and NGC 5457.

In analysing galaxy spectra to study the role of thermal absorption, we were primarily interested in galaxies with Hubble types corresponding to spiral and irregular objects, galaxies that are known to have recent star formation. In order to have a uniform sample of such galaxies, we chose the following selection criteria:

• included in Yun et al. (2001), supplemented with the above seven additional galaxies;

• radio flux density > 50 mJy at 1.4 GHz;

• located in the northern hemisphere (Dec > 0°);

• morphological type T ≥ 0 (according to the HyperLeda database) to avoid elliptical galaxies;

• not dominated by an active galactic nucleus (AGN), for example, excluding NGC 1275 and NGC 4258, but taking into consideration M 51 and M 81 (with a low-luminosity LINER/Seyfert nucleus).

2.2. MSSS survey

We used the High Band Antenna (HBA) part of MSSS to measure low-frequency flux densities of our sample galaxies. Images from MSSS are available as sets of mosaics representing sky images of 10 ° ×10° in size. Each set consists of eight narrowband images generated with 2 MHz bandwidth each, and with central frequencies ranging from 120 to 157 MHz, as well as one broadband image obtained by averaging images from all bands.

We used the preliminary MSSS mosaics from the first internal version of the MSSS source catalogue (e.g. see details in Kokotanekov et al. 2017). The typical resolution of the mosaics is 3′ and the root mean square (rms) noise level varies between about 15 and 25 mJy beam⁻¹. Differences in noise levels may result from fluctuations of the Milky Way foreground emission, distributions of strong sources in mosaics, imperfect calibration and ionospheric weather conditions during observations. The catalogue of automatically detected sources (the MSSS source catalogue; Heald et al. 2015) provides us with flux densities from the individual bands. The MSSS images, and thus the catalogued flux densities, were corrected as part of the MSSS analysis to mitigate the well-known LOFAR flux calibration transfer issues using the “bootstrap” method described by Hardcastle et al. (2016); the residual flux calibration error should not exceed 10%. A polynomial function was used with a Levenberg–Marquardt minimisation algorithm in order to determine the best-fit parameters of the multi-band spectrum in the logarithmic flux density-frequency space. The catalogue provided us with the parameters of the fit as well as the interpolated flux density of each galaxy at 150 MHz used by us in this work.

2.3. Flux density measurements

Apart from the MSSS source catalogue, we also used our own flux density measurements based on the MSSS image mosaics; for these measurements we used the Common Astronomy Software Applications package (CASA; McMullin et al. 2007). We applied a polygon mask to outline galactic emission in the averaged mosaic. This was then used to integrate flux densities in the eight individual mosaics of the various bands. The uncertainties of the flux densities were calculated as the local rms sensitivities σ_rms measured in the maps multiplied by the square root of the flux integration area Ω_s in units of the beam area Ω_b:

(1)

Finally, we employed a Levenberg–Marquardt minimisation method to fit the multi-band power-law spectra to our measurements and to obtain the interpolated flux densities at 150 MHz S₁₅₀ together with their uncertainties σ_L − M. We added an independent uncertainty of 5% of the total flux density, stemming from the uncertainty of the absolute flux density scale and from calibration errors:

(2)

Uncertainties of the flux densities from the MSSS source catalogue were calculated in a similar way.

The flux densities from CASA were generally similar to those from the MSSS catalogue, although we noticed that in some cases the CASA flux densities were slightly smaller, but still in agreement within the uncertainties. We attributed this to different ways of masking the area during the flux density integration: for the MSSS catalogue, each source comprises one or more Gaussian components identified within the boundaries of an automatically determined PyBDSF “island mask”, whereas in CASA we manually specified polygon regions. The two distinct approaches caused slightly different pixel areas to be considered, with generally minimal discrepancy. If the flux densities from both methods were consistent, we used the catalogue flux density for further analysis. If, however, the flux densities differed significantly, for example due to a complicated source structure or due to a confusing background source close in projection to the target galaxy, we used the CASA measurements instead and applied a mask to omit background sources when integrating the flux. Due to the resolution of MSSS (3′) not being significantly finer than the typical angular sizes of the galaxies in our sample (mean value of ), we were not able to deblend and remove background sources that are spatially coincident with a given galaxy. However, from our experience with observations of individual galaxies at various frequencies, the typical contribution from background sources is small, below 10% of the integrated flux (see e.g. Chyży et al. 2007).

Table A.1 presents the galaxies from our original sample that were detected in the MSSS survey and have reliable flux densities at 150 MHz. All galaxies have redshifts less than 0.04 and constitute our sample of 129 sources. Some of these objects have been identified as galaxy pairs, unresolved in the MSSS mosaics. In these cases, the measured flux density refers to the system as a whole. Table A.1 also provides the galaxy morphology (Hubble T-type), inclination angle of the galaxy disk (90° corresponds to an edge-on object), the distance, and supplemental information concerning the flux density measurements and the constructed global spectra characterized by four different kinds of flags:

• interaction flag (according to the NED database): 1 – single source; 2 – luminous infrared galaxy (LIRG); 3 – strongly interacting galaxy but not LIRG;

• flux density flag: 1 – flux density for a single source with an undisturbed disk; 2 – flux density for a single source but with unclear or amorphous disk plane, morphologically irregular dwarf, or starburst galaxy; 3 – flux density for a double or a triple source;

• spectral flag (see Sect. 3.2): 1 – MSSS flux density fits well (within σ₁₅₀) the spectrum constructed from the literature data; 2 – MSSS flux density is within 1 − 2σ₁₅₀ of the general spectral trend, or the spectral trend is estimated from only 2–3 flux density measurements; 3 – there is a lack of reliable/sufficient data to construct the global spectrum;

• flux density measurement method flag: this indicates whether the presented 150 MHz flux densities come from: 1 – the MSSS source catalogue; or 2 – our CASA measurements from the MSSS image mosaics.

3. Results

3.1. Low-frequency radio-FIR correlation

In order to test the sample quality and identify possible objects with a substantial AGN contribution, we studied the radio-FIR relation for all 129 sources using galaxy luminosities at 150 MHz (see Table A.1 for the flux densities and distances we used) and at 60 μm (Yun et al. 2001). Furthermore the sample was divided into three subsamples, depending on the NVSS flux density at 1.4 GHz (Yun et al. 2001): bright sources with S_NVSS >  143 mJy; medium bright sources with 87 mJy < S_NVSS <  143 mJy; and faint sources with S_NVSS <  87 mJy. Each subsample consisted of 43 objects.

The low-frequency radio-FIR relation that we obtain is similar to the high-frequency one (Fig. 1), which confirms the high quality of the MSSS data and corroborates the existence of a well-defined relation between the FIR emission of heated dust grains and the synchrotron emission of CR electrons, which dominates the radio emission of galaxies at low frequencies. There are no distinct outliers in this relation since we removed them from the sample as described in Sect. 2.1. The spread of data points in the diagram of MSSS and NVSS flux densities is even smaller (Fig. 1, right panel), which suggests that a significant part of the dispersion in the radio-FIR relation comes either from the uncertainties of the IRAS data or from the intrinsic spread in the relation. Either way, our radio data do not appear to be the limiting factor.

Fig. 1. Radio-FIR relation for the low-frequency MSSS (left) and high-frequency NVSS data (middle). The red, blue and green symbols represent faint, medium bright, and bright NVSS galaxies, respectively. The MSSS–NVSS luminosity diagram is shown on the right. The solid line is a bisector fit, while the dashed lines represent simple X vs. Y and Y vs. X linear regressions in the logarithmic space of the parameters.

Using the bisector method we estimated a slope in the radio-FIR power-law relation for the MSSS flux densities at 150 MHz as 1.12 ± 0.04. For the same set of galaxies we obtained the slope for the NVSS data at 1.4 GHz of 1.09 ± 0.03. Our results are consistent with those found (once AGN had been removed) by Hardcastle et al. (2016), Calistro Rivera et al. (2017), Gürkan et al. (2018), and Magnelli et al. (2015), who investigated the low-frequency radio-FIR relationships in individual deep fields. The wide sky coverage of MSSS enables us to explore this relationship in more detail in the local Universe (at z ≈ 0), which will be presented in a forthcoming paper.

One would expect that the effects of absorption by dense thermal gas could be particularly evident in galaxies with high SFRs, which are also more luminous in the infrared band (Condon et al. 1991). Powerful infrared galaxies would thus become weaker in the low radio frequency range. Our galaxies do not show such behaviour, as can be seen in Fig. 1, which may indicate that the effect of thermal absorption in our sample is relatively weak.

3.2. Deriving galaxy spectra

Our analysis of the spectral characteristics of radio emission of galaxies is based on their global spectra. We carried out an extensive search of the literature for galaxy flux densities at different frequencies, using mainly information from the NED and SIMBAD databases. We avoided old measurements, for example early interferometric arrays before 1980. For 31 galaxies we found data in the GaLactic and Extragalactic All-sky Murchison Widefield Array (GLEAM; Hurley-Walker et al. 2017) survey, which provides a catalogue of integrated fluxes from 20 narrow-band images as well as a fitted 200 MHz integrated flux density. For 19 galaxies from our sample, we did not find literature data of sufficient quality to derive spectra. Such objects, marked in Table A.1 by spectral flag 3, were excluded from further analysis. We also dropped double and multiple galaxy systems marked by the flux density flag 3 (Sect. 2.3). The resulting sample therefore consists of 106 galaxies.

Some examples of spectral trends, integrated flux density measurements from the MSSS survey, and published data are presented in Fig. 2. The MSSS and GLEAM measurements are consistent with each other and follow the general spectral trends found in the published data. Our main interpretation of these diagrams is that the spectra tend to slightly flatten at low frequencies (e.g. NGC 1569 and NGC 4102). However, there are also spectra that appear to be quite straight (e.g. NGC 3646 and NGC 5936).

Fig. 2. Examples of radio spectra for a subset of MSSS galaxies. The flux densities from eight individual MSSS spectral bands are in green, the interpolated MSSS flux densities at 150 MHz are in blue, the GLEAM flux densities from 20 individual spectral bands (if available) are in red, the GLEAM catalogue fitted flux densities at 200 MHz are in orange, and other published measurements are in black. The pink and black lines represent the low- and high-frequency power-law fits, respectively. Highlighted regions show the 68% confidence bands of the fits.

3.3. Spectral flattening

The spectral curvatures seen in Fig. 2 are small and usually not limited to just the lowest frequencies (e.g. NGC 972 and NGC 4102). Therefore, in order to evaluate spectral slopes, we constructed both low- and high-frequency spectral indices. The low-frequency index (α_low) was obtained by fitting a power-law to the interpolated MSSS flux density at 150 MHz and to literature data from 50 MHz up to 1.5 GHz. When the GLEAM data were available we used the fitted flux density at 200 MHz, and not the correlated data from individual sub-bands, which is similar to the approach with the MSSS data. The high-frequency index (α_high) between 1.3 GHz and 5 GHz was calculated in a similar way. We did not consider frequencies higher than this so as to avoid the effects of spectral flattening due to the increasing thermal component: for example for a sample of galaxies, the average thermal fraction given by Niklas et al. (1997) is 8 ± 1% at 1 GHz, and by Tabatabaei et al. (2017) is 10% at 1 GHz and 23% at 5 GHz. Both the indices derived for all our galaxies are listed in Table A.1 and some examples of the power-law spectral fits together with their 95% confidence bands are presented in Fig. 2.

As can be seen in Fig. 3, there is a systematic difference between α_low and α_high, which reflects the prevalence of galaxies that have flatter spectra at low frequencies than they have at high frequencies. The different spectral indices may result either from flattening of the low-frequency spectra or from steepening of the high-frequency ones; both have different underlying physical processes. The distributions of low- and high-frequency spectral indices are presented in Fig. 4. The distribution of the low-frequency index is wider and is shifted to smaller absolute values (indicating flatter spectra of galaxies). Accordingly, the high-frequency spectra can be observed as either steep or flat, while the low-frequency spectra are relatively flat. The two-sample Anderson–Darling, Kolmogoroff–Smirnov tests, and the Kruskal–Wallis rank test indicate that the hypothesis that the values of α_high and α_low are derived from the same population is highly unlikely since the p-values¹ in all tests are much smaller than 1%.

Fig. 3. Radio low-frequency αlow and high-frequency αhigh spectral indices for our sample of 106 galaxies. The straight line corresponds to identical spectral index values, i.e. to simple power-law spectra.

Fig. 4. Distributions of the radio low-frequency αlow and high-frequency αhigh spectral indices for the MSSS galaxies (in blue and pink respectively). Overlapping parts of the histograms are in violet.

The median value for α_low, with the uncertainty obtained using the bootstrap method as the standard deviation, is −0.57 ± 0.01. Analogously, for the high-frequency part of the spectrum, the spectral index α_high is −0.77 ± 0.03. This high-frequency estimate is very close to the spectral index of −0.74 ± 0.12, as estimated by Gioia et al. (1982) between 408 MHz and 10.7 GHz for a sample of 57 galaxies, which confirms the statistical consistency of both samples. It is also close to the 1.4–10 GHz spectral index of −0.79 ± 0.15 estimated for a sample of 36 galaxies by Tabatabaei et al. (2017). We can simply express the observed spectral curvature for each galaxy as the difference between the low- and high-frequency spectral indices Δα = α_low − α_high. The estimated median value of this parameter for the whole sample of galaxies turns out to be relatively small but positive: 0.18 ± 0.02, which underpins the low-frequency spectral flattening. This value is statistically significant because a non-parametric sign test allows us to reject a null hypothesis that the median of Δα is equal to zero (i.e. the p-value appears much smaller than 1%). In a similar way, we estimated a mean value of Δα as 0.20 ± 0.02, which according to a parametric Student’s t-test is also statistically significant (the p-value is again much lower than 1%).

We also analysed galaxy spectra by fitting two models to their flux densities S_ν between 50 MHz and 5 GHz: a power-law spectrum of the form

(3)

and a curved spectrum,

(4)

which are the simplest versions of polynomial models, having the lowest number of free parameters. In the model fitting we used the Marquardt–Levenberg algorithm to find the lowest reduced chi-square parameter.

In Table 1 we show the best-fit parameters A₁ and A₂ from the power-law and curved models, respectively, as well as the obtained goodness of fit for all galaxies from Fig. 2. We also compared the models by the F-test, using the ratio of the residual sum of squares from the power-law and curved models, scaled by corresponding degrees of freedom. If the p-value obtained was small (≤5%) we concluded that the curved model was statistically significantly better than the simple power-law model. Otherwise, we inferred that there is no convincing evidence to support the curved model and considered the power-law model as a satisfactory one (see the last column of Table 1). We found that only about 35% of galaxies had simple spectra described by the power-law model, while the other objects had curved spectra. For galaxies with curved spectra, the fitted values of A₂ were always negative. This corresponds to flatter spectra of galaxies at lower frequencies which is consistent with the statistical differences we found between distributions of α_low and α_high parameters (Fig. 3), and their median values. Moreover, the model with a curved spectrum was always better fitted than the power-law model for galaxies with higher values of Δα.

Therefore, our analysis shows that the spectral flattening towards low frequencies is small in nearby galaxies, yet is statistically significant. A similar conclusion was derived by Israel & Mahoney (1990), who noticed that flux densities at 57.5 MHz are systematically lower than expected from an extrapolation from measurements at higher frequencies (Sect. 1), and, for example, more recently by Calistro Rivera et al. (2017) using LOFAR observations of the Boötes deep field.

In Fig. 5 we present the two-point spectral index between 150 MHz and 1.4 GHz derived from the MSSS and NVSS data against the MSSS flux densities, and plot theoretical lines indicating different flux density limits at 1.4 GHz. The distribution of data points shows a distinct cut corresponding to 50 mJy at 1.4 GHz. This is understandable, since this was actually the limit applied as one of the selection criteria for our sample (Sect. 2.1). Although the scarcity of steep spectrum galaxies among the weakest sources (log(S_MSSS/mJy)< 2.5) results from a selection bias, it cannot account for the curvature observed in the global galaxy spectra.

Fig. 5. Two-point spectral index between 150 MHz and 1.4 GHz, as estimated from MSSS and NVSS data, for the sample of 106 galaxies. The dashed, solid, and dotted lines correspond to galaxy flux density limits of 25, 50, and 100 mJy at 1.4 GHz, respectively.

3.4. Spectra versus galaxy inclination

We noticed in Sect. 3.2 that spectra of our galaxies are slightly curved, which could indicate some kind of underlying physical process. One possibility could be internal thermal absorption of nonthermal emission (Sect. 1). We expect strong absorption effects to be at work for a well-mixed thermal and nonthermal plasma when galactic lines of sight are long, for example in edge-on galaxies (see Sect. 4.1). Also in the case when many localised H II regions are situated along the line of sight, the synchrotron emission from all regions along the line of sight have to pass through them, and their number, and hence their absorption, increases with higher inclination. To examine these predictions, we constructed a diagram to show the difference between low- and high-frequency spectral indices against the inclination angle i (Fig. 6). For strong absorption, we would find galaxies in the top-right corner of this diagram. However, there are no objects in this area and no distinct spectral dependence on galaxy tilt can be observed. The Kendall correlation coefficient between Δα and i is just −0.05, which confirms this finding. A similar conclusion could be drawn from the two-point spectral index versus the inclination angle (Fig. 7).

Fig. 6. Difference between low- and high-frequency spectral indices, Δα, vs. inclination angle i for the MSSS sample of 106 galaxies (i = 0° corresponds to a face-on galaxy). Galaxies with well-defined disks (93 objects) are indicated by red circles.

Fig. 7. Two-point spectral index between 150 MHz and 1.4 GHz, determined from MSSS and NVSS data, vs. the inclination angle i for the sample of 106 galaxies. Galaxies with well-defined disks (93 objects) are additionally marked by red circles. The horizontal line represent the median value of the spectral index of −0.56.

One can also consider a situation where the absorbing thermal gas is situated completely outside the synchrotron medium as a foreground layer while still within the galaxy. In that case, the synchrotron emission would be free-free absorbed, irrespective of the galaxy inclination, and this could potentially explain our results. However, such a configuration is rather unlikely as the thermal (H II) gas has a smaller spatial extent than the nonthermal emission, as in NGC 6946 (Tabatabaei et al. 2013), or in NGC 4254 (Chyży et al. 2007). Accordingly, we conclude that the flattening observed at low frequencies (Sect. 3.3) is not due to free-free absorption. These results contradict the claim of Israel & Mahoney (1990), who used a similar diagram to show that galaxy spectra are flatter for highly inclined objects. As a consequence, we do not see any need for a special low-temperature ionised gas postulated by those authors. We also considered that some galaxies in our sample may not have a well-determined value of viewing angle as their disks are not particularly regular. Such objects could have somehow affected probing the absorption origin of spectral flattening. Therefore, in another approach we excluded mergers (like UGC 8696), amorphous compact galaxy pairs (like NGC 5929/30), and dwarf-irregular galaxies (like NGC 4449) from our sample. These objects are indicated in Table A.1 by the flux density flag 2. The remaining 93 objects with well-defined disks are marked red in Fig. 6. It can be seen again that the distribution of Δα does not reveal any systematic dependence on the viewing angle, which is further confirmed by a low value (−0.04) of the Kendall correlation coefficient.

Apart from attempting to account for spectral flattening in terms of orientation, we also investigated its possible relationship with the morphological type of galaxies. The resulting diagram of the Δα as a function of the Hubble type T is shown in Fig. 8. No dependence of flattening on the galaxy’s morphology is to be seen; although there is a small number of flatter low-frequency spectra (smaller α_low) for late-type spirals, spectra such as these are also found for early-type objects (Fig. 9). Relatively flat spectra for dwarf galaxies have already been noticed by Klein et al. (2018) and interpreted as low CR confinement in low-mass galaxies. Considerably lower synchrotron components, and therefore systematically weaker magnetic fields, have recently also been found in late-type spiral galaxies, suggesting that a similar process could be at work in those objects as well (Chyży et al. 2017). Our results indicate that, apart from these trends that depend on galaxy properties generally associated with different morphological types, the spectral curvature (Δα) does not correlate with galaxy morphology.

Fig. 8. Difference between low- and high-frequency spectral indices, Δα, vs. morphological type T for the sample of 106 galaxies. Galaxies with well-defined disks (93 objects) are indicated by red circles.

Fig. 9. Two-point spectral index between 150 MHz and 1.4 GHz, computed from the MSSS and NVSS data, vs. morphological type T for the sample of 106 galaxies. Galaxies with well-defined disks (93 objects) are indicated by red circles.

4. Discussion

Our examination of the spectra of 106 galaxies (Sect. 3.4) has revealed gentle flattenings at low radio frequencies, which, contrary to the claims of Israel & Mahoney (1990), do not seem to be caused by thermal gas absorption. Such a conclusion is well supported by the results of Marvil et al. (2015), although they used different selection criteria, including, for instance, elliptical galaxies with radio emission possibly influenced by AGNs. This approach could lead to various effects in their sample, such as a positive spectral index. In order to avoid further complication, we deliberately left out such types of objects from our analysis (Sect. 2.1).

If the curved spectra observed in our sample of galaxies do not originate from absorption, we are facing a contradiction with earlier studies. Strong thermal absorption seems to be the only reasonable explanation for the distinct drop of low-frequency radio emission in the central part of M 82 (Wills et al. 1997; Varenius et al. 2015) and in other starburst galaxies (e.g. NGC 253 Kapińska et al. 2017). Also, the turnovers found in the spectra of some regions in the centre of the Milky Way (Roy & Pramesh Rao 2006) are well accounted for by absorption effects. So how can we possibly reconcile these results with our analysis of a large sample of galaxies?

In order to answer this question, we decided to construct a simple model of galaxy emission allowing us to analyse the effects of absorption with a realistic treatment of projection effects. The modelled radio emission originates from different regions of a galaxy that we assume to be composed of compact H II regions, supernovae, and diffuse plasma in both disk and halo. These regions emit and absorb synchrotron and free-free radiation of various amounts, which can be estimated from observational data. We did not take into account the Razin-Tsytovich effect or synchrotron self-absorption, which could affect the galaxy spectra at extremely low frequencies below 10 MHz (Fleishman & Tokarev 1995).

4.1. Numerical model of radio emission

In our numerical modelling we considered the 3D shapes of individual components of galaxies (such as the core and thin and thick disks), which were positioned on a 3D grid. We set the parameters describing thermal emission and absorption in each grid element, chose the galaxy inclination angle, and solved the radiative transfer equation along the line of sight at various frequencies. On the path through the source towards the observer, we distinguished two solutions for radiative transfer in the cell with index n, depending on the content of this and the preceding cells:

(5)

where

(6)

is the optical thickness of thermal gas in the cell with index n, T_e is the thermal electron temperature, is the emission measure, N_en is the thermal electron density in the cell, s is the cell size, B_n is the synchrotron intensity, and Ω is the solid angle of the cell. The first solution in Eq. (5) corresponds to a cell filled with well-mixed synchrotron- and thermal-emitting gas. The radio emission (S_n − 1) from the previous cell along the line of sight is absorbed by thermal gas in the cell with index n. The thermal and synchrotron gas components of this cell contribute with their emission minus the thermally absorbed part. The second solution applies to a cell with a solely synchrotron-emitting gas (e.g. in the galaxy halo). By solving the radiative transfer equation for all cells in the grid for a particular viewing angle and at different frequencies, we obtained synthetic maps of radio emission at various frequencies. We then integrated the flux density in the maps to construct the modelled global galaxy spectra for a variety of inclination angles.

First, we considered a simple model for a starburst galaxy represented by a single cylinder with a well-mixed thermal and nonthermal plasma. Our 3D model well reproduces the results of the analytical 1D modelling by Condon et al. (1991). Galaxies with higher thermal gas densities, but with the same thermal fraction (fixed at 10% at 1.4 GHz), have a higher turnover frequency in their integrated spectra (Fig. 10, left panel). In our modelling, we were also able to simulate what such a starburst would look like at different viewing angles. Less inclined galaxies are stronger radio emitters (less thermally absorbed), and have spectral turnovers at lower frequencies (Fig. 10, right panel).

Fig. 10. Left: effect of the thermal electron density (determined by the EM) on low-frequency spectra of starburst galaxies seen edge-on. The thermal and nonthermal emitting gases are fully mixed with a constant 10% thermal fraction at 1.4 GHz (to be compared with the analytical model of Condon et al. (1991). Right: similar model applied for a galaxy seen at a different viewing angle with Ne = 80 cm−3.

Because the starburst model of Condon et al. (1991) did not include a synchrotron halo, which is apparently present in such objects (e.g. Adebahr et al. 2013; Varenius et al. 2016), we constructed another model with synchrotron emission coming also from beyond the thermally emitting volume. In this model, the low-frequency spectra of galaxies strongly depended on the size of the synchrotron halo, measured by the ratio b of the volumes of thermal to non-thermal-emitting regions (Fig. 11, left panel). Value b = 1 corresponds to a region radiating both thermally and nonthermally, whereas b = 0.1 means that the central part of the galaxy, containing both thermal and synchrotron gas, accounts for only 10% of the entire synchrotron halo volume. In this model, the thermal fraction at 1.4 GHz was kept constant at 10% independent of b. We found that the resulting spectra depend strongly on b (Fig. 11, left panel). Furthermore, we found that even for a fixed halo size at the low value of b = 0.1, the spectra can still be easily modified by changing either the level of synchrotron intensity in the halo (Fig. 11, right panel) or the value of the nonthermal spectral index.

Fig. 11. Left: model of a starburst galaxy with a synchrotron halo of varying size, as measured by the ratio b (see text for details). The thermal fraction of the disk emission is assumed to be 10% at 1.4 GHz for all models. Right: similar model with the same synchrotron halo of b = 0.1 but different synchrotron intensity.

The above examples show that the shape of integrated spectra at low frequencies depends strongly on both the specifics of geometry of the radio emitting regions and the parameters of the thermal and nonthermal emissions. Adding a halo component to a simple starburst region introduces a further potential ambiguity in modelling if only integrated spectra are considered. Therefore, it is not possible to fully interpret integrated spectra without detailed information on the distribution of thermal and nonthermal radio emissions throughout the galaxy.

Therefore, in the following sections, we carefully investigate those galaxies from our sample (Sect. 2.1) for which such details are available. Firstly we model M 51, a face-on spiral galaxy with relatively low star formation activity, and secondly we model M 82, a starburst galaxy seen nearly edge-on. With our models, we are able to modify inclination angles and compare the constructed spectra to observed ones. This analysis enables us to draw general conclusions on the role of thermal absorption as well as on the origin of the curved spectra that are observed in nearby galaxies.

4.2. Modelling M 51-like galaxies

We modelled M 51 as an example of a non-starbursting galaxy with a low inclination angle (i ≈ 20°). The thermal emission coming from the ionised gas was represented by two vertical components: the thin and thick disk. Just as in the case of the Milky Way, we described the distribution of the thermal electron density (N_e) in these components with exponential functions with vertical scale heights of about 100 pc and 1 kpc, respectively (cf. Cordes 2004). The radial dependence of N_e was approximated by inner and outer exponential disks. In the inner disk, we determined the radial profile of N_e from Hα-derived emission measure maps by Greenawalt et al. (1998), while for the outer disk we assumed a radial scale length of 10 kpc (Gutiérrez & Beckman 2010). The density of free electrons in the plane of the outer disk is approximated by 0.02 cm⁻³, like in the models of the Milky Way (Ferrière 2001). Therefore, the 3D distribution of thermal electrons was modelled as:

(7)

where A^a is a constant and the superscript a used in this and following equations signifies the modelling of M 51-like galaxies. The value of A^a was determined so that the integrated radio thermal emission corresponds to a thermal fraction of 28% at 14.7 GHz, estimated for M 51 by Klein et al. (1984). This equation was used to find the EM_n in each model cell and then τ_n according to Eq. (6), where we assumed T_e ≈ 10⁴ K everywhere.

The approximate properties of the nonthermal emission throughout the galaxy were derived by the iteration method outlined below. The multifrequency observations of M 51 revealed that the radial profiles of the radio intensity and the spectral index between 151 MHz and 1.4 GHz vary significantly within the galaxy (Mulcahy et al. 2014). Furthermore, the local spectra in the central part of M 51 are relatively flat, but become steeper closer to the galaxy edges, as expected for spectral ageing of CR electrons by synchrotron and inverse Compton radiation. There is also evidence for diffusion of CRs from the star-forming regions into the inter-arm regions and outer parts of the galaxy (Mulcahy et al. 2016). Accordingly, we modelled the unabsorbed nonthermal intensity B_n as two exponential vertical disks with a radial variation that we described using four different continuous functions:

(8)

where the consecutive equations are for different sections of radial distance from the galactic centre: r <  0.5 kpc, 0.5 ≤ r <  1.7 kpc, 1.7 ≤ r <  10 kpc, and 10 ≤ r ≤ 16 kpc respectively. The C^a parameter is a scale factor, , , are scale-lengths of separate exponential profiles describing radial dependence of nonthermal emission in corresponding sections. In the first section (r <  0.5 kpc), we kept the emission constant along the radius, because the central part of M 51 is not resolved in the observational profiles (Mulcahy et al. 2016); furthermore, an exponential profile would lead to a strong centrally concentrated source, which is not expected, due to CRs diffusion. The parameters and denote exponential scale heights of the vertical thin and thick nonthermal disks, respectively.

In order to reproduce the radial profile of the spectral index of M 51 (Mulcahy et al. 2014), the scale heights and lengths of the modelled radio emission profiles have to be frequency-dependent. Such a frequency dependence can be underlined by the various CR propagation processes (Mulcahy et al. 2014; Krause et al. 2018). We assumed the power-law functions: (ν),(ν) ∝ (ν/ν₁)^Da and ,  ∝ (ν/ν₂)^Ea. The values of , , and at ν₁ = 0.15 GHz were estimated as: 1.1 kpc, 5.3 kpc, and 2.1 kpc, respectively (Mulcahy et al. 2014). The values of and at ν₂ = 4.85 GHz we approximated as 0.3 kpc and 1.8 kpc, respectively (Krause et al. 2018). We notice that the scale factor C^a should also depend on frequency, C^a = f(ν). Initially, we assumed f in the form of a power-law function. However, we found that in order to obtain a better fit of the global spectrum it was necessary to introduce a third-degree polynomial: .

After assuming some initial values for Da, Ea, and C^a in Eq. (8), the above values of , , ,  , and values of parameters describing the thermal emission (Eq. (7)), we solved the transfer equation (Eq. (5)) for various frequencies, and constructed model maps of radio emission. Next, we integrated flux densities using these maps, derived the modelled global spectrum, and compared it with the observed spectrum of M 51. We tried out a number of models with different values of model parameters. The reduced chi-square value was used as a measure of goodness of the fit. We also qualitatively compared the model results with the radial profile of the radio emission of M 51 and the spectral index profile from Mulcahy et al. (2014). As the final stage, we changed the model parameters by a small value (±10%), identifying the model that gave the best fit. The parameters of the best-fitting model resulting from this procedure are given in Table 2.

The obtained model (Fig. 12) reproduced well the observed radial intensity and spectral index profiles from Mulcahy et al. (2014). The break in the radial intensity profile at r ≈ 2 kpc defines the transition from the inner to the outer star-forming disk. Another break at r ≈ 10 kpc corresponds to the virtual disappearance of CR sources at larger galactocentric radii where, as suggested by Mulcahy et al. (2014), the outer edge of the star-forming disk is located. The resulting global spectrum for the modelled galaxy also matches well the literature data (Fig. 13, left panel).

Fig. 12. Left: modelled radial profiles of total radio emission of an M 51-like galaxy seen face-on at 30, 150, and 1400 MHz. Middle: model vertical profiles of total radio emission of an M 51-like galaxy seen edge-on at 30, 150, and 1400 MHz. Right: local synchrotron spectra of the galaxy at different distances along the galaxy major axis from the centre.

Fig. 13. Left: resulting global spectrum of an M 51-like galaxy from our 3D model as seen face-on (red solid line) with the synchrotron component without absorption (orange dotted line) and thermal free-free emission (black dashed line). Black solid circles denote literature data used during the modelling and the interpolated flux density at 150 MHz from MSSS survey. For the sake of completeness, data from the individual MSSS sub-bands (blue circles) are also shown. Middle: global spectra for different inclination angles of 0, 30, 60, and 90 degrees, respectively. Right: local spectrum of an area surrounding the galaxy centre for different inclination angles.

We calculated the extent to which the model parameters affect the fitted spectrum and the goodness of the model fit by changing the obtained best-fit values by +20% and −20%. The resulting chi-square values are presented in Table 2. The influence is also illustrated in Fig. 14 for several selected parameters. The ranges of the fitted model spectra are narrow, comparable to the data spread. The results were similar for all the other parameters.

Fig. 14. Ranges of the modelled global spectrum of M 51-like galaxy as affected by increasing and decreasing best-fit model parameters , , and by ±20%. Black solid circles denote literature data used during the modelling and the interpolated flux density at 150 MHz from the MSSS survey. For the sake of completeness, the data from individual MSSS sub-bands are also shown (brown rectangles).

Using the obtained model parameters to best fit the M 51 observations, we then solved the transfer equation for a number of viewing angles and frequencies. The results are shown in the form of synthetic maps of the radio continuum emission in Fig. 15. Our modelling revealed no significant absorption effects at 150 MHz, that is, the frequency of the MSSS survey, independently of viewing angle. This result corresponds well to the observed global spectrum of the edge-on galaxy NGC 891 recently obtained by Mulcahy et al. (2018), in which the effect of thermal absorption is also barely visible. Strongly reduced emission due to absorption was only found in our maps at 30 MHz for highly inclined galaxies (i ≥ 60°).

Fig. 15. Synthetic radio emission maps of a nonstarburst M 51-like galaxy; from top to bottom, rows show results at 30 MHz, 150 MHz, and 1400 MHz; from left to right, columns show results for inclination angles of 0, 30, 60, and 90 degrees. The colour scale is in mJy beam−1.

The effects of thermal absorption on the global spectrum are very limited; they can only be seen as a weak flattening below 20 MHz (Fig. 13, middle panel). This is due to the relatively small scale-height of thermal electron distribution in M 51 (100 pc of the thin disk and 1 kpc in the thick disk, see Eq. (7)) as compared to the synchrotron scale-height (0.3 kpc of the thin disk and 1.8 kpc of the thick disk, see Table 2). Moreover, thermal emission has scale-heights that are half the size of the mentioned scales of thermal electrons due to the direct proportionality of thermal emission to . A more quantitative analysis of absorption is presented in Table 3, which contains the calculated optical depths at different radial distances from the centre of the modelled face-on galaxy. We note that only in the centre of the galaxy is the optical depth significant (τ = 0.032) giving a reduction of synchrotron emission at 150 MHz by only (1 − exp(−0.3)) × 100 = 3%. Our model predicts that at a much lower frequency of 30 MHz much stronger absorption, by up to 30% (for τ = 0.36), is to be observed, but only in the central region up to about 1 kpc radius (Table 2).

Therefore, we expect no spectral turnovers due to absorption in the global spectra of galaxies of this type at any observation frequency available from the Earth. Our modelling properly accounts for the results of our statistical analysis of the observed spectra of nearby galaxies (Sect. 3.4): the statistically weak effect of thermal absorption on the global galaxy spectra was suggested by the radio-FIR relation and confirmed by the spectral index versus inclination angle diagram (Fig. 6). Our analysis therefore indicates that for typical late-type spiral galaxies with only moderate stellar activity, such as M 51, any observed curvature visible in the global spectrum at frequencies above 20 MHz is the result of curved local synchrotron spectra rather than of internal thermal absorption. Such curved nonthermal spectra can be explained by the combined effects of CR electron energy losses and CR transport.

As a case in point, observations of the Milky Way spectrum at very low frequencies (e.g. Brown 1973; Ellis 1982) show a turnover at about 3 MHz, interpreted as the effect of thermal absorption. This is different to our modelling of M 51-like galaxies, where no turnover is found down to 1 MHz. However, it should be noted that the measurements of the Milky Way spectrum are by necessity made from the inside, while external galaxies are seen from the outside. It is likely that due to our location within the plane of the Galaxy’s disk, we can observe efficient absorption of radio emission by the thermally ionised gas in the thin disk. If the Milky Way does not differ too much from M 51, our modelling would therefore indicate that there would be no turnover in the Galaxy’s spectrum as seen from outside.

We also noticed that there could be another effect associated with the Milky Way. It can be expected that its thermal gas absorbs the incoming emission of external galaxies, which could be seen as spectral turnovers at about 3 MHz for objects outside the Galactic plane. Very-low-frequency observations of galaxies distributed over the celestial sphere and studies of their low-frequency spectra would also contribute substantially to our knowledge of the Milky Way’s ISM and facilitate modelling of its detailed structure.

From our modelling, we also obtained local spectra of the central parts of M 51-like galaxies viewed at a different inclination angle (Fig. 13, right panel). The spectra of these compact regions are different from the modelled integrated galaxy spectra (Fig. 13, middle panel), because the shapes of local spectra highly depend on the galaxy tilt. In the case of edge-on galaxies, they show distinct turnovers at frequencies of about 100 − 200 MHz due to strong absorption in the thin disk. The observed local spectra in the star-forming regions in the disk of NGC 891 seem to be indeed clearly affected by thermal absorption, showing low nonthermal spectral index, below the injection value (Mulcahy et al. 2018).

For the less inclined galaxies (i <  60°), turnovers become less pronounced, and at very low frequencies (especially below 10 − 20 MHz) a rise of spectra is observed due to unabsorbed synchrotron emission in the thick galactic disk. Our modelling therefore indicates that observations of compact regions in galaxies at high resolution and very low frequencies would enable us to constrain models of their constituents as well as properties of the ISM. Unfortunately, such observations would require radio instruments located in space or on the Moon.

4.3. Modelling of M 82-like galaxies and relation to sample galaxies

We also modelled M 82, a starburst galaxy visible at a high inclination angle (i ≈ 80°). It shows an extremely complex ISM, due to strong gravitational interaction with M 81. Numerous young massive star clusters and supernovae in the central starburst region drive an Hα- and X-ray-bright superwind (Westmoquette et al. 2009). In the circumnuclear starburst region (with a radius of about 250 pc) the emission measures in the compact sources are typically about 10⁵ cm⁻⁶ pc (Varenius et al. 2015; Wills et al. 1997). Pervasive diffuse Hα emission at low levels is found throughout the halo of the galaxy up to about 1 kpc height (Shopbell & Bland-Hawthorn 1998). Therefore, in our model of M 82, we approximated the 3D distribution of thermal electrons in the form of thin and thick exponential disks and estimated their characteristic scales by fitting a double-exponential function to vertical profiles of Hα emission obtained from the Hα map of M 82 from the SINGS survey (Kennicutt et al. 2003), which gave the following formula:

(9)

where A^b and B^b are scale factors and the superscript b denotes, in this and the following equations, the modelling of M 82-like galaxies. The values of these factors were determined so that the integrated radio thermal emission should correspond to a thermal fraction of 3% at 1.4 GHz, as estimated for M 82 by Basu et al. (2012). The peak of the emission measure in the central part of the disk, resulting from the modelled distribution of N_e given by Eq. (9), is 2.3 × 10⁶ cm⁻⁶ pc.

In the nonthermal part of the model, we introduced a thick disk with the scale height similar to the scale length. Such a synchrotron halo extending up to a height of about (≈3.56 kpc) is seen in the observations of M 82 at 330 MHz and 1.4 GHz (Adebahr et al. 2013). In our model, we did not include details such as the asymmetries between the northern and southern parts of the galaxy, which could be the result of faster CR electron energy losses or slower wind speeds in the southern outflow. Neither did we take into account the emission from individual supernovae and star clusters. Therefore, the nonthermal radio emission is represented in the model in the following way:

(10)

where C^b and F^b are scale factors, and denote exponential scale-lengths of inner and outer synchrotron disks, and and are scale-heights of exponential thin and thick disks, respectively. Similar constituents of M 82 were involved in the modelling of Yoast-Hull et al. (2013), but in a simplified cylindrical geometry.

We employed a similar fitting procedure as for M 51, comparing modelled and observed radio emission profiles along the major and minor axes, using the data of Adebahr et al. (2013) as reference. We assumed power-law functions: (ν),(ν) ∝ (ν/ν₁)^Db and ,  ∝ (ν/ν₁)^Eb, where ν₁ refers to 0.33 GHz. The values obtained for the scale lengths R₁ and R₂ were 0.4 and 0.8 kpc, respectively, and for the scale heights Z₁ and Z₂ we found values of 0.13 and 0.6 kpc, respectively (see Table 4).

As for M 51 (Sect. 4.2), the integrated spectra obtained from the range of M 82 models were compared with the observed spectrum by measuring the goodness of fit with the chi-square value. The best-fit spectrum is shown in Fig. 16 (left panel) and the corresponding values of model parameters are given in Table 4. We also present how the goodness of fit was affected when model parameters were changed by ±20%. In Fig. 17, we illustrate this influence for three selected parameters. The range of the obtained modelled spectra is narrow (close to the data spread), and its shape confirms spectral flattening at low frequencies.

Fig. 16. Global spectra from our 3D model of an M 82-like galaxy. Left: global spectrum of a highly inclined galaxy (red solid line) with the synchrotron component without absorption (orange dotted line) and thermal free-free emission with absorption (black dashed line). Black solid circles denote literature data used during modelling, including the interpolated flux density at 150 MHz from the MSSS survey. For the sake of completeness, data from individual MSSS sub-bands (blue circles) and flux densities determined during the commissioning phase of the 42-element Allen Telescope Array (Williams & Bower 2010) (orange triangles) are also shown. Middle: comparison of global spectra for edge-on (orange dotted line) and face-on galaxies (red solid line). Right: local spectra from the galaxy centre for different inclination angles.

Fig. 17. Ranges of the modelled global spectrum of M 82-like galaxy as affected by increasing and decreasing best-fit model parameters , , and by ±20%. Black solid circles denote literature data used during modelling including the interpolated flux density at 150 MHz from the MSSS survey. For the sake of completeness data from the individual MSSS sub-bands (brown rectangles) and flux densities determined during commissioning phase of the 42-element Allen Telescope Array (Williams & Bower 2010) (orange triangles) are also shown.

Some modelled maps that we obtained by solving the radiative transfer Eq. (5) for various frequencies and viewing angles are presented in Fig. 18. The effect of free-free absorption can already be seen in the maps at 150 MHz as a significant drop of the emission in the thin star-forming disk. This region corresponds to a modelled emission measure of several 10⁵ cm⁻⁶ pc, which well reproduces the results of Wills et al. (1997). The LOFAR high-resolution observations of M 82 by Varenius et al. (2015) reveal a similar decrease of emission in the galactic disk. While in our model maps the drop at 150 MHz is relatively localised, the maps at 30 MHz predict a much larger depression all around the galactic centre (Fig. 18, top-right panel). The existence of such a depression could be verified by future observations of starburst galaxies with the LOFAR low-band antenna (LBA).

Fig. 18. Synthetic radio emission maps of a starburst M 82-like galaxy; from top to bottom, rows show results at 30 MHz, 150 MHz, and 1400 MHz; from left to right, columns show results for inclination angles of 0, 30, 60, and 90 degrees. The colour scale is in mJy beam−1.

Unlike typical nearby galaxies like M 51, starbursting objects like M 82 show strong absorption effects that appear in the integrated spectra for all viewing angles (Fig. 16, left and middle panels). Towards lower frequencies, the modelled integrated spectrum of M 82 first flattens but then slightly rises again below frequencies of approximately 10 MHz due to its synchrotron halo. Depending on the viewing angle of the modelled sources, a spectral flattening due to absorption is observed between 200 and 400 MHz. As can be inferred from Fig. 19, the global spectrum is mainly shaped by the emission from the inner part of the galaxy, at least to a height of about 500 pc. At heights in excess of about 700 pc, the profiles show no indication of spectral flattening.

Fig. 19. Modelled vertical profiles of the total radio emission of an M 82-like galaxy seen edge-on at 30, 150, and 1400 MHz.

Modelling M 82 as a single starburst region involving mixed thermal and nonthermal gas components, as in the geometrically simplified starburst model presented in Fig. 10, cannot account for the observational data at frequencies below 200 MHz (Fig. 16, left panel and Fig. 17). One needs to add a synchrotron envelope with layers in front of the disk that are not thermally absorbed to increase the integrated flux densities at the lowest frequencies (see Yoast-Hull et al. 2013). Such an effect was achieved even in our geometrically simplified model after introducing a halo, as presented in Fig. 11.

The local spectra of the modelled compact starburst core show strong turnovers that move from ≈1 GHz for edge-on galaxies to ≈500 MHz for face-on objects. For M 51-like galaxies (Sect. 4.2), due to the lower emission measure in the galactic centre, the turnover appears at lower frequencies (of about 200 MHz) and for the edge-on objects only (Fig. 13, right panel).

The SFR of M 82 is relatively high as compared to most of the galaxies in our sample, while not being particularly extreme among starburst galaxies. Accordingly, we can predict even stronger absorption effects in the integrated spectra of violently star-forming galaxies as well as of more distant galaxies. The local spectra of compact star-forming regions within these objects should reveal turnovers at even higher frequencies (around or even exceeding 1 GHz). Starburst galaxies are in fact rare in the local Universe, which can explain why we found just very weak spectral curvatures in our sample of galaxies (Sect. 3.2). However, we note that Calistro Rivera et al. (2017) searched for evolution of spectral curvature in star-forming galaxies out to high redshift (which in their flux-limited sample corresponds to luminosity) and did not find any evidence for it, albeit with much sparser sampling of the radio spectra than we present here. Clear evidence of thermal absorption in spectra of starburst galaxies was recently reported by Galvin et al. (2018).

Our modelling of both M 51 and M 82 required the introduction of locally curved nonthermal spectra, which could be the effect of spectral ageing of the CR electrons due to for example synchrotron losses growing stronger at high CR energies. Such losses become most pronounced at high radio frequencies and are likely to be best observed in galactic outskirts due to CR transport (Mulcahy et al. 2014). A difference between the spectral indices at low and high frequencies could also be influenced by other CR electron energy losses: nonthermal bremsstrahlung, and ionisation, which is most effective at low energies (Murphy 2009). For typical field strengths of 10 μG and atomic hydrogen densities of 1 cm⁻³, the transition energy above which synchrotron cooling begins to dominate over bremsstrahlung and ionisation cooling is at about 1 GeV, which corresponds to a spectral break at about 200 MHz. For stronger fields (e.g. 15 μG) or less dense ISM, synchrotron cooling can dominate over CR electron losses down to frequencies of about 100 MHz.

On the other hand, Klein et al. (2018) claim that the integrated spectra of most spiral and dwarf galaxies can be fitted with a break or an exponential decline in their high-frequency part (> 1 GHz). Such features could be due either to synchrotron and inverse-Compton losses accompanied by advective transport in a galactic wind or due to a cutoff in recent times of the CR production rate. These results are different from those of Niklas & Beck (1997) and Tabatabaei et al. (2017), who, in spite of not considering breaks in the synchrotron spectrum in their models, were still able to fit most galaxy spectra quite well. Our modelling shows that in general the validity of the approach used in all the aforementioned studies, which are based entirely on the integrated spectra of galaxies, is very limited. As local synchrotron spectra are different in different galaxy regions (as can be seen e.g. in M 51; Fig. 12, right panel), the assumption of a constant break frequency or thermal fraction in the entire confinement volume is not correct. Different break/flattening frequencies due to locally different CR losses or thermal absorption are smeared out after adding the emission in constructing the global spectrum (see Lisenfeld & Völk 2000; Basu et al. 2015). Spatially resolved radio observations of galaxies and modelling of spatial composition of the ISM components are advisable in further studies of radiation processes and CR transport as well as for a proper interpretation of the global spectra of galaxies. For most of our sample galaxies, such high-resolution observations at low frequencies are not yet available; however, they will soon become available from the ongoing LOFAR Two-metre Sky Survey (LoTSS; Shimwell et al. 2017).

5. Summary and conclusions

We compiled a sample of 129 nearby star-forming galaxies from the LOFAR MSSS survey and measured their integrated flux densities at 150 MHz (Table A.1). Combining them with flux densities at various frequencies from the literature, we constructed global radio spectra for 106 objects. A statistical analysis of the data led us to the following conclusions.

• The low-frequency radio-FIR relation resembles the high-frequency one, which confirms (i) the good quality of the MSSS data, (ii) the proper removal of AGN-dominated galaxies from the sample, and (iii) the tight relation between dust emission and the almost pure synchrotron radiation at low radio frequencies. It may also indicate that the effect of thermal absorption in the sample is relatively weak at 150 MHz.

• In general, the global spectra of galaxies in our sample are statistically flatter at lower frequencies, or equivalently, steeper at higher frequencies (Fig. 4). The sample median value of the low-frequency spectral index α_low, measured between ≈50 MHz and 1.5 GHz, is −0.57 ± 0.01, which is significantly larger than the corresponding value of the high-frequency spectral index α_high = −0.77 ± 0.03, calculated between 1.3 and 5 GHz. For each galaxy we estimated the spectral curvature as the difference Δα = α_low − α_high and found that its median value is 0.18 ± 0.02, indicating a small but statistically significant spectral flattening at low frequencies.

• There is no tendency for highly inclined galaxies to have more flattened low-frequency spectra (Fig. 6). Hence the observed flattening is not due to thermal absorption, which contradicts the suggestion by Israel & Mahoney (1990). We show that the flattening does not depend on galaxy morphology either.

Interpretation of these results and discussion of processes to shape the low-frequency spectra were performed with a numerical model of radio emission involving absorption and projection effects. Our modelling of M 51-like galaxies, which represent low-SFR objects in the local Universe, suggests that local depressions of emission due to thermal absorption can only be seen in spatially resolved radio maps at frequencies below ≈30 MHz. We show that thermal absorption can influence the global spectrum at frequencies only below 20 MHz. A weak flattening observed in some galaxy spectra above this frequency (cf. Fig. 2) can be accounted for by the curvature in the synchrotron spectrum due to for example CR electron energy losses and diffusion of CRs from their regions of acceleration. We also predict a possible influence of the Milky Way foreground gas on external galaxy spectra at frequencies below 10 MHz.

Our modelling of M 51-like galaxies indicates that local radio spectra of central parts of galaxies may differ much from their global spectra, and may even show turnovers due to thermal absorption in edge-on objects. Such turnovers could appear at frequencies of about 100–200 MHz (Fig. 13) and high-resolution observations are necessary to reveal them. Our models also show that even galaxies with simple, power-law-like global spectra can locally show significantly curved spectra, especially far away from the galactic centres. Therefore, integrated spectra alone cannot be properly interpreted without supplementary data on the properties of the local ISM within the galaxies.

Moreover, we modelled galaxies with higher SFRs based on observations of M 82. Compared to M 51, the modelled M 82-like galaxies show much stronger absorption effects irrespective of the galaxy tilt. At 150 MHz, strong absorption was found mainly along the thin disk, while at 30 MHz our model predicts a more extensive depression of emission around the galaxy centre (Fig. 18). Towards low frequencies, the integrated spectra of modelled starbursts first flatten and then slightly increase again at about 10 MHz due to presence of the synchrotron halo. This could explain why simple modelling based on just one mixed region of thermal and nonthermal gas is not able to reproduce the observed spectra of M 82.

In the modelled local spectra of the central starbursts strong turnovers were found below about 1 GHz for galaxies viewed edge-on and below about 500 MHz for galaxies viewed face-on (Fig. 16). Therefore one can expect much stronger free-free absorption effects and significant changes in the radio spectra of LIRGs and more distant galaxies with high SFRs. Since starburst galaxies are not frequent in the local Universe we observed mostly weak spectral curvatures in our sample.

In this work we primarily discussed the impact of free-free absorption on the low-frequency spectra of galaxies. The forthcoming, more sensitive LoTSS HBA survey (Shimwell et al. 2017), and particularly the planned LBA survey at ≈50 MHz, as well as Murchison Widefield Array data covering a wide frequency range (e.g. Hurley-Walker et al. 2017) are highly desirable to distinguish between a number of processes shaping the spectra in galaxies at frequencies below 100 MHz (Sect. 1) as well as to verify our conjectures. Broad-band polarization observations spanning low and high frequencies would also help to constrain the thermal content of galaxies.

Acknowledgments

WJ acknowledges support by the Polish National Science Centre grant No. 2013/09/N/ST9/02511. MJH acknowledges support from the UK Science and Technology Facilities Council [ST/M001008/1]. GG acknowledges the CSIRO OCE Postdoctoral Fellowship. AOC and AMS gratefully acknowledge support from the European Research Council (ERC) under grant ERC-2012-StG-307215 LODESTONE. RJvW acknowledges support from the ERC Advanced Investigator programme NewClusters 321271. RPB acknowledges funding from the ERC under the European Union’s Horizon 2020 research and innovation programme (ERC Starting Grant agreement 715051 “Spiders”). We acknowledge the help of the Academic Computer Centre CYFRONET AGH in Krakow, Poland, where our numerical analyses were partly performed with the use of the computing cluster Prometheus. We acknowledge the use of the HyperLeda (http://leda.univ-lyon1.fr), NED (http://nedwww.ipac.caltech.edu), and SINBAD (http://simbad.u-strasbg.fr) databases. We thank Fatemeh Tabatabaei and an anonymous referee for helpful comments on the manuscript. This paper is based (in part) on data obtained with the International LOFAR Telescope (ILT). LOFAR (van Haarlem et al. 2013) is the Low Frequency Array designed and constructed by ASTRON. It has observing, data processing, and data storage facilities in several countries, that are owned by various parties (each with their own funding sources), and that are collectively operated by the ILT foundation under a joint scientific policy. The ILT resources have benefited from the following recent major funding sources: CNRS-INSU, Observatoire de Paris and Université d’Orléans, France; BMBF, MIWF-NRW, MPG, Germany; Science Foundation Ireland (SFI), Department of Business, Enterprise and Innovation (DBEI), Ireland; NWO, The Netherlands; The Science and Technology Facilities Council, UK, Ministry of Science and Higher Education, Poland.

References

Appendix A Additional table

All Tables

All Figures

Fig. 1. Radio-FIR relation for the low-frequency MSSS (left) and high-frequency NVSS data (middle). The red, blue and green symbols represent faint, medium bright, and bright NVSS galaxies, respectively. The MSSS–NVSS luminosity diagram is shown on the right. The solid line is a bisector fit, while the dashed lines represent simple X vs. Y and Y vs. X linear regressions in the logarithmic space of the parameters.

In the text

Fig. 2. Examples of radio spectra for a subset of MSSS galaxies. The flux densities from eight individual MSSS spectral bands are in green, the interpolated MSSS flux densities at 150 MHz are in blue, the GLEAM flux densities from 20 individual spectral bands (if available) are in red, the GLEAM catalogue fitted flux densities at 200 MHz are in orange, and other published measurements are in black. The pink and black lines represent the low- and high-frequency power-law fits, respectively. Highlighted regions show the 68% confidence bands of the fits.

In the text

	Fig. 3. Radio low-frequency αlow and high-frequency αhigh spectral indices for our sample of 106 galaxies. The straight line corresponds to identical spectral index values, i.e. to simple power-law spectra.
In the text

	Fig. 4. Distributions of the radio low-frequency αlow and high-frequency αhigh spectral indices for the MSSS galaxies (in blue and pink respectively). Overlapping parts of the histograms are in violet.
In the text

	Fig. 5. Two-point spectral index between 150 MHz and 1.4 GHz, as estimated from MSSS and NVSS data, for the sample of 106 galaxies. The dashed, solid, and dotted lines correspond to galaxy flux density limits of 25, 50, and 100 mJy at 1.4 GHz, respectively.
In the text

	Fig. 6. Difference between low- and high-frequency spectral indices, Δα, vs. inclination angle i for the MSSS sample of 106 galaxies (i = 0° corresponds to a face-on galaxy). Galaxies with well-defined disks (93 objects) are indicated by red circles.
In the text

	Fig. 7. Two-point spectral index between 150 MHz and 1.4 GHz, determined from MSSS and NVSS data, vs. the inclination angle i for the sample of 106 galaxies. Galaxies with well-defined disks (93 objects) are additionally marked by red circles. The horizontal line represent the median value of the spectral index of −0.56.
In the text

	Fig. 8. Difference between low- and high-frequency spectral indices, Δα, vs. morphological type T for the sample of 106 galaxies. Galaxies with well-defined disks (93 objects) are indicated by red circles.
In the text

	Fig. 9. Two-point spectral index between 150 MHz and 1.4 GHz, computed from the MSSS and NVSS data, vs. morphological type T for the sample of 106 galaxies. Galaxies with well-defined disks (93 objects) are indicated by red circles.
In the text

Fig. 10. Left: effect of the thermal electron density (determined by the EM) on low-frequency spectra of starburst galaxies seen edge-on. The thermal and nonthermal emitting gases are fully mixed with a constant 10% thermal fraction at 1.4 GHz (to be compared with the analytical model of Condon et al. (1991). Right: similar model applied for a galaxy seen at a different viewing angle with Ne = 80 cm−3.

In the text

	Fig. 11. Left: model of a starburst galaxy with a synchrotron halo of varying size, as measured by the ratio b (see text for details). The thermal fraction of the disk emission is assumed to be 10% at 1.4 GHz for all models. Right: similar model with the same synchrotron halo of b = 0.1 but different synchrotron intensity.
In the text

	Fig. 12. Left: modelled radial profiles of total radio emission of an M 51-like galaxy seen face-on at 30, 150, and 1400 MHz. Middle: model vertical profiles of total radio emission of an M 51-like galaxy seen edge-on at 30, 150, and 1400 MHz. Right: local synchrotron spectra of the galaxy at different distances along the galaxy major axis from the centre.
In the text

Fig. 13. Left: resulting global spectrum of an M 51-like galaxy from our 3D model as seen face-on (red solid line) with the synchrotron component without absorption (orange dotted line) and thermal free-free emission (black dashed line). Black solid circles denote literature data used during the modelling and the interpolated flux density at 150 MHz from MSSS survey. For the sake of completeness, data from the individual MSSS sub-bands (blue circles) are also shown. Middle: global spectra for different inclination angles of 0, 30, 60, and 90 degrees, respectively. Right: local spectrum of an area surrounding the galaxy centre for different inclination angles.

In the text

	Fig. 14. Ranges of the modelled global spectrum of M 51-like galaxy as affected by increasing and decreasing best-fit model parameters , , and by ±20%. Black solid circles denote literature data used during the modelling and the interpolated flux density at 150 MHz from the MSSS survey. For the sake of completeness, the data from individual MSSS sub-bands are also shown (brown rectangles).
In the text

	Fig. 15. Synthetic radio emission maps of a nonstarburst M 51-like galaxy; from top to bottom, rows show results at 30 MHz, 150 MHz, and 1400 MHz; from left to right, columns show results for inclination angles of 0, 30, 60, and 90 degrees. The colour scale is in mJy beam−1.
In the text

Fig. 16. Global spectra from our 3D model of an M 82-like galaxy. Left: global spectrum of a highly inclined galaxy (red solid line) with the synchrotron component without absorption (orange dotted line) and thermal free-free emission with absorption (black dashed line). Black solid circles denote literature data used during modelling, including the interpolated flux density at 150 MHz from the MSSS survey. For the sake of completeness, data from individual MSSS sub-bands (blue circles) and flux densities determined during the commissioning phase of the 42-element Allen Telescope Array (Williams & Bower 2010) (orange triangles) are also shown. Middle: comparison of global spectra for edge-on (orange dotted line) and face-on galaxies (red solid line). Right: local spectra from the galaxy centre for different inclination angles.

In the text

Fig. 17. Ranges of the modelled global spectrum of M 82-like galaxy as affected by increasing and decreasing best-fit model parameters , , and by ±20%. Black solid circles denote literature data used during modelling including the interpolated flux density at 150 MHz from the MSSS survey. For the sake of completeness data from the individual MSSS sub-bands (brown rectangles) and flux densities determined during commissioning phase of the 42-element Allen Telescope Array (Williams & Bower 2010) (orange triangles) are also shown.

In the text

	Fig. 18. Synthetic radio emission maps of a starburst M 82-like galaxy; from top to bottom, rows show results at 30 MHz, 150 MHz, and 1400 MHz; from left to right, columns show results for inclination angles of 0, 30, 60, and 90 degrees. The colour scale is in mJy beam−1.
In the text

	Fig. 19. Modelled vertical profiles of the total radio emission of an M 82-like galaxy seen edge-on at 30, 150, and 1400 MHz.
In the text