An Ultralong-wavelength Sky Model with Absorption Effect

The extragalactic background is assumed to be largely isotropic, and a summary of some past results is given in Guzmán et al. (2011). Two methods are commonly used to extract it from the observed sky maps (Kogut et al. 2011; Seiffert et al. 2011). In the first approach the Galactic radiation is described by a plane-parallel model, so the total radiation is

$$\begin{eqnarray}&&T(\nu )={T}_{{\rm{E}}}^{\mathrm{iso}}(\nu )+{T}_{{\rm{G}}}(\nu )\times \csc | b| ,\end{eqnarray} \tag{ 11 }$$

where both ${T}_{{\rm{E}}}^{\mathrm{iso}}$ and T_G depend only on frequency. In the second approach, the Galactic radiation is derived from a tracer, usually the Galactic [C ii] emission:

$$\begin{eqnarray}&&T(\nu )=\sum _{i=1}^{2}{T}_{{\rm{E}},i}^{\mathrm{iso}}(\nu )+{a}_{i}(\nu )\times \sqrt{{I}_{\mathrm{CII}}},\end{eqnarray} \tag{ 12 }$$

where I_CII is from the CORE/FIRAS measurements. In practice, all coefficients are fitted from observed radio maps and these two methods give consistent results (Kogut et al. 2011; Dowell & Taylor 2018). However, whether the extracted isotropic background really originates from extragalactic sources is still being debated. Seiffert et al. (2011) extracted a background that is larger than the integrated radio emission of external galaxies (Gervasi et al. 2008). The excess could be from unknown extragalactic source populations, or could be unaccounted-for Galactic foreground (see also Subrahmanyan & Cowsik 2013).

Here we assume that the isotropic background is from extragalactic sources. We adopt the following simple model of isotropic extragalactic background:

$$\begin{eqnarray}&&{T}_{{\rm{E}}}^{\mathrm{iso}}(\nu )=1.2{\left(\displaystyle \frac{\nu }{1\mathrm{GHz}}\right)}^{-2.58}\,[{\rm{K}}].\end{eqnarray} \tag{ 13 }$$

This is the fit obtained by the ARCADE-2 balloon experiment, which observed the absolute brightness temperature of the sky between 3 and 90 GHz (Seiffert et al. 2011), minus the cosmic microwave background. This is in good agreement with the fit obtained from the Long Wavelength Array (LWA) (Dowell & Taylor 2018) below 100 MHz. Strictly speaking, "background" is the cumulative radiation from unresolved sources, so it depends on the sensitivity and angular resolution of the telescope. For simplicity we directly use the background fit given by the ARCADE-2 experiment, since in our model we do not specify any particular instrument.

In addition to absorption by the Milky Way, the extragalactic radio background would also be absorbed by the ISM in the host galaxies. Moreover, the spectrum of relativistic electrons could become softer at lower energy, resulting in a softer extragalactic radio background (Protheroe & Biermann 1996). In a recent estimate for the extragalactic radio background, Niţu et al. (2021) found the extragalactic radio background deviates from power law below ∼1 MHz. In Section 3.3 we investigate the case where the extragalactic radio background is not in power-law form.

Free–free emission also contributes to the low-frequency sky brightness and should be subtracted from the observed maps for the purpose of assessing Galactic synchrotron radiation. A free–free template could be constructed by using the Hα emission line as a tracer. After correcting the absorption and scattering, the intrinsic Hα flux could be converted into free–free flux (Dickinson et al. 2003; Lian et al. 2020). Moreover, it can be constructed from multifrequency microwave observations. In this work, we use EM and T_e maps constructed from Planck multifrequency observations (Planck Collaboration et al. 2016b) ⁶ to derive the free–free map:

$$\begin{eqnarray}\begin{array}{rcl}{g}_{\mathrm{ff}}(\nu ) & = & \mathrm{ln}\left\{\exp \left[5.960-\frac{\sqrt{3}}{\pi }\mathrm{ln}\left({\nu }_{9}\,{T}_{4}^{-1.5}\right)\right]+\exp (1)\right\}\\ {\tau }_{\mathrm{ff}}(\nu ) & = & 0.05468\,{T}_{{\rm{e}}}^{-1.5}\,{\nu }_{9}^{-2}\,\mathrm{EM}\,{g}_{\mathrm{ff}}\\ {T}_{\mathrm{ff}}(\nu ) & = & {10}^{6}\,{T}_{{\rm{e}}}[1-\exp (-{\tau }_{\mathrm{ff}})],\end{array}\end{eqnarray} \tag{ 14 }$$

where ν₉ is the frequency in units of 10⁹ Hz and T₄ is the electron temperature in units of 10⁴ K (Draine 2011). In the top panel of Figure 2, we show the 408 MHz Galactic free–free emission map constructed from the Planck observations. The bottom panel shows the Galactic free–free spectrum. Clearly, most of the free–free emission is from the Galactic plane.

Zoom In Zoom Out Reset image size

Comparing the above free–free map with observations at 408 MHz (see next subsection), we find that the free–free emission contributes 4% of the mean observed sky brightness. At ∣b∣ > 10°, the free–free fraction decreases to 1%. If we smooth both maps with a resolution of 10°, then the rms contributed by the free–free emission is 15% of the total, while at ∣b∣ > 10° the rms fraction decreases to 6%. This is consistent with Dickinson et al. (2003).

We first derive the synchrotron spectral index from the observed sky maps at frequencies where absorption is assumed to be negligible. At low frequencies, there are full-sky maps at 408 MHz (Haslam et al. 1974, 1981, 1982; Remazeilles et al. 2015) and 1.42 GHz (Reich 1982; Reich & Reich 1986; Reich et al. 2001). Here we use the free–free subtracted Haslam 408 MHz map as the base, assuming that in the absence of absorption, the brightness temperature of the Galactic synchrotron radiation follows a power-law form (Platania et al. 1998; Kogut 2012; McKinley et al. 2018),

$$\begin{eqnarray}&&{T}_{{\rm{G}}}(\nu )={T}_{{\rm{G}}}({\nu }_{* }){\left(\displaystyle \frac{\nu }{{\nu }_{* }}\right)}^{{\beta }_{{\rm{G}}}},\end{eqnarray} \tag{ 15 }$$

where ν_* = 408 MHz. We determine β_G by minimizing

$$\begin{eqnarray}\begin{array}{rcl}{\chi }^{2} & = & \displaystyle \sum _{{\nu }_{j}}\displaystyle \sum _{i}\displaystyle \frac{1}{{T}_{i}^{2}({\nu }_{j})}\\ & & \times \ {\left[{T}_{{\rm{E}}}({\nu }_{j})+{T}_{{\rm{G}},i}({\nu }_{* }){\left(\displaystyle \frac{{\nu }_{j}}{{\nu }_{* }}\right)}^{{\beta }_{{\rm{G}}}}+{T}_{\mathrm{ff},i}(\nu )-{T}_{i}({\nu }_{j})\right]}^{2},\end{array}\end{eqnarray} \tag{ 16 }$$

where T_i (ν_j ) is the observed brightness temperature of the ith pixel in the sky map at frequency ν_j , and T_G,i (ν_*) is the Galactic radiation for the ith pixel at ν_*. The summation is performed for all available pixels at all observed frequencies of our selected data.

Below several hundred megahertz, there are observed sky maps, each of which covers a part of the sky. In determining the spectral index, in addition to the Haslam 408 MHz map, we use a selection of observations that have large sky coverages: the 45 MHz sky map in combination with two observations (Guzmán et al. 2011), and the 35 MHz, ⁷ 38 MHz, 40 MHz, 50 MHz, 60 MHz, 70 MHz, 74 MHz, and 80 MHz maps observed by the LWA (Dowell et al. 2017). The Guzman 45 MHz map covers 96% of the sky, with resolutions of 46 × 24 and 36 × 36 for the southern and northern sky, respectively. It combines two surveys performed by devices of similar characteristics. All of the LWA maps cover a total of about 80% of the sky, with resolutions ranging from 48 × 45 to 21 × 20. As the different maps have different angular resolutions, we smooth all the maps with a beam with FWHM = 5° when deriving the power-law index. For all the maps we subtract the free–free emission given in Section 2.2. Our fiducial model assumes that the spectral index is a constant, but we will also investigate a model with a frequency-dependent spectral index in Section 3.3 and a model with direction dependence in the spectral index in Section 3.4.

Using the above method, we derive a spectral index β_G = −2.51 for the Galactic synchrotron radiation. This is actually quite similar to the spectral index of the extragalactic background, β_E = −2.58. In Figure 3 we plot the 1 MHz sky map that is extrapolated from the free–free subtracted Haslam 408 MHz map using β_G = −2.51 for the Galactic component and β_E = −2.58 for the extragalactic component. As this is a simple power-law extrapolation, its structure is the same as that of the Haslam 408 MHz map, except for the different temperature scale. This map could be compared with the absorption-included maps presented later in this paper.

Zoom In Zoom Out Reset image size

We tentatively adopt an axisymmetric form for the Galactic emissivity . In cylindrical coordinates, the emissivity can be written as

$$\begin{eqnarray}&&\epsilon (\nu ,R,Z)=A{\left(\displaystyle \frac{R+{r}_{1}}{{R}_{0}}\right)}^{\alpha }{e}^{-R/{R}_{0}}{e}^{-| Z/{Z}_{0}{| }^{\gamma }}{\left(\displaystyle \frac{\nu }{{\nu }_{* }}\right)}^{{\beta }_{{\rm{G}}}}.\end{eqnarray} \tag{ 17 }$$

Here r₁ is a small cutoff radius to avoid singularity at the Galactic center, and we take r₁ = 0.1 kpc. The five frequency-independent free parameters—A, R₀, α, Z₀, and γ—will be obtained by fitting the free–free subtracted Haslam 408 MHz map. The parameter β_G is generally a function of frequency and spatial position. Below we take the case of a constant β_G as our fiducial model, but also discuss the cases of a frequency-dependent β_G and a direction-dependent β_G.

Before deriving the emissivity parameters, we mask the Loop I and North Polar Spur (NPS) regions. They are generally believed to be nearby objects (see, e.g., Wolleben 2007), which do not suffer from absorption on the Galaxy scale; hence one should not use these to derive the parameters of the Galactic emission and absorption model. We fit the emissivity parameters from the free–free subtracted Haslam 408 MHz map and the fitted parameters are listed in Table 1.

Most free electrons in the ISM are either in the warm ionized medium (WIM), with typical density of ∼0.01–0.1 cm⁻³ and typical temperature of ∼10⁴ K (Gaensler et al. 2008; de Avillez et al. 2012), or in the hot ionized medium, with typical density of ∼10⁻³ cm⁻³ and typical temperature of ∼10⁵–10⁶ K (Ferrière 2001). The free–free absorption is dominated by the WIM. In addition, there are some dense H ii regions with typical density of ∼10² cm⁻³ and typical temperature of several thousand kelvins (Hindson et al. 2016). They are almost opaque to low-frequency radio radiation, and hence good for separating contributions to the synchrotron radiation from different distances along lines of sight (Su et al. 2017, 2018; Polderman et al. 2019). Moreover, around the dense H ii regions, there could be some extended H ii region envelopes (EHEs), as inferred from observations (Anantharamaiah 1986; Kassim 1989). They have lower densities (∼0.5–10 cm⁻³) but larger sizes (∼0.05–0.2 kpc) than the classical H ii regions. These could also potentially contribute much to the absorption.

There are a variety of models for the Galactic electron distribution (e.g., Gaensler et al. 2008; Gómez et al. 2001; Taylor & Cordes 1993; Cordes & Lazio 2002, 2003; Yao et al. 2017 (YMW16); Schnitzeler 2012). In the present work, we mainly focus on the general effect of the absorption, so we adopt the NE2001 model (Cordes & Lazio 2002, 2003) for the electron density, which is one of the most well-accepted models. In NE2001, the Galactic free electrons have five components, i.e., smooth components including a thick disk, a thin disk, and five spiral arms; the Galactic center component; the local ISM; a list of known dense clumps; and a list of voids. The dense clumps are mainly H ii regions around massive OB stars or supernova remnants (SNRs). The voids are low-density regions found between the Sun and some pulsars. The smooth components contain lots of free-electron clouds and their volume filling factor is η. The cloud-by-cloud density fluctuations are described by the parameter

$$\begin{eqnarray}&&\zeta =\displaystyle \frac{\left\langle {n}_{{\rm{c}}}^{2}\right\rangle }{{n}_{\mathrm{smooth}}^{2}},\end{eqnarray} \tag{ 18 }$$

where n_c is the free-electron density of the cloud, $\left\langle \right\rangle$ denotes the cloud-by-cloud average, and n_smooth is the density of the smooth component. Inside each cloud, the free electrons also have small-scale structures below the largest scale L₀. It is assumed that the power spectrum of such electron density fluctuations follows a power-law distribution with index −11/3 as in the classical Kolmogorov turbulence model (Cordes & Lazio 2002). Inside a cloud the small-scale free-electron fluctuations are described by the fractional variance

$$\begin{eqnarray}&&{\omega }^{2}=\frac{\bar{{(\delta {n}_{{\rm{e}}})}^{2}}}{{n}_{{\rm{c}}}^{2}},\end{eqnarray} \tag{ 19 }$$

where the bar represents the average inside the cloud for all scales below L₀. The EM is then

$$\begin{eqnarray}&&\mathrm{EM}=\int {\eta }^{-1}\zeta \left\langle {\omega }^{2}\right\rangle {n}_{\mathrm{smooth}}^{2}{ds}\end{eqnarray} \tag{ 20 }$$

$$\begin{eqnarray}&&=\int {L}_{0}^{2/3}{{Fn}}_{\mathrm{smooth}}^{2}{ds},\end{eqnarray} \tag{ 21 }$$

where the fluctuation parameter F is defined as

$$\begin{eqnarray}&&F={L}_{0}^{-2/3}{\omega }^{2}\zeta {\eta }^{-1}.\end{eqnarray} \tag{ 22 }$$

Hence, given the smooth density distribution by NE2001, the fluctuation parameter determines the strength of absorption/emission related to the density squared. The fluctuation parameter can be derived from the observed scattering measure of Galactic pulsars (Cordes et al. 1985; Cordes & Lazio 1991; Cordes et al. 1991; Taylor & Cordes 1993) and extragalactic active galactic nuclei and gamma-ray bursts (Cordes & Lazio 2003). In NE2001, the fluctuation parameter is given for each of the three smooth components. The default values for the thick disk, thin disk, and spiral arms are F₁ = 0.18, F₂ = 120, and F_a = 5, respectively. F₁ is the most relevant one for our purpose, while F₂ and F_a are for components on the Galactic plane.

The free–free absorption also depends on the electron temperature (see Equation (2)). The temperature of the WIM is generally ∼8000 K (Gaensler et al. 2008). Throughout this paper we adopt a constant T_e = 8000 K for all Galactic electron components. The result is actually not sensitive to this. In Figure 4 we plot the optical depth of the free–free absorption toward extragalactic sources, i.e., τ(ν, s_G) in Equation (10) at 1 MHz, based on the NE2001 model. As expected, the optical depth is large near the Galactic plane, especially in the direction of the Galactic center, while outside the Galactic plane the optical depth is small. At the Galactic plane at l ∼ 90° and at l ∼ 240°, the optical depth is small because there are gaps between the two spiral arms.

Zoom In Zoom Out Reset image size

In Figure 5, we plot the absorption-free sky map at 1 MHz, obtained by integration of our fitted emissivity Equation (17) without considering absorption. This is a symmetric distribution, as it is derived from a model with cylindrical symmetry. If we mask out the Loop I and NPS regions, where the bright nearby feature dominates, the relative difference with the actual map is at the level of ∼16%.

Zoom In Zoom Out Reset image size

With the cylindrical emissivity model and the NE2001 free-electron distribution, one can construct maps with free–free absorption. In Figure 6, we plot the sky map with absorption using the NE2001 fiducial parameters. As one can see in these figures, at the lower frequencies the absorption becomes very significant. In particular, the absorption is strongest along the Galactic plane, and as a result the Galactic disk is darker than high Galactic latitude regions at the lower frequencies. This is in stark contrast with observations at higher frequencies and with the model map without absorption (see Figure 3), where the Galactic disk is the brightest part of the sky.

Zoom In Zoom Out Reset image size

How do these maps compare with the observation? The data at this low-frequency range are scarce, but there are some global spectrum measurements, i.e., the average spectrum for the whole or a large part of the sky. In Figure 7 we plot the global spectrum derived from averaging the brightness of the maps of our model, together with some observational data points. The data were taken by the IMP-6 satellite (filled circles; Brown 1973) and the RAE-2 satellite (upward and downward triangles; Cane 1979). Besides the instrumental differences, the differences between the measurements may be due partly to the different sky areas being averaged over for the two satellites (Keshet et al. 2004). The models are plotted as curves, including the absorption-free case (dotted–dashed line), the extragalactic source (dots), the fiducial NE2001 model (F₁ = 0.18, cyan dashed line), and a model with enhanced fluctuations (F₁ = 3.0, orange solid line; see discussion below). Compared with the data, we find that the sky brightness below ∼3 MHz is overpredicted when using the NE2001 fiducial fluctuation parameter for the thick disk (F₁ = 0.18), as shown by the cyan dashed curve. This has already been noticed by several earlier works (e.g., Peterson & Webber 2002; Webber et al. 2008).

Some works (Gaensler et al. 2008; Sun et al. 2008) show that the scale height of the thick disk could be twice the default value used in NE2001. This would in principle enhance the free–free absorption. However, adopting a larger scale height of 1.8 kpc, as compared to the default 0.97 kpc, results in a decrement of the mean sky brightness by only 11% at 1 MHz, because a large fraction of the sky radiation is from high Galactic latitude regions that are not influenced much by the increase of disk thickness.

Now there are some large uncertainties in the data itself, but the downturn of the spectrum below ∼3 MHz seems robust. If we take the data points below 3 MHz at face value, a possible solution of this discrepancy is that there may be small-scale structures not modeled in the current electron distribution model that enhance the absorption. In order to reach consistent results with the low-frequency observations at ≲3 MHz, one can enhance the absorption by adopting larger fluctuation parameters for the Galactic diffuse components or alternatively by adopting a softer spectral index at lower frequencies. Since the thick disk has the largest influence on the global mean sky brightness, we test the model in which the thick disk has a larger fluctuation parameter of F₁ = 3.0, while keeping other parameters the same as the NE2001 fiducial values. Note that the fiducial value of this fluctuation parameter is F₁ = 0.18 in NE2001, so the optical depth is enhanced by more than 15 times. The result is the solid curve in Figure 7, which now agrees with the observational data below 3 MHz. The corresponding maps are shown in Figure 8, for frequencies 10, 3, and 1 MHz.

Zoom In Zoom Out Reset image size

The maps in Figures 6 and 8 begin to show differences at 3 MHz, and the two 1 MHz maps are quite different. For the 1 MHz map, in Figure 8 one can see huge shadows of the spiral arms and the Galactic center component in the sky map (the blue–cyan regions at the leftmost and rightmost parts and near the center). Note that the spiral arms not only have higher WIM density than the gap between them, but also are more crowded with discrete H ii regions (Paladini et al. 2004) and their extended envelopes. They produce especially prominent large-scale shadows even at higher frequencies, as has been recognized much earlier (Bridle 1969; Roger et al. 1999; Guzmán et al. 2011). On the other hand, for regions near the Galactic pole, and at some voids, the absorption is much weaker. Such directions would be recommended for detecting extragalactic sources. There are some dark spots, most of which are near the Galactic plane. They are H ii regions around massive stars or SNRs, with enhanced absorption. For example near the Galactic plane at l ∼ 270°, there is the Gum Nebula. It is an ancient SNR with a distance between 200 pc and 500 pc (Woermann et al. 2001). So far, we have not taken into account the radiation from these sources themselves, so the brightness at their locations may be underestimated in such models.

Zoom In Zoom Out Reset image size

The IMP-6 satellite has detected maximum radiation from directions near the Galactic poles at frequencies between 0.13 and 2.6 MHz (Brown 1973), and the modulation index (defined as the ratio between the maximum flux minus the minimum flux and the maximum flux plus the minimum flux) is ∼15% at 1 MHz. The satellite has a poor resolution of ∼100° however. We confirmed that if we smooth our 1 MHz map by a window with FWHM ∼100°, the modulation index is ∼20%. Qualitatively, IMP-6 is in line with our prediction that at ultralong wavelengths, the Galactic plane becomes darker because of the absorption, leaving brighter regions at high Galactic latitudes. However, our smoothed map shows a darker Galactic plane and brighter Galactic poles also significantly modulated by the projected spiral arms, so the direction of the maximum-flux pole is shifted from the Galactic poles. It is not clear if the anisotropy seen by IMP-6 has other causes.

Because of the free–free absorption, the sky radiation observed at ultralong wavelengths would be dominated by nearby sources. As the frequency increases, more and more contribution will be gained from further sources. Therefore, by combining observations at various frequencies, ultralong-wavelength observations provide a 3D tomography tool for measuring the electron distribution and emissivity distribution in the Milky Way. Moreover, with the synchrotron emissivity distribution, one can further derive the cosmic-ray electron spectrum and the magnetic field distribution if gamma-ray observations are analyzed jointly (see Nord et al. 2006; Polderman et al. 2020). In Figure 9, we show maps of the critical distance d₅₀, within which about 50% of the radiation we observe is emitted. Here we adopt a thick-disk fluctuation parameter of F₁ = 3.0, and from left to right, we show the results for 10, 3, and 1 MHz. For each line of sight, by comparing the observed specific intensities at different frequencies, one can derive the contribution from different distances.

Zoom In Zoom Out Reset image size

We also show the optical depth from the Z-axis and Y-axis perspectives in Figure 10. On the Galactic plane, the shape of the local optical depth distribution is irregular, with obvious anisotropy, as well as raylike features due to the nearby ISM distribution. On larger scales, the electron density increases rapidly toward the inner Galaxy region, forming a fanlike structure. From the side view, we can see there are high-transparency cones above and below the Galactic plane. The overall structure is what we would expect: there is more opacity toward the inner region of the Galaxy, but the local ISM could have significant impact.

The model discussed above is based on an emissivity model with cylindrical symmetry with respect to the Galactic center. For practical usage, it is desirable to have mock sky maps that have nonsymmetric brightness temperature variations. We therefore consider a model that includes such variations, given by

$$\begin{eqnarray}&&\begin{array}{l}{\rm{\Delta }}{T}_{{\rm{G}}}(\nu ,l,b)={\displaystyle \int }_{0}^{{s}_{{\rm{G}}}}{\rm{\Delta }}\epsilon (\nu ,R,Z,\phi ,s){e}^{-\tau (l,b)}{ds}.\end{array}\end{eqnarray} \tag{ 23 }$$

The emissivity fluctuations Δ(ν, R, Z, ϕ, s) are of course very poorly known at present, so we adopt an approximate approach as shown below. The cylindrical absorption-free sky brightness temperature is

$$\begin{eqnarray}&&{T}_{{\rm{G}}}^{\mathrm{af}}(\nu ,l,b)={\int }_{0}^{{s}_{{\rm{G}}}}\epsilon (\nu ,R,Z){ds},\end{eqnarray} \tag{ 24 }$$

and the absorption-included sky brightness temperature is

$$\begin{eqnarray}&&\begin{array}{l}{T}_{{\rm{G}}}(\nu ,l,b)={\displaystyle \int }_{0}^{{s}_{{\rm{G}}}}\epsilon (\nu ,R,Z){e}^{-\tau (l,b)}{ds}.\end{array}\end{eqnarray} \tag{ 25 }$$

A "transmittance" for each line of sight is defined as

$$\begin{eqnarray}&&{R}_{1}(\nu ,l,b)=\displaystyle \frac{{T}_{{\rm{G}}}(\nu ,l,b)}{{T}_{{\rm{G}}}^{\mathrm{af}}(\nu ,l,b)}.\end{eqnarray} \tag{ 26 }$$

We then assume that the nonsymmetric temperature variations have the same transmittance as R₁ defined above; therefore, instead of looking for a Δ in Equation (23), we directly derive the nonsymmetric variations by

$$\begin{eqnarray}&&{\rm{\Delta }}{T}_{{\rm{G}}}(\nu ,l,b)={R}_{1}{\rm{\Delta }}{T}_{{\rm{G}}}^{\mathrm{af}}(\nu ,l,b),\end{eqnarray} \tag{ 27 }$$

where ${\rm{\Delta }}{T}_{{\rm{G}}}^{\mathrm{af}}(\nu ,l,b)$ is the difference between the extrapolated sky map (using Equation (15)) and the cylindrical ${T}_{{\rm{G}}}^{\mathrm{af}}$. Note that for some large-scale features on the sky, this may not be a good approximation. For example Loop I is more likely to be a nearby object (Dickinson 2018); therefore R₁ should be higher than the average. We leave this for future work.

Zoom In Zoom Out Reset image size

We add these fluctuations to the cylindrical symmetry model sky maps and the results are shown in Figure 11. In Figure 12, we also show the evolution of Galactic absorption in frequencies of 10 to 0.1 MHz. We compare our generated map at 10 MHz with the observed map by Caswell (1976) at the same frequency. In the common sky region, the relative mean deviation is 18.6%. We also compare our generated map at 22 MHz with the observed one by Roger et al. (1999)—the mean relative difference is 19.2%.

Zoom In Zoom Out Reset image size

So far we have only considered models with constant spectral indices for both the extragalactic background and Galactic synchrotron radiation. However, the spectral index may vary with frequency. For example, the energy spectrum of interstellar cosmic-ray electrons derived from the combination of PAMELA observations at near-Earth space and Voyager observations at distant space (≳100 au) can be described by a broken power-law form, with a broken point at ∼0.1–1 GeV (Potgieter & Nndanganeni 2013; Bisschoff et al. 2019). Assuming a typical interstellar magnetic field ∼5 μG, then the broken point corresponds to synchrotron critical frequency ∼1–100 MHz. Inspired by this, in this section we investigate a model in which the synchrotron spectral index depends on frequency.

The frequency dependence of the spectral index is not yet confirmed in current observations at ≳10–22 MHz. To proceed, we tentatively assume the following frequency dependence for the spectral index, for both Galactic radiation and the extragalactic background:

$$\begin{eqnarray}&&\beta ={\beta }_{0}+{\beta }_{1}\exp \left(-\displaystyle \frac{\nu }{{\nu }_{1}}\right),\end{eqnarray} \tag{ 28 }$$

where β₀ = −2.58 and −2.51 is the constant spectral index in our fiducial model for the extragalactic background and Galactic radiation, respectively, and ν₁ is a critical frequency—for ν ≫ ν₁, β ≈ β₀; for ν ≪ ν₁, β ≈ β₀ + β₁. For simplicity we adopt the same β₁ = 0.7 and ν₁ = 1.0 MHz for both Galactic radiation and the extragalactic background.

The absorption-included sky maps with a frequency-dependent spectral index are shown in Figure 13. Here we adopt F₁ = 0.18, as the default of NE2001. Comparing them with Figure 8, we find that although the mean sky brightness is similar, the morphology is rather different, especially at lower frequencies. Ultralong-wavelength observations would provide necessary information to distinguish between these two models.

Zoom In Zoom Out Reset image size

With this frequency-dependent spectral index, the tension on the fluctuation parameter discussed in Section 3.1 can be mitigated; see Figure 14, where we plot the mean sky brightness in this model. Even for the fiducial fluctuation parameter F₁ = 0.18 in NE2001, the spectrum is consistent with the ultralong-wavelength observations.

Zoom In Zoom Out Reset image size

The spectral index could also have spatial variations. However, at low frequency only the Haslam 408 MHz map covers the full sky. To derive the full-sky spectral index map, we first generate a spectral index map by combining the Haslam 408 MHz map with the LWA maps, and another one by combining the 408 MHz map with the Guzman 45 MHz map. These two maps are incomplete around the south and north celestial poles, respectively, but for sky regions common to both maps, the two maps show similar structures, and the difference is about a 13.6% standard deviation. We produce a combined map by connecting these two spectral index maps by a smooth function. The spectral index at decl. δ is

$$\begin{eqnarray}\begin{array}{rcl}{\beta }_{{\rm{G}}}(\delta ) & = & \displaystyle \frac{1}{2}\left[1+\mathrm{erf}\left(A\displaystyle \frac{\delta -{\delta }_{0}}{| {\delta }_{0}| }\right)\right]{\beta }_{{\rm{G}}1}(\delta )\\ & & +\displaystyle \frac{1}{2}\left[1-\mathrm{erf}\left(A\displaystyle \frac{\delta -{\delta }_{0}}{| {\delta }_{0}| }\right)\right]{\beta }_{{\rm{G}}2}(\delta ),\end{array}\end{eqnarray} \tag{ 29 }$$

where β_G1 is the spectral index derived by the Haslam 408 MHz map and LWA maps, while β_G2 is the spectral index derived by the Haslam 408 MHz map and Guzman 45 MHz map. When δ ≫ δ₀, we have β ≈ β_G1, and when δ ≪ δ₀, β ≈ β_G2. We take δ₀ = −15° and A = 1.7.

All maps used here are smoothed by a Gaussian beam with FWHM = 5°. The resulting spectral index map is shown in Figure 15. At the Galactic plane the spectrum slope is shallower than the higher Galactic latitude, and the southern sky is also shallower than the northern sky. Note here the spectral index is for the frequency range of a few tens of megahertz to 408 MHz, where the free–free absorption is still small except for dense H ii regions, so it reflects the intrinsic emission spectrum.

Zoom In Zoom Out Reset image size

We then model the emissivity as

$$\begin{eqnarray}&&\epsilon (\nu ,R,Z)=A{\left(\displaystyle \frac{R+{r}_{1}}{{R}_{0}}\right)}^{\alpha }{e}^{-R/{R}_{0}}{e}^{-| Z/{Z}_{0}{| }^{\gamma }}{\left(\displaystyle \frac{\nu }{{\nu }_{* }}\right)}^{{\beta }_{{\rm{G}}}(l,b)}.\end{eqnarray} \tag{ 30 }$$

It is formally the same as Equation (17), except that now β_G is allowed to vary with direction (direction-dependent), and the value is taken from the spectral index map obtained above. Using the emissivity given above, and adding the free–free absorption and small-scale fluctuations (here we also adopt F₁ = 3.0), we finally obtain the map shown in Figure 16. Compared with that in Figure 11, the overall brightness is higher in these maps, because the pixels with steeper spectral indices become brighter when extrapolated to the lower frequencies.

Zoom In Zoom Out Reset image size

The mean sky brightness temperature is also plotted in Figure 14. In this model, for the same absorption strength (F₁ parameter) the sky brightness temperature is higher than the constant spectral index model, at frequencies below ∼3 MHz. As we see in Figure 15, the spectrum slope at higher Galactic latitude is steeper than the lower Galactic latitude, and as a result, the absorption-free brightness at high Galactic latitude is higher than that in the constant spectral index model. In addition, the free–free absorption at higher Galactic latitude is weaker. Consequently, the ultralong-wavelength sky is brighter in the direction-dependent spectral index model.

As noted above, the NE2001 model with its fiducial parameters may predict inadequate absorption to produce the observed downturn in the global radio spectrum. To overcome this difficulty we considered adopting a larger fluctuation parameter, i.e., by assuming there are more small-scale fluctuations in the WIM. We now consider another, more physical solution to this problem.

NE2001 includes 175 known H ii objects (SNRs or H ii regions) with enhanced electron densities ranging from 0.01 to 40 cm⁻³. Most of them have density ∼0.5 cm⁻³ at the center, and are modeled with Gaussian profiles of width ∼ 0.01 kpc. Of course this is not a complete list; there should be many more SNRs, probably between 1000 and 10,000 (Berkhuijsen 1984), and more H ii regions (the Wide-field Infrared Survey Explorer mission has already detected more than 8000; see Anderson et al. 2014). H ii regions could be small and dense (compact or classical H ii regions), or large and diffuse (diffuse H ii regions; e.g., Lockman et al. 1996; Anderson et al. 2018). Around the classical H ii regions there could be EHEs, with sizes extending to ∼0.05–0.2 kpc, density ∼0.5–10 cm⁻³, and temperature 3000–8000 K, as found by low-frequency radio recombination lines (Anantharamaiah 1985a, 1985b, 1986; Roshi & Anantharamaiah 2001) or absorption of the radio continuum toward SNRs (Kassim 1989; Lacey et al. 2001).

We now consider the contribution to free–free absorption if more absorbers are included, to see if they can account for the downturn of the observed global radio spectrum at ≲3–5 MHz. Although the classical H ii region is opaque to low-frequency radio flux, it is difficult for them to significantly reduce the global spectrum because their volume filling factor is too small, so here we focus only on EHEs.

We assume the spatial distribution of EHEs in the Milky Way has a form similar to that of SNRs (Strong & Moskalenko 1998):

$$\begin{eqnarray}&&q(R,Z)={q}_{0}{\left(\displaystyle \frac{R}{{R}_{\odot }}\right)}^{{\eta }_{q}}\exp \left(-\xi \displaystyle \frac{R-{R}_{\odot }}{{R}_{\odot }}\right)\exp \left(-\displaystyle \frac{| Z| }{{Z}_{s}}\right),\end{eqnarray} \tag{ 31 }$$

where Z_s = 0.2 kpc, η = 1.25, and ξ = 3.56 as constrained by the Fermi Galactic gamma-ray observations (Trotta et al. 2011). The normalization factor q₀ is set by

$$\begin{eqnarray}&&{N}_{\mathrm{EHE}}=\int {dZ}\int {dR}\,2\pi {Rq}(R,Z),\end{eqnarray} \tag{ 32 }$$

where N_EHE is the total number of EHEs in the Milky Way, and is treated as a free parameter here.

We make random realizations of the EHEs that follow the spatial distribution of Equation (31), and further assume that their size follows a log-normal distribution, with central value $\left\langle {r}_{\mathrm{EHE}}\right\rangle =0.1\,\mathrm{kpc}$ and variance ${\sigma }_{\mathrm{log}{r}_{\mathrm{EHE}}}=0.5$, and is limited in the range 0.05–0.2 kpc. For simplicity, we model them with sharp edges, and within the edge the electron density is 1.0 cm⁻³, to be consistent with the typical EHEs. We adopt the same density for all mock EHEs. Regarding the total number, we adopt N_EHE = 20,000. We artificially remove one EHE that occasionally encompasses our Sun; it will not influence the statistics of the sample.

In the top panel of Figure 17 we plot the free–free absorption optical depth to extragalactic sources, purely by our mock EHEs, at 1 MHz. We then replace the H ii clump list in NE2001 with these mock EHEs, and then recalculate the sky map. For all other parameters we still use the fiducial values in NE2001, because we want to see if the mock EHEs could account for the downturn on the global radio spectrum without increasing the WIM absorption.

Zoom In Zoom Out Reset image size

In the middle panel of Figure 17 we plot the 1 MHz sky map, and in the bottom panel we plot the global spectrum. For both of them the absorption from the WIM and mock EHEs is involved. In the bottom panel we also plot the global spectrum if only the absorption from mock EHEs is involved. Clearly, the mock EHEs cover ∼80% of the sky area and can absorb the ultralong-wavelength signal efficiently. However, when we look through the shape of the global spectrum shown in the bottom panel of Figure 17, we find that the predicted shape does not agree with the observations well.

In our random sample, all EHEs are opaque to radio signals at 1 MHz and they occupy ∼80% of the sky area. Therefore ∼20% of the flux leaks between gaps in the EHEs and can be detected unabsorbed. That is to say, the observed global spectrum is

$$\begin{eqnarray}&&I(\nu )\approx (1-{f}_{\mathrm{sky}}){I}_{0}(\nu )+{f}_{\mathrm{sky}}{I}_{0}(\nu )\exp [-\tau (\nu )],\end{eqnarray} \tag{ 33 }$$

where I₀ is the global spectrum without EHE absorption, and f_sky is the sky coverage of the EHEs (not the volume filling factor, so it could be close to 1). Because here f_sky ∼ 1, 1 − f_sky ≪ 1. When ν ≫ 1 MHz, the optical depth is still not that large. As a result the second term of the right-hand side (RHS) dominates the observed global spectrum. When ν ∼ 1 MHz, the EHEs are opaque enough so that ${I}_{0}(\nu )\exp [-\tau (\nu )]\sim 0$ and the first term of the RHS dominates. Therefore I(ν) shows a complex trend with decreasing frequency, as shown by the bottom panel of Figure 17. The WIM absorption would change the I₀(ν), so it also changes the behavior of I(ν), but not significantly; see the difference between the curves with and without the WIM in the bottom panel of Figure 17.

We further check that, to produce the downturn in such a model, we must have one low-density EHE encompassing the Sun. If our Sun is inside an EHE with a typical density of ∼0.1 cm⁻³, we can get the global spectrum to turn over at ∼3–5 MHz and be consistent with observations. However, we note that it is known that the Sun resides in a region known as the Local Bubble, which has much lower density than the abovementioned EHE (Frisch et al. 2011). So if there are many EHEs they could indeed affect the global radio spectrum significantly, but to produce a downturn consistent with observations, while also being consistent with the observed local density around the Sun, is still difficult.

We developed a Python code package named Ultralong-wavelength Sky Model with Absorption (ULSA) to generate the sky map at very low frequencies, taking into account the free–free absorption effect by free electrons in both discrete H ii clumps and the WIM. The code first sets up the model for the Galactic emissivity distribution by fitting available observational data at higher frequencies, and the free-electron distribution, then generates the absorbed full-sky map. The code files are organized as follows (all under the ULSA directory), with brief explanations in italics:

• 1.  
  spectral_index_fitting/ Spectral indices, three models are provided. spectral_index_constant.py spectral_index_frequency_dependent.py spectral_index_direction_dependent.py
• 2.  
  emissivity_fitting/ Set up emissivity model. produce_data_for_fitting.py Preprocess the data for emissivity fitting, removing nearby structures "NPS" and "Loop I." fit_emissivity_params.py Fit the emissivity parameters.
• 3.  
  NE2001/* The Galactic electron model.
• 4.  
  sky_map/ Make the sky map. produce_absorbed_sky_map.py

The computations are done in the following three steps:

• 1.  
  Spectral Index Parameter Fitting. As described in the main text, we have developed three models with different spectral indices: (i) a constant spectral index, (ii) a frequency-dependent spectral index, and (iii) a direction-dependent spectral index. The user may choose one of these. In this step the code produces the spectral index parameters for use in the next step. Once this is done, the sky map can be extrapolated from 408 MHz to lower frequencies.
• 2.  
  Galactic Emissivity Model Construction. We then derive the emissivity model parameters using the code in the emissivity_fitting/ directory. First, we preprocess the data by removing local structures such as the NPS and Loop I from the 408 MHz map; then the emissivity model parameters are obtained by fitting this map.
• 3.  
  Low-frequency Sky Map Generation. Finally, the sky map is generated with the emissivity model constructed above and the free–free absorption effect is computed using the NE2001 model.

To run the code, one can simply call the absorption_JRZ function, which serves as the drive routine, with appropriate choices of various parameters:

• >>> from ULSA.sky_map.produce_absorbed_sky_map import absorption_JRZ
• >>> f = absorption_JRZ(v, nside, index_type, distance, using_raw_diffuse, using_default_params,
• critical_dis,output_absorp_free_skymap).mpi()

The parameters in the command are as follows:

• 1.  
  v (float): the frequency in megahertz for the output map.
• 2.  
  nside (int): the Healpix NSIDE value for the sky map.
• 3.  
  index_type (str): the spectral index modeling option, one of constant_index, freq_dependent_index, and direction_dependent_index.
• 4.  
  distance (kpc): maximum integration distance along line of sight; the default is set to 50 kpc.
• 5.  
  using_raw_diffuse (bool): if False, the data will be smoothed by the Gaussian kernel; otherwise the raw data will be used.
• 6.  
  v_file_dir (dict): a dictionary structure used to specify additional input map data. Specifies the frequency of the map by the dictionary key, and the relative path of the map data as the dictionary value. The input sky map file should be in HDF5 format, such as {XX:"/dir/xxx.hdf5"}. If None, the spectral index is calculated with the existing data.
• 7.  
  using_default_params (bool): if True, it uses the default spectral index value in the code, or recalculates the spectral index value if False.
• 8.  
  input_spectral_index (array): one can specify the spectral index value by putting in an array containing the spectral index map in the direction-dependent case, or an array containing one element for the constant and frequency-dependent cases.
• 9.  
  params_408 (list): the emissivity model parameters (A, R₀, α, Z₀, and γ) obtained by fitting the Haslam 408 MHz sky map. If this parameter is omitted, the values given in Table 1 will be used as defaults. One can also specify these parameters directly by putting in the values, or force the code to refit by setting it to [0.,0.,0.,0.,0.].
• 10.  
  critical_dis (bool): if True, it calculates the half-brightness distance; otherwise this is not calculated.
• 11.  
  output_absorp_free_skymap (bool): if True, it produces an absorption-free sky map at frequency v as well.

The function returns a 2D matrix, the first column is the Healpix pixel number, and the second column gives the corresponding sky map. One can then plot this map using healpy.mollview. The default coordinate system is set to the Galactic coordinate system, but one can transform it to the equatorial coordinate system by calling the Healpix function f.change_coord(m,["G", "C"]). See the healpy documentation for more information. Additionally, some results are automatically rendered as an HDF5 file. The map is saved in a file with a name of the form freqMHz_sky_map_with_absorption.hdf5, where freq is the frequency in megahertz. This file contains two keys that are respectively named "spectral_index" for the spectral index map and "data" for the sky map at the given frequency.

An absorption-free sky map is computed as an intermediate result though not saved automatically, but if the parameter output_absorp_free_skymap is set as True, it will be saved in a file named in the format freqMHz_absor p_free_skymap.hdf5. Moreover, the half-brightness distance will be computed and saved as freqMHz_critical_dist.hdf5 if the parameter critical_dis is set as True.

In the example shown below, sky maps are produced at 1, 2, ... 10 MHz with the constant spectral index model.

All the observation maps used in this work are under the directory obs_sky_data/; one can also add new map data (e.g., at additional frequencies) there. The map file should be in the HDF5 format, with the first column indicating the HEALPIX pixel number, and the second column the sky temperature values in equatorial coordinates. Most observed data in their original form are in these coordinates; if not, one must convert them into equatorial coordinates. To use such new data, first, put the new data under obs_sky_data/; then, when calling the function, register the new data by a dictionary whose key is the frequency of this new input map, and the value is the path relative to obs_sky_data/. For example, if there is a new map at 22 MHz, one first creates a file 22MHz_sky_map.hdf5, putting it under obs_sky_data/22MHz/22MHz_sky_map.hdf5; then, when calling the absorption_JRZ function, provide v_file_dir={22:''/22MHz/22MHz_sky_map.hdf5''}. The code will then add the new data when calculating the spectral index.