© The Authors 2024

1. Introduction

Understanding how and when the Universe underwent reionisation requires knowledge of the nature of the sources that provide the ionising UV photons. Although quasars are very bright in the UV, their number density at z ≥ 4 is likely too low to contribute significantly to the reionising budget of the Universe (see discussions in Kulkarni et al. 2019; Giallongo et al. 2015; Grazian et al. 2020). However, at high redshift, low-mass star-forming galaxies dominate the luminosity (mass) functions (e.g., Ouchi et al. 2009; Bouwens et al. 2015; Grazian et al. 2015; Song et al. 2016) which can be significant contributors to reionising the Universe (Robertson et al. 2010; Ishigaki et al. 2018; Finkelstein et al. 2019). Bright galaxies could also contribute to reionisation, as discussed by several authors in previous works (see e.g., Sharma et al. 2016; Naidu et al. 2020; Marques-Chaves et al. 2021, 2022a).

The escape of Lyman continuum (LyC) photons from galaxies before the Universe was completely reionised cannot be directly studied since the LyC photons will be absorbed by the intergalactic medium (IGM) (Inoue et al. 2014). Thus, we have to rely on local (low redshift) analogs of high-redshift galaxies to measure the escape of LyC photons where the IGM is mostly transparent to ionizing photons. Moreover, low redshift galaxies are easier to study in detail. Typically these local analogs (z ∼ 0.3) are extremely young starbursts with low dust attenuation, characterized by the slope of the non-ionizing UV continumm, β ≤ −2 given by f_λ ∝ λ^β (e.g., Chisholm et al. 2022). They also display a high star-formation rate density (Σ_SFR ≥ 0.1 M_⊙ yr⁻¹ kpc⁻²) and high ionization ratio ([OIII]/[OII]; O₃₂ ≥ 3), along with several other extreme emission line properties (e.g., Cardamone et al. 2009; Yang et al. 2017; Izotov et al. 2011, 2021a). These galaxies can efficiently produce ionising photons (Schaerer et al. 2016) and also leak a large amount of them (Izotov et al. 2016b,c,2018a,b, 2021b); thus, their counterparts at high-z can significantly contribute to reionisation. Recent James Webb Space Telescope (JWST) observations have confirmed that these are excellent local analogs of z > 6 galaxies (i.e., during the epoch of reionisation) since they show several similar extreme emission line properties (e.g., Schaerer et al. 2022; Rhoads et al. 2023; Matthee et al. 2023; Cameron et al. 2023; Mascia et al. 2023).

A systematic study of low redshift LyC emitters (LCEs, often also referred to as LyC leakers) was lacking until recently. The high success rate in finding LCEs in low-z analogs has motivated the Low-z Lyman Continuum Survey (Flury et al. 2022a, LzLCS, PI: Jaskot HST Project ID: 15626). This has enabled us to measure the escape fraction of LyC photons ($f_{\mathrm{esc}}^{\mathrm{LyC}}$) in a statistically large sample of galaxies. Here, $f_{\mathrm{esc}}^{\mathrm{LyC}}$ is defined as the ratio between the observed and the intrinsic ionizing fluxes and it is crucial parameter to constrain various models of reionization. LzLCS is a systematic survey of 66 star-forming galaxies at z ∼ 0.3 using Hubble Space Telescope (HST) Cosmic Origins Spectrograph (COS) to measure their LyC. They are selected from the Sloan Digital Sky Survey (SDSS) and Galaxy Evolution Explorer (GALEX), with a range of O₃₂, star-formation rate densities (Σ_SFR) and UV-slopes β. The LzLCS has found strong (> 2σ) LyC detection in 35 galaxies with a range of LyC escape fractions ($f_{\mathrm{esc}}^{\mathrm{LyC}}$ up to ∼50%), and has been combined with previous LyC measurements to provide a sample of 89 low-z galaxies with LyC measurements. Such a sample can then be used to study the dependence of various UV/optical-based physical properties on $f_{\mathrm{esc}}^{\mathrm{LyC}}$ (Flury et al. 2022b; Saldana-Lopez et al. 2022; Marques-Chaves et al. 2022b; Chisholm et al. 2022; Wang et al. 2021; Xu et al. 2023).

The exact physical mechanisms which lead to the leakage of LyC photons are still a matter of debate. LyC photons can either escape through low HI column density holes carved in the interstellar medium (ISM; Gazagnes et al. 2018; Saldana-Lopez et al. 2022) or from a density-bounded nebula (e.g., Nakajima & Ouchi 2014; Zackrisson et al. 2013; Bremer & Dayal 2023). Simulations show that strong supernova (SN) driven feedback can play an important role in the escape of LyC photons (Kimm & Cen 2014; Paardekooper et al. 2015; Kimm et al. 2017; Trebitsch et al. 2017; Ma et al. 2020; Barrow et al. 2020). Kakiichi & Gronke (2021) also show the effect of turbulent ISM on the escape of LyC photons. Recently Katz et al. (2023) showed that escape of LyC photons can occur from different modes of star formation. Ma et al. (2020) used zoom-in cosmological simulations to find that most LCEs are associated with a feedback-driven kiloparsec-scale superbubble and host very young stellar populations (≲10 Myr). Observationally, this can be seen, for example, in the Sunburst arc, where the regions that are leaking exhibit a very young stellar population and outflow signatures, as opposed to the non-leaking regions (Mainali et al. 2022; Kim et al. 2023). Simulations have also shown rapid fluctuations in $f_{\mathrm{esc}}^{\mathrm{LyC}}$ (on 10 Myr timescales) due to the stochastic and bursty nature of star formation (e.g., Ma et al. 2015; Trebitsch et al. 2017; Kimm et al. 2017; Barrow et al. 2020; Endsley et al. 2024). All these effects complicate the physical understanding and predictability of $f_{\mathrm{esc}}^{\mathrm{LyC}}$ using various diagnostics.

Radio observations provide a complimentary view into these galaxies. For example, the radio continuum (RC) is a dust-free tracer of the total star-formation rate (SFR) in galaxies (Yun & Carilli 2002; Murphy et al. 2011). For normal galaxies, the RC at 1.4 GHz and other mid-gigahertz (mid-GHz) frequencies (3, 6, and 10 GHz) has been extensively used to measure the SFR of galaxies both nearby and at high redshifts (e.g., Kennicutt & Evans 2012; Tabatabaei et al. 2017). Modelling the radio-spectral energy distribution (radio-SED), especially across several widely separated frequencies, enables us to study several unique physical parameters (Tabatabaei et al. 2017); for example, the supernova rate, the role of cosmic rays (CRs), emission measure based on the free-free absorption optical depth (Condon 1992; Hunt et al. 2004), and the equipartition magnetic field (Beck & Krause 2005; Beck 2015).

The RC emission at GHz frequencies in normal galaxies arises due to a combination of non-thermal synchrotron emission, which approximately follows a power-law with a slope of α_nth ∼ −0.8 and a thermal free-free component with a power law slope of −0.1 (Condon et al. 1991). In more complex scenarios the radio-SED shows a deviation from a simple power law behaviour. For example, in the case of dense starbursts in ultraluminous infrared galaxies (ULIRGS), blue-compact dwarfs (BCDs), and metal-poor dwarf galaxies, it has been observed that the radio-SED shows a turnover below ∼1 − 2 GHz (e.g., Condon et al. 1991; Hunt et al. 2004). Such a turnover in dense starbursts can arise from free-free absorption (FFA) due to large optical depths. Moreover, the RC spectrum can also show multiple FFA components (e.g., Clemens et al. 2010; Galvin et al. 2018). In addition at mid-to-high frequencies, they can show a much steeper spectrum than the canonical slope of −0.7, due to a steeper non-thermal spectral index (Galvin et al. 2018). Deviations in the non-thermal spectral index (α_nth) are also observed in local BCDs, showing both steeper and flatter α_nth than normal galaxies, possibly due to the inhomogeneous nature of star-formation in these galaxies (e.g., Klein et al. 1984, 1991; Deeg et al. 1993). Steep α_nth values are associated with a general lack of CR confinement (Skillman & Klein 1988; Klein et al. 1991). The exact physical mechanisms for this steepening are not well known, but it can be due to synchrotron ageing particularly in a strong magnetic field, inverse-Compton (IC) losses or even the escape of CRs (see Klein et al. 1991; Ramya et al. 2011, for a discussion). Such deviations have systematic effects on the radio-SED and can be used to distinguish different physical processes (see e.g., Klein et al. 2018; Clemens et al. 2010; Galvin et al. 2018). Overall the radio spectral index can also be used as a proxy for the age of starburst (Cannon & Skillman 2004; Hirashita & Hunt 2006). In young, metal-poor, and dense starbursts (e.g., SBS 0335−052; Hunt et al. 2004), the non-thermal emission is associated with an ensemble of compact SN remnants (SNRs) that are expanding in a dense ISM (e.g., Hunt et al. 2004; Hirashita & Hunt 2006) instead of diffuse CRs in normal galaxies.

Of particular interest is the non-thermal luminosity which (along with the non-thermal spectral index) can be used to estimate the Type II SN rate (Condon 1992). This can help us in understanding the role of SN feedback and, thus, also indirectly the contribution of SN to the escape of ionizing photons. The radio spectral index and the deviations between radio and other UV/optical/IR tracers are also used to understand changes in the recent star-formation history (e.g., Cannon & Skillman 2004; Bressan et al. 2002; Hirashita & Hunt 2006; Arango-Toro et al. 2023). There is no systematic RC study of LCEs to date. Past RC studies of local LCEs on Haro 11 (Hayes et al. 2007) and Tol 1247−232 (Puschnig et al. 2017) found a lower radio-based SFR compared to the SFR from other tracers. Furthermore, RC studies of sources most closely resembling LCEs, for instance green peas (GPs) and blueberries (BBs) also have suppressed radio-based SFR compared to other SFR indicators (Chakraborti et al. 2012; Sebastian & Bait 2019). It is speculated that such a suppression could be because galaxies with young stellar populations do not have any SN activity yet; thus, they would have a lower non-thermal emission component than normally expected (Hayes et al. 2007; Chakraborti et al. 2012; Sebastian & Bait 2019). Moreover, the loss of CRs either via diffusion or outflows might also be playing a role in causing the deficit (Sebastian & Bait 2019). Such a deficit has also been observed in local analogs of Lyman-break galaxies (LBGs) and for z ∼ 3 LBGs based on a stacking analysis (Carilli et al. 2008; Greis et al. 2017). The exact causes for this deficit are not yet completely investigated due to their unknown radio-SED. Several of the factors mentioned above for powerful starbursts are also valid for GPs and BBs which can lead to these variations.

In this study, for the very first time, we present RC observations of low-z star-forming galaxies with LyC measurements, which allows us to examine if and how their RC emission relates to the escape of ionizing photons. For instance, the radio thermal fraction and spectral index depend on the past SN activity and star-formation history (SFH) of galaxies. Importantly, the radio spectral index can be used as a proxy for the age of the starburst (Cannon & Skillman 2004; Hirashita & Hunt 2006). The dependence of these radio properties on $f_{\mathrm{esc}}^{\mathrm{LyC}}$ can tell us, for example, about the role of SN-driven feedback and SFH of galaxies on the escape of LyC photons. To do so, we target z ∼ 0.3 − 0.4 galaxies from the LzLCS sample. Our observations were conducted using the NRAO¹ Karl G. Jansky Very Large Array (VLA) for a sample of 53 LzLCS sources at 6 (C-band), 3 (S-band), and (a subsample) at 1.5 (L-band) GHz.

Our paper is organized as follows. In Section 2 we describe our observation setup, data analysis and detection, and non-detection statistics. In Section 3, we present the radio spectral index (across 3−6 GHz; ${\alpha }_{6\mathrm{GHz}}^{3\mathrm{GHz}}$) and thermal fraction distribution, and radio-SED for our sample of detections. Our main result is described in Section 3.5 where we show a dependence between ${\alpha }_{6\mathrm{GHz}}^{3\mathrm{GHz}}$ and $f_{\mathrm{esc}}^{\mathrm{LyC}}$. In Section 4, we discuss the various physical reasons that can give rise to such a relation. In Section 5, we summarize our main results and conclude. Throughout this paper, we assume H₀ = 70 km s⁻¹ Mpc⁻¹, Ω_m = 0.3, and Ω_Λ = 0.7, which is same as that used in Flury et al. (2022a).

2. Data

In this section, we describe our sample selection and VLA observations at multiple observing bands, namely: the C (4−8 GHz), S (2−4 GHz), and L (1−2 GHz) with the B-configuration (PI: Borthakur; Project Code: 21B-111). We also describe our calibration and imaging strategy.

2.1. Sample selection

Our current sample consists of 53 LzLCS sources observed with the VLA. This sample is a subset of the complete 66 galaxies from LzLCS. The LzLCS sample spans a wide range in $f_{\mathrm{esc}}^{\mathrm{LyC}}$ and can hence be used to do a statistical study of LCEs and its dependence on other physical properties (Flury et al. 2022b; Saldana-Lopez et al. 2022; Marques-Chaves et al. 2022b; Xu et al. 2023). Our 53 LzLCS sources with VLA observations were selected on the basis of their expected flux densities at L-band being above 40 μJy at 5-σ. This corresponds to a SFR of 6.6 M_⊙ yr⁻¹ for the median redshift (z ∼ 0.3) of our sample based on the standard radio–SFR relation (Murphy et al. 2011). The various physical properties of the LzLCS sources studied here are taken from Flury et al. (2022a) and Saldana-Lopez et al. (2022).

2.2. Observation setup

We used the default VLA setup for our C-, S-, and L- band continuum observations. For the C-band (4−8 GHz), we used a three-bit sampler with a total bandwidth of 4 GHz. For the S-band (2−4 GHz) and L-band (1−2 GHz), we used the eight-bit sampler with a total bandwidth of 2 and 1 GHz respectively. For the C-band setup, we split the entire bandwidth into two basebands each having 16 × 128 MHz subbands. Each subband was further split into 64 channels. Similarly for the S-band setup, we split the entire bandwidth into two basebands each having 8 × 128 MHz subbands where each subband was further split into 64 channels. And for the L-band setup, we split the entire bandwidth into two basebands each having 8 × 64 MHz sub-bands where each sub-band was further split into 64 channels. All our observations were conducted in the B-configuration with a few in the BtoBnA, BnA, and BnAtoA move time, depending on the scheduling constraints. We were thus able to typically achieve an angular resolution of 1.6″, 3.2″ and 6.9″ for C-, S-, and L-bands, respectively.

We followed the standard observing strategy used for continuum observations. We typically spent a minute on the initial setup and then 3−5 minutes on the flux density scale calibrator. We then observed the phase calibrator for a minute and the target for 5−6 minutes. We then cycled between the phase calibrator and the target. We had a total on-source time of 22 minutes for the C- and S-band observations, as well as a total on-source time of 75 minutes for the L-band observations. The C- and S-band observations of a given target source were combined in a single ∼1 h long scheduling block (SB), while their L-band observations were in separate SB with a length of ∼1 h and 40 min. Tables A.1 and A.2 summarizes the S-, C- and L-band observing details of our sample. In total, we were able to target 53 sources in both the S- and C-bands, and a subset of 17 in the L-band.

2.3. Data analysis: Calibration and imaging

We used the Common Astronomy Software Applications (CASA) data processing software (McMullin et al. 2007; CASA Team 2022) for all our data analysis. Each SB containing the multi-band C- and S-band data was initially flagged and calibrated using the VLA calibration pipeline v6.2.1. The calibrated flux density scale and phase calibrators for each SB were manually inspected for any leftover RFI and poorly performing antennas. These were then flagged and the datasets were re-run using the latest CASA pipeline v6.4.1. For some datasets, several rounds of the pipeline were run, after having performed flagging in each round.

The final calibrated target data were then extracted and any leftover RFI in some spectral windows was flagged using the task RFLAG. We then imaged these data using CASA task TCLEAN, separately by splitting the target in C- and S-bands. We use the multi-term multi-frequency synthesis (mtmfs) algorithm (Rau & Cornwell 2011) with NTERMS = 2, and a ROBUST = 0.5. We then manually made elliptical regions around each of our target sources using CASA task IMVIEW. We then estimated the flux densities of these sources using 2D Gaussian fit using the task IMFIT. We used this same region to estimate the noise in the images from a emission free area. For the one resolved source in our sample (J095700+235709), we estimated the fluxes in different bands and the corresponding errors using PyBDSF (Mohan & Rafferty 2015). The estimated flux densities and noise (σ) for our sources in C-, S-, and L-bands are presented in Table 1. For sources with a signal-to-noise ratio (S/N) below 3, we quote the 3σ upper limits in the flux density. In Table 2 we present a summary of the number of detections and non-detections in different bands. In all further analysis we add a 5% calibration error (e.g., Shao et al. 2022) in quadrature to the noise presented in Table 1.

2.4. VLA images

In Fig. 1, upper panel (lower panel), we overlay the S- and C-band radio contours on SDSS r-band images in the background for J124423+021540 (J095700+235709). In Appendix B, we show the C- and S-band VLA images for the rest of the detection sample. We note that except for J095700+235709, shown in Fig. 1, all the sources in our sample are unresolved in both radio and optical wavelengths.

Fig. 1. VLA images of two sources in our sample, J124423+021540 (top panels) and J095700+235709 (bottom panels), shown in white contours with the SDSS r-band image (Blanton et al. 2017) in the background. On the left (right) panels we show the VLA S-band (C-band) image contours. The lowest contours in each panels corresponds to 3σ level and they rise in powers of $\sqrt{2}$. σ values correspond to ∼8 and ∼5 μJy beam−1 for S- and C-band respectively.

3. Results

3.1. Radio-SED at mid GHz frequencies

The RC in the simplest case is modeled by a single power law given by²: S​_ν = Aν^α, where α is the power-law index, also termed the total spectral index, ν is the frequency, and A is a proportionality constant. At frequencies below 30 GHz, the radio-SED for normal star-forming galaxies can be modeled as a combination of (i) free-free thermal emission from star-forming regions having a relatively flat spectral index of −0.1 and (ii) non-thermal emission having a steeper spectral index (α_nth) arising due to synchrotron emission from relativistic electrons (see e.g., Condon 1992). Thus, the RC SED in such a case can be written with respect to a reference frequency ν_° following Tabatabaei et al. (2017), as:

$$\begin{array}{c}\hfill S_{\nu }=S_{\nu }^{\mathrm{th}}+S_{\nu }^{\mathrm{nth}}=A_1{\left(\frac{\nu }{{\nu }_{°}}\right)}^{-0.1}+A_2{\left(\frac{\nu }{{\nu }_{°}}\right)}^{{\alpha }_{\mathrm{nth}}}.\end{array}$$(1)

We can then define the thermal fraction at the reference frequency $f_{\text{th}}=S{\text{}}_{{\nu }_{°}}^{\text{th}}/S{\text{}}_{{\nu }_{°}}$.

Under the assumption that the thermal emission is optically thin and no escape of LyC photons, $S_{\nu }^{\mathrm{th}}$ is directly related to the intrinsic ionizing photon rate and, thus, to the SFR and the electron temperature (see e.g., Eqs. (10) and (11) in Murphy et al. 2011). For normal star-forming galaxies, the non-thermal emission is dominated by CRs which were accelerated by supernova and diffuse in the ISM. Thus, the non-thermal emission is directly related to the supernova rate and α_nth (see Eq. (13) in Murphy et al. 2011). The relative contribution of the non-thermal component in normal star-forming galaxies increases as we go to lower frequencies (30 GHz and below). We show the radio-SED, for a normal star-forming galaxy with f_th = 0.1 (10%) at 1 GHz and α_nth = −0.8 in Fig. 2 as a dashed blue line. The dotted blue line shows the corresponding thermal emission component. The combined effect of the steep non-thermal and flat thermal component leads to an apparent spectral index between C- and S-bands (${\alpha }_{6\mathrm{GHz}}^{3\mathrm{GHz}}$) equal to −0.64. In the following, we refer to such a radio-SED and its apparent spectral index between C- and S-bands as standard and/or canonical, representing normal star-forming galaxies. We use this reference value of ${\alpha }_{6\mathrm{GHz}}^{3\mathrm{GHz}}$ to compare with the spectral index of LzLCS sources.

Fig. 2. Typical radio-SED with a mixture of thermal and non-thermal emission (Eq. (1)) in blue line. The corresponding thermal and non-thermal emission components are shown in blue dashed and dotted lines. We also show the effect of FFA (Eq. (2)) with a turnover at 1 GHz in orange. And the cutoff and curved model (from Klein et al. 2018) with a high-frequency break at 6 GHz in green and red respectively. The grey-shaded area shows the region of the SED explored in this work.

For young galaxies (ages < 10 Myr), it is possible that the f_th is higher than 0.1 due to a lack of enough supernovae and/or flat CR spectrum (≥0.5) and thus closer to that found in SNRs (e.g., Hunt et al. 2004, 2005). Moreover, in dense starbursts as discussed before, FFA has an effect on both the thermal and non-thermal emission and leads to a turnover in the radio spectra below a characteristic frequency (Condon 1992). This is termed the turnover frequency (ν_t) and is the frequency at which the free-free optical depth (τ_ν) reaches unity. The FFA τ_ν is proportional to the emission measure (EM), ${\tau }_{\nu }\propto n_{\text{e}}^{2}L$ (pc cm⁻⁶), where n_e is the electron number density and L is the path length along the line of sight. Thus ν_t is also directly related to EM and increases with increasing EM. Following Clemens et al. (2010) and Galvin et al. (2018), we can parameterize τ_ν using ν_t, such that τ_ν = (ν/ν_t)^−2.1. We can thus modify the simple radio-SED from Eq. (1) under FFA in the simplest case as:

$$\begin{array}{c}\hfill S_{\nu }=(1-e^{-{\tau }_{\nu }})[A_1+A_2{\left(\frac{\nu }{{\nu }_{°}}\right)}^{0.1+{\alpha }_{\mathrm{nth}}}]{\left(\frac{\nu }{{\nu }_{°}}\right)}^2.\end{array}$$(2)

In Fig. 2, we show the effect of FFA with ν_t at 1 GHz as orange line. In more complex scenarios, there can be multiple FFA components leading to even more complex features in the radio-SED (e.g., Clemens et al. 2010; Galvin et al. 2018).

Loss of CRs particularly due to synchrotron and inverse-Compton (IC) mechanisms can lead to curvatures in the spectrum towards high frequencies (Klein et al. 2018). We show such a curved spectrum in Fig. 2 with the red line. Finally, in the scenario of a single injection of CRs there can be a sharp cutoff at higher frequencies in the CR energy spectrum leading to a steepening at frequencies above a characteristic break frequency (cutoff model in Klein et al. 2018). We show such a spectrum with a cuttoff with the green line in Fig. 2. We refer to Klein et al. (2018) for a detailed discussion on various modifications to the standard radio-SED. The effect of steepening of the spectrum due to these mechanisms will show an effect on the spectral index in the frequency range explored in this work between (2−8 GHz), although we would need high-frequency (up to 30 GHz) observations to distinguish between different scenarios.

Figure 3 shows the observed RC SED for three LzLCS sources. In all three panels, the blue dots show the observed flux densities, the gray line corresponds to a simple radio-SED Eq. (1) normalized to the 3 GHz observed flux density with f_th = 0.1 at 1 GHz and α_nth = −0.8 which is typical of a star-forming galaxy. In the left panel is the radio-SED for J144010+461937 which has a radio-SED close to a star-forming galaxy. In Appendix C we can see that several sources in our sample show similar SEDs.

Fig. 3. Observed radio continuum SEDs at three frequencies: 1.5, 3, and 6 GHz for three sources in our sample. The error bars represent 1σ errors on the flux measurement. The gray line in each panel is an overlay (not a fit) for a simple radio-SED (Eq. (1)) normalized to the 3 GHz observed flux density with fth = 0.1 at 1 GHz and αnth = −0.8 commonly seen in star-forming galaxies. The left panel shows an example of a radio-SED which is close to the standard radio-SED (Eq. (1)). The middle panel shows an example that has a steep spectrum above 3 GHz and is close to the standard radio-SED below 3 GHz. The right panel shows the case of a broken SED.

The middle panel shows the radio-SED for J095700+235709 where we see a steepening above 3 GHz. Several sources shown in Appendix C are very steep at GHz frequencies. Moreover, the radio-SED of J082652+182052 and J091113+183108 shown in Appendix C show evidence for FFA based on the 3-σ upper limits on the L-band flux densities. The right panel of Fig. 3 exhibits the radio-SED of J083440+480541 which is an example of a broken SED above and below 3 GHz. J081409+211459 also shows a similar feature (Appendix C). Such a SED can arise due to a mix of multiple star-forming regions, one with an FFA component at high (∼3 GHz), and the other with a steep spectrum. Multiple ULIRGs in the Clemens et al. (2010) sample show such a radio-SED. Interestingly, we also find that these two sources are outliers when compared to their LyC properties which we discuss in Sect. 4. Appendix C shows the radio-SED for our entire sample of C- and S-band detections. Modelling the radio-SED to derive the relevant physical parameters has several free parameters and requires a flux density measurement across a large frequency range in order to remove the various degeneracies in parameter estimation (see e.g., Galvin et al. 2018). Future observations of the radio-SED at lower frequencies (< 1 GHz, e.g., with the Giant Meterewave Radio Telescope and Low-Frequency Array) should help in constraining some of these physical mechanisms.

3.2. Total spectral index distribution

We calculated the total spectral indices using the slope of the observed flux densities at 6 and 3 GHz (${\alpha }_{6\mathrm{GHz}}^{3\mathrm{GHz}}$) for our sources. We propagated the 1-σ flux uncertainties on the flux measurements to estimate the uncertainties on the spectral indices. For sources with non-detections in only one of the S- and C-band, we estimated the corresponding lower/upper limit on the spectral indices using the 3σ upper limit on the flux density. In Table 1 we present the ${\alpha }_{6\mathrm{GHz}}^{3\mathrm{GHz}}$ values for our sample. From our sample of 24 C- and S-band detections, we could not derive reliable ${\alpha }_{6\mathrm{GHz}}^{3\mathrm{GHz}}$ for three sources due to unreliable flux density measurements. Here, the fluxes are unreliable due to the presence of strong sidelobes arising from very strong sources in the field. We use the rest of the sample of 21 sources with reliable ${\alpha }_{6\mathrm{GHz}}^{3\mathrm{GHz}}$ for all further analysis involving spectra indices. Overall in this sample 13 sources show steep (${\alpha }_{6\mathrm{GHz}}^{3\mathrm{GHz}}<-0.64$) and 8 show flat radio spectra (${\alpha }_{6\mathrm{GHz}}^{3\mathrm{GHz}}\ge -0.64$) between 2−8 GHz. For 7 of these sources we have L-band observations. Of them, 3 sources show a turnover at frequencies below 3 GHz, which is evidence for a FFA component between 1−3 GHz. And 2 sources show a broken radio-SED where the spectral indices between 1−3 GHz (L- and S-bands) is steeper than between 3−8 GHz (C- and S-bands). The remaining one source shows a typical radio-SED of a normal star-forming galaxy.

In Fig. 4, we show the distribution of ${\alpha }_{6\mathrm{GHz}}^{3\mathrm{GHz}}$ in the blue histogram. The mean ${\alpha }_{6\mathrm{GHz}}^{3\mathrm{GHz}}$ for our sample is −0.92 with a scatter of 0.71. The total spectral index also varies in a very broad range from −2.18 to 0.22. On the contrary, the mean total spectral index of normal star-forming galaxies from the KINGFISHER galaxies was found to be −0.79 with a scatter of 0.15 (Tabatabaei et al. 2017), albiet in a wider frequency range of 1−10 GHz. Thus LzLCS galaxies overall have a steeper mean radio spectra with a larger scatter than normal star-forming galaxies.

Fig. 4. Distribution of the total spectral index of LzLCS sources which spans a wide range from flat to steep. The dashed line shows the value of ${\alpha }_{6\mathrm{GHz}}^{3\mathrm{GHz}}$ expected in normal star-forming galaxies with a fth = 0.1 at 1 GHz and αnth = −0.8, ${\alpha }_{6\mathrm{GHz}}^{3\mathrm{GHz}}=-0.64$.

3.3. The RC–SFR relation for LzLCS sources

The RC of normal star-forming galaxies at GHz frequencies can be used to determine their SFR (Kennicutt & Evans 2012; Murphy et al. 2011; Tabatabaei et al. 2017). In this section, we compare the observed RC of LzLCS galaxies to those predicted by such a calibration. The standard method used to calibrate the RC–SFR relation relies on the tight radio (at 1.4 GHz)-Far-infrared (FIR) correlation (de Jong et al. 1985; Helou et al. 1985). However, since our observations are made at 3−6 GHz we need to scale to the 1.4 GHz luminosity assuming a spectral index. Since LzLCS sources show a large variety in their radio-SED, such an assumption may not be valid across such a large frequency range.

An alternate approach as discussed in Murphy et al. (2011) relates the thermal RC (RC_th) luminosity at a given frequency, ν, to the SFR, as follows:

$$\begin{array}{c}\hfill {\mathrm{SFR}}_{\nu }^{\mathrm{th}}=4.6\times {10}^{-28}{\left(\frac{T_{\mathrm{e}}}{{10}^4\mathrm{K}}\right)}^{-0.45}{\left(\frac{\nu }{\mathrm{GHz}}\right)}^{0.1}\left(\frac{L_{\nu }^{\mathrm{th}}}{\mathrm{ergs}/\mathrm{s}/\mathrm{Hz}}\right),\end{array}$$(3)

where T_e is the electron temperature. This is because both the SFR and RC_th emission are proportional to the ionizing photon rate in an optically thin plasma. Moreover, the non-thermal RC (RC_nth) luminosity at a given frequency, ν, can also be related to the SFR. This is because the SFR in a galaxy is related to the SN rate and so is the synchrotron emission at GHz frequencies. Therefore owing to a calibration between the SN rate and non-thermal luminosity (Condon & Yin 1990) we can derive a RC_nth–SFR calibration (see Eq. (15) in Murphy et al. 2011). Here, it is assumed that the calibration between SN rate and non-thermal relation is the same as that observed in our Galaxy. These two relations can then be combined to relate the total RC luminosity and the SFR as follows:

$$\begin{array}{cc}\hfill {\mathrm{SFR}}_{\nu }=& {10}^{-27}[2.18{\left(\frac{T_{\mathrm{e}}}{{10}^4\mathrm{K}}\right)}^{0.45}{\left(\frac{\nu }{\mathrm{GHz}}\right)}^{-0.1}\hfill \\ \hfill & {+15.1{\left(\frac{\nu }{\mathrm{GHz}}\right)}^{{\alpha }_{\mathrm{nth}}}]}^{-1}\left(\frac{L_{\nu }}{\mathrm{ergs}/\mathrm{s}/\mathrm{Hz}}\right).\hfill \end{array}$$(4)

Using this relation is more convenient than the radio-FIR correlation to compare galaxies with a non-standard α_nth and thermal fraction.

In Fig. 5, we show the RC luminosity at 3 and 6 GHz for LzLCS sources against the extinction corrected UV-based SFR in the left panel and the Hβ-based SFR in the right panel. For converting the observed S-band (and C-band) flux density of LzLCS sources to the RC luminosity at 3 GHz (and 6 GHz) for the sample with S-band (and C-band) detections we used the observed total ${\alpha }_{6\mathrm{GHz}}^{3\mathrm{GHz}}$, and for the non-detections we assumed a total ${\alpha }_{6\mathrm{GHz}}^{3\mathrm{GHz}}=-0.64$, as that of normal galaxies. As such there is a relatively weak dependence of the spectral index on the luminosities as the K-corrections are small since our sources are at a relatively low redshift of ∼0.3. The blue points in both panels are for the C- (and S-band) detection samples and the orange points are the 3σ upper limits for the sample with S-band (and C-band) non-detections.

Fig. 5. LzLCS sources on the RC–SFR relation at 3 GHz top panels and 6 GHz in the bottom panels. The left (right) panel shows the UV-based (Hβ-based) SFRs for the LzLCS sources. In all panels, the blue points represent detections in the corresponding band and the orange points represent the 3σ upper limits for non-detections. The solid black line is plotted for the standard RC–SFR relation with αnth = −0.8. The dotted, dashed and dotted dashed lines are for different values of αnth from flat (−0.3) to very steep (−2), respectively. The red solid line represents the RCth–SFR relation from Eq. (3). In all these cases, we have used a fixed value of Te = 10 000 K.

In Fig. 5, the solid black lines show the total RC–SFR calibration from Eq. (4) in all the panels. Here, we assumed a standard α_nth = −0.8 and T_e = 10 000 K. We also show how this relation changes purely if the α_nth steepens to −1.15 and −2. Steep values of α_nth ∼ −2.1 are observed in nearby BCDs (Klein et al. 1991), and this is the reason we chose this range. We note that the relation has the same slope, but it systematically moves down as the α_nth steepens. This happens because, for a fixed frequency, a steeper non-thermal spectrum has a lower contribution from the non-thermal emission. On the contrary, a flatter α_nth moves the relation upwards. In the case where there is no non-thermal emission namely, purely thermal radio spectra, we use the RC_th–SFR relation from Eq. (3). The red solid lines in both panels show this limit. In principle, galaxies can also be below this line if the thermal emission is reduced due to FFA around 3 GHz. Such FFA at GHz frequencies is observed in compact starbursts with high EM (e.g., Hunt et al. 2004). Most of the points in Fig. 5 are above the red line due to the contribution of non-thermal emission in LzLCS sources. We compare the RC luminosity with the SFR derived from two different tracers, UV-based and Hβ-based. Thermal radio emission is due to the Coulomb interaction of free thermal electrons in an ionized gas and Hβ emission arises due to the recombination of electrons in the same region. And the UV emission is simple starlight arising from massive OB stars. The timescale of star-formation traced by the Hβ line is of the order of 10 Myr and is much shorter than the star-formation timescale traced by the UV emission which is ∼100 Myr (Kennicutt & Evans 2012).

Overall, we observe that although LzLCS sources show a correlation between their RC luminosity and SFR, there is a lot of scatter from the standard RC–SFR relation, with several sources much above and below this relation. It is important to note that several points are above the RC_th–SFR relation, thus suggesting that LzLCS sources have non-thermal emission which suggests past SN activity.

We next investigate if the offset between the observed and expected RC luminosity can be explained by variations in the observed ${\alpha }_{6\mathrm{GHz}}^{3\mathrm{GHz}}$ of these sources. In Fig. 6, we show the ratio between the observed RC luminosity and that expected from the RC–SFR relation based on Eq. (4) (termed as the “SFR offset” here on) with respect to the total ${\alpha }_{6\mathrm{GHz}}^{3\mathrm{GHz}}$ for only the C- and S-band detection sample. To estimate the RC luminosities here we have used the observed spectral index values as done previously. We study the RC-offset for both the S- and C-bands and for SFR measured using UV and Hβ emission. We notice that the Hβ-based SFR offset systematically increases for steep spectrum sources (top right panel) for the S-band, whereas the flat spectrum sources are closer to that expected. The SFR offsets with UV-based SFRs at C-band are higher for flat spectrum sources (bottom left panel). Such a deviation can lead to values above the expected RC–SFR relation in Fig. 5. The SFR offsets with UV-based SFRs for S-band (top left panel) and Hβ-based SFRs for C-band (bottom right panel) do not show any such offsets although they have a large scatter. We currently do not completely understand the origin of such systematic and large offsets. Overall this suggests that accurate SFR estimation using different tracers is difficult for LzLCS sources and particularly so with radio-based SFR tracers.

Fig. 6. Ratio between the observed RC luminosity and that expected from the RC–SFR relation from a standard radio-SED with αnth = −0.8, fth = 0.1 at 1 GHz, and Te = 104 K, at 3 GHz and 6 GHz (top and bottom panels respectively) vs. the spectral index ${\alpha }_{6\mathrm{GHz}}^{3\mathrm{GHz}}$. The left (right) panel shows the ratio with the expected RC for the UV-based (Hβ-based) SFR. The dashed line shows a ratio of 1 for reference. Here, we have used only the S- and C-band detection sample.

Finally, an important reason for the scatter in the RC–SFR relation can be that for galaxies underdoing starbursts like the LzLCS sources, the SFR traced by different methods are in general not consistent due to the stochastic/bursty nature of star formation (e.g., Izotov et al. 2016a; Atek et al. 2022). This inconsistency is also evident between the UV-based and Hβ-based SFRs. Finally, for the leakers (weak and strong) the Hβ-based SFR can underestimate the true SFR, since not all the ionizing photons are absorbed by the star-forming regions. Most of the leakers in our sample are weak (< 10%) and the two strong leakers have $f_{\mathrm{esc}}^{\mathrm{LyC}}$ below 30%. Thus this effect is marginal compared to the large SFR offsets we observe. We also see that the RC is significantly below the RC–SFR relation in several sources. This can happen due to a combination of a steep α_nth or a relatively low non-thermal fraction. In a few sources, their RC luminosity is fairly close to that expected from the RC_th–SFR relation, and are thus possibly thermally dominated systems. Such a deficit has also been observed in a sample of green-peas (Chakraborti et al. 2012), blueberries (Sebastian & Bait 2019) and in HII galaxies (Rosa-González et al. 2007). The effect of FFA on both thermal and non-thermal emissions can also be a further cause of such a deficit. A few sources are above the black solid line, where it is possible that their α_nth is flatter than the canonical value. It is also possible that the calibration between the SN rate-SFR observed in our Galaxy is not valid for LzLCS sources.

3.4. Comparison between radio-based SFR with other SFR tracers

We compare the radio-based SFR derived at 6 GHz (SFR${}_{6\mathrm{GHz}}^{\mathrm{RC}}$) with the SFR_UV and SFR_Hβ for our LzLCS sample in Fig. 7. We use Eq. (4) to derive the SFR^RC where we assumed a α_nth of −0.8. We find that the SFR^RC shows a systematic offset from the SFR${}_{6\mathrm{GHz}}^{\mathrm{RC}}$–SFR_UV (SFR_Hβ) relation for local normal star-forming galaxies from Tabatabaei et al. (2017) similar to that observed in the bottom-left panel of Fig. 5. The median value of this offset from SFR_UV (SFR_Hβ) is −0.59 dex (−0.26 dex). The only difference between the SFR calibrations from Tabatabaei et al. (2017) is that it uses Hα based SFR as opposed to Hβ line used in our work. Tabatabaei et al. (2017) found a slightly higher value of α_nth of −0.97, in their KINGFISHER sample. Assuming this value to calculate the SFR^RC decreases the offsets by a small amount to −0.52 dex (−0.18 dex) for SFR_UV (SFR_Hβ).

Fig. 7. Comparison between SFR${}_{6\mathrm{GHz}}^{\mathrm{RC}}$ and SFRUV (circles) and SFRHβ (triangles). Each of the points is color coded by their apparent ${\alpha }_{6\mathrm{GHz}}^{3\mathrm{GHz}}$. We over plot the relation between SFR${}_{6\mathrm{GHz}}^{\mathrm{RC}}$ and other tracers from Tabatabaei et al. (2017) in blue (SFRUV calibration) and orange (SFRHβ calibration). The blue and orange shaded region shows the scatter in these relations.

We also color-code the points in Fig. 7 with ${\alpha }_{6\mathrm{GHz}}^{3\mathrm{GHz}}$ to study it dependence of the SFR offset. As discussed in Sect. 3.3 we observe that the ratio of the observed and expected RC luminosity shows a systematic offset with ${\alpha }_{6\mathrm{GHz}}^{3\mathrm{GHz}}$. We hence fit the logarithm of SFR_UV (SFR_Hβ) with a linear combination of log₁₀ SFR${}_{6\mathrm{GHz}}^{\mathrm{RC}}$ and ${\alpha }_{6\mathrm{GHz}}^{3\mathrm{GHz}}$ as follows:

$$\begin{array}{c}\hfill {log}_{10}{\mathrm{SFR}}_{\mathrm{UV}(\mathrm{H}\beta )}=a+b\times {log}_{10}{\mathrm{SFR}}_{\mathrm{RC}}^{\nu }+c\times {\alpha }_{6\mathrm{GHz}}^{3\mathrm{GHz}}.\end{array}$$(5)

The details on the fitting procedure can be found in Appendix E. We show this fit between log₁₀ SFR_UV and log₁₀ SFR${}_{6\mathrm{GHz}}^{\mathrm{RC}}$ and ${\alpha }_{6\mathrm{GHz}}^{3\mathrm{GHz}}$ in Fig. 8 where the residual is 0.21. This is comparable to the typical scatter found between different SFR tracers from Tabatabaei et al. (2017). We provide the fit parameters at different radio frequencies in Table 3. In the future, with observations spanning a wide frequency range from 1−10 GHz, we can use the Mid RC (MRC)-SFR calibration suggested by (Tabatabaei et al. 2017) to further improve the SFR^RC calibration.

Fig. 8. Combined fit between the linear combination of log10 SFR${}_{6\mathrm{GHz}}^{\mathrm{RC}}$ and ${\alpha }_{6\mathrm{GHz}}^{3\mathrm{GHz}}$, and log10 SFRUV. The blue line shows the x = y line. Note: the residual in the fit is 0.21, which is similar to the typical scatter between different SFR tracers from Tabatabaei et al. (2017).

3.5. Relations between LyC escape and RC emission

In this section, we for the first time study the dependence of the Lyman continuum escape fraction, $f_{\mathrm{esc}}^{\mathrm{LyC}}$, on the total spectral index and thermal fraction empirically. Figure 9 shows the dependence of $f_{\mathrm{esc}}^{\mathrm{LyC}}$ on the total spectral index measured in the two bands, C- and S. For this and all further analysis we use only the sources with reliable ${\alpha }_{6\mathrm{GHz}}^{3\mathrm{GHz}}$ for our analysis. The sources with LyC detections are shown in blue points and the non-detections are shown as upper limits (1-σ) in orange points.

Fig. 9. Dependence of $f_{\mathrm{esc}}^{\mathrm{LyC}}$ on the total spectral index measured across C- (6 GHz) and S- (3 GHz) bands. The left panel shows the dependence on $f_{\mathrm{esc}}^{\mathrm{LyC}}$ measured using the UV-fit method. And the right panel shows the dependence on $f_{\mathrm{esc}}^{\mathrm{LyC}}$ measured using the Hβ method Flury et al. (2022a). The blue points in both panels show LyC detections and the orange points show LyC non-detections (1σ upper limits). We also show two known nearby leakers from the literature, Haro 11 (brown star) and Tol 1247−232 (black star). In both panels the horizontal lines demarcate the regions for strong leakers ($f_{\mathrm{esc}}^{\mathrm{LyC}}>0.1$), weak leakers ($0.01<f_{\mathrm{esc}}^{\mathrm{LyC}}<0.1$), and non-leakers ($f_{\mathrm{esc}}^{\mathrm{LyC}}<0.01$). The vertical line shows the region above and below a spectral index of a standard radio-SED (${\alpha }_{6\mathrm{GHz}}^{3\mathrm{GHz}}=-0.64$) demarcating the region for flat and steep spectrum sources. See text for details.

The left and right panels shows the radio spectral index dependence of $f_{\mathrm{esc}}^{\mathrm{LyC}}$ measured using two different approaches, named the “UV-fit” and “Hβ” methods (see Flury et al. 2022a, for details). In both methods, the escape fraction is calculated as the ratio between the “observed” (escaping) and the “intrinsic” (produced) ionizing fluxes, both measured at a narrow wavelength interval just below the Lyman edge (912 Å). Both former approaches differ according to the way the intrinsic ionizing flux is estimated. In the UV-fit method, the intrinsic LyC flux was predicted via detailed SED-fitting to non-ionizing part of the COS/G140L spectra of LzLCS galaxies (see Saldana-Lopez et al. 2022). In the Hβ method, the intrinsic LyC flux density was estimated from the measured equivalent-width of the Hβ line from SDSS spectra, as a tracer of the age of the most massive stellar populations (see also Izotov et al. 2016c). In instantaneous-burst models (Leitherer et al. 1999, used consistently for both former approaches), the age of the burst is inherently related to the production of LyC photons.

The horizontal dashed lines show the region for strong leakers ($f_{\mathrm{esc}}^{\mathrm{LyC}}>0.1$), weak leakers ($0.01<f_{\mathrm{esc}}^{\mathrm{LyC}}<0.1$), and non-leakers ($f_{\mathrm{esc}}^{\mathrm{LyC}}<0.01$). The vertical dashed line shows the value of ${\alpha }_{6\mathrm{GHz}}^{3\mathrm{GHz}}$ for a standard radio-SED from Eq. (1) with a α_nth = −0.8 and f_th = 0.1 at 1 GHz. This line splits the figure into regions with flat spectrum (${\alpha }_{6\mathrm{GHz}}^{3\mathrm{GHz}}>-0.64$) and steep spectrum sources (${\alpha }_{6\mathrm{GHz}}^{3\mathrm{GHz}}<-0.64$). We also show two known nearby weak leakers, Haro 11 and Tol 1247−232 on this plot. We used the literature value of $f_{\mathrm{esc}}^{\mathrm{LyC}}$ from Leitet et al. (2011) and Leitherer et al. (2016) for Haro 11 and Tol 1247−232 respectively. The radio spectral index for Haro 11 was measured between 3.0−8.46 GHz. We used Schmitt et al. (2006) for the 8.46 GHz flux density and the Very Large Array Sky Survey (VLASS) catalogue (Gordon et al. 2021) for the 3.0 GHz flux density. The radio spectral index for Tol 1247−232 was measured between 5−8 GHz using VLA observations from Rosa-González et al. (2007). The radio spectral index is measured in a very similar range frequency range as our C- and S-bands, enabling us to use them for comparison purposes.

It is interesting to find a positive correlation between $f_{\mathrm{esc}}^{\mathrm{LyC}}$, particularly measured using the UV-fit and ${\alpha }_{6\mathrm{GHz}}^{3\mathrm{GHz}}$. The correlation appears to be marginally weak for the $f_{\mathrm{esc}}^{\mathrm{LyC}}$ measured using Hβ, however, that seems to be largely driven by the difference in the upper limits on the LyC non-detections between the two methods. This indicates that strong leakers show overall flatter radio spectra, and non-leakers have steeper ${\alpha }_{6\mathrm{GHz}}^{3\mathrm{GHz}}$, while weak leakers typically have ${\alpha }_{6\mathrm{GHz}}^{3\mathrm{GHz}}$ close to the canonical value (${\alpha }_{6\mathrm{GHz}}^{3\mathrm{GHz}}=-0.64$). Interestingly, Haro 11 and Tol 1247−232 are also within our ${\alpha }_{6\mathrm{GHz}}^{3\mathrm{GHz}}-f_{\mathrm{esc}}^{\mathrm{LyC}}$ relation. We also find two outliers (J155945+403325 and J083440+480541) from the ${\alpha }_{6\mathrm{GHz}}^{3\mathrm{GHz}}-f_{\mathrm{esc}}^{\mathrm{LyC}}$ relation (in both the plots) and J081409+211459 which is clearly an outlier in the left-panel of Fig. 9 and possibly in the right panel. These sources show a flat spectrum but much lower values of $f_{\mathrm{esc}}^{\mathrm{LyC}}$ than other LzLCS sources. Interestingly, J083440+480541 and J081409+211459 have L-band observations and show a steeper spectrum at low frequencies. However, in other sources with L-band observations we see either a lower L-band flux density or that consistent with a standard radio-SED. Finally in Table 4 we quantify these correlations using the Kendal τ correlation coefficient. Here we take into account the upper limits on the $f_{\mathrm{esc}}^{\mathrm{LyC}}$ values following Akritas & Siebert (1996) which is the procedure following in Flury et al. (2022b). We find that ${\alpha }_{6\mathrm{GHz}}^{3\mathrm{GHz}}$ correlates with both $f_{\mathrm{esc}}^{\mathrm{LyC}}$ UV-fit and Hβ method (strong τ values) particularly, after removing the two outliers. However the significance of this correlation is less compared to other diagnostics discussed in Flury et al. (2022b). It is important to note that our current sample is quite small compared to the full sample studied in Flury et al. (2022b) and also lacks sources at high LyC escape. In the future, we aim to increase the sample size with radio observations of LyC leakers, and particularly with L-band data. This should give us a better picture of the dependance of $f_{\mathrm{esc}}^{\mathrm{LyC}}$ and ${\alpha }_{6\mathrm{GHz}}^{3\mathrm{GHz}}$ and on the properties of outliers.

3.6. Relation between β₁₅₅₀, O₃₂, and Σ_SFR and the RC emission

The escape of LyC photons from galaxies is known to be correlated with several other physical parameters from the literature (e.g., Izotov et al. 2021b; Flury et al. 2022b; Saldana-Lopez et al. 2022; Chisholm et al. 2022; Marques-Chaves et al. 2022b; Xu et al. 2023), particularly the UV slope β₁₅₅₀, the O₃₂ emission line ratio, and SFR surface density Σ_SFR. In particular, galaxies with high $f_{\mathrm{esc}}^{\mathrm{LyC}}$ are expected to have steeper values of β₁₅₅₀, higher O₃₂, and a higher Σ_SFR and vice versa. Thus the correlation between ${\alpha }_{6\mathrm{GHz}}^{3\mathrm{GHz}}$ and $f_{\mathrm{esc}}^{\mathrm{LyC}}$ observed here leads us to explore if the ${\alpha }_{6\mathrm{GHz}}^{3\mathrm{GHz}}$ of LzLCS sources also depends on these parameters with similar trends. Moreover, the LzLCS sample is also selected based on a range of β, O₃₂, and Σ_SFR. In Fig. 10, we show the dependence of β₁₅₅₀, O₃₂, and Σ_SFR on ${\alpha }_{6\mathrm{GHz}}^{3\mathrm{GHz}}$. We see that ${\alpha }_{6\mathrm{GHz}}^{3\mathrm{GHz}}$ and β₁₅₅₀ are generally anti-correlated, so the β₁₅₅₀ slope steepens as ${\alpha }_{6\mathrm{GHz}}^{3\mathrm{GHz}}$ flattens. The presence of flatter β₁₅₅₀ slopes (i.e., redder UV spectra) in sources with steep ${\alpha }_{6\mathrm{GHz}}^{3\mathrm{GHz}}$ can be associated with the declining phase of the starburst (see Sect. 4.3 for more discussion) and also to higher dust attenuation (e.g., Saldana-Lopez et al. 2022; Chisholm et al. 2022). We also notice that ${\alpha }_{6\mathrm{GHz}}^{3\mathrm{GHz}}$ overall flattens with increasing O₃₂ and Σ_SFR which is both in line with that expected with their dependence on $f_{\mathrm{esc}}^{\mathrm{LyC}}$. However, it is difficult to associate a direct causal relation between ${\alpha }_{6\mathrm{GHz}}^{3\mathrm{GHz}}$ and these parameters as $f_{\mathrm{esc}}^{\mathrm{LyC}}$ depends on multiple physical parameters (e.g., Flury et al. 2022b). We report the Kendall τ values between these quantities in Table 4. We find a strong correlation between these quantities but with a low significance (low p-values). We discuss the possible physical reasons driving a correlation between ${\alpha }_{6\mathrm{GHz}}^{3\mathrm{GHz}}$ and $f_{\mathrm{esc}}^{\mathrm{LyC}}$ in the next section.

Fig. 10. Dependence of β1550, O32 and ΣSFR on the spectral index ${\alpha }_{6\mathrm{GHz}}^{3\mathrm{GHz}}$ for the LzLCS sample. The vertical line demarcates the flat (${\alpha }_{6\mathrm{GHz}}^{3\mathrm{GHz}}>-0.64$) and steep spectrum (${\alpha }_{6\mathrm{GHz}}^{3\mathrm{GHz}}<-0.64$) sources.

4. Discussion

4.1. Estimated thermal fraction of LzLCS sources

We estimated the thermal flux density at two different observed frequencies, 3 GHz and 6 GHz, using Eq. (3). We use the SFR (derived from the dust-corrected Hβ line intensity) and electron temperature derived using [OIII] λλ4363; 4959, 5007 (see Flury et al. 2022a, for details). We further correct this thermal flux density by a factor of ($1-f_{\mathrm{esc}}^{\mathrm{LyC}}$) to account for the escaping LyC photons. We use this thermal flux density and the corresponding observed flux density at 3 GHz (6 GHz) to derive the thermal fraction at these two frequencies.

The mean value of the thermal fraction at 3 GHz (6 GHz) is 0.21 ± 0.11(0.35 ± 0.14). The mean thermal fraction at 6 GHz for nearby star-forming galaxies from the KINGFISHER sample was found to be 0.23 ± 0.13 (Tabatabaei et al. 2017) which is slightly higher although within the uncertainties. Moreover, the range of the thermal fractions of our sample is similar to that of the KINGFISHER sample. Figure 11 shows the distribution of the thermal fraction at 3 GHz in blue and 6 GHz in the orange histogram. The fact that the radio emission from these galaxies is not dominated by thermal emission provides evidence for a non-zero contribution from synchrotron emission in these galaxies. Interestingly, this provides the first direct evidence of SN activity in these galaxies.

Fig. 11. Distribution of the estimated thermal fraction for LzLCS sources. The blue (orange) histograms show the thermal fraction at 3 GHz (6 GHz). The solid (dashed) line shows the expected thermal fraction at 3 GHz (6 GHz) for normal star-forming galaxies with a fth = 0.1 at 1 GHz and αnth = −0.8.

Next, we investigate if SN activity could play a role in facilitating LyC escape as expected from various simulations (e.g., Trebitsch et al. 2017; Ma et al. 2020). In Fig. 12 we study the dependence of $f_{\mathrm{esc}}^{\mathrm{LyC}}$ on the thermal fraction which informs us about increased SN activity or the lack of it. We do not observe any correlation, rather most of our sources including the strongest leakers have thermal fractions close to that expected for normal star-forming galaxies (vertical dashed line). Most of the non-leakers instead are non-thermal dominated thus suggesting older episodes of star formation. A similar lack of dependence is observed between the thermal fraction derived at 6 GHz and $f_{\mathrm{esc}}^{\mathrm{LyC}}$.

Fig. 12. Dependence of $f_{\mathrm{esc}}^{\mathrm{LyC}}$ measured using the UV-fit on the estimated thermal fraction at 3 GHz. The horizontal dashed lines demarcate the region for strong, weak and non-leakers. The vertical dashed line shows the expected thermal fraction at 3 GHz for normal star-forming galaxies with a fth = 0.1 at 1 GHz and αnth = −0.8.

A few sources in our sample are very compact starbursts where it is common to have a significant FFA component. Therefore, it is possible that the true thermal fractions can be even higher than our estimates. A complete radio-SED fitting, particularly with even higher frequency data, is required to accurately constrain the thermal fractions.

4.2. Considering why strong leakers have a flat radio spectrum

As seen in Fig. 9, the observed spectral index across the C- and S-bands, ${\alpha }_{6\mathrm{GHz}}^{3\mathrm{GHz}}$, flattens with increasing $f_{\mathrm{esc}}^{\mathrm{LyC}}$. This hints towards a flattening of ${\alpha }_{6\mathrm{GHz}}^{3\mathrm{GHz}}$ due to the young stellar ages of high $f_{\mathrm{esc}}^{\mathrm{LyC}}$ galaxies ($f_{\mathrm{esc}}^{\mathrm{LyC}}>0.1$). Young stellar systems (ages ≤5 Myr) would have a high thermal fraction due to a lack of significant SN activity (and thus nonthermal emission). This will lead to a gradual flattening of ${\alpha }_{6\mathrm{GHz}}^{3\mathrm{GHz}}$ with decreasing stellar age. In Fig. 13, we show the dependence of ${\alpha }_{6\mathrm{GHz}}^{3\mathrm{GHz}}$ on the estimated thermal fraction at 3 GHz (top panel) and 6 GHz (bottom panel) for our LzLCS sample. We quantitatively show the effect of increasing the f_th, defined at a reference frequency of 3 GHz (top panel) and 6 GHz (bottom panel), on the observed ${\alpha }_{6\mathrm{GHz}}^{3\mathrm{GHz}}$. This is calculated using the simple model for the radio-SED consisting of a thermal and a non-thermal component as described in Eq. (1). The blue line shows the variation for a canonical value of α_nth = −0.8. We also show the effect of steepening α_nth to up to −1.5 in different coloured lines. The dashed horizontal line shows the value of ${\alpha }_{6\mathrm{GHz}}^{3\mathrm{GHz}}=-0.64$ which is expected for a standard f_th = 0.1 at 1 GHz, and α_nth = −0.8. The dashed vertical line shows the value of the expected thermal fraction at 3 and 6 GHz for such a model.

Fig. 13. Effect of increasing fth at a reference frequency of 3 GHz (top panel) and 6 GHz (bottom panel) on the observed spectral index ${\alpha }_{6\mathrm{GHz}}^{3\mathrm{GHz}}$. Here we consider the simplest radio-SED composed of a thermal and non-thermal component as shown in Eq. (1). The different coloured lines show the effect of increasing the αnth from the canonical value of −0.8 to up to −1.5. In both the panels the dashed horizontal lines show the value of ${\alpha }_{6\mathrm{GHz}}^{3\mathrm{GHz}}=-0.64$ for fth = 0.1, and αnth = −0.8 from Eq. (1). And the dashed vertical line shows the expected value of the thermal fraction at 3 and 6 GHz for such a radio spectrum. We also show our LzLCS observations on this plot in blue points.

As expected in all the models with both flat and steep values of α_nth, ${\alpha }_{6\mathrm{GHz}}^{3\mathrm{GHz}}$ flattens as we increase f_th. Almost half of our sample lies within the region bracketed by these models. Thus their variation in the spectral indices in principle can be associated with their different values of f_th and α_nth. However, very flat spectrum sources (${\alpha }_{6\mathrm{GHz}}^{3\mathrm{GHz}}>-0.5$) would need other mechanisms to further flatten the spectrum, for instance, FFA as discussed below. Similarly very steep spectrum sources (${\alpha }_{6\mathrm{GHz}}^{3\mathrm{GHz}}<-1.5$) would either need unusually steep α_nth, or other mechanisms to steepen α_nth, e.g., in a break model or CR escape which we discuss in Sect. 4.3.

Also, FFA has a strong effect on flattening the radio-SED. In normal galaxies, the effect of FFA is usually observed at low frequencies below ∼30 MHz (Chyży et al. 2018). For dense starbursts, FFA can be important at relatively higher frequencies around ∼700 MHz for starburst galaxies like in M 82 (see e.g., Chyży et al. 2018) and close to GHz for ULIRGS (e.g., Clemens et al. 2010; Galvin et al. 2018). However, for low metallicity BCDs, for example in SBS 0335−052, which host a dense super star cluster, Hunt et al. (2004) observed a turnover frequency between 2−3 GHz, corresponding to an EM of ∼8 × 10⁷ pc cm⁻⁶. In Fig. 14, we show the effect of increasing ν_t, from 10 MHz to ∼5 GHz, on the observed ${\alpha }_{6\mathrm{GHz}}^{3\mathrm{GHz}}$. We also show the effect of increasing the f_th from 0.1 to 1.0 on ${\alpha }_{6\mathrm{GHz}}^{3\mathrm{GHz}}$ shown in solid lines. In dashed lines, we show the same set of curves with the effect of steepening the α_nth to −2.0. We notice that in order to flatten ${\alpha }_{6\mathrm{GHz}}^{3\mathrm{GHz}}>0.0$, in all the models we need a ν_t ≳ 1.5 GHz, which corresponds to an EM ∼1.2 × 10⁷ pc cm⁻⁶ assuming a z ∼ 0.3 to convert the observed ν_t to the rest-frame and the T_e of 10⁴ K. This is also in line with the correlation seen between $f_{\mathrm{esc}}^{\mathrm{LyC}}$ and Σ_SFR (e.g., Izotov et al. 2016b; Naidu et al. 2020; Flury et al. 2022b). Hence, increasing Σ_SFR (and thus the EM) can indirectly lead to an increase in ν_t and thus flatten ${\alpha }_{6\mathrm{GHz}}^{3\mathrm{GHz}}$. This can become particularly important for the high $f_{\mathrm{esc}}^{\mathrm{LyC}}$ which has the highest Σ_SFR. Finally, for normal star-forming galaxies in the KINGFISHER sample, it is observed that α_nth flattens with increasing Σ_SFR (Tabatabaei et al. 2017, 2022). Such an effect is also observed in the arms and inner regions within a galaxy where the gas density is higher and the α_nth can flatten to even higher than −0.5 (Basu et al. 2015). More importantly in the case of IZw 18, a nearby low-metallicity BCD, Hunt et al. (2005) finds a rather flat α_nth ≈ −0.5. This is associated with the CR energy index being close to the injection spectrum during SN acceleration (Tabatabaei et al. 2017; Basu et al. 2015) which is expected in galaxies with high SFRs. Alternatively, the high SFR can lead to an amplification of turbulent magnetic fields (e.g., Gressel et al. 2008; Gent et al. 2023) which can scatter high-energy CRs before they lose their energy to synchrotron losses (Tabatabaei et al. 2017). Thus leading to a flat non-thermal spectral index. Our LzLCS sources are several orders of magnitude denser than KINGFISHER sample leading to the possibility that they too have relatively flat α_nth which can overall flatten ${\alpha }_{6\mathrm{GHz}}^{3\mathrm{GHz}}$.

Fig. 14. Effect of increasing the turnover frequency (νt) on ${\alpha }_{6\mathrm{GHz}}^{3\mathrm{GHz}}$. The solid lines are for a canonical αnth = −0.8, with the blue line corresponding to a fth = 0.1. The other lines shows the effect of increasing the fth. The dashed lines are for the same range of fth but with a more steep αnth = −2. The dashed horizontal line shows a value of ${\alpha }_{6\mathrm{GHz}}^{3\mathrm{GHz}}=0$. This shows that we need FFA with a νt of ≳1.5 GHz to flatten ${\alpha }_{6\mathrm{GHz}}^{3\mathrm{GHz}}>0$.

We conclude that since high $f_{\mathrm{esc}}^{\mathrm{LyC}}$ galaxies have a relatively young age of the starburst and relatively high Σ_SFR they can have a high thermal fraction, EM, and possibly a flatter α_nth. These effects can altogether flatten the ${\alpha }_{6\mathrm{GHz}}^{3\mathrm{GHz}}$ which can explain the observed trend for strong leakers in Fig. 9. Finally, a flat spectrum source is not always a strong leaker source. We have observed a few outliers in the ${\alpha }_{6\mathrm{GHz}}^{3\mathrm{GHz}}-f_{\mathrm{esc}}^{\mathrm{LyC}}$ relation. This indicates flat spectrum is possibly necessary but not a sufficient condition for identifying strong leakers.

4.3. Considering why non-leakers have a steep radio spectrum

It is quite remarkable that the correlation between ${\alpha }_{6\mathrm{GHz}}^{3\mathrm{GHz}}$ and $f_{\mathrm{esc}}^{\mathrm{LyC}}$ extends all the way towards non-leakers, and in particular, non-leakers have very steep ${\alpha }_{6\mathrm{GHz}}^{3\mathrm{GHz}}$ (≤1). Several processes related to either the complete loss of CRs or energy-dependent losses are usually invoked to explain the steepening of ${\alpha }_{6\mathrm{GHz}}^{3\mathrm{GHz}}$, for example, a break, cutoff, or curvature at a higher frequency in the radio spectrum (e.g., see the models listed in Klein et al. 2018). A break in the spectrum is produced in an open-box model where the CRs can escape the halo of the galaxy (Lisenfeld et al. 2004). This break occurs due to a high convective speed due to which the CRs can leave the halo (e.g., Lisenfeld et al. 2004; Klein et al. 2018). The break steepens α_nth by 0.5 above a characteristic frequency termed as the break frequency (Klein et al. 2018). The break frequency is related to the convection speed and size of the halo. Thus, an original α_nth ∼ −0.8 can easily steepen to ∼ − 1.3 above the break frequency. A cutoff in a spectrum occurs when there is an abrupt change in the CR energy spectrum arising due to a relatively recent starburst a few Myr ago. This happens in a so-called “single-injection scenario”, where most of the CRs were produced in a recent starburst. The high-energy CRs will lose their energy faster than the low-energy ones, and it will produce a strong break above some CR energy scale (Schlickeiser 1984; Klein et al. 2018). This CR break energy is related to the break frequency and depends inversely on the time of CR injection (see Eq. (4) in Lisenfeld et al. 2004).

Breaks in the radio spectrum were observed in NGC 1569, a nearby dwarf starburst, at around 8 GHz (Israel & de Bruyn 1988), and several other nearby dwarf galaxies at GHz frequencies ranging from 1−12 GHz (Deeg et al. 1993; Klein et al. 2018). Israel & de Bruyn (1988) applied the single-injection CR acceleration model to NGC 1569 and estimated the injection time to ∼5 Myr ago. Lisenfeld et al. (2004) has interpreted the same break in NGC 1569 due to a single-injection CR model (also termed as the cutoff model) or due to CR outflows, however, they favor the latter scenario. Deeg et al. (1993) interpret several of their convex radio spectra in BCDs due to a break from a single-injection CR model.

If non-leakers in the LzLCS sample are at the “the low-phase of a bursty SF episode”, namely, those that recently stopped star formation or exhibit a declining SFR over time, we expect to see to first-order differences in their SFRs derived using UV and Hβ luminosities, since the two trace SFR over long (∼100 Myr) and short timescales (∼10 Myr), respectively (e.g., Kennicutt & Evans 2012). In Fig. 15 we compare the ratio SFR_UV/SFR_Hβ with the EW(Hβ) for our sample. We notice that sources with the lowest EW(Hβ) also have the highest ratio of SFR_UV/SFR_Hβ, which is naturally explained by relative age differences, with younger objects having higher EW(Hβ). To demonstrate this, we overplot predictions from the synthesis models of Schaerer (2003) for metallicities bracketing those of the LzLCS and both for instantaneous burst models and constant SFR. Qualitatively the observed trend of increasing SFR_UV/SFR_Hβ with decreasing EW(Hβ), corresponding to a temporal evolution, is reproduced, and models of instantaneous burst (or other declining SF histories) also naturally lead to SFR_UV/SFR_Hβ > 1.

Fig. 15. Ratio between the SFR derived using UV and Hβ luminosity as a function of the EW(Hβ). We note that typically galaxies with lower EW show a higher UV-based SFR compared to Hβ based SFR. We also overplot the stellar population synthesis models from Schaerer (2003) for a single burst model in solid lines, and continuous star-formation in dashed lines. The different colours blue, orange, green, and red represent the different metallicities corresponding to Z = 0.02 (solar), 0.008, 0.004, and 0.001, respectively.

In Fig. 16, we compare this SFR ratio with the observed ${\alpha }_{6\mathrm{GHz}}^{3\mathrm{GHz}}$ for our LzLCS sample. We note that the sources with the steepest spectra (thus non-leakers), show the highest ratio of SFR_UV/SFR_Hβ. In such a cutoff model we can also estimate the CR injection timescale Lisenfeld et al. (using Eq. (4) from 2004) for a break frequency close to 6 GHz and with magnetic fields ranging from 10−20 μG. We find that in such a model for the entire range of magnetic fields, a break can occur below 6 GHz within 10 Myr. This timescale is consistent with the timescale on which we expect to see differences in the UV-based and Hβ-based SFR, and further supports our idea that the steepening of our LzLCS sources is because of their declining SFR with time.

Fig. 16. Ratio between the SFR derived using UV and Hβ luminosity vs. the ${\alpha }_{6\mathrm{GHz}}^{3\mathrm{GHz}}$ for our LzLCS sample. The horizontal line represents a ratio of 1. And the vertical line shows the ${\alpha }_{6\mathrm{GHz}}^{3\mathrm{GHz}}$ for a canonical radio spectrum. We notice that galaxies with steep spectra have a higher UV-based SFR compared to Hβ based SFR (and also lower EW(Hβ) as seen in Fig. 15).

Further, Heesen et al. (2022) has shown a steepening of the radio spectral index with galaxy stellar mass for nearby galaxies in the LOw Frequency ARray (LOFAR) survey. Their study which is at lower frequencies (144−1440 MHz) than ours, suggests that more massive galaxies have experienced a stronger synchrotron energy loss leading to a steeper radio spectra. In Fig. 17 we observe a similar trend for our LzLCS sources. In Table 4 we also find a good correlation between stellar mass and ${\alpha }_{6\mathrm{GHz}}^{3\mathrm{GHz}}$. Overall, it is difficult to distinguish which galaxy physical parameter drives the steepening of ${\alpha }_{6\mathrm{GHz}}^{3\mathrm{GHz}}$. However, physically in both the scenarios, a significant loss of high-energy CRs can lead to the steepening of ${\alpha }_{6\mathrm{GHz}}^{3\mathrm{GHz}}$.

Fig. 17. Dependance of ${\alpha }_{6\mathrm{GHz}}^{3\mathrm{GHz}}$ on the stellar mass of LzLCS sources. We notice that ${\alpha }_{6\mathrm{GHz}}^{3\mathrm{GHz}}$ steepens with increasing stellar mass. The dashed line shows the ${\alpha }_{6\mathrm{GHz}}^{3\mathrm{GHz}}$ for a canonical radio spectrum.

Alternatively, it is possible that the steepening is caused due to CR escape from advection or diffusion. In the intriguing case of the starburst dwarf galaxy, IC 10, there is evidence of SN feedback in the form of a non-thermal superbubble (Yang & Skillman 1993). The spectral index is steep (close to −1.0) and is found to be consistent with CR escape due to advection or diffusion (Heesen et al. 2011, 2018). Strong outflows are also observed in compact starbursts (e.g., Weiner et al. 2009; Heckman et al. 2011; Chisholm et al. 2017). Given the typically small sizes of these galaxies, the advection/diffusion timescales are expected to less than 1 Myr (Sebastian & Bait 2019). However, it is not directly very clear why the advection/diffusion driven CR escape should be more prominent in non-leakers as opposed to leakers. Farcy et al. (2022) showed that the presence of CRs in simulated galaxies tend to smooth out the ISM thus filling it with denser gas and hence reducing $f_{\mathrm{esc}}^{\mathrm{LyC}}$. In the future, the inclusion of high and low-frequency radio data should better constrain the radio-SED and help us in distinguishing between the break or cutoff model for the steepening of ${\alpha }_{6\mathrm{GHz}}^{3\mathrm{GHz}}$.

4.4. Overall picture relating LyC escape and RC variations

In summary, since the radio spectral index is a proxy for the age of the starburst (Cannon & Skillman 2004; Hirashita & Hunt 2006), we may possibly observe a relation between ${\alpha }_{6\mathrm{GHz}}^{3\mathrm{GHz}}$ and $f_{\mathrm{esc}}^{\mathrm{LyC}}$. Using a timescale argument we propose the following evolutionary picture. After the onset of a starburst (< 5 Myr) typically in dense ISM conditions, the radio spectrum will be flat (${\alpha }_{6\mathrm{GHz}}^{3\mathrm{GHz}}>-0.65$) due to (i) thermally dominated radio emission as most of the young stars haven’t yet undergone supernova (ii) the less dominant non-thermal emission which would be arising when some of the first supernovae start to explode will accelerate fresh CRs that have a flat non-thermal spectrum (0.65) and (iii) the presence of dense ISM conditions can lead to a flat radio spectrum due to the strong effect of FFA at GHz frequencies. These physical conditions which lead to a flat ${\alpha }_{6\mathrm{GHz}}^{3\mathrm{GHz}}$ are also when the $f_{\mathrm{esc}}^{\mathrm{LyC}}$ is highest. At later times (5−10 Myr) when the starburst is declining along with the $f_{\mathrm{esc}}^{\mathrm{LyC}}$, high-energy CRs also lose their energies (in ∼10 Myr) which results in the steepening of ${\alpha }_{6\mathrm{GHz}}^{3\mathrm{GHz}}$ (≤0.65). Consequently, the galaxies will quickly move down along the ${\alpha }_{6\mathrm{GHz}}^{3\mathrm{GHz}}-f_{\mathrm{esc}}^{\mathrm{LyC}}$ relation in Fig. 9 in around 10 Myr. If there are any further bursts of star formation in these galaxies they can move along this sequence. Alternatively, the steepening can also be caused due to CR escape. For galaxies showing continuous star formation over 100 Myr timescale, the CRs would have diffused in the galaxy and have a non-thermal spectral index close to the canonical value of −0.8 (Condon 1992). Such galaxies would then have a total ${\alpha }_{6\mathrm{GHz}}^{3\mathrm{GHz}}$ close to −0.65 and $f_{\mathrm{esc}}^{\mathrm{LyC}}<1$%.

Finally, the fact that some properties of the radio continuum correlate with the LyC escape fraction could indicate that the escape of ionizing photons is not highly directional, namely, which happens in a relatively isotropic manner. This is due to the fact that the observational determination of LyC escape ($f_{\mathrm{esc}}^{\mathrm{LyC}}$) measures a directional quantity, namely the escape of ionizing photons along the line-of-sight towards the observer, whereas RC emission is a priori a volumetric probe and orientation-independent. Therefore the finding of a correlation between $f_{\mathrm{esc}}^{\mathrm{LyC}}$ and RC continuum properties would indicate that $f_{\mathrm{esc}}^{\mathrm{LyC}}$ is also orientation-independent, to first order. This is consistent with previous results inferred from correlations between global galaxy properties and $f_{\mathrm{esc}}^{\mathrm{LyC}}$ (Wang et al. 2019). Now in reality the relation between ${\alpha }_{6\mathrm{GHz}}^{3\mathrm{GHz}}$ and $f_{\mathrm{esc}}^{\mathrm{LyC}}$, for instance, shows a considerable scatter, which might partly be driven by some orientation-dependence. Better statistics, and a more fundamental understanding of the origin of these correlations will therefore be needed to properly address the question of non-isotropic LyC escape. Other volumetric measures at other wavelengths also exist. For example, the nebular emission lines in the optical are thought to be essentially optically thin, and $f_{\mathrm{esc}}^{\mathrm{LyC}}$ is found to correlate with O₃₂, which could again indicate relatively isotropic LyC emission. However, the scatter in these relations is also large (see Izotov et al. 2021b; Flury et al. 2022b), and there is no unique and direct causal connection between LyC escape and O₃₂, which a priori depends on several parameters, such as the ionization parameter, metallicity, age of stellar population(s), etc., and $f_{\mathrm{esc}}^{\mathrm{LyC}}$ (Nakajima & Ouchi 2014; Ramambason et al. 2020). Possible orientation and geometrical effects of LyC escape are so far difficult to establish observationally from the available low-z galaxy samples.

5. Conclusions and future work

For the first time, we studied the radio continuum (RC) properties of galaxies with LyC measurements based on VLA C- and S-band observations for a sample of 53 galaxies. We also obtained L-band observations for 17 of these sources. These galaxies were selected from the recent Low-z Lyman Continuum Survey (LzLCS; Flury et al. 2022a,b), which contains a large sample of low-z (z ∼ 0.3 − 0.4) star-forming galaxies with LyC measurements. We detected the RC at both C- and S-bands for about half of our sample, where we were then able to measure the total spectral index (${\alpha }_{6\mathrm{GHz}}^{3\mathrm{GHz}}$) and thermal fraction.

Our study highlights an interesting connection between the radio spectral index and the escape of LyC photons. Overall, we associate this with the rapidly changing SFH (on a few Myr timescale) due to the onset of starburst and its decline in LzLCS galaxies. These changes in SFHs are traced by the CR energy spectrum (and, thus, by the radio spectral index) and they also have an imprint on the LyC escape fraction. Our main results can be summarized as follows.

1. The RC spectral index between C- and S-bands (${\alpha }_{6\mathrm{GHz}}^{3\mathrm{GHz}}$) for LzLCS sources shows a large variation from that of typical star-forming galaxies, with sources showing a very flat (≥0.1) to very steep spectral index (≤1).

2. We find that LzLCS galaxies have a normal to significant contribution of the non-thermal component to the total radio emission. This provides the first direct evidence of SN activity in these galaxies and suggests that SN can contribute to the escape of LyC photons.

3. We found a positive correlation between ${\alpha }_{6\mathrm{GHz}}^{3\mathrm{GHz}}$ and the escape fraction of Lyman continuum photons, $f_{\mathrm{esc}}^{\mathrm{LyC}}$; thus, strong leakers have flat ${\alpha }_{6\mathrm{GHz}}^{3\mathrm{GHz}}$ and non-leakers exhibit a very steep ${\alpha }_{6\mathrm{GHz}}^{3\mathrm{GHz}}$.

4. We argue that a combination of young ages (high f_th), FFA with a turnover at GHz frequencies and CRs with a flat spectrum (young) can lead to the flattening of ${\alpha }_{6\mathrm{GHz}}^{3\mathrm{GHz}}$. Such physical conditions typical of young and dense starbursts are also seen when galaxies are expected to exhibit the highest $f_{\mathrm{esc}}^{\mathrm{LyC}}$.

5. Steep spectrum sources are associated with low $f_{\mathrm{esc}}^{\mathrm{LyC}}$. They show evidence of declining SFRs on short timescales and with higher stellar masses, both of which may lead to strong synchrotron losses and a steeper radio spectrum as a result.

6. Our work highlights the evolutionary nature of $f_{\mathrm{esc}}^{\mathrm{LyC}}$, on short timescales (∼10 Myr) and its association with SNs, the CR energy spectrum, and magnetic fields related to the escape of LyC photons.

7. The RC–SFR relation of LzLCS sources show a large scatter. Using RC as an SFR tracer for such galaxies is challenging when only a single or narrow frequency range is available.

8. The SFR^RC shows a systematic offset from the SFR${}_{6\mathrm{GHz}}^{\mathrm{RC}}$–SFR_UV (SFR_Hβ) relation for local normal star-forming galaxies from Tabatabaei et al. (2017). We calibrated the SFR_UV (and SFR_Hβ) for a combination of SFR^RC and ${\alpha }_{6\mathrm{GHz}}^{3\mathrm{GHz}}$ values, which can be used to estimate their SFRs. Information on the radio-SED will be crucial to robustly calibrate the radio-based SFRs in such galaxies. Since these galaxies are analogs of high-z galaxies (z > 3), our study gives important insights into such galaxies. This is particularly relevant in the context of the upcoming highly sensitive radio facilities, such as SKA and ngVLA, aimed at detecting star-forming galaxies at high redshifts.

In our current study, using just the C- and S-band data, we cannot robustly estimate all of the physical parameters we are looking to obtain. In the future, we plan to conduct low-frequency (< 1.5 GHz) and high-frequency (> 8 GHz) follow-up observations of LzLCS sources to better study their radio-SEDs. This will allow us to fit complex radio-SED models, which will be useful in accurately estimating the non-thermal spectral index and the thermal fraction. This will also help us estimate some of the physical parameters that are lacking; for example, the SN rate, CR energy spectrum, and magnetic fields, and their relation to LyC escape. These approaches will improve our understanding of the role of SN feedback in the escape of LyC photons. Finally, our current radio follow-up of LzLCS has missed a number of strong leakers. This will be the subject of a follow-up study in the future.

Acknowledgments

We thank the referee for their comments which helped improve the draft. O. Bait (OB) is supported by the AstroSignals Sinergia Project funded by the Swiss National Science Foundation. OB would like to thank Jennifer Schober and Aayush Saxena for useful comments on the manuscript. OB would also like to thank Drew Medlin and the NRAO science helpdesk team for their extensive help with running the VLA pipeline and with the VLA data reduction. We would like to thank Brad Koplitz for his help with setting up some of the VLA scheduling blocks. ASL acknowledges support from Knut and Alice Wallenberg Foundation unding for the Sloan Digital Sky Survey IV has been provided by the Alfred P. Sloan Foundation, the U.S. Department of Energy Office of Science, and the Participating Institutions. SDSS-IV acknowledges support and resources from the Center for High Performance Computing at the University of Utah. The SDSS website is www.sdss4.org. SDSS-IV is managed by the Astrophysical Research Consortium for the Participating Institutions of the SDSS Collaboration including the Brazilian Participation Group, the Carnegie Institution for Science, Carnegie Mellon University, Center for Astrophysics | Harvard & Smithsonian, the Chilean Participation Group, the French Participation Group, Instituto de Astrofísica de Canarias, The Johns Hopkins University, Kavli Institute for the Physics and Mathematics of the Universe (IPMU)/University of Tokyo, the Korean Participation Group, Lawrence Berkeley National Laboratory, Leibniz Institut für Astrophysik Potsdam (AIP), Max-Planck-Institut für Astronomie (MPIA Heidelberg), Max-Planck-Institut für Astrophysik (MPA Garching), Max-Planck-Institut für Extraterrestrische Physik (MPE), National Astronomical Observatories of China, New Mexico State University, New York University, University of Notre Dame, Observatário Nacional/MCTI, The Ohio State University, Pennsylvania State University, Shanghai Astronomical Observatory, United Kingdom Participation Group, Universidad Nacional Autónoma de México, University of Arizona, University of Colorado Boulder, University of Oxford, University of Portsmouth, University of Utah, University of Virginia, University of Washington, University of Wisconsin, Vanderbilt University, and Yale University.

References

Appendix A: Observation details of VLA S-, C-, and L-band sources in our sample

In Table A.1, we describe the details of the C- and S-band observations for our sources. And Table A.2 shows the details of the L-band observations for a sub-sample of these sources.

Appendix B: VLA S- and C-band images for the detections in our sample

In Figs. B.1 and B.2, we show the VLA C-band (in orange) and S-band (in blue) contours overlaid on each other for our detection sample.

Fig. B.1. RC images for all the LzLCS detections in S-band (blue contours) and C-band (orange contours). The black crosshair represents the optical centre. The synthesized beam is shown blue (orange) filled ellipticals for S-band (C-band). The contours start from 3-σ, and successively increase in powers of $\sqrt{2}$. The dashed contours represent the -3-σ contour.

Fig. B.2. Continued from Fig. B.1.

Appendix C: Radio-SED for our LzLCS sample at GHz frequencies

In Figs. C.1 and C.2, we present the radio-SED across C-, S-, and also L-band (for those with C, S-band detections, and L-band observations).

Fig. C.1. Same as Fig. 3 for the remaining LzLCS sample with C- and S-bands detections. The red arrows represent the 3-σ upper limits for L-band non-detections in our sample.

Fig. C.2. Continued from Fig. C.1.

Appendix D: Possible contamination of the radio flux of J130128+510451 from a neighbouring AGN jet

In Fig. D.1, we show the large-scale VLA S-band image around J130128+510451. We notice that a neighbouring radio-loud source (J130126+510500: RA: 13h01m25.908s; DEC: +51d05m00.66s) shows a strong jet-like emission along the north-west direction. It possibly has diffuse extensions along the south-east direction. This latter extension coincidentally passes directly through our target, thus potentially contaminating the radio flux of our target. However, deeper radio data is required to confirm this diffuse south-east jet-like feature and its association to J130126+510500.

Fig. D.1. Large-scale VLA S-band image around the target, J130128+510451. The radio emission from the target seems to be contaminated by the jet emission from a nearby radio-loud AGN. The radio contour levels start from 3.5σ and rise in powers of 1.2. We also added three more contours at 0.1, 1.0, and 2.0 mJy/beam to show the bright core emission from the neighbouring source, J130126+510500.

Appendix E: Details on the SFR^RC, ${\alpha }_{6\mathrm{GHz}}^{3\mathrm{GHz}}$ fit to SFR from other tracers

We perform a linear fit between log SFR^RC, ${\alpha }_{6\mathrm{GHz}}^{3\mathrm{GHz}}$, and the log SFR measured using the UV flux and Hβ line flux. Here, the fit is performed separately for SFR${}_{3\mathrm{GHz}}^{\mathrm{RC}}$ and SFR${}_{6\mathrm{GHz}}^{\mathrm{RC}}$, to SFR_UV and SFR_Hβ. Since all three variables have error bars on them, we use the bootstrap method to take into account the uncertainties. We draw 100 random samples for each of the three variables assuming a Gaussian distribution with the 1-σ errors on each variable representing the width of the Gaussian. This leads to a million combinations of variable pairs (a pair consists of three variables, two SFRs, and ${\alpha }_{6\mathrm{GHz}}^{3\mathrm{GHz}}$). We performed a linear fit on each of these pairs using Eq. 5. The best-fit parameters and the uncertainty are estimated by taking the mean and standard deviation of each of the fit parameters from the entire sample. The residual in the fit is the mean of the residual from the entire sample. Table 3 provides the fit parameters, corresponding errors, and the residual from the fit.

All Tables

All Figures

Fig. 1. VLA images of two sources in our sample, J124423+021540 (top panels) and J095700+235709 (bottom panels), shown in white contours with the SDSS r-band image (Blanton et al. 2017) in the background. On the left (right) panels we show the VLA S-band (C-band) image contours. The lowest contours in each panels corresponds to 3σ level and they rise in powers of $\sqrt{2}$. σ values correspond to ∼8 and ∼5 μJy beam−1 for S- and C-band respectively.

In the text

Fig. 2. Typical radio-SED with a mixture of thermal and non-thermal emission (Eq. (1)) in blue line. The corresponding thermal and non-thermal emission components are shown in blue dashed and dotted lines. We also show the effect of FFA (Eq. (2)) with a turnover at 1 GHz in orange. And the cutoff and curved model (from Klein et al. 2018) with a high-frequency break at 6 GHz in green and red respectively. The grey-shaded area shows the region of the SED explored in this work.

In the text

Fig. 3. Observed radio continuum SEDs at three frequencies: 1.5, 3, and 6 GHz for three sources in our sample. The error bars represent 1σ errors on the flux measurement. The gray line in each panel is an overlay (not a fit) for a simple radio-SED (Eq. (1)) normalized to the 3 GHz observed flux density with fth = 0.1 at 1 GHz and αnth = −0.8 commonly seen in star-forming galaxies. The left panel shows an example of a radio-SED which is close to the standard radio-SED (Eq. (1)). The middle panel shows an example that has a steep spectrum above 3 GHz and is close to the standard radio-SED below 3 GHz. The right panel shows the case of a broken SED.

In the text

	Fig. 4. Distribution of the total spectral index of LzLCS sources which spans a wide range from flat to steep. The dashed line shows the value of ${\alpha }_{6\mathrm{GHz}}^{3\mathrm{GHz}}$ expected in normal star-forming galaxies with a fth = 0.1 at 1 GHz and αnth = −0.8, ${\alpha }_{6\mathrm{GHz}}^{3\mathrm{GHz}}=-0.64$.
In the text

Fig. 5. LzLCS sources on the RC–SFR relation at 3 GHz top panels and 6 GHz in the bottom panels. The left (right) panel shows the UV-based (Hβ-based) SFRs for the LzLCS sources. In all panels, the blue points represent detections in the corresponding band and the orange points represent the 3σ upper limits for non-detections. The solid black line is plotted for the standard RC–SFR relation with αnth = −0.8. The dotted, dashed and dotted dashed lines are for different values of αnth from flat (−0.3) to very steep (−2), respectively. The red solid line represents the RCth–SFR relation from Eq. (3). In all these cases, we have used a fixed value of Te = 10 000 K.

In the text

Fig. 6. Ratio between the observed RC luminosity and that expected from the RC–SFR relation from a standard radio-SED with αnth = −0.8, fth = 0.1 at 1 GHz, and Te = 104 K, at 3 GHz and 6 GHz (top and bottom panels respectively) vs. the spectral index ${\alpha }_{6\mathrm{GHz}}^{3\mathrm{GHz}}$. The left (right) panel shows the ratio with the expected RC for the UV-based (Hβ-based) SFR. The dashed line shows a ratio of 1 for reference. Here, we have used only the S- and C-band detection sample.

In the text

Fig. 7. Comparison between SFR${}_{6\mathrm{GHz}}^{\mathrm{RC}}$ and SFRUV (circles) and SFRHβ (triangles). Each of the points is color coded by their apparent ${\alpha }_{6\mathrm{GHz}}^{3\mathrm{GHz}}$. We over plot the relation between SFR${}_{6\mathrm{GHz}}^{\mathrm{RC}}$ and other tracers from Tabatabaei et al. (2017) in blue (SFRUV calibration) and orange (SFRHβ calibration). The blue and orange shaded region shows the scatter in these relations.

In the text

	Fig. 8. Combined fit between the linear combination of log10 SFR${}_{6\mathrm{GHz}}^{\mathrm{RC}}$ and ${\alpha }_{6\mathrm{GHz}}^{3\mathrm{GHz}}$, and log10 SFRUV. The blue line shows the x = y line. Note: the residual in the fit is 0.21, which is similar to the typical scatter between different SFR tracers from Tabatabaei et al. (2017).
In the text

Fig. 9. Dependence of $f_{\mathrm{esc}}^{\mathrm{LyC}}$ on the total spectral index measured across C- (6 GHz) and S- (3 GHz) bands. The left panel shows the dependence on $f_{\mathrm{esc}}^{\mathrm{LyC}}$ measured using the UV-fit method. And the right panel shows the dependence on $f_{\mathrm{esc}}^{\mathrm{LyC}}$ measured using the Hβ method Flury et al. (2022a). The blue points in both panels show LyC detections and the orange points show LyC non-detections (1σ upper limits). We also show two known nearby leakers from the literature, Haro 11 (brown star) and Tol 1247−232 (black star). In both panels the horizontal lines demarcate the regions for strong leakers ($f_{\mathrm{esc}}^{\mathrm{LyC}}>0.1$), weak leakers ($0.01<f_{\mathrm{esc}}^{\mathrm{LyC}}<0.1$), and non-leakers ($f_{\mathrm{esc}}^{\mathrm{LyC}}<0.01$). The vertical line shows the region above and below a spectral index of a standard radio-SED (${\alpha }_{6\mathrm{GHz}}^{3\mathrm{GHz}}=-0.64$) demarcating the region for flat and steep spectrum sources. See text for details.

In the text

	Fig. 10. Dependence of β1550, O32 and ΣSFR on the spectral index ${\alpha }_{6\mathrm{GHz}}^{3\mathrm{GHz}}$ for the LzLCS sample. The vertical line demarcates the flat (${\alpha }_{6\mathrm{GHz}}^{3\mathrm{GHz}}>-0.64$) and steep spectrum (${\alpha }_{6\mathrm{GHz}}^{3\mathrm{GHz}}<-0.64$) sources.
In the text

	Fig. 11. Distribution of the estimated thermal fraction for LzLCS sources. The blue (orange) histograms show the thermal fraction at 3 GHz (6 GHz). The solid (dashed) line shows the expected thermal fraction at 3 GHz (6 GHz) for normal star-forming galaxies with a fth = 0.1 at 1 GHz and αnth = −0.8.
In the text

	Fig. 12. Dependence of $f_{\mathrm{esc}}^{\mathrm{LyC}}$ measured using the UV-fit on the estimated thermal fraction at 3 GHz. The horizontal dashed lines demarcate the region for strong, weak and non-leakers. The vertical dashed line shows the expected thermal fraction at 3 GHz for normal star-forming galaxies with a fth = 0.1 at 1 GHz and αnth = −0.8.
In the text

Fig. 13. Effect of increasing fth at a reference frequency of 3 GHz (top panel) and 6 GHz (bottom panel) on the observed spectral index ${\alpha }_{6\mathrm{GHz}}^{3\mathrm{GHz}}$. Here we consider the simplest radio-SED composed of a thermal and non-thermal component as shown in Eq. (1). The different coloured lines show the effect of increasing the αnth from the canonical value of −0.8 to up to −1.5. In both the panels the dashed horizontal lines show the value of ${\alpha }_{6\mathrm{GHz}}^{3\mathrm{GHz}}=-0.64$ for fth = 0.1, and αnth = −0.8 from Eq. (1). And the dashed vertical line shows the expected value of the thermal fraction at 3 and 6 GHz for such a radio spectrum. We also show our LzLCS observations on this plot in blue points.

In the text

Fig. 14. Effect of increasing the turnover frequency (νt) on ${\alpha }_{6\mathrm{GHz}}^{3\mathrm{GHz}}$. The solid lines are for a canonical αnth = −0.8, with the blue line corresponding to a fth = 0.1. The other lines shows the effect of increasing the fth. The dashed lines are for the same range of fth but with a more steep αnth = −2. The dashed horizontal line shows a value of ${\alpha }_{6\mathrm{GHz}}^{3\mathrm{GHz}}=0$. This shows that we need FFA with a νt of ≳1.5 GHz to flatten ${\alpha }_{6\mathrm{GHz}}^{3\mathrm{GHz}}>0$.

In the text

Fig. 15. Ratio between the SFR derived using UV and Hβ luminosity as a function of the EW(Hβ). We note that typically galaxies with lower EW show a higher UV-based SFR compared to Hβ based SFR. We also overplot the stellar population synthesis models from Schaerer (2003) for a single burst model in solid lines, and continuous star-formation in dashed lines. The different colours blue, orange, green, and red represent the different metallicities corresponding to Z = 0.02 (solar), 0.008, 0.004, and 0.001, respectively.

In the text

Fig. 16. Ratio between the SFR derived using UV and Hβ luminosity vs. the ${\alpha }_{6\mathrm{GHz}}^{3\mathrm{GHz}}$ for our LzLCS sample. The horizontal line represents a ratio of 1. And the vertical line shows the ${\alpha }_{6\mathrm{GHz}}^{3\mathrm{GHz}}$ for a canonical radio spectrum. We notice that galaxies with steep spectra have a higher UV-based SFR compared to Hβ based SFR (and also lower EW(Hβ) as seen in Fig. 15).

In the text

	Fig. 17. Dependance of ${\alpha }_{6\mathrm{GHz}}^{3\mathrm{GHz}}$ on the stellar mass of LzLCS sources. We notice that ${\alpha }_{6\mathrm{GHz}}^{3\mathrm{GHz}}$ steepens with increasing stellar mass. The dashed line shows the ${\alpha }_{6\mathrm{GHz}}^{3\mathrm{GHz}}$ for a canonical radio spectrum.
In the text

	Fig. B.1. RC images for all the LzLCS detections in S-band (blue contours) and C-band (orange contours). The black crosshair represents the optical centre. The synthesized beam is shown blue (orange) filled ellipticals for S-band (C-band). The contours start from 3-σ, and successively increase in powers of $\sqrt{2}$. The dashed contours represent the -3-σ contour.
In the text

	Fig. B.2. Continued from Fig. B.1.
In the text

	Fig. C.1. Same as Fig. 3 for the remaining LzLCS sample with C- and S-bands detections. The red arrows represent the 3-σ upper limits for L-band non-detections in our sample.
In the text

	Fig. C.2. Continued from Fig. C.1.
In the text

	Fig. D.1. Large-scale VLA S-band image around the target, J130128+510451. The radio emission from the target seems to be contaminated by the jet emission from a nearby radio-loud AGN. The radio contour levels start from 3.5σ and rise in powers of 1.2. We also added three more contours at 0.1, 1.0, and 2.0 mJy/beam to show the bright core emission from the neighbouring source, J130126+510500.
In the text