Evidence for Cold Plasma in Planetary Nebulae From Radio Observations With the LOw Frequency ARray (LOFAR)

We present observations of planetary nebulae with the LOw Frequency ARray (LOFAR) between 120 and 168 MHz. The images show thermal free–free emission from the nebular shells. We have determined the electron temperatures for spatially resolved, optically thick nebulae. These temperatures are 20%–60% lower than those estimated from collisionally excited optical emission lines. This strongly supports the existence of a cold plasma component, which co-exists with hot plasma in planetary nebulae. This cold plasma does not contribute to the collisionally excited lines, but does contribute to recombination lines and radio flux. Neither of the plasma components are spatially resolved in our images, although we infer that the cold plasma extends to the outer radii of planetary nebulae. However, more cold plasma appears to exist at smaller radii. The presence of cold plasma should be taken into account in modeling of radio emission of planetary nebulae. Modelling of radio emission usually uses electron temperatures calculated from collisionally excited optical and/or infrared lines. This may lead to an underestimate of the ionized mass and an overestimate of the extinction correction from planetary nebulae when derived from the radio flux alone. The correction improves the consistency of extinction derived from the radio fluxes when compared to estimates from the Balmer decrement flux ratios.


INTRODUCTION
Planetary nebulae (PNe) are detectable at a broad range of wavelengths, from X-rays up to radio frequencies.Continuum radio emission originates from thermal free-free emission of ionized elements.It traces all of the ionized ejecta in PNe.Radio emission is not affected by interstellar or circumstellar extinction caused by dust.
Radio observations constrain the physical parameters of astrophysical plasma.In particular, optically thick free-free emission allows the electron temperature to be determined from the Rayleigh-Jeans law.Optically thick free-free emission decreases very quickly as frequency squared making it difficult to detect.However, the LOw Frequency ARray (LOFAR van Haarlem et al. 2013) provides enough sensitivity and spatial resolution to image optically thick radio emission of PNe.
The electron temperature is one of the most important parameters in studying PNe.It governs the energy balance and is a very important parameter when assessing the chemical composition of PNe (Stasińska 2002).Electron temperatures can be measured from the flux ratios of collisionally excited lines (CEL) (Kaler 1986).
These depend on the electron density and rely on the accuracy of the determination of transition probabilities and collision strengths.CELs are suppressed by collisional de-excitation when the critical density is exceeded in plasma.They are also weak at low electron temperatures.Thus, electron temperatures measured from CELs are weighted towards hot regions that do not exceed the critical density, and are less sensitive to dense and cool plasma.
Electron temperatures can be alternatively derived from recombination lines (RLs) or from free-bound emission (e.g.Balmer jump).RLs are in general much fainter than CELs, and therefore more difficult to measure.However, electron temperatures derived from RLs are systematically lower than from CELs.Moreover, the abundances derived from RLs are higher than the abundances obtained from CELs (Peimbert 1971;Stasińska & Szczerba 2001;Zhang et al. 2004;Wesson et al. 2005).The difference between these two determinations is referred to as an abundance discrepancy factor.Peimbert (1971) attributed this discrepancy to temperature fluctuations in PNe.If large temperature fluctuations occur, then the temperature measured from CELs is overestimated and RLs appear stronger.However, photoionization models failed to reproduce temperature fluctuations sufficiently large enough to account for the temperature and abundance discrepancy (Kingdon & Ferland 1995).Liu et al. (2000) showed that the inclusion of highdensity hydrogen-deficient plasma can explain the RL and CEL temperature and abundance discrepancies.This has been subsequently confirmed by Tsamis et al. (2004), Zhang et al. (2005), and Wesson et al. (2005).A further insight came from Corradi et al. (2015), who linked the abundance discrepancy with binarity of their central stars.García-Rojas et al. (2019) present the most recent review of the topic.
Optically thick free-free radio emission provides an alternative method to assess electron temperatures in PNe.Brightness temperature is simply equal to electron temperature for optically thick thermal radiation.Unlike other methods, this determination does not depend on electron density.The only assumption is Maxwellian distribution of electrons in nebular plasma, which is most likely fulfilled (Draine & Kreisch 2018).In this paper we report electron temperatures for a sample of PNe using LOFAR observations of optically thick freefree emission.

RADIO EMISSION FROM PLANETARY NEBULAE
In the case of an ionized nebula with constant electron density and temperature T e (hereinafter referred to as an homogeneous nebula or a homogeneous model) which covers a solid angle of Ω, the free-free flux density is given by in the Rayleigh-Jeans approximation.The optical depth is given as τ ν = 0.0544 × T −1.5 e ν −2 g ff (ν, T e ) EM (Olnon 1975).EM denotes emission measure, k -Boltzmann constant, ν -frequency, and c -the speed of light.van Hoof et al. (2014) computed non-relativistic Gauntt factors g ff for a wide range of frequencies.
The brightness temperature depends on the surface brightness of the object The brightness temperature approaches electron temperature in an optically thick case.It is often assumed, that PNe with T B determined from Equation 2 lower than some arbitrary value (e.g. 1 kK in Ruffle et al. (2004), 3 kK in Stasińska et al. (1992)) are optically thin, which is well justified in the case of homogeneous nebula for a typical electron temperature of 10 kK.
However, the study of the radio spectral indices (SI) F (5 GHz)/F (1.4 GHz) and brightness temperatures by Phillips (2007) revealed that the majority of PNe show an excess of the F (5 GHz)/F (1.4 GHz) ratios over the value predicted by homogeneous model.Phillips (2007) has attributed this excess to the existence of strong radial density gradients in the nebulae.In such cases, nebulae become partially optically thick over a wide range of frequencies (Wright & Barlow 1975).Siódmiak & Tylenda (2001) attempted to explain the F (5 GHz)/F (1.4 GHz) index excess using an alternative approach.They used two components instead of one in Equation 1.One of the components covers only a fraction of the solid angle ξΩ and has an optical thickness of τ ν .The other component covers the rest of the solid angle (1 − ξΩ) and has an optical thickness of ητ ν :  (Masson 1990;Aaquist & Kwok 1996) has a higher F (5 GHz)/F (1.4 GHz) index compared to the homogeneous model.Other studies have shown that the prolate ellipsoidal shell model provides a better fit to the observed surface brightness distribution of PNe than the homogeneous model.However, ellipsoidal shells would have to be enormously elongated to account for the high excesses observed in some PNe.Equation 1 gives a satisfactory fit to most of the PNe using electron temperatures derived from CELs given that Ω is smaller than the observed size of the nebula, i.e. when the bulk of the emission comes from a fraction of the solid angle.This is equivalent to Equation 3 for η = 0 and 0 < ξ < 1.It is impossible to find a single value of EM which would allow fitting the optically thin and optically thick parts of the spectrum simultaneously for ξ = 1 in most cases.With higher EM Hajduk et al. (2018) could reproduce the turnover frequency, but overestimated the optically thin flux.Lower values of EM allowed us to fit the optically thin part of the spectrum, but shifted the turnover to lower frequencies than observed.

OBSERVATIONS AND DATA ANALYSIS
LOFAR is a radio interferometer which consists of 52 stations distributed in Europe.The Netherlands host 24 core and 14 remote stations operating at the shortest baselines.The remaining 14 stations are located in other countries and provide the longest baselines.Each single station consists of a set of low-band (LBA) and highband (HBA) antennas observing in 30-80 and 110-240 frequency ranges, respectively (van Haarlem et al. 2013).We used the radio continuum 120-168 MHz images (central frequency of 144 MHz) of PNe collected by the LOFAR Two-Metre Survey (LoTSS) (Shimwell et al. 2019).The survey uses only the data from core and remote stations.The collected visibilities are processed with direction-dependent calibration (van Weeren et al. 2016).The clean algorithm is replaced with a spectraldependent deconvolution algorithm, which improves dynamic range of the obtained images (Tasse et al. 2018).de Gasperin et al. (2019) presents the calibration strategy and examples.
The survey provides low-and high-resolution images with the full width at half maximum of the restoring beam being 20 and 6 arcsec, respectively.We assumed the absolute flux density scale accuracy of 10% (Shimwell et al. 2019).The median positional accuracy of the high-resolution images is 0.2 arcsec, though it may range from 0.1 arcsec to 4.8 arcsec for individual fields.LoTSS fields reach a flux accuracy of 100 − 500 µJy beam −1 .
Good sampling of the uv plane by short baselines provides LOFAR with an excellent sensitivity to extended emission.Some examples are presented in Shimwell et al. (2019).An upgraded pipeline improved the reduction of extended emission and removed artifacts which were present in the preliminary LoTSS release (Shimwell et al. 2017).With the shortest baseline of about 80 m the largest angular scale of LoTSS reaches 40 arc mins (Savini et al. 2018).
The LoTSS observations and data processing are still ongoing.We included observations which were pro-cessed before April 2021 1 .This largely overlapped with the upcoming LoTSS-DR2 (T.Shimwell et al., 2021, in prep.).LoTSS-DR2 includes overlapping fields which are mosaiced to produce the final survey images.We also included additional pointings which have not been yet mosaiced.Their quality will improve in the future after LoTSS completes observations and produces final mosaics.
We selected 165 PNe in the observed part of sky using the SIMBAD database (Wenger et al. 2000) and the catalogue by Parker et al. (2016).Out of them, 30 were detected.Table 1 presents the nebular sizes and flux densities of these PNe at 144 MHz.The fluxes and diameters of compact PNe Θ d were measured with Gaussian deconvolution using CASA (McMullin et al. 2007).We multiplied the deconvolved diameters by correction factors to account for more realistic surface brightness distribution than a simple Gaussian (van Hoof 2000).We applied the correction factors computed for disk geometry for optically thick PNe.This choice is justified by a flat surface brightness profile of a spherically symmetric nebular model at optically thick 20 cm (equivalent to frequency of about 1400 MHz) computed by van Hoof (2000).Flat surface brightness profile represents a circular, constant surface brightness disk.
The correction factors were not applied for the large PNe and for unresolved PNe.For well-resolved PNe we fitted an ellipse to the emission which exceeded background by 3σ, which is marked with a thick line in Figures 2 -5.To measure the flux density, we integrated the emission within this area.

Spectral fitting
We combined our new 144 MHz flux densities with flux densities collected at different frequencies with other surveys, which are listed in Hajduk et al. (2018).We fitted the spectra with Equation 3 for η = 0 using the derived sizes (Figures 2 -5).Only two of the three unknown parameters: ξ, EM , and T e could be fitted independently.EM is parametrized in the optical depth term.We used electron temperatures derived from CELs by Kaler (1986), leaving ξ and EM as free parameters.Using electron temperature derived from CELs or fixing it to an arbitrary value (e.g. 10 4 K) is a common practice in fitting radio spectral energy distribution (SED) of PNe (Pazderska et al. 2009;Hajduk et al. 2018;Bojičić et al. 2021).This reduces the number of unknown variables in the fit to two.However -as we will show later -cool plasma component may also contribute to radio emission and bias the results.Optical depth does not strongly depend on the assumed temperature and can be robustly derived from the fit.Lower temperature would lead to higher ξ and lower EM .
We fitted the spectra to check if PNe are optically thick at 144 MHz.In such a case electron temperature could be determined from the brightness temperature.Our fits also confirm that PN radio SEDs are consistent with free-free emission.The LoTSS source density is 770 per square degree (Shimwell et al. 2019).The probability of finding one confusing source closer than 6 arcsec in the sample of 165 objects is about 30%.Some PNe show background sources nearby, but far enough to be separated.
Table 3 lists upper flux limits measured from the maps for 135 undetected PNe.The median root-mean square (rms) is about 700 µJy/beam (Table 3).A 3 rms upper limit of 2.1 mJy/beam at 144 MHz corresponds to the brightness temperature of about 4800 K for a 6 arc sec beam.Out of 135 undetected PNe 32 objects have been detected at 1.4 GHz with the NRAO VLA Sky Survey (NVSS) by Condon & Kaplan (1998).The sensitivity of NVSS is about 450 µJy/beam.We measured upper limits for spectral indices between 1.4 GHz and 144 MHz (Table 4), assuming the PN sizes of ≤ 6 arcsec.All but two were between optically thin (SI = −0.1)and optically thick (SI = 2) limit.This indicates that most of the PNe detected in NVSS show optically thick effects at 144 MHz.However, PNe with higher optical thickness are brighter in radio and more likely to be detected.

Electron temperatures
The nebular images and spectra are shown in Figures 2  through 5. We converted the intensity scale in the images from F 144 MHz /beam to T B /beam.The converted images map electron temperature for resolved and optically thick PNe.However, for PNe more compact than the instrument beam, the peak flux and brightness temperature are dilluted by the squared ratio of the source size to the beam size (Θ/Θ beam ) 2 .
We determined average electron temperatures for well resolved PNe which are optically thick at 144 MHz.For this purpose, we substituted the measured size and the integrated 144 MHz flux to equation 2 with T B = T e .The flux and diameter uncertainties propagate to the calculated T e error.The derived temperatures are presented in Table 2 along with electron temperatures from the literature obtained using alternative methods.The electron temperatures derived from the 144 MHz images do not exceed 9.6 kK.Table 2 shows, that our derived temperatures are 20% to 60% lower than the temperatures derived from CELs of [O iii] and [N ii].They are also lower from the temperatures derived from the Balmer jump, although they agree within 1 σ in two cases (NGC 6543 and NGC 6826).The mean temperature determined from the 144 MHz images is about 7.0 kK, which is about 35% lower than the 10.7 kK mean temperature derived for [O iii].The low T e derived from 144 MHz optically thick emission results from the presence of the cold plasma component, which is observed in RLs (Liu et al. 2000).Radio flux is strongly affected by the coldest and most dense regions, even if they contain only a small fraction of the total ionized mass in PNe.The optical thickness of plasma at radio wavelengths is approximately proportional to T −1.35 e .Thus, low electron-temperature plasma has much higher opacity from hot plasma and can become a strong opacity source for low-frequency radio emission.
We modelled a radio spectrum from an ionized nebula filled with the plasma with T e of 10 kK, a typical value for PNe, with 10% of the volume filled with randomly distributed cool plasma component with T e of 2 kK.The electron temperature averaged over volume is thus 9.2 kK.Analysis of RLs confirms, that the cold component can indeed have electron temperature as low as T e ≈ 1 kK (Corradi et al. 2015).The resulting free-free radio continuum spectrum is compared to the spectrum of the homogeneous model with a temperature of 10 kK (Equation 1) in Figure 1, upper panel.
An inclusion of cold plasma increases flux emitted at optically thin high frequencies with respect to the homogeneous model (Figure 1).This is because optically thin cold plasma emits much more flux than hot plasma.The turnover is shifted to higher frequency.The inclusion of low temperature plasma reduces the brightness temperature in the optically thick part of the spectrum compared to that expected from an homogeneous model by as much as 50% at 144 MHz.Hence, in this scenario the electron temperature determined from an optically thick radio image would be around 50% lower than the temperature of the hot component.This is consistent with our observations, with T e (144 MHz) lower by 20-60% than T e ([O iii]), which represents the hot component.
The model with two plasma components is compared with the homogeneous model with lower electron temperature of T e = 5 kK (Figure 1, lower panel).The spectra appear quite similar.The turnover is approximately at the same frequency.Both models emit similar amounts of optically thick and optically thin flux.However, the optical depth effects become visible already before the turnover frequency in the two component model.This would produce an excess of the 5 GHz to 1.4 GHz flux ratio with respect to a homogeneous model.
The distribution, temperature, chemical composition, and density of cold plasma may vary from one PN to another.In every case cold plasma adds more flux to the optically thin part of the spectrum compared to the homogeneous model assuming T e derived from CELs, and suppresses optically thick flux, if it extends to the outer boundary of the PN.T e should be treated as a free parameter in modeling of radio spectra of PNe, but unfortunately it is not independent from ξ and EM .

Comparison of radio and optical emission
We compared radio continuum images of PNe with the Hubble Space Telescope (HST) and the Isaac Newton Telescope (INT), collected in most cases by the INT Photometric H-Alpha Survey (Drew et al. 2005, IPHAS) Hα images in Figures 6, 7, and 8.We convolved Hα images with 6 or 20 arcsec Gaussian in order to match the resolution of radio and optical images.Hα emission is optically thin in PNe.The brightness distribution in Hα should be similar to optically thin at 144 MHz emission, as both of them depend primarily on the emission measure.However, cold and hydrogen deficient regions should stand out in the 144 MHz images.On the other hand, the 144 MHz surface brightness distribution for optically thick PNe is not expected to correlate with the Hα image since it depends on the local electron temperature close to the outer radius of a PN.
PNe K 3-17, IC 4593, NGC 6543, NGC 6826, and NGC 7027 are optically thick in radio.IC 4593 (Figure 6) is not very well resolved and does not allow for a detailed comparison of optical and radio surface brightness distribution.NGC 7027 (Figure 6) is slightly more resolved.It shows a maximum of the radio emission in the north western part of the nebula.NGC 6826, NGC 6543 (Figure 6), and NGC 40 (Figure 8) have similar size in radio and Hα (Table 1).Low brightness temperatures of these three PNe indicate that cold plasma exists near the outer radius of the shell.Otherwise their brightness temperature would reflect the temperature of the hot plasma component, close to the [O iii] or [N ii] temperatures (Tab.2).We were able to measure temperature fluctuations (Peimbert 1967), since these three PNe are well resolved and optically thick.The fluctuations clearly exceed the 3σ noise in the 144 MHz images (Figure 6).The point-to-point temperature fluctuations measured in the plane of sky in the central part of the disk are t 2 = 0.006 ± 0.001 in NGC 6826, t 2 = 0.004 ± 0.001 for NGC 6543, and t 2 = 0.013 ± 0.006 for NGC 40.t 2 is defined as the standard deviation of the temperature distribution computed close to the neb-   ular center, so that it would not be affected by drop of the flux at the edges.We subtracted the background root mean square error (rms).In order to compute the background rms we first measured the rms in the whole image.Then we flagged all the regions exceeding 3 rms and we computed a new value of the rms in the unflagged image.We repeated the process if there were still regions exceeding 3 rms in the image.The temperature fluctuations in NGC 6543 agree with the 0.004 estimate by Wesson & Liu (2004) from optical imaging.Temperature fluctuations which could explain the abundance discrepancy would need to be one order of magnitude bigger.However, the temperature fluctuations which could be responsible for the dichotomy of abundance determination (Peimbert 1971) may exist on lower spatial scales.NGC 40 shows patchy structure in the 144 MHz image.The optical and radio images are not correlated.NGC 40 a born-again candidate (Toalá et al. 2019).In such a case, a new, hydrogen-free ejecta can be mixed with the previously ejected hydrogen-rich envelope and cause significant inhomogenities of the chemical composition and electron temperature within the nebula.
K 3-17 shows a weak trace of the bipolar structure in radio (Figure 7).The waist of the hourglass nebula is optically thick in radio, whereas the bipolar structure is optically thin.
The radio images of NGC 6720 (Figure 6) and NGC 2371 (Figure 6) resemble the Hα images.However, their radio spectra indicate considerable optical thickness at 144 MHz.The 144 MHz flux is dominated by the brightest regions in these two nebulae.This is confirmed by non-uniform brightness distributions observed in the 144 MHz images and small values of ξ derived in the fit of the radio spectra.The small value of ξ indicates that most of the emission comes from a small fraction of the nebulae, while the remaining part remains optically thin.Both nebulae have significantly smaller τ 144 MHz than PNe which are fully optically thick.
NGC 6720 is very well resolved in the LOFAR image.Both 144 MHz and Hα images show an oval ring.The brightest part, reaching a temperature of 5000 K at maximum, is optically thick.The center of the nebula and the part of the ring close to the long axis, which contribute less to the radio flux, remain optically thin.O'Dell et al. ( 2013) modelled the optical image of NGC 6720 with a triaxial ellipsoid seen nearly pole-on.The projected ring is brighter on its shorter axis.The maximum of the optical emission along the shorter axis is close to the outer edge of the ring.The maximum of radio emission is shifted toward the center of the nebula with respect to Hα emission.The reason for this can be that more cold plasma exists closer to the central star.As a result, the 144 MHz opacity increases toward the central star, so the inner part of the ring is brighter.
Figure 9 compares the 70 µm and 144 MHz continuum images of NGC 6720.The maximum of the 144 MHz emission is closer to the central star than 70 µm emission.It appears, that cold plasma is not associated with dust emitting at 70 µm.
NGC 2371 is a bipolar nebula.The brightest, barrellike structure contains a collection of knots.Two pair of brightest knots in the Hα image lie at the position angle of 60 degree in the NW and SE direction from the central star, although they are not perfectly aligned with the central star (Figure 6).The 144 MHz image of NGC 2371 does not trace the Hα emission in detail.In particular, the maxima of 144 MHz emission are not centered on the brightest knots in the Hα image, but located on fainter knots closer to the central star.This is more clearly seen when comparing a full resolution optical image with the radio emission (Fig. 10).
Gómez-González et al. ( 2020) used spatially resolved spectroscopy to study the electron temperature in NGC 2371.The brightest clump in the nebula, located in the NE direction from the central star (designed by them as A7, Fig. 10) has a temperature of 13.8 kK.It is classified as a low-ionization knot (Gonçalves et al. 2001).For comparison, the neighbouring region A6, which is the brightest region in the 144 MHz map, has significantly higher temperature of 18 kK.The A7 clump should stand out at 144 MHz since it is brighter in optical and cooler.However, it is at least 2 times fainter in radio than the A6 clump.This suggests, that the low-ionization knot A7 does not contain cold plasma or contains less cold plasma than the A6 clump.

DISCUSSION AND CONCLUSIONS
We observed 144 MHz free-free radio emission in a sample of PNe using LOFAR.Optically thick emission allows for relatively straightforward measurement of local electron temperature.The observations confirm a presence of significant amount of cold plasma, which was first proposed from the study of RLs and CELs.Cold and hot components remain spatially unresolved in the 144 MHz images.However, cold plasma has much higher opacity compared to hot plasma.In the result, the determined electron temperatures are significantly weighted toward cold plasma component.Thus, the previous approaches which assumed a homogeneous model  of PNe with the temperature derived from CEL or an arbitrary value of 10 kK are incorrect.
Different studies assumed homogeneous models with the temperature corresponding to the hot plasma component.In particular, ionized mass determination relies on the assumption that PNe are optically thin at 5 GHz and their plasma temperature is similar to T e derived from CELs (Buckley & Schneider 1995).This approximation remains valid as soon as the electron temperatures used in the computation are reduced by 40% on average.Taking this into account it would scale down the ionized masses of PNe derived from optically thin radio emission.Lower T e would also result in lower diameters derived from the radio SED fit compared to observed diameters (Hajduk et al. 2018;Bojičić et al. 2021).Stasińska et al. (1992) used an homogeneous model for extinction determination from the ratio of the optically thin radio to hydrogen flux.They used electron temperatures from Kaler (1986) or derived using his formulae.The radio flux was used to determine the dereddened Hβ 0 flux.This dereddened flux was compared with the observed Hβ flux to derive the extinction C rad .Another extinction determination C opt comes from the observed Hα to Hβ ratio.Stasińska et al. (1992) showed that C opt is systematically larger by a factor of about 1.2 than C rad for the PNe in the direction to the Galactic center.Ruffle et al. (2004) postulated steeper extinction law toward the Galactic center to explain the difference between radio and optical extinction, which was later confirmed by Hajduk & Zijlstra (2012).Finally, Pottasch & Bernard-Salas (2013) suggested that the 5 GHz emission in PNe is not optically thin, which, however, contradicted most of other studies.
The ratio of the optically thin radio to Hβ flux depends on electron temperature ∼ T 0.53 e (Pottasch 1984).If PNe were modeled with the temperature lower by 40%, the dereddened Hβ fluxes would be higher by a factor of 1.3.This would decrease C rad by 0.12 and improve the consistency of both extinction determinations, though not yet fully explain it.
The 144 MHz images allowed us to constrain the spatial distribution of cold plasma.It extends to the outer radii of the nebulae.However, the partially optically thick image of NGC 6720 shows, that cool plasma is more abundant toward the center of the nebula.It is noteworthy, that the abundance discrepancy factor also increases toward the center of the PN (Garnett & Dinerstein 2001).The 70 µm image has also different brightness distribution from radio emission, which suggests that dusty regions observed at 70 µm do not harbour cold plasma.Another example, in which radio emisson is more concentrated in the inner ring of the nebula than dust emission is the Helix nebula (Planck Collaboration et al. 2015).
Low-ionization structures are often present in PNe hosting binary central stars (Miszalski et al. 2009), and they are one of the candidates to explain the abundance discrepancy (Corradi et al. 2015).If low-ionization structures in NGC 2371 contained cold plasma, they would stand out in the 144 MHz radio image.Instead, they are fainter.This could be explained if they did not contain cold plasma, or contained significantly less cold plasma than other regions of the nebula.Observations of a larger sample when LoTSS is completed will allow for more deep and statistically important study of the cold plasma component in PNe.Multi-frequency analysis of PNe radio continuum images may allow us to better constrain the spatial distribution of cold plasma in PNe.
We will continue to study low-frequency radio emission of PNe using more complete data from the LoTSS survey.The number of observed PNe will increase rapidly when the survey improves completeness at low Galactic latitudes.The low-frequency survey LoLSS will observe at 42 − 66 MHz, but with reduced spatial resolution of 15 arcsec (de Gasperin et al. 2021) and sensitivity (1 mJy/beam) compared to LoTSS.Further advance will be made with the Square Kilometer Array (SKA) (Umana et al. 2015), which will carry out an extremely sensitive radio continuum survey at 1.4 GHz.Multi-frequency images will allow us to obtain accurate spectral index maps of PNe and possibly model the spatial distribution of cold plasma in PNe.

Figure 1 .
Figure 1.Top: Comparison of a computed radio spectrum of a PN with uniform electron temperature of 10 kK (blue dashed line) and a PN with an inclusion of cold plasma with Te = 2 kK randomly distributed filling 10% of the radius (red solid line).Bottom: the same as above, with the electron temperature of the uniform model is 5 kK.

Figure 2 .Figure 3 .
Figure 2. Left: Images of PNe at 144 MHz.We converted the flux intensity scale to brightness temperature, represented by the colorbar.The contours levels are spaced by 3 σ.White circle marks the size of the beam.Right: observed and fitted radio spectra of PNe.Lower panel shows the difference between the fit and the observed fluxes.

Figure 6 .
Figure 6.Radio contours and optical HST images convolved with a 6 arcsec gaussian in the F656N filter.The contours are separated by 3σ.North is at the top, east is to the left.The size of the restoring beam of 6 arcsec is marked in the bottom left corner.

Figure 7 .
Figure 7. Radio contours and optical IPHAS images in the Hα filter.The contours are separated by 3σ for most of PNe.For K 3-17 the contours are plotted at the level of 2, 3, 4.5, and 6 σ which better show the weak radio emission from the bipolar lobes.PM 1-305 is not centered in the image and is affected by two background stars in the Hα image.The size of the restoring beam of 6 arcsec for high-resolution images or 20 arcsec for low-resolution images is shown in the bottom left corner.North is at the top, east is to the left.

Figure 8 .
Figure 8. Radio contours and optical INT images in the Hα filter.The contours are separated by 3σ.The size of the restoring beam of 6 arcsec is shown in the bottom left corner.

Figure 9 .
Figure 9.Comparison of the 144 MHz radio continuum (contours) and Herschel 70 µm image (background) of NGC 6720.The Herschel image is taken from van Hoof et al. (2010).The size of the box is 100 arcsec height and 120 arcsec width.The 144 MHz beam is marked with a circle.

Figure 10 .
Figure 10.Comparison of the 144 MHz radio continuum (contours) and Hα HST image (background) of NGC 6371.The regions mark the extraction boxes used in Gómez-González et al. (2020).

Table 1 .
Frew et al. (2016), deconvolved diameters Θ d , corrected diameters Θ, and optical diameters taken fromFrew et al. (2016)of PNe detected in the LoTSS survey.The deconvolved and corrected diameters are not given for unresolved PNe.Large PNe were not fitted with Gaussian and their diameters Θ d refer to the size which exceeded 3σ (see text).

Table 2 .
Comparison of mean electron temperatures derived from radio observations optically thick PNe at 144 MHz, Balmer jump, helium lines, [O iii] and [N ii] CELs, and RL of O ii in Kelvin.He i λ7281/λ5876 c [O iii] d,e [N ii] d,e O ii f This work b Zhang et al. (2004) c Zhang et al. (2005) d Kaler (1986) e Kaler et al. (1996) f McNabb et al. (2013) a Figure 5. Images and spectra of PNe -continued.

Table 3 .
Upper limits for nebular flux densities at 144 MHz and corresponding brightness temperatures.