© ESO, 2013

1. Introduction

M 82 is a nearby prototypical starburst galaxy with a distance of ~3.5 Mpc (Jacobs et al. 2009) in the M 81 group of galaxies. It hosts a complicated distribution of different gas phases in both its halo and violently starbursting core. Material from the core region is accelerated into the halo via a galactic superwind and dragged out of the galaxy by tidal interaction. A prominent large-scale biconical outflow has been observed in Hα(Devine & Bally 1999; Ohyama et al. 2002) and in X-rays (Strickland & Heckman 2009). Dust, neutral gas phases, and even molecular material are mixed in the ionised phase and seem to be swept along by this superwind (e.g. CO: Walter et al. 2002; H₂: Veilleux et al. 2009; dust: Ohyama et al. 2002; Leeuw & Robson 2009). Large-scale molecular and atomic streamers are anchored at the edges of the galactic disk and therefore show a completely different distribution (Walter et al. 2002; Yun et al. 1993b). These are interpreted to be remnants of the interaction with M 81 approximately 100 Myr ago (Yun et al. 1993a).

While these kinematic parameters are confirmed by a large number of observations and studies, the contribution of the magnetic field to the outflow is still a matter of debate. For normal spiral galaxies, the mean field αω-dynamo in combination with galactic winds driven by supernova explosions in the spiral arms is the favoured explanation for the observed magnetic field patterns (Beck et al. 1996; Krause 2009). However, the role of the magnetic field in starburst galaxies is completely unclear. It can be a passive one, where the magnetic field is simply dragged along the velocity gradient of the wind, or an active one, where the magnetic field influences the gas motions in the outflow and it becomes a wind driven by cosmic rays and its associated field. The total magnetic field energy in normal spiral galaxies contributes only marginally to the energy budget, but the high particle densities and possible high magnetic field strengths in starbursts can increase it to a value where it needs to be accounted for. A detailed study of these parameters is needed to constrain the inputs for simulations to understand material circulation and the observed morphology of galactic halos. This has a significant influence on models for the origin of magnetic fields and the magnetisation of the intergalactic medium in the early Universe, where starbursts and interactions are known to be more common. Starbursting dwarf galaxies in the early Universe can explain the magnetisation and enrichment of the intergalactic medium without a primordial seed field (Kronberg et al. 1999).

Current instruments lack the sensitivity as well as the resolution for a study of such objects in the early Universe, so that an extrapolation from nearby observable targets is needed. Due to its proximity, large apparent size of 112   ×   43, and the available data at nearly all wavelengths, M 82 is the ideal target for such an analysis. In addition, its interaction scenario makes it possible to examine the influence of the group medium on a galactic superwind scenario hosting a starburst.

Yun et al. (1994) and Chynoweth et al. (2008) found a large amount of H i filling the whole group as well as three distinct spurs, which connect the main galaxies of the group: M 81, M 82, and NGC 3077. Additionally, M 82 is well known for its large fraction of molecular gas in the central region (Aladro et al. 2011; Naylor et al. 2010). Observations by Greaves et al. (2000) show a giant magnetic molecular bubble that was possibly blown out by the superwind. But how far this gas distribution influences the evolution of the whole group and whether it is reaccreted to their host galaxies is still not clear. An extensive star formation in the recently formed spurs could have enriched the intragroup medium with metals and magnetic fields (de Mello et al. 2008).

In how far the superwind in M 82 is constantly fuelling the halo with fresh material is a matter of discussion. Evolution scenarios reach from an increased star-formation period lasting more than 10 Myr to the occurrence of several episodes of increased star formation over the last Myr. A detailed near- and mid-infrared spectroscopy by Förster Schreiber et al. (2003) suggests two successive starburst periods about 10 and 5 Myr ago. This would explain the enrichment of the halo with ionised polycyclic aromatic hydrocarbons (PAH⁺) observed by Kaneda et al. (2010), which would otherwise have been destroyed by the starburst of the violent radiation. Additionally extended emission has been found ~11 kpc to the north of the galaxy in Hα, X-rays, and UV (Devine & Bally 1999; Lehnert et al. 1999; Hoopes et al. 2005), which is spatially not connected to the main halo of M 82. This might be a collision of the hot superwind with a preexisting neutral cloud (Lehnert et al. 1999) and a remnant of a previously stronger star-bursting activity that was blowing material even further into the halo. The detection of hard γ-ray emission by the VERITAS Collaboration et al. (2009) further verifies a closed environment scenario since it favours the idea of M 82 being a proton calorimeter (Lacki et al. 2011). Additionally, the high particle and energy densities in the core region lead to a scenario where even charged particles accelerated by supernova explosions would not be able to escape this region and could not build up a synchrotron halo, assuming the recent activity of star formation and particle densities.

Previous studies by Seaquist & Odegard (1991) and Reuter et al. (1992) found emission, which was mostly confined to the core region of the galaxy, and discussed the possibilities for reacceleration and/or dynamo processes (Reuter et al. 1994) inside the turbulent medium of M 82. However, these studies could not sufficiently explain the diffuse emission reaching up into the halo. Improved data reduction and instrumental techniques, together with more precise measurements over the whole wavelength regime, now allow us to study this prototypical starburst galaxy in increased detail.

Since synchrotron emission from gyrating electrons is the best tracer for magnetic fields, a complementary study at several different radio wavelengths between λ3 cm and λ92 cm is needed to get a full picture of the loss processes and the magnetic field morphology of this galaxy and its surrounding halo. The high flux difference between the highly luminous core and the very diffuse halo structure, especially at longer wavelengths, makes such a study challenging for the dynamic range requirements of today’s instruments and data reduction techniques. Therefore data from the Very Large Array (VLA) at λ3 cm and λ6 cm and data from the Westerbork Synthesis Radio Telescope (WSRT) at λ22 cm and λ92 cm were used for the following study. Both telescopes are well known for their exceptionally good properties at these wavelengths.

The goal of this paper is to investigate the outflow mechanism and to describe a revised scenario for the built-up of the observed radio halo. As a first step we investigate the magnetic field strength in M 82 using an energy equipartition ansatz. These results can then be used to calculate the different loss timescales and designate the dominant physical processes influencing the wind in the two different regions. Arguments for the existing halo are discussed, and a scenario with a low filling factor of the ionised as well as magnetised medium in the core region is described. These results and the comparison with observations at other wavelengths are used to understand the outflow mechanism of this archetypal galaxy and draw conclusions for the early Universe.

This paper is organised as follows: Sect. 2 shows the reduced data and explains the special data reduction techniques used to reach the high dynamic ranges needed for this study. In Sect. 3 the total intensity maps and vertical scaleheights are presented. Section 4 shows the spectral index maps. Estimates for the magnetic field strength are given in Sect. 5, and a detailed discussion about the different cosmic ray losses is presented in Sect. 6. We summarise our outcomes in Sect. 7 and discuss them in Sect. 8.

2. Observations and data reduction

2.1. VLA observations

The λ3 cm and λ6 cm data were obtained from the archive (already published by Seaquist & Odegard (1991) and Reuter et al. 1992), combined and calibration improved for our studies. We chose only the data from the compact D-configuration to match the resolution to our WSRT data and to restore most of the diffuse halo structure of M 82.

All observations used two subbands at central frequencies of 8.465 GHz and 8.415 GHz and 4.885 GHz and 4.835 GHz, respectively. Each subband had a bandwidth of 50 MHz, resulting in a total bandwidth of 100 MHz for each of the observed wavelengths. The observations were accompanied by calibrator observations of 3C 138, 3C 286, 3C 48, 1044+719, and 0836+710 to determine the fluxes, phase corrections, polarisation leakages and polarisation angles over the time of the observations. Data reduction and flagging were done using the Common Astronomy Software Applications package (CASA). Initially calibrator gains were applied independently to each subband and observation using the flux density scale given by Perley & Butler (2013). Then self-calibration was performed in the same way, calibrating four times on phase only, followed by one simultaneous amplitude and phase calibration. For each iteration a model was derived using the Clark-Clean algorithm (Clark 1980). Inside each self-calibration loop, the clean iterations were increased to include more of the diffuse flux in the model and the time interval for solutions decreased until the integration time was reached. The data were then cleaned down to a 3σ-noise level and baseline-based corrections were applied using a long solution interval of ten minutes for the total power images to increase the dynamic range. Integrated fluxes for the total power images before and after the baseline-based corrections were applied and differed by less than 1.5%. Not enough flux was available in the Stokes Q and U images, so that the baseline-based corrections were not used for the production of the polarisation images. For each clean step, masks were used from previous data reductions iteratively, including more sources inside the field. The resulting final images of each observed wavelength were combined using an inverse square weighting, followed by primary beam correction.

This subband-based calibration strategy lowers the frequency dependence of the gain solutions as well as the flux variations of sources inside the field over the complete bandwidth. Care has to be taken for observations with low flux intensities since this method lowers the signal-to-noise for the determination of the gain solutions. This can lead to an image quality that is worse than calibrating using the whole data.

2.2. Westerbork observations

The WSRT archive was inspected for archival datasets of M 82. One dataset at 18 cm, two at 22 cm, and one at 92 cm were found. The 18 cm dataset and one of the 22 cm datasets were recorded in all four cross-correlations. Therefore these data contained polarisation information. For the other, only the XX- and YY-correlations were available. From this data, only a total power image of the 18 cm dataset has been published so far by Braun et al. (2007).

All datasets were first read into the Astronomical Image Processing System (AIPS) to apply the system temperatures. Since AIPS can normally only handle circular polarised feeds, the data headers were changed to circular polarisation during this step and changed back to linear polarisation later. The data were converted into the CASA MS-format and inspected using RFI GUI, the graphical frontend to Radio Frequency Interference (RFI) console (Offringa et al. 2010). A special flagging strategy was developed for each wavelength and source (calibrators/M 82). Since the bandpass for the WSRT is very steep at the subband edges, a preliminary bandpass was applied to the data which was not used for calibration during later stages of the data reduction. For this, the calibrator observation least affected by RFI was used for the appropriate wavelength. This procedure helped the RFI console to distinguish more easily between real RFI and bandpass effects. After this initial flagging, each dataset was again checked for leftover RFI, which was then flagged manually.

Calibrator gains were determined from the unpolarised calibrator, and polarisation leakage and angle terms were determined for the polarised calibrator. This was done for each channel separately. The model parameters for each calibrator were based on a fourth-order polynomial with coefficients A,B,C, and D from the VLA calibrator manual¹. Additionally, the models of the polarised calibrators were given a fractional polarisation p, an intrinsic polarisation angle Φ₀, and a rotation measure RM. This information was gathered from the newest Expanded Very Large Array (EVLA) measurements by R. Perley (priv. comm.). According to the equations $\begin{array}{ccc}\mathit{I}& \mathrm{=}& {\mathrm{10}}^{\mathit{A}\mathrm{+}\mathit{B}\log \mathit{\nu }\mathrm{+}\mathit{C}\mathrm{(}\log \mathit{\nu }{\mathrm{)}}^{\mathrm{2}}\mathrm{+}\mathit{D}\mathrm{(}\log \mathit{\nu }{\mathrm{)}}^{\mathrm{3}}}\\ \mathit{Q}& \mathrm{=}& \mathit{pI}\cos \left[\mathrm{2}{\left({\mathrm{\Phi }}_{\mathrm{0}}\mathrm{+}\mathit{RM}\frac{\mathit{c}}{\mathit{\nu }}\right)}^{\mathrm{2}}\right]\\ \mathit{U}& \mathrm{=}& \mathit{pI}\sin \left[\mathrm{2}{\left({\mathrm{\Phi }}_{\mathrm{0}}\mathrm{+}\mathit{RM}\frac{\mathit{c}}{\mathit{\nu }}\right)}^{\mathrm{2}}\right]\\ \mathit{V}& \mathrm{=}& \mathrm{0}\mathit{,}\end{array}$where I, Q, U, and V are the flux densities of the Stokes parameters, c is the speed of light, ν is the frequency of the channel, and using the relations for the linear feeds of the WSRT $\mathit{XX}\mathrm{=}\mathit{I}\mathrm{+}\mathit{Q} \mathit{XY}\mathrm{=}\mathit{U}\mathrm{+}\mathit{V} \mathit{YX}\mathrm{=}\mathit{U}\mathrm{-}\mathit{V} \mathit{YY}\mathrm{=}\mathit{I}\mathrm{-}\mathit{Q,}$(5)the fluxes for all four cross-correlations could easily be calculated for each individual channel. A summary of the parameters used is listed in Table 3.

This channel-based derivation of the calibrator gains minimises the effects of the 17 MHz-ripple in the pointing centre of the beam and reduces the instrumental polarisation significantly, which can vary by up to 1.5% across the frequency band (Brentjens 2008). Figure 1 shows the difference between assuming a constant gain across one subband and the channel-based calibration performed here.

After this initial calibration step, each dataset was self-calibrated with different strategies, which are described in the next subsections.

2.2.1. The 18 cm and 22 cm datasets

The first 22 cm dataset contained a single 12 h run with a total bandwidth of 80 MHz split into eight subbands centered at 1450, 1424, 1419, 1398, 1374, 1349, 1326, and 1303 MHz. Each subband was covered by 128 channels in two cross-correlations. The integration time was 30 s, which was sufficient to sample the ionospheric fluctuations. The observed calibrators for this dataset are 3C 147 and 3C 295. The second dataset contained a complete 12 h run with frequency switching between 18 cm and 22 cm every five minutes, yielding an integration time of ~360 min for each frequency band. The subbands were centred at 1753, 1737, 1721, 1705, 1689, 1673, 1657, and 1641 MHz for the 18 cm band and at 1422, 1406, 1390, 1374, 1358, 1342, 1326, and 1310 MHz for the 22 cm band. Each subband was covered by 64 channels in all four cross-correlations with an integration time of one minute. The observed calibrators for these two datasets were 3C 138 and CTD93.

After applying the calibrator gains with the strategy described in Sect. 2.2, both datasets were exported to MIRIAD for faster imaging and cleaning. Due to the extremely luminous extended core of M 82, a special self-calibration technique, described in Sect. 2.2.2, had to be used to reach a dynamic range of ~50 000. For each step in this data reduction involving cleaning, masks were used from previous data reductions to include all visible flux inside the field.

Because of the better data quality of the observations with 80 MHz bandwidth, the shorter integration time, and the similar subband frequencies, clean models from the former dataset at appropriate frequencies were used to start the self-calibration process for the dataset with 160 MHz bandwidth at 22 cm. This led to 13 usable images, which were then combined with an inverse square weighting and primary beam corrected.

The 18 cm dataset was only used for polarisation purposes and therefore no total power images were produced.

2.2.2. Calibration technique for strong extended sources

The east-west alignment of the WSRT and the observation mode with only one calibrator observation in the beginning and one at the end make self-calibration the most important step in the data reduction of high dynamic-range images. One has to make sure that especially the first models during the self-calibration process do not include spurious features, which would then be transported into images and models for later stages. These features are normally easily visible in the images since they are biased towards calibrator observations and array layout. For the strong extended source M 82, this means that artificial clean components, which elongate the source in the east-west direction, would be used.

Unresolved sources are normally very good calibrators since their total emission is covered by a clean beam. This makes them much better sources for self-calibration than strong extended sources because their clean component model is only one delta function and not multiple delta functions, resembling the complex diffuse emission in M 82-type sources. The clean components for these sources would resemble the artificial elongation in the east-west direction. Therefore, the best choice would be to construct a (u,v)-dataset of the M 82 field with a dominating flux contribution from point sources. This was achieved by using the following technique:

• 1.

  The channel-based initial calibration(Sect. 2.2) was applied to the dataset. Thisminimises the contribution of the 17 MHz ripple inthe central part of the field.

• 2.

  We applied a normal self-calibration approach for the entire field using only baselines ≥4 kλ for the first self-calibration loop with a long solution interval. This parameter was reduced as was the solution interval during further loops.

• 3.

  We subtracted M 82 from the dataset using the gains and the clean model from the last self-calibration loop. Afterwards the inverted gains from the same loop were applied, resulting in a dataset without M 82 and only the gains from the initial calibration that was applied.

• 4.

  A self-calibration on this dataset was applied using only the point sources, yielding an average gain over the entire field. After two rounds of self-calibration with long solution intervals and (u,v)-ranges ≥4 kλ, the gains were transferred to the original initially calibrated dataset containing all point sources and M 82. Then self-calibration was applied to this dataset in the normal manner, reducing the solution interval and including smaller (u,v)-ranges to improve the model.

This technique enabled us to use the high signal-to-noise ratio of the strong extended source without artificially including clean components from the east-west biasing. A comparison of image quality before and after the described technique was applied is shown in Fig. 2.

2.2.3. Polarisation calibration

Due to high rotation measures already detected by Reuter et al. (1994), the derivation of a model for all four Stokes parameters for each single subband is not sufficient. The Stokes Q and U fluxes vary significantly over the bandwidth of one subband. This would result in a wrong model for a significant number of channels in each subband, shifting flux from Q to U or vice versa and producing artificially false rotation measure components. Therefore, self-calibration for each single channel was done independently, as described by de Bruyn & Brentjens (2005). The high flux of M 82 was sufficient to reach the signal-to-noise required for a satisfactory dynamic range in the Stokes Q and U parameters.

The whole calibration was done in MIRIAD using a script and including all (u, v)-ranges in the self-calibration. Masks from the total power data reduction were used to include all diffuse emission in the model. The final self-calibration loop used an amplitude and phase calibration with a solution interval of one minute. The final images were then cleaned to two times the noise level and primary beam corrected. After that, all images were inspected by eye for significant artefacts due to calibration errors and leftover RFI. At the time of the observation, the sun was at a maximum of its activity, which led to significant signals in the polarisation images far off the pointing centre. Since these signals are heavily time- and frequency-dependent, a subtraction of the signal was not possible, so that these images had to be discarded. However, images at frequencies with a lower contribution of the signal from the sun were usable. Additionally, any images suffering from a significant elongation of the beam as a result of time-based flagging of the data of an east-west interferometer like the WSRT were not used. Images with a noise floor three times higher than the average were also excluded. This yielded only 396 (187 for the λ18 cm and 209 for the λ22 cm observations) out of 1024 final images for each polarisation product.

2.2.4. The 92 cm dataset

The 92 cm dataset contained a single 12 h exposure with a total bandwidth of 20 MHz. Each of the eight subbands was sampled by 128 channels in two cross-correlations. The central subband frequencies were set to 321, 325, 332, 336, 372, 376, 380, and 385 MHz. The observed calibrators were 3C 147 and 3C 295. The integration time was 30 s, which is sufficient to sample the ionospheric fluctuations during the night. Unfortunately, some of the observation was carried out during the day and at sunset, when the interference from the sun is especially high on the short baselines. Subtracting a clean model of the sun after the initial calibration was not successful, so that several hours on the shorter baselines had to be flagged.

We then exported the data to MIRIAD and carried out self-calibration for each subband independently, using masks during the clean step to include all diffuse emission in the field. The small bandwidth of each single subband of only 2.5 MHz made it possible to use a normal Clark-clean algorithm rather than a multi-frequency deconvolution (Sault & Wieringa 1994) to reach a dynamic range of ~10 000. The subband images were then combined with an inverse square weighting and primary beam corrected.

3. Total power emission

3.1. Morphology

The maps showing the total power emission are presented in Figs. 3–6 at wavelengths of 3 cm, 6 cm, 22 cm, and 92 cm, respectively. The morphology of the total power emission depends strongly on the observed wavelength. While the maps at 3 cm and 6 cm show only emission in the central disk and some extraplanar emission to the north and south of the very central region of M 82, the maps at 22 cm and 92 cm show the structure of the radio halo and the disk emission. Overall, the halo emission in the southern part seems to be less extended than in the northern part, but the whole halo emission correlates with the Hα- as well as the PAH⁺-morphology (Kaneda et al. 2010).

The most obvious difference between the northern and southern parts of the halo are two extensions accompanied by a suppression of emission in between them in the northern halo. While the 92 cm map shows basically the same features towards the north and south, with a smaller extent towards the south, the 22 cm morphology is completely different in the southern halo. Here the morphology is dominated by an extension south of the starbursting, which lacks emission in the southwestern part.

While the extent of the southern extensions at λ92 cm is about the same, it is different in the northern halo at λ22 cm and at λ92 cm. Here the northwestern spur seems to be more extended than the northeastern one. McDonald et al. (2002) have already mentioned a difference in the distribution of compact radio sources in the core region between the eastern and western part. The western part seems to be more densely populated by supernova remnants than by H ii regions; This situation is vice versa in the eastern part. This could mean that the western part is further evolved and has ejected the material into the halo earlier than the eastern part, where a higher supernova activity is expected within the next Myr.

Remarkable is the morphology of the radio emission along the disk of M 82. Observations by Seaquist & Odegard (1991) and Reuter et al. (1992) already showed this extent of the emission. Therefore it is most likely that the radio morphology of M 82 consists of a superposition of a disk and a biconical outflow.

While the lowest contour levels for the 3 cm, 6 cm, and 22 cm maps have similar noise levels, the large halo emission is only visible in the 22 cm and 92 cm map. One could argue that this might be partly due to lower sensitivities of the λ3 cm and λ6 cm maps to extended emission because of missing spacings. Comparing our integrated flux of 3.23 Jy ± 0.16 Jy with the single dish observations of Gregory & Condon (1991) at λ6 cm, who measured a flux of 3.96 Jy ± 0.60 Jy, may support this statement. The surrounding point sources were not subtracted for the single-dish measurement due to the lacking resolution, but adding up the flux of all point sources in our interferometer image within the beam of the single-dish telescope yields only ~20 mJy.

Hence, even if this is less than the observed difference of the two measurements, it brings the values in agreement within their errors. The largest observable scale of the VLA at λ3 cm is 3′ and 4′ at λ6 cm in comparison to the WSRT with 21′ at λ22 cm and 81′ at λ92 cm. This is similar in size to the observed structures at λ3 cm and λ6 cm respectively and larger than the observed structures at λ22 cm and λ92 cm. However, the rapid decrease of our emission at λ3 cm and λ6 cm in direction away from the galaxy is stronger than expected from missing spacings alone.

Soida et al. (2011) also used the VLA D-array configuration for their studies of NGC 5775 and could recover emission with a size of 35 at λ3 cm and λ6 cm. This is two times larger in comparison to our structure size of 1.5′ ~ 2′. A combination of their VLA measurements with single-dish data from the Effelsberg telescope recovered only 3% more flux at λ3 cm and no increase was recognised at λ6 cm. Therefore we conclude that missing spacings do not significantly influence the integrated fluxes or the overall morphology.

3.2. Total intensity distribution

To isolate the emission of the galaxy from contributions of background sources, all point-like sources exceeding three times the rms noise in the images were fitted with Gaussians in the image plane and subtracted. This method was successful for all sources excluding the double source to the southeast of M 82 due to its complicated structure and the fact that it is embedded in the diffuse emission. Therefore this source was left in the images; for parts of the data analysis where its flux contributes, the data points were not included in the analysis.

3.2.1. Radial total intensity profile

The resulting total intensity profiles along the major axis are shown in Fig. 7 for all four observed wavelengths. The individual values were obtained by regridding the total power maps to a Sloan Digital Sky Survey g′ image of the galaxy and integrating the flux over each pixel column. Since we are only interested in the overall shape of the profile and the relative differences between the wavelengths, the absolute integrated fluxes are irrelevant.

While the emission in the central part is dominated by the star-forming nucleus of M 82, the 22 cm and 92 cm observations show significant emission in the disk of the galaxy. Remarkable is a slight shift of the emission maximum to the west, which is visible in the higher frequency images due to their better resolution. This seems to be a consequence of the higher star-forming activity in the western part of the central region of the galaxy (Muxlow et al. 1994) and might even cause the difference between the extent of the spurs in the northern region. The most extended feature reaching up to ~4.5 kpc into the halo is positioned directly north of this region. Here the cosmic ray gas may be able to escape the galaxy more easily since the increased star formation simply blew away the dense medium, thereby providing a channel towards which the gas can escape from the central starburst.

3.3. Vertical scaleheights

The observed radio emission of an inclined galaxy is a mixture of the halo emission and the projected disk emission contributing to the halo emission (Dumke et al. 1995). In order to simulate the latter, we compressed the intensity distribution along the major axis by the inclination of the galaxy. The resulting distribution was then convolved with the half power beam width (HPBW) of the observation to account for beam-smearing effects. The full width at half maximum (FWHM) of a Gaussian fitted to the major axis profile of this image is then the effective beam size. This led to values of 16′′ (265 pc) at λ3 cm, 19′′ (330 pc) at λ6 cm, 21′′ (350 pc) at λ22 cm, and 50′′ (860 pc) at λ92 cm. All emission outside this distance from the major axis can be attributed to extraplanar emission.

To analyse the propagation of cosmic rays and their losses, a parameter which is independent of the total flux is needed. Therefore the vertical scaleheights were used. This approach has already been performed for other objects (Dumke & Krause 1998; Heesen et al. 2009; Soida et al. 2011) and makes the values comparable to each other.

According to Dumke et al. (1995) the fitted function in case of a Gaussian distribution is ${\mathit{W}}_{\mathrm{Gauss}}\mathrm{\left(}\mathit{z}\mathrm{\right)}\mathrm{=}\frac{{\mathit{w}}_{\mathrm{0}}{\mathit{z}}_{\mathrm{0}}}{\sqrt{\mathrm{2}{\mathit{\sigma }}^{\mathrm{2}}\mathrm{+}{\mathit{z}}_{\mathrm{0}}^{\mathrm{2}}}}\exp \left(\mathrm{-}\frac{{\mathit{z}}^{\mathrm{2}}}{\mathrm{2}{\mathit{\sigma }}^{\mathrm{2}}\mathrm{+}{\mathit{z}}_{\mathrm{0}}^{\mathrm{2}}}\right)\mathit{,}$(6)with σ being the effective beam size, w₀ the maximum of the distribution, and z₀ the scaleheight. In case of an exponential distribution, this becomes $\begin{array}{ccc}{\mathit{W}}_{\exp }\mathrm{\left(}\mathit{z}\mathrm{\right)}& \mathrm{=}& \frac{{\mathit{w}}_{\mathrm{0}}}{\mathrm{2}}\exp \left(\mathrm{-}\frac{{\mathit{z}}^{\mathrm{2}}}{\mathrm{2}{\mathit{\sigma }}^{\mathrm{2}}}\right)\\ & & \mathrm{\times }[\exp {\left(\frac{{\mathit{\sigma }}^{\mathrm{2}}\mathrm{-}\mathit{z}{\mathit{z}}_{\mathrm{0}}}{\sqrt{\mathrm{2}}\mathit{\sigma }{\mathit{z}}_{\mathrm{0}}}\right)}^{\mathrm{2}}\mathrm{erfc}\left(\frac{{\mathit{\sigma }}^{\mathrm{2}}\mathrm{-}\mathit{z}{\mathit{z}}_{\mathrm{0}}}{\sqrt{\mathrm{2}}\mathit{\sigma }{\mathit{z}}_{\mathrm{0}}}\right)\\ & & \mathrm{+}\exp {\left(\frac{{\mathit{\sigma }}^{\mathrm{2}}\mathrm{+}\mathit{z}{\mathit{z}}_{\mathrm{0}}}{\sqrt{\mathrm{2}}\mathit{\sigma }{\mathit{z}}_{\mathrm{0}}}\right)}^{\mathrm{2}}\mathrm{erfc}\left(\frac{{\mathit{\sigma }}^{\mathrm{2}}\mathrm{+}\mathit{z}{\mathit{z}}_{\mathrm{0}}}{\sqrt{\mathrm{2}}\mathit{\sigma }{\mathit{z}}_{\mathrm{0}}}\right)]\end{array}$(7)with the complementary error function $\mathrm{erfc}\hspace{0.17em}\mathit{x}\mathrm{=}\mathrm{1}\mathrm{-}\mathrm{erf}\hspace{0.17em}\mathit{x}\mathrm{=}\frac{\mathrm{2}}{\sqrt{\mathit{\pi }}}\underset{\mathit{x}}{\overset{\mathrm{\infty }}{\mathrm{\int }}}\exp \mathrm{(}\mathrm{-}{\mathit{r}}^{\mathrm{2}}\mathrm{)}\hspace{0.17em}\mathrm{d}\mathit{r}$(8)Integrations over rectangular areas with a width of 60′′ parallel to the galaxy’s minor axis were used as an input to a Gaussian, a one-component exponential fit, and a two-component one. This divides the galaxy into three strips each north and south of the major axis for the λ22 cm and λ92 cm data and one strip each for the λ3 cm and λ6 cm data, which resemble a good database for the spurs in the east and west and the outflow cone in the centre. This integration technique ensures that the fits do not suffer from small-scale fluctuations in the images.

We used the fitting routines by Dumke et al. (1995), which take the calculated effective beam size into account. Our measured values and the three different types of fit functions are presented in Figs. 8 to 11, while the deprojected scaleheights are presented in Table 5.

Data points where the background radio source in the southwest contributed to the diffuse emission were not used for the fitting routines. Errors σ_S for the data points were calculated using ${\mathit{\sigma }}_{\mathrm{S}}\mathrm{=}\sqrt{{\mathit{\sigma }}_{\mathrm{rms}}^{\mathrm{2}}\mathrm{+}{\mathit{\sigma }}_{\mathrm{B}}^{\mathrm{2}}}\mathit{,}$(9)where σ_rms is the error due to the noise in the image, which was estimated by selecting an emission-free area close to the centre of the image and calculating the standard deviation, and ${\mathit{\sigma }}_{\mathrm{B}}\mathrm{=}{\mathit{\sigma }}_{\mathrm{rms}}\sqrt{\mathit{n}}$ the calculated baselevel error in the image, with n being the number of integrated beam areas. Additionally, the standard deviation inside each integration area was calculated, and in cases where it was higher than the calculated baselevel error, σ_B was replaced by this value. This is especially critical in areas with strong gradients like the core region, where the actual fluctuations of the source over the area is dominating the error for the integration.

3.3.1. Scaleheights at λ3 cm and λ6 cm

Since the observations at λ3 cm and λ6 cm only show emission in the central part of the galaxy, the analysis of the scaleheights is mostly limited to the disk and inner halo emission. Figures 8 and 9 clearly show that the one-component and two-component exponential fits fit much better than the Gaussian fit.

The absolute values for the λ3 cm data and λ6 cm data do not differ a lot from each other. This is a clear sign for a homogeneous halo property at these wavelengths for this part of the galaxy. Therefore an average value of the scaleheights can be used by simply calculating the mean. This yields average scaleheights of 20 ± 5 pc for the thin disk and 140 ± 30 pc for the thick disk at λ3 cm and λ6 cm.

The influence of missing spacings on these values is assumed to be negligible since the fits are weighted by the relative error of the individual data points. This error normally scales with the signal-to-noise ratio, which is much higher in the central part than in the diffuse outer regions, where missing spacings might influence the integrated flux values.

3.3.2. Scaleheights at λ22 cm and λ92 cm

Looking at the emission profiles at longer wavelengths with high sensitivity shows more extended disk and halo emission. The fits in Figs. 10 and 11 clearly show that there is no function which fits all the data points, but for most cases the two-component exponential fit is the best solution. The exception is the southwestern part of the galaxy, where the emission drops very rapidly in z-direction. This seems to be a clear sign of high losses of the cosmic ray electrons in this area. No successful fit could be found for the northwestern part at λ22 cm.

Remarkable is the difference in the scaleheights between the northern and southern part of the galaxy. Using again the mean of the λ22 cm values yields 690    ±    50 pc for the northern thick disk and 390    ±    70 pc for the southern one, and for the λ92 cm emission 820    ±    170 pc and 450    ±    20 pc respectively. Since most of the galaxies in the M 81/M 82 group core are located in the south of M 82, the group medium is expected to be denser towards this direction. The reduced scaleheights in the southern part are a good argument for higher losses of cosmic rays towards the group centre. Kaneda et al. (2010) detected a stronger radiation field only in the southeastern part of the galaxy, which favours higher inverse Compton losses in this area. However, this cannot explain the overall reduced extent of the southern halo and its reduced scaleheights. A slower southern outflow would be a possible explanation, but lower starburst activity and/or density in the southern central part would be needed to explain the reduced kinetic energy input. A motion of M 82 towards the group centre and a compression of the cosmic ray gas, as has been seen in clusters of galaxies, cannot be confirmed by simulations (Yun et al. 1993a), but is also not fully ruled out due to the uncertainties in the possible evolution scenarios.

Since it is most likely that secondary electrons are produced in the core of M 82 and its surrounding region up to several hundred parsecs (see Sect. 6 for a detailed discussion), a one-component exponential fit was also used to model loss processes that are independent of the core region. No strong difference is seen between the scaleheights at λ22 cm and λ92 cm. Averaging the scaleheights for the thick disk of both wavelengths results in 610 ± 230 pc. All values for the thick disk in Table 5 are still smaller by a factor of 2−4 than the ones from normal spiral galaxies (Dumke & Krause 1998; Krause 2011).

4. Spectral index distribution

4.1. Morphology

Spectral index maps were produced between λ3 cm and λ6 cm (Fig. 12), λ6 cm and λ22 cm (Fig. 13), and λ22 cm and λ92 cm (Fig. 14) using a 5σ cutoff level for the total power images. To check the contribution of different (u,v) coverages for the different wavelengths, we produced each spectral index map once by limiting the used baselines to match the shortest and longest baselines of each other and once by using all data and convolving the final maps with the same beam after image restoration and primary beam correction. No significant difference was visible, so that the latter maps were used due to their higher sensitivity. Additionally, error images were calculated, leading to errors of σ_α = 0.02 for the core regions and up σ_α = 0.2 for regions close to the cutoffs.

Thermal emission is known to contribute less than 10% of the total flux density of M 82 at λ6 cm and does not even dominate the spectrum at frequencies up to 25 GHz (Seaquist & Odegard 1991; Klein et al. 1988), which implies a negligible thermal contribution to the measurements at λ22 and λ92 cm.

All maps show an overall flat spectral index between α =  − 0.6 and α =  − 0.8 for the central part of the galaxy, which decreases towards steeper values of α =  − 1.5 for the outer parts of the disk and halo emission. The outermost values in the maps between λ3 cm and λ6 cm and between λ6 cm and λ22 cm might be influenced by missing spacings. Therefore an overall steepening of the spectral index towards the halo cannot be neglected. The average spectral index gradient in the halo between λ22 cm and λ92 cm is relatively flat with an index of α =  − 1.2 and does not steepen with increasing distance from the midplane (Fig. 14). Two remarkable regions, one towards the northwest and one to the south, are visible between wavelengths of λ3 cm and λ22 cm, where the spectral index first steepens and then rises with increasing distance from the disk. Maps from Seaquist & Odegard (1991), who used different data for λ22 cm and λ92 cm, independently show the same characteristics at these locations.

Looking into these two regions at spectral indices derived from longer wavelengths does not show this behaviour to be due to the smaller beam size, but following the same direction deeper into the halo shows spectral indices of α =  − 2.5. While the northern region is located directly in the outflow cone, which indicates high inverse Compton, adiabatic, and/or synchrotron losses in this region, the southern region can be associated with higher inverse Compton losses towards a recently formed OB-association due to the stronger radiation field.

4.2. Disentangling the core and halo emission

The loss processes of the cosmic ray electrons are expected to be different in the high-density region of the core and in the halo. To analyse this, the flux at each wavelength was integrated once over the whole galaxy, once in the starburst region with a typical radius of 450 pc (Leeuw & Robson 2009), and once over the halo (Fig. 15). An additional measurement could be acquired from Braun et al. (2007) for the core region, but the dynamic range of their maps was not sufficient to extract a measurement for the halo.

The spectral index α is defined as $\mathit{S}\mathrm{=}{\mathit{S}}_{\mathrm{0}}{\left(\frac{\mathit{\nu }}{{\mathit{\nu }}_{\mathrm{0}}}\right)}^{\mathit{\alpha }}\mathrm{·}$(10)The slope of the fits to the intensities S in Jy at the wavelengths ν in GHz resembles the spectral index. Looking at the data points in a log-log-diagram (Fig. 15) shows that only the halo component can be fitted with this function, yielding α_h =  − 1.20 ± 0.01. Using only this function for the total emission resembles the average of the spectral index over the whole galaxy with α_a =  − 0.52 ± 0.04, but this is strongly influenced by the core region, where Eq. (10) does not provide an adequate fit.

A large amount of flux in the core region of M 82 has been resolved into ~100 independent discrete sources, which are a mixture of supernova remnants and ionised regions around newly forming star clusters (Fenech et al. 2010; McDonald et al. 2002; Wills et al. 1997). The prior authors found indications that particles in these dense regions are suffering from free-free absorption processes inside these ionised regions and in the sources itself, which reduces the observed flux for low frequencies. Under the assumption that the measured flux for the resolution of the presented observations is dominated by these dense objects in the core region, the opacity τ of the medium can be included in Eq. (10), resulting in $\mathit{S}\mathrm{=}{\mathit{S}}_{\mathrm{0}}{\left(\frac{\mathit{\nu }}{{\mathit{\nu }}_{\mathrm{0}}}\right)}^{\mathit{\alpha }}{\mathrm{e}}^{\mathrm{-}\mathit{\tau }}$(11)in case of free-free absorption of sources inside an ionised medium and $\mathit{S}\mathrm{=}{\mathit{S}}_{\mathrm{0}}{\left(\frac{\mathit{\nu }}{{\mathit{\nu }}_{\mathrm{0}}}\right)}^{\mathrm{2}}\mathrm{(}\mathrm{1}\mathrm{-}{\mathrm{e}}^{\mathrm{-}\mathit{\tau }}\mathrm{)}$(12)in case of free-free absorption in the sources itself with $\mathit{\tau }\mathrm{=}\frac{\mathrm{8.2}\mathrm{\times }{\mathrm{10}}^{-2}{\mathit{\nu }}^{-2.1}\mathit{EM}}{{\mathit{T}}_{\mathrm{e}}^{\mathrm{1.35}}}\mathit{,}$(13)where EM is the emission measure of the region in pc cm^-6 and T_e the electron temperature of the warm medium in K (Wills et al. 1997). Assuming T_e = 10⁴   K, the emission measure as well as the spectral index can be calculated by fitting Eqs. (11) and (12) to the data. Only the fit for Eq. (11) converged sufficiently, which suggests that sources surrounded by dense ionised star-forming regions are the origin for this low-frequency turnover. This results in a spectral index of α_t =  − 0.67 ± 0.01 for the total emission and of α_c =  − 0.62 ± 0.01 in the core region, with emission measures of EM_t = 1.12    ±    0.02 × 10⁵   pc   cm^-6 and EM_c = 3.16 ± 0.1 × 10⁵   pc   cm^-6 respectively. Using the calculated values results in an opacity of the core region of τ_c = 0.92 at λ92 cm.

Since the emission measure is defined as $\mathit{EM}\mathrm{=}{\mathrm{\int }}_{\mathrm{0}}^{{\mathit{s}}_{\mathrm{0}}}{\mathit{n}}_{\mathrm{e}}^{\mathrm{2}}\mathrm{d}\mathit{s,}$(14)where n_e is the electron density of the ionised medium and s₀ the line of sight through the medium, n_e can easily be calculated. Using s₀ = 900 pc results in an electron density of n_e = 18.7   cm^-3. Assuming that nearly all electrons are confined inside the star-forming regions and supernova remnants, the filling factor f_e can be calculated. Using a density of n_e = 1000   cm^-3 from Westmoquette et al. (2009) for the star-forming regions leads to a filling factor of f_e = 1.9%, which is confirmed by theoretical predictions from Yoast-Hull et al. (2013). While these calculations suffer from a contribution to the integrated area of the disk emission in the front and on the back of the galaxy as well as the unknown morphology in the core region, they are good enough for an order of magnitude estimate.

5. Magnetic field strength

Assuming equipartition between the magnetic field energy and the energy density of the cosmic rays, one way to determine the total magnetic field strength is using the revised equipartition formula by Beck & Krause (2005). Due to the complicated structure of M 82, the high dynamic range inside the galaxy in total intensity, and the complex environment it is embedded in, we used a pixel-based method to determine the magnetic field strength and cosmic ray losses. This ensured that the steep gradient over the galaxy of the total power emission was accounted for and the error for the calculation of the magnetic field strength due to the integration over the line of sight minimised. The λ22 cm fluxes were used for the entire calculation because of their superior dynamic range over the λ92 cm map and smaller influence of free-free absorption processes in the core region.

A cylindrical geometry was assumed for the synchrotron-emitting regions, for which the projected profile is a Gaussian with constant values along the minor axis and a Gaussian profile along the major axis. We calculated the FWHM for these geometries by convolving the 22 cm map with the 92 cm beam to match the resolution to the spectral index map. This made the result comparable. An integration of the emission along the minor axis and using a Gaussian fit to the resulting profile yielded a FWHM of 930 pc. Moreover, the assumed non-thermal spectral index was used from the halo fit in Fig. 15 since the equipartition between the magnetic field and cosmic ray energy is assumed to take several Myr. This leads to reliable magnetic field strengths only in the halo region. To estimate the magnetic field strength in the core region, we used a spectral index of α =  −0.6 from the spectral index map between λ6 cm and λ22 cm (Fig. 13), which is less affected in this region by ionisation losses. The corresponding fluxes and intensities were used from the λ22 cm map.

The pixel-based technique of calculating the magnetic field strength is only applicable for the total intensity emission as an estimate of the total magnetic field strength. The polarisation and its implication for the ordered magnetic field in M 82 will be discussed in Adebahr et al. (in prep., hereafter Paper II). The number density ratio of protons to electrons was assumed to be constant at a value of K₀ = 100. This assumption may not be fulfilled due to the possible high fraction of protons suffering less energy losses due to their higher mass than the electrons on their path away from the disk. In contrast to that, the production of secondary electrons via pion decay would lead to a decrease of K₀. It is very likely that both effects are viable and could compensate each other. The uncertainty here scales only approximately with the fourth root of the assumed ratio, so that even for values between K₀ = 10 to K₀ = 1000 the error is less than ±30%.

The map resulting from these calculations is presented in Fig. 16. Averaging the magnetic field strength over the whole galaxy results in a mean equipartition field strength of B_eq ≈ 35   μG for the assumed geometry. This is high but still close to the values for spiral galaxies, which were calculated as having 9 ± 3   μG from a sample of 74 galaxies (Niklas et al. 1995). For galaxies with higher star formation rates like M51 the average field strengths are 20   μG (Fletcher et al. 2011). Looking into the central region of M 82 gives an even higher magnetic field strength of B_eq ≈ 98   μG. While this may be uncertain by an order of magnitude for the above-mentioned reasons, it is of the same order as the value derived by Klein et al. (1988) (B_eq ≈ 50   μG). Lacki & Beck (2013) suggested a revised equipartition formula for starburst galaxies, which includes pion decay processes and the generation of secondary electrons. They calculated a magnetic field strength of B_eq = 250   μG, but this value was derived by using the integrated flux over the whole galaxy, not distinguishing between the core and halo emission, and using a constant scalelength of l = 500 pc. As a result, it is higher.

The calculated average halo field strength for M 82 is B_eq ≈ 24   μG. Some rising gradient is visible towards the edges of the map. Since the resulting field strength is proportional to (1/l)^α, the assumed geometry influences the values. A more spherical geometry would result in a more homogeneous magnetic field gradient over the halo. Since the geometry of the synchrotron-emitting medium is mostly unknown, the errors are quite huge. We stress, however, that the cylindrical geometry seems to be the preferred one for M 82 due to its correlation to other wavelengths like Hα(Ohyama et al. 2002) or X-ray (Strickland & Heckman 2009). These Hα observations also showed two different components of the outflow cone, one with a small opening angle and a cylindrical geometry and the other with a larger opening angle and a plume-like geometry. A magnetic field morphology following the latter one would reduce the magnetic field strength in the mentioned regions.

6. Cosmic ray electron losses

Cosmic ray electrons, which are assumed to be produced in star-forming regions, mainly suffer energy losses from five different processes: synchrotron radiation, the inverse Compton (IC) effect, non-thermal bremsstrahlung, ionisation, and adiabatic losses. In order to get an estimate of their contribution to the total losses, the different timescales were calculated for the derived values for the core and the halo region using the equations given by Longair (2011). All derived values are listed in Table 6.

Electrons gyrating in a magnetic field B   (μG) of a given strength and radiating at a frequency ν(Hz) have the energy E(GeV) according to the equation $\mathit{E}\mathrm{\approx }\mathrm{2.5}\mathrm{\times }{\mathrm{10}}^{-4}{\mathit{\nu }}^{\frac{\mathrm{1}}{\mathrm{2}}}{\mathit{B}}^{\mathrm{-}\frac{\mathrm{1}}{\mathrm{2}}}\hspace{0.17em}\mathrm{GeV}\mathit{.}$(15)The synchrotron lifetime τ_syn can then be calculated using ${\mathit{\tau }}_{\mathrm{syn}}\mathrm{\approx }\mathrm{8.35}\mathrm{\times }{\mathrm{10}}^{\mathrm{9}}{\mathit{E}}^{-1}{\mathit{B}}^{-2}\hspace{0.17em}\mathrm{yrs}\mathit{.}$(16)To calculate the inverse Compton lifetime τ_IC, an energy density U_ph   ⁽erg   cm^-3) is needed: ${\mathit{\tau }}_{\mathrm{IC}}\mathrm{\approx }\mathrm{3.55}\mathrm{\times }{\mathrm{10}}^{-4}{\mathit{E}}^{-1}{{\mathit{U}}_{\mathrm{ph}}}^{-1}\hspace{0.17em}\mathrm{yrs}\mathit{.}$(17)This can be calculated using the radio-far-infrared (FIR) correlation from Helou et al. (1988): ${\mathit{L}}_{\mathrm{FIR}}\mathrm{=}\mathrm{1.26}\mathrm{\times }{\mathrm{10}}^{-14}\mathrm{·}\mathrm{·}\mathrm{·}\mathrm{(}\mathrm{2.58}\hspace{0.17em}{\mathit{S}}_{\mathrm{60}}\mathrm{+}{\mathit{S}}_{\mathrm{100}}\mathrm{)}\hspace{0.17em}{\mathrm{erg}\mathrm{cm}}^{-3}\mathit{,}$(18)with L_FIR being the FIR luminosity in W m^-2 and S_60/100 the FIR luminosity in Jy at wavelengths of 60   μm and 100   μm. Using a value of S₆₀ = 630 Jy and S₁₀₀ = 529 Jy for the central 1′ radius and S₆₀ = 1700 Jy and S₁₀₀ = 1840 Jy for a radius of 4′ representing the whole galaxy (Kaneda et al. 2010) leads to an energy density of U_ph,c = 4.43 × 10^-11   erg   cm^-3 for the core and U_ph,h = 5.57 × 10^-12   erg   cm^-3 for the halo.

To estimate the non-thermal bremsstrahlung lifetime τ_brems and the ionisation lifetime τ_ion, the ISM particle density n   (cm^-3) is needed: ${\mathit{\tau }}_{\mathrm{brems}}\mathrm{\approx }\mathrm{4}\mathrm{\times }{\mathrm{10}}^{\mathrm{7}}{\mathit{n}}^{-1}\hspace{0.17em}\mathrm{yrs}$(19)${\mathit{\tau }}_{\mathrm{ion}}\mathrm{\approx }\mathrm{9.5}\mathrm{\times }{\mathrm{10}}^{\mathrm{6}}{\mathit{n}}^{-1}\mathit{E}\hspace{0.17em}\mathrm{yrs}$(20)Using n_c = 250   cm^-3 for the core region from Weiß et al. (2001) gives a lower limit because Naylor et al. (2010) have reported particle densities of up to several 10⁴   cm^-3. This would reduce τ_brems and τ_ion drastically, making them the dominant loss processes for the core region. The measured HI mass of 0.75 × 10⁹  M_⊙ from Yun et al. (1994) is used for the halo. Assuming a spherical volume with a radius of 5 kpc leads to a particle density in the halo of n_h = 0.058   cm^-3.

The adiabatic losses are only dependent on the derivative of the velocity in direction of the minor axis: ${\mathit{\tau }}_{\mathrm{ad}}\mathrm{=}\mathrm{3}{\left(\frac{\mathrm{d}\mathit{v}}{\mathrm{d}\mathit{z}}\right)}^{-1}\hspace{0.17em}\mathrm{yrs}\mathit{.}$(21)The values in the literature differ between velocities of v = 1200   km   s^-1(Reuter et al. 1992) for a non-reaccelerated cosmic ray transport and v = 4000   km   s^-1(Seaquist & Odegard 1991) for one with reacceleration. However, with the close correlation between the radio continuum morphology and the kinematic studies using charged particles, we can assume that the magnetic field is frozen into the ionised medium and is moving at the same speed. Therefore we can take values of v = 600   km   s^-1 at a height above the disk of z = 450 pc from Bland & Tully (1988), which is very close to the value derived for starbursting galaxies by Martin (1999) of v = 500   km   s^-1. The cosmic ray bulk speed is an addition of the outflow speed of the charged particles and the Alfvén speed v_A. The latter can be calculated using the standard textbook formula ${\mathit{v}}_{\mathrm{A}}\mathrm{=}\mathrm{2.18}\hspace{0.17em}{\mathit{n}}_{\mathrm{e}}^{\mathrm{-}\frac{\mathrm{1}}{\mathrm{2}}}\mathit{B}\hspace{0.17em}{\mathrm{km}\hspace{0.17em}\mathrm{s}}^{-1}\mathit{.}$(22)This results in a value of v_A = 10   km   s^-1 using n_e = 30   cm^-3 from McKeith et al. (1995), showing that in this environment the Alfvén speed contributes negligibly to the bulk speed. Even using the value found in the cap region by Devine & Bally (1999) of n_e = 0.1   cm^-3 as a lower limit would result in an Alfvén speed of v_A = 165   km   s^-1, which is still much lower than the outflow speed.

The adiabatic losses are still 1−2 orders lower than the ionisation and bremsstrahlung losses, which dominate the loss processes in the core region. This supports our interpretation of the lower fluxes in this region at longer wavelengths.

The VERITAS Collaboration et al. (2009) detected γ-ray emission from M 82 with a luminosity of L_1   GeV ≈ 1.9 × 10⁴⁰   erg   s^-1 at 1 GeV. Scaling this value down to 100 MeV using their spectral index of α ≈ 2.2 results in L_{100   MeV} ≈ 8.3 × 10³⁹   erg   s^-1. This value is of the same order as the predicted value of L_{100   MeV} = 3 × 10³⁹   erg   s^-1 using L_{100   MeV} = 2 × 10^-5  L_IR(Thompson et al. 2006) in case of pion decay. This would produce secondary electrons in the star-forming regions with the same spectrum as the protons and additionally flatten the spectrum towards lower frequencies. From the particle density n, one can calculate the pion loss timescale (Mannheim & Schlickeiser 1994): ${\mathit{\tau }}_{\mathit{\pi }}\mathrm{\approx }\mathrm{2}\mathrm{\times }{\mathrm{10}}^{\mathrm{5}}{\left(\frac{\mathit{n}}{\mathrm{250}\hspace{0.17em}{\mathrm{cm}}^{-3}}\right)}^{-1}\hspace{0.17em}\mathrm{yrs}$(23)If this is shorter than τ_esc, the starburst region is a proton calorimeter and pion decay has to be taken into account. This is the case for the calculated values. But care must be taken here with the estimated values for n since no filling factor was included in Eq. (23). It is very unlikely that these secondary electrons are uniformly produced in the whole core region due to the high clumpiness of the neutral as well as ionised gas.

We can now calculate if the cosmic rays can escape the galaxy using the cooling timescale τ_cool, which can be estimated for a continuous injection of cosmic rays by combining all the timescales: ${\mathit{\tau }}_{\mathrm{cool}}^{-1}\mathrm{=}{\mathit{\tau }}_{\mathrm{syn}}^{-1}\mathrm{+}{\mathit{\tau }}_{\mathrm{IC}}^{-1}\mathrm{+}{\mathit{\tau }}_{\mathrm{brems}}^{-1}\mathrm{+}{\mathit{\tau }}_{\mathrm{ion}}^{-1}\mathrm{+}{\mathit{\tau }}_{\mathrm{ad}}^{-1}\mathrm{+}{\mathit{\tau }}_{\mathit{\pi }}^{-1}\mathit{.}$(24)This timescale can then be compared to the escape time for the cosmic rays due to advection with a galactic wind, assuming a constant acceleration of the galactic wind from the core region into the halo: ${\mathit{\tau }}_{\mathrm{esc}}\mathrm{\approx }\mathrm{1.02}\hspace{0.17em}\mathit{h}{\mathit{v}}^{-1}\hspace{0.17em}\mathrm{yrs}\mathit{,}$(25)where h is the scaleheight of the disk, which was calculated to h = 70  pc from the average of the first component of the two-component fit to the scaleheights for the core region and to h = 670 pc from the exponential fit to the halo (Sect. 3.3.2), and v the bulk speed of the outflowing cosmic ray gas.

Comparing now the different timescales shows that synchrotron losses are not the dominant loss process for the cosmic ray electrons in the core region, even assuming energy equipartition. Since τ_cool is an order of magnitude smaller than τ_esc, the cosmic ray electrons cannot escape the recent starburst directly from the starbursting region into the halo within a normal galactic wind. Even more important is that τ_cool in the halo is of the same order as τ_esc, making the built-up of a kiloparsec size halo difficult.

This shows that there is no easy explanation for the existence of the large radio continuum halo around M 82 and that more complicated processes like secondary electrons originating from pion decay and multiple starbursting episodes have to be taken into account.

6.1. Building up the halo

Using a simple hydrostatic equilibrium approach from Thompson et al. (2006), the magnetic field strength for the core region can be calculated with ${\mathit{B}}_{\mathrm{HE}}\mathrm{=}\sqrt{\frac{\mathrm{8}{\mathit{\pi }}^{\mathrm{2}}\mathrm{G}}{\mathit{\eta }{\mathit{f}}_{\mathrm{g}}^{-1}}}{\mathrm{\Sigma }}_{\mathrm{g}}\mathit{,}$(26)where G is the gravitational constant, η the ratio for which the gravitational and magnetic energy are in equilibrium, f_g the complete gas fraction, and Σ_g the surface density. Assuming a value of η = 1 to set an upper limit for the magnetic field strength, a gas fraction of ${\mathit{f}}_{\mathrm{g}}\mathrm{=}\frac{\mathrm{1}}{\mathrm{3}}$(Förster Schreiber et al. 2001), and a surface density of 0.85 g cm^-2(Kennicutt 1998) results in an average magnetic field strength of B_HE = 1.13 mG. This would reduce τ_syn to 2.36 × 10⁴ yrs and increase τ_IC to 2.89 × 10⁷ yrs, making synchrotron losses the dominant loss process in this region. Using these values and an average bulk speed for the cosmic rays of v_c = 1000   km   s^-1 in the core region results in a travel distance of s_c = 24 pc within the loss timescale, which is too short to build up the observed extended synchrotron halo.

A possible explanation results from a scenario with a filling factor of the magnetic field similar to the ionised H ii medium of f_e = 1.9%. Using our value of B_eq = 98   μG for the core region results in an average magnetic field strength of B ≈ 5 mG for the individual H ii regions, which yields a synchrotron loss timescale of τ_syn = 2.42 × 10³ yrs. Using the same value for the bulk speed as above gives a travel distance of s_c = 2.5 pc. Muxlow et al. (1994) found an average size of 2.4 ± 1.0  pc for the compact sources in M 82. This value is by more than one order of magnitude smaller than the one for compact sources in the Large Magellanic Cloud or the Milky Way, and therefore magnetic field strength in the mG regime is a reasonable value. It was interpreted as being smaller due to its young population, but it can also be smaller due to the higher pressure of the medium surrounding those complexes and constraining their size. This would additionally support the arguments for a low filling factor of the ionised medium and the magnetic field.

Figure 17 shows the dependency of the cosmic ray energy vs. their free path lengths in the case of an exploding supernova inside an H ii region using the above-mentioned field strength and velocity at frequencies between 100 MHz and 10 GHz. A comparison with the average compact source size (dashed line) shows that only electrons injected at frequencies ν < 1.4 GHz can escape the H ii regions. Since synchrotron emission from the halo is only visible at λ22 cm and λ92 cm, this explanation describes the observations well.

The pressure to keep these compact regions confined would then originate from the hot X-ray gas as well as from the molecular and neutral gas filling the cavity in between them with a filling factor close to 100%. The gas density and photon pressure is up to 10⁴ times higher than those in the Milky Way on larger scales, resulting in a much higher filling factor for these gas phases. A simple estimate using pressure equilibrium between the ionised gas in the compact regions with a temperature of 10⁴ K and the X-ray gas with an assumed temperature of 10⁷ K leads to an energy density of 50 eV cm^-3. This is on the same order as the value derived by Klein et al. (1988) (60 eV cm^-3), and very similar to our calculated photon field strength of 28 eV cm^-3.

As we have shown, electrons produced by supernovae, in contrast to their accompanying protons, lose a large part of their injected energy over the first parsecs of their travel. As a result, the ratio between the energy densities in protons and electrons increases towards the halo. Additionally, the hot ionised gas filling the whole galaxy is charged and couples the gyrating electrons to its kinematics together with the magnetic field as soon as they leave the high magnetic field in the ionised star-forming complexes. This explains the strong correlation between the Hα, PAH⁺, and radio continuum morphologies and accounts for the observed indications concerning pion decay and secondary electrons.

A simple estimate for the travel distance in the halo with a cosmic ray bulk speed of v = 600   km   s^-1 and the cooling timescale of τ_cool,h = 1.61 × 10⁶  yrs results in a travel distance of s = 990 pc. This is still by a factor of four lower than the observed 4 kpc extent of the synchrotron halo. A radial decrease of the dust temperature in the halo directing outwards was observed by Roussel et al. (2010), which indicates lower values for the radiation field. This increases the travel distance by decreasing the inverse Compton losses, but not the adiabatic or synchrotron losses, which dominate the loss timescales in the halo and limit the travel distances to the estimated value. A dominance of synchrotron losses in the halo would additionally produce a steepening of the spectral index in this region, which is not observed. A scenario which agrees with the observations is the built-up of the halo within multiple starbursting periods, where the magnetic field and the bulk speed vary, allowing the fuelling of the halo with fresh material.

6.2. Multiple starbursting periods

Förster Schreiber et al. (2003) describe a scenario of at least two successive starburst epochs for M 82 about 10 and 5 Myr ago, in addition to the recent one, where the material was released into the halo and accelerated by shocks and the high pressure from the core region (Kaneda et al. 2010). This would also explain the observed kpc extent of the synchrotron halo. Within periods where the magnetic field declines due to less star formation and adiabatic losses are less efficient due to a lower bulk speed, the cosmic rays can leave the inner region of the galaxy and build up a halo of the observed size. This would need variations of the loss timescales with at least a factor of four over the time of a starburst period to reach heights above the disk of several kpcs. Lisenfeld et al. (1996) showed that variations in the far-infrared and radio luminosity of one order of magnitude are reasonable for young, compact starbursts like M 82.

Such an evolutionary scenario would explain not only the extent of the radio halo but also the abundance of PAH⁺ in the halo, which would, if expelled into the halo from the current starburst, have been destroyed due to the strong radiation from the radiation field (Bendo et al. 2008). Furthermore, the hot gas cap visible far to the north in Fig. 16 might be a remnant of an older starburst phase at the interface of the intergalactic medium.

The difference between the extents of the synchrotron halo at λ22/λ92 cm and λ3/λ6 cm can now be explained. The recently injected higher energetic cosmic rays would not reach the halo due to their high synchrotron losses in the compact H ii regions in the core. In this case the cosmic rays could be expelled into the halo due to their increased escape timescale, which is not possible during a starburst period such as M 82 is encountering now.

Possibilities for particle reacceleration and diffusion were discussed by Seaquist & Odegard (1991). While reacceleration by small-scale turbulence would lead to depolarisation inside the beam, large-scale shocks originating from an epoch of starburst activity can still reaccelerate the electrons in the halo and lead to observable polarisation degrees as seen by Reuter et al. (1994). On the other hand, the recent high-resolution Hα observations show small shocks inside the outflow cone, suggesting that at least on these scales the cosmic ray electrons reside inside a small-scale turbulent field. This cannot rule out that on larger scales the outflow may still be polarised. Lerche & Schlickeiser (1980) suggested a diffusion-dominated model, which would result in a flat spectral index gradient for the halo. This would agree with the observations and therefore may be the preferred one. We refer to Paper II for a detailed analysis of the polarisation data.

7. Summary

A collection of datasets from the VLA at λ3 cm and λ6 cm and from the WSRT at λ22 cm and λ92 cm was reduced to carry out a multifrequency analysis of the starburst galaxy M 82. The physical parameters in the core region as well as its influence on the halo were investigated and connected to the evolution scenario of the galaxy.

Advanced data calibration techniques for radio continuum data are needed to achieve dynamic ranges sufficient for the scientific analysis of strong extended sources like M 82. With the newly proposed technique in this publication, we were able to reach dynamic ranges of ~50 000 for λ22 cm observations and even ~10 000 for λ92 cm.

These high dynamic range images allowed us to look into the morphology of the halo at these long wavelengths. While the extent of the halo at λ3 cm and λ6 cm is limited to the inner 1 kpc radius from the core, diffuse emission up to 4 kpc radius shows up at λ22 cm and λ92 cm. A large difference in the extent of the northern and southern outflow is recognised. While the northern emission reaches up to 4−5 kpc into the halo, the southern one can only be traced up to 2.5 kpc. This asymmetry is also visible in the X-ray observations (Fig. 16). A possible explanation would be a motion of M 82 towards the group centre in the south, causing ram pressure effects as in clusters of galaxies, but such an evolutionary scenario contradicts the results by Yun et al. (1993a).

The northern outflow cone shows a lack of emission at long wavelengths in contrast to the two spurs surrounding it as well as steep spectral indices of up to −2.5 in its central part. This indicates that the dominant loss processes for the outflow cone in the north are either adiabatic or synchrotron losses. Numerical studies by Miyake et al. (2007) confirm this in the case of NGC 253, where a similar but not as well-pronounced morphology was recognized by Heesen et al. (2009).

A fitting to the scaleheights led to very different results than in spiral galaxies. While for spiral galaxies typically a thin and a thick disk are found with scaleheights of 0.3  kpc and 1.8  kpc respectively (Krause 2011), the values for M 82 were lower by a factor of ~3. Additionally, an overall small gradient for the spectral index was found for the halo between λ22 cm and λ92 cm with α =  −1.2−1.5, while an increase is expected for galaxies where inverse Compton and synchrotron losses dominate the cosmic ray losses. Using an energy equipartition ansatz between the magnetic field energy and the cosmic ray particles results in an average field strength of 24   μG for the halo and 98   μG for the core region.

An integration over the emission of the core and the halo showed that free-free absorption lowers the flux of the core at wavelengths from λ22  cm on, while the halo emission is completely free of this feature, at least up to λ92 cm. This indicates a high contribution of bremsstrahlung losses in the core region. From a fit to the data of the core, an average electron density of 19   cm^-3 could be obtained, leading to a filling factor of ~1.9% for the ionised hydrogen. These fits are limited to just five datapoints and therefore have to be taken with care. However, additional observations between λ22 cm and λ92 cm with the WSRT and even at lower wavelengths with LOFAR will give better constraints on this.

A calculation of the loss timescales shows that at the current stage the ionisation and bremsstrahlung losses are dominating the expansion of the cosmic rays in the core region, which confirms the previous statement. Additionally, the recent γ-ray detections support the thesis of M 82 encountering pion decay in the core region; this indicates very high particle densities, but lacks the resolution to determine the exact origin of these particles within the galaxy.

To account for the above observations, we pointed out a scenario where H ii regions and supernova remnants dominate the radio flux in the core region. Cosmic rays are generated inside these small regions, where the medium is so dense that secondary electron production via pion decay is possible and cosmic ray electrons suffer high losses due to strong magnetic fields. Using the same filling factor for the magnetic field as for the ionised medium, the magnetic field strength inside these compact regions was estimated to 5.16 mG. This makes synchrotron losses the dominant loss process inside these regions and the free path lengths for these particles at ν ≥ 1.4  GHz are smaller than the size of these regions. This explains the difference of the extent of the synchrotron halo between frequencies of λ3 cm/λ6 cm and λ22 cm/λ92 cm. Higher energetic particles lose their energy faster than lower energetic ones and are not propagated out into the halo. Since the pion decay occurs mostly inside these dense regions, it cannot help but produce electrons less affected by these high synchrotron losses.

In contrast to this, heavier charged particles like protons or molecules are less affected by the magnetic field and therefore can escape from these small regions more easily. The pressure to constrain these regions to sizes of several pc can be generated by the hot ionised X-ray and the neutral hydrogen and molecular gas filling the whole core region of the galaxy. As soon as charged particles reach the lower density regions surrounding the star-forming regions, they become coupled to the kinematic flow of the lower density X-ray gas and may leave the core region, which explains the similarity in the morphologies of Hα, PAH⁺, and radio.

The loss timescale in the halo is of the same order as the escape timescale. Using a simple scenario where cosmic rays are propagated from the core region into the halo results in free path lengths of only 990 pc, which is too low to build up a 4 kpc size halo. This gives evidence for a scenario where M 82 experienced several starburst periods in the last 5−10 Myr, in which the FIR- and radio-luminosity varied by one order of magnitude and it was possible to eject cosmic rays into the halo. The small gradient in the spectral index in the halo indicates that different propagation effects than in normal spiral galaxies are dominant. It cannot explain a halo where a continuous propagation of cosmic rays suffering from inverse Compton and synchrotron losses is fuelling the halo. Instead we propose that the synchrotron halo was produced by multiple ejections of cosmic rays and the observed spectrum is a superposition of several populations of those.

8. Conclusions

The described scenario could imply that the magnetisation of the intergalactic medium by starbursts is limited to intervals shorter than the actual time of the starburst. This is because losses for the cosmic rays are too high to let cosmic ray electrons escape the starbursting galaxy as soon as a strong magnetic field and high densities have been reached. Only heavier charged particles might still escape the galaxy due to their longer loss timescales. This would have a significant influence on the proton-to-electron ratio in the intergalactic medium and significantly influence theoretical predictions of the propagation of cosmic rays and the strength of magnetic fields in the intergalactic space. In simulations of three interacting galaxies Kotarba et al. (2011) showed, that the whole group medium can become magnetised over a timescale of several hundred Myr.

This evolution scenario and the discovered high clumpiness of the medium in M 82 even has an influence on predictions for the early Universe. Due to the overall increased star-formation rate and number of starburst galaxies, the release of magnetic fields into the intergalactic medium may not be as easy as predicted by current simulations. Additionally, the high-pressure environments in such galaxies might sufficiently suppress the expansion of supernova shells. Hence the dominant mechanism for material leaving the originating region are large-scale shocks from only several supernova explosions inside an H ii region. To verify our evolution scenario further, it is necessary to have numerical magnetohydrodynamic simulations which use the discussed loss processes together with the constraints given here and in several other publications.

With LOFAR becoming more and more scientifically usable and complementary observations from the WSRT and VLA with wavelengths of up to λ92 cm available for a lot of nearby galaxies, free-free absorption at long wavelengths is a viable technique to estimate thermal electron density and the filling factor of the ionised medium of star-forming spiral and starbursting galaxies. The integrated spectrum over whole galaxies does not sufficiently account for the difference between the complicated physics in the star-forming regions and the halo and interarm areas. Due to their high resolution, the mentioned observations are able to differentiate between these regions and provide new input for theoretical predictions and numerical simulations.

---

Acknowledgments

The authors wish to express thanks to Rainer Beck, John S. Gallagher, and V. Heesen for useful scientific discussions and to Mark Westmoquette for providing the Hα map and additional input to the publication. We want to express our gratitude to Carlos Sotomayor from the Astronomisches Institut der Ruhr-Universität Bochum for his help in developing scripts for the WSRT data reduction. The work at the AIRUB is supported by the DFG through FOR 1048. The Westerbork Synthesis Radio Telescope is operated by ASTRON (Netherlands Foundation for Research in Astronomy) with support from the Netherlands Foundation for Scientific Research (NWO). The National Radio Astronomy Observatory is a facility of the National Science Foundation, operated under cooperative agreement by Associated Universities, Inc.

References

All Tables

All Figures