Low-frequency turnover star-forming galaxies I: Radio continuum observations and global properties

Abstract

There is growing evidence that the broadband radio spectral energy distributions (SEDs) of star-forming galaxies (SFGs) contain a wealth of complex physics. In this paper we aim to determine the physical emission and loss processes causing radio SED curvature and steepening to see what observed global astrophysical properties, if any, are correlated with radio SED complexity. To do this, we have acquired radio continuum data between 70 MHz and 17 GHz for a sample of 19 southern local ($z \lt 0.04$) SFGs. Of this sample 11 are selected to contain low-frequency ($ \lt $300 MHz) turnovers (LFTOs) in their SEDs and eight are control galaxies with similar global properties. We model the radio SEDs for our sample using a Bayesian framework whereby radio emission (synchrotron and free-free) and absorption or loss processes are included modularly. We find that without the inclusion of higher frequency data ($ \gt $17 GHz) single synchrotron power-law based models are always preferred for our sample; however, additional processes including free-free absorption (FFA) and synchrotron losses are often required to accurately model radio SED complexity in SFGs. The fitted synchrotron spectral indices range from $-0.45$ to $-1.07$ and are strongly anticorrelated with stellar mass suggesting that synchrotron losses are the dominant mechanism acting to steepen the spectral index in larger/more massive nearby SFGs. We find that LFTOs in the radio SED are independent from the inclination of SFGs; however, higher inclination galaxies tend to have steeper fitted spectral indices indicating losses to diffusion of cosmic ray electrons into the galactic halo. Four of five of the merging systems in our SFG sample have elevated specific star formation rates and flatter fitted spectral indices with unconstrained LFTOs. Lastly, we find no significant separation in global properties between SFGs with or without modelled LFTOs. Overall these results suggest that LFTOs are likely caused by a combination of FFA and ionisation losses in individual recent starburst regions with specific orientations and interstellar medium properties that, when averaged over the entire galaxy, do not correlate with global astrophysical properties.

1. Introduction

Galaxy spectral energy distributions (SEDs) contain crucial information about the astrophysical processes occurring within them, with varying contributions from gas, dust, stars, and, active galactic nuclei (AGN) at different wavelengths (Walcher et al. 2011; Conroy 2013). Different emission mechanisms dominate at distinct wavelength regimes and vary over time allowing us to measure certain astrophysical quantities and estimate their histories (e.g. Thorne et al. 2023). For example, stellar masses can be measured in the optical and/or near-infrared regime where stellar emission is dominant (Taylor 2011). Star formation rates (SFRs) can be measured instantaneously at ultraviolet (UV) wavelengths due to the emission being dominated by short-lived high mass OB stars (HMS; $M_{\odot} \geq 8 \, {\rm M}_{\odot}$ ) (Kennicutt & Evans 2012). However, optical and UV photons are absorbed by intervening dust which is heated, reradiating this energy at infrared (IR) wavelengths (Draine 2003). Thus corrections, or a combination of IR and UV emission, are often required to provide accurate estimates of the SFR (Bell et al. 2005; Kennicutt et al. 2009; Davies et al. 2016; Delvecchio et al. 2021). Radio continuum emission, however, is impervious to the effects of dust attenuation due to its long wavelength and can be used as a direct probe into the star-formation activity of ‘normal’ SFGs (galaxies without AGN; Condon 1992).

The shape of a galaxy’s radio SED can provide further insight into the astrophysical processes occurring within but has been often simplified and underutilised due to lack of observations and poor spectral sampling (often including only two photometric points over four orders of magnitude in radio frequency). The spectral index ( $\alpha$ , where flux density, $S_{\nu} \propto \nu^{\alpha}$ ) of a galaxy’s radio emission can be used to explore radiation laws and examine the interplay between the heating and cooling mechanisms of the interstellar medium (ISM), magnetic fields and relativistic cosmic rays (Murphy et al. 2008; Lacki & Thompson 2010; Roth et al. 2023; Thorne et al. 2023). In galaxies that contain an AGN, the accretion activity and timescales, as well as how AGN jets impact the ISM of the host galaxy, can also be explored through the spectral index and variability over time (Fabian 2012; Slob et al. 2022; Ross et al. 2023). More widely studied is the direct proportionality between the radio luminosity and SFR in normal SFGs (Condon 1992; Bell 2003; Davies et al. 2017; Heesen et al. 2022).

Radio continuum emission from normal SFGs is produced as a result of two different emission processes of which both are dependent on HMS formation. Thermal free-free emission is produced in the HII envelopes around HMS’s that are ionised by UV flux and provides a direct, near-instantaneous measure of SFR. Thermal emission therefore strongly depends on the number of ionising UV photons with only a weak dependence on electron temperature (see Eq. 2; Condon 1992). Non-thermal synchrotron emission is produced as relativistic cosmic ray electrons (CREs), accelerated by the shocks produced after Type II and Type Ib supernova, interact with the large scale magnetic fields within a galaxy. Synchrotron emission typically dominates below 30 GHz and is delayed by at least $\approx$ 10 Myr (Condon 1992). The CREs that produce the synchrotron emission lose energy as they propagate away from star-forming regions with lifetimes of $\sim$ 50–300 Myr and scale lengths up to $\approx$ 7 kpc dependent on their injection energy and magnetic field strength (Condon 1992; Murphy et al. 2008; Heesen et al. 2023; Heesen et al. 2024). Synchrotron emission has a characteristic power-law emission spectrum ( $\alpha \approx$ -0.8) whilst free-free emission has almost a flat spectral index ( $\alpha \approx$ -0.1). These differences in emission spectrum can be used to separate their contribution to a galaxy’s radio SED.

Current radio continuum SFR indicators, mostly related to the 1.4 GHz luminosity, are typically calibrated against far-infrared (FIR) or total-infrared (TIR) measurements using the FIR/TIR to radio correlation (FRC/TRC; Condon 1992; Yun, Reddy, & Condon 2001; Bell 2003; Molnár et al. 2021; Heesen et al. 2022). The FRC/TRC has been shown to have a tight linear relationship over many orders of magnitude in FIR/TIR and radio luminosities for SFGs (Condon 1992; Yun et al. 2001; Bell 2003; Molnár et al. 2021). This strong correlation arises as a result of the common origin of HMS formation for both radio and IR emission and is often interpreted using calorimetric models (Voelk 1989; Lacki & Thompson 2010). These calorimetric models assume that galaxies are optically thick to UV emission such that all UV emission is absorbed and re-emitted as FIR radiation (Holwerda et al. 2005) and that CREs release all their energy as synchrotron emission before they leave the galaxy. Current research, however, suggests a ‘conspiracy’ which maintains the FRC/TRC in non-calorimetric models. For example, lower mass galaxies with UV continuum leakage are also small enough such that their CRE diffusive escape time is less than than the synchrotron cooling time (Niklas, Klein, & Wielebinski 1997; Murphy et al. 2009; Lacki & Thompson 2010; Basu et al. 2015; Heesen et al. 2022). Recent findings by Cook et al. (2024) have also shown the SFR history is an important factor when considering the FRC/TRC for calibration of the radio SFR. Synchrotron emission relies on core-collapse supernova from HMS to accelerate the CREs such that the emission is delayed compared to IR emission which is partially responsible for the non-linearity at low masses and scatter in the relationship. Current research suggests the FRC/TRC for SFGs decreases with redshift (Sargent et al. 2010; Magnelli et al. 2015; Delhaize et al. 2017; Delvecchio et al. 2021; Molnár et al. 2021); however, this is primarily thought to be a selection effect of sampling bias towards higher mass/luminosity galaxies at higher redshift.

Radio continuum emission can also impacted by a range of loss and absorption processes that leads to the increased complexity in their radio SEDs. At low frequencies free-free absorption (FFA) and ionisation losses are thought to contribute to low-frequency turnovers (LFTOs) recently observed in the radio SEDs of some ultra-luminous/luminous IR galaxies (ULIRGS/LIRGS Clemens et al. 2010; Galvin et al. 2018; Dey et al. 2022; Dey et al. 2024). FFA is dependent on the electron density along our line of sight. As such, high density or heavily obscured starburst regions are expected to absorb the radio emission towards lower frequencies. There is little to no evidence of a relationship between FFA and galactic inclination (Hummel 1991; Chyży et al. 2018), which therefore suggests that the dense ISM within individual star-forming regions is the primary absorber rather than the lower density ISM within the galactic disc. Low energy CREs can also lose energy due to the ionisation of atomic and molecular hydrogen and the energy loss is directly proportional to the number density of neutral atoms and molecules (Longair 2011). Ionisation losses act to flatten the low frequency spectral index ( $\Delta\alpha \leq$ 0.5; Basu et al. 2015) if the ionisation loss timescale is less than the synchrotron loss timescale ( $\approx$ 100 Myr; Murphy et al. 2009; Heesen et al. 2023) which can occur in regions of a relatively low total gas mass surface density ( $\Sigma_{gas}$ = 2.5–5 M $_\odot\,$ pc^-2; Basu et al. 2015).

At higher frequencies synchrotron cooling and inverse-Compton (IC) losses become more dominant with their combined effect gradually steepening the synchrotron spectral index by $\Delta\alpha = -0.5$ (see Eq. 12; Condon 1992). IC losses occur as the synchrotron emitting CREs are scattered by photons and is dependent on the photon energy density, both from IR ISM emission and the cosmic microwave background (CMB), CRE energy and redshift (in the case of CMB photons; Lacki & Thompson 2010; Roth et al. 2023). Synchrotron cooling losses, however, are primarily dependent on a galaxy’s magnetic field and size such that the synchrotron electrons lose energy as they propagate away from their injection site before escaping into the galactic halo (Murphy et al. 2008; Lacki & Thompson 2010; Heesen et al. 2024; Roth et al. 2024). Therefore higher frequency spectral steepening in nearby SFGs is expected to be more prevalent in both larger galaxies and extremely compact starbursts with high gas surface densities. Bremsstrahlung, adiabatic and diffusion losses also play a factor in reducing the overall radio emission observed but they do not impact the curvature of synchrotron spectral index as they are generally subdominant and their CRE loss timescales are almost energy independent (Roth et al. 2023).

Somewhat surprisingly, the combined effect of these emission, loss, and absorption processes result in an observed SFG spectral index of $\alpha \approx -0.7$ at the typically surveyed 1.4 GHz radio frequency. As all-sky radio continuum surveys such as GLEAM (Hurley-Walker 2017) and the LOw-Frequency ARray (LOFAR) Two-metre Sky Survey (LoTSS) (Shimwell et al. 2017) have begun to fill out the low frequency regime (50–300 MHz) and targeted observations cover individual sources or regions above 1.4 GHz, the formerly hidden complexity of the radio SED is being revealed. Recent work has begun to uncover and explain this radio SED complexity with the inclusion of FFA components (Clemens et al. 2010; Galvin et al. 2018; Dey et al. 2022; Dey et al. 2024); however, the samples are limited to unresolved LIRGS/ULIRGS and do not comprehensively explore other loss processes and their possible relationship with galactic astrophysical properties. With our improved spectral sampling and a sample of more nearby, resolvable lower SFR SFGs, one of the primary goals of this paper is to explore how the shapes of radio SEDs and their parameters are related to a galaxy’s global astrophysical properties. This will allow us to infer the dominant physical processes, including emission processes and cooling in the ISM, for a diverse sample of SFGs.

In this study we select a sample of twenty nearby $(z \lt 0.04)$ SFGs including twelve with LFTOs and eight controls to investigate whether their global astrophysical properties differ. We perform detailed radio SED modelling by constructing a series of increasingly complex modular, radio emission models including loss and absorption processes. We then compare the model results to a number of global properties to investigate how they are connected to the radio SED features observed and the emission and loss processes occurring within SFGs. In Sections 2 and 3, we describe the sample selection and data acquisition. Section 4 details the SED model construction, fitting, and selection. Sections 5 and 6 present the results and discussion of these results. Lastly Section 7 presents our conclusions. Throughout this paper we assume a Hubble constant of 70 km s^-1 Mpc^-1 (h = 0.70), and matter and cosmological constant density parameters of $\Omega_{M}$ = 0.3 and $\Omega_{\Lambda}$ = 0.7.

2. Sample and ancillary data

2.1. Sample selection

We initially select SFGs from the GLEAM-6dFGS catalogue (Franzen et al. 2021) of $z\leq 0.1$ sources selected at 200 MHz (1590 total sources). Galaxies are categorised as SF based on their optical spectrum containing H $\alpha$ and H $\beta$ emission lines typical of HII regions (427 SFGs or 26.9%). The candidate LFTO galaxies are selected based on the measured spectral index between 76–227 MHz ( $\alpha_\textrm{L} \geq -0.2$ , 15 sources) which are then visually inspected and confirmed to exhibit flattening/peaks in their radio SED. Visual inspection is necessary due to the large uncertainties in flux density for some GLEAM photometry in a few sources. GLEAM sources with any negative measured sub-band flux densities were not fitted with spectral indices but are also considered upon visual inspection (54 sources with no measured $\alpha_\textrm{L}$ , 27 of which have $\alpha^\textrm{GLEAM}_\textrm{low} \geq -0.2$ ). The candidate LFTO galaxies were then further selected with declination < $-14^{\circ}$ to minimise elongation of the synthesised beam during Australian Telescope Compact Array (ATCA; Frater, Brooks, & Whiteoak 1992; Wilson et al. 2011) observations and to maintain comparable beamshapes to Widefield Infrared Survey Explorer (WISE; Wright et al. 2010) observations. Using the high-frequency spectral index ( $\alpha_\textrm{H}$ ) between 227–843/1 400 MHz (depending on whether Sydney University Molonglo Sky Survey (SUMSS; Mauch et al. 2003) or NRAO VLA Sky Survey (NVSS; Condon et al. 1998) observations are available) we estimated 9.5 GHz flux density and performed a flux cut of $S_\nu=~$ 8 mJy to select brighter SFGs allowing them to be observed by ATCA and ensure sensitivity within the requested observing time. Further LFTO sample restriction was then performed by limiting $z \leq 0.04$ . Lastly visual inspection of the RACS 888 MHz and NVSS/SUMSS emission for the remaining LFTO candidates was performed to remove sources which were confused in the large ( $\sim+2'$ ) GLEAM beam which left a final sample of 11 LFTO SFGs.

A control sample of eight sources is selected first by limiting the low-frequency spectral index to steeper values ( $\alpha_\textrm{L} \leq -0.5$ ). The distribution of their redshifts, K-band absolute magnitude and estimated 9.5 GHz flux density are then limited in the same fashion as to the LFTO sample and their distributions matched. Control galaxies are then chosen to be within 10° of LFTO galaxies to allow the use of the same phase calibrator for multiple sources and reduce overheads during ATCA observations. We then perform the same visual inspection to remove GLEAM confused sources. The complete SFG sample including common parameters and a flag for whether they are LFTO or control galaxies are given in Table 1.

Table 1. SFG sample properties.

Note: Column (1): GLEAM source ID. Column (2): Common ID from NED. Column (3): SFG sample membership. Column (4-5): Right ascension and declination (J2000) of the source from the WXSC (Jarrett et al. 2019). Column (6): Redshift from GLEAM-6dFGS as in Franzen et al. (2021). Column (7): Luminosity distance from the WXSC (Jarrett et al. 2019). Column (8): K-band absolute magnitude from Franzen et al. (2021). Column (9): Stellar mass from the WXSC calculated using the light and colours method (Jarrett et al. 2023). Column (10): Estimated source inclination calculated using K-band light axis ratios from the Two Micron All-Sky Survey (Jarrett et al. 2000). Column (11): Lyon-Meudon Extragalactic Database (LEDA; Paturel et al. 2003) morphology.

^a K-band light axis ratios were unavailable so light axis ratios from the Morphological catalogue of galaxies (Vorontsov-Vel’Yaminov & Arkhipova 1974) were used.

^b K-band light axis ratios were unavailable so light axis ratios from the ESO-Uppsala galaxies catalogue (Lauberts & Valentijn 1989) were used.

^cMerging systems; inclination estimates are unreliable.

2.2. Radio

We make use of the radio data included in the GLEAM-6dFGS SFG catalogue (Franzen et al. 2021) which includes 20 flux measurements in the GLEAM frequency bands (72–232 MHz) with channel widths of 7.68 MHz. A declination dependent SUMSS (Mauch et al. 2003, 843 MHz;) or NVSS (Condon et al. 1998, 1.4 GHz;) measurement is also included. The GLEAM-6dFGS SFG catalogue (Franzen et al. 2021) also includes a measurement of the GLEAM spectral index as measured in Hurley-Walker (2017) by linear least squares fitting a single power law spectrum to sources with no missing or negative flux density values. We measure the GLEAM spectral index for all sources after removal of the negative flux values using the same method as Hurley-Walker (2017). This results in three LFTO sources with measured GLEAM spectral indices $\alpha^\textrm{GLEAM}_\textrm{ low} \leq$ -0.2 primarily due to large GLEAM flux density uncertainties at low frequencies. The new $\alpha^\textrm{GLEAM}_\textrm{low}$ values for the remaining sources agree with the values in Hurley-Walker (2017) and new values are marked in Table A1. The GLEAM to SUMSS/NVSS spectral index is also updated using RACS-mid flux densities to improve consistency and labelled $\alpha^\textrm{RACS}_\textrm{mid}$ .

We create new 31-MHz bandwidth flux densities by combining the flux densities in five sets of four adjacent sub-bands. These are determined by using a weighted average:

(1) \begin{equation} S_{\nu} = \frac{\Sigma_{i}\;S_{\nu,i}/\sigma_{i}^{2}}{\Sigma_{i}\;1/\sigma_{i}^{2}},\end{equation}

where S $_{\nu,i}$ are the flux measurements of the individual 7.68 MHz sub-bands with uncertainties $\sigma_{i}$ . This is done to increase model fitting accuracy and convergence times as a number of individual sub-bands contain significant outliers including negative flux densities. It also avoids the need for the SED fitting to deal with the correlated noise in each set of four adjacent sub-bands which were deconvolved at the same time.

Flux densities and images from the Rapid Australian SKA Pathfinder (ASKAP) Continuum Survey (RACS; McConnell et al. 2016), a survey performed by the ASKAP currently covering dec. $ \lt +30^{\circ}$ at 888 MHz (Hale et al. 2021) and 1.37 GHz (Duchesne et al. 2024), are also utilised. Lastly we collated any supplementary radio data in the NASA Extragalactic Database (NED ^a ). These flux values are presented in Table A1.

2.3. Infrared

FIR flux density measurements at 60 and 100 $\mu$ m from Infrared Astronomy Satellite (IRAS Neugebauer et al. 1984) are collated from the IRAS Faint Source Catalogue (Moshir et al. 1990).

Although the GLEAM-6dFGS catalogue also includes WISE source catalogue measurements (Wright et al. 2010) at 3.4, 4.6, 12, and 23 $\mu$ m (WISE band W1, W2, W3, and W4), due to the morphological complexity of many of the sources in our sample we choose to instead use values from the WISE eXtended Source Catalogue (WXSC Jarrett et al. 2019). The WXSC more accurately encapsulates the emission from extended sources and removes contamination from foreground or background sources. For blended or merging systems the WXSC flux density and Vega magnitude values are measured for the entire system such that they are comparable to the radio continuum observations which have a larger beam. These WXSC flux density values are used to derive the mid-IR based global parameters including the mid-IR SFR (SFR $_{mircor}$ ; Cluver et al. 2024), stellar mass (Jarrett et al. 2023), and specific SFR (sSFR $_{mircor}$ ). SFRs are from Cluver et al. (2024) that updates the calibration in Cluver et al. (2017) which related WISE W3 and W4 to L $_{TIR}$ (total infrared luminosity). SFRs are determined using an invariance-weighting of SFR $_{12}$ and SFR $_{23}$ (WISE band W3 and band W4 based SFR) and include a correction to account for star formation in low mass galaxies with relatively low dust content (which makes use of the relationship between dust density and UV to IR emission). The WXSC also provides the luminosity distance D $_L$ for each source. The IR flux density measurements shown in Table B1 and global derived parameters in Table 8.

2.4. Ancillary data

We make use of optical background images to present the stellar distribution of the galaxies in our sample. These optical images are collated from the DataCentral Data Aggregation Service (Miszalski et al. 2022) ^b and are primarily g-band images. Images from the Panoramic Survey Telescope and Rapid Response System (PanSTARRS) data release 2 (PS1-DR2; Flewelling et al. 2020), the Dark Energy Survey (DES; Abbott et al. 2021), the DESI Legacy Survey (Dey et al. 2019), SkyMapper Southern Survey (Onken et al. 2019) and VISTA Hemisphere Survey (McMahon et al. 2013, where we used K-band where g-band images are not available) are presented in Fig. A1. Estimated source inclinations were calculated using K-band light axis ratios from the Two Micron All-Sky Survey (Jarrett et al. 2000).

3. ATCA observations and data processing

3.1. ATCA observation details

ATCA observations were performed across three 24 h periods and two 14 h periods (project ID C3483) using four different ATCA configurations and are summarised in Table 2. Different configurations are chosen to obtain adequate (u,v) coverage and total integration time for accurate flux measurements for SED modelling as well as high ( $\sim$ 6”) resolution observations at 5.5 and 9.5 GHz.

Table 2. ATCA observation details.

Note: Column (1): universal time date. Column (2): ATCA configuration. Column(3): Number of antennas used to observe. Column (4): Integration time per u-v scan on source. Column (5): Average total integration time on source. Column (6-7): Central frequencies of ATCA observations. Column (8): Bandwidths at each central frequency. Column (9): The largest angular scale probed at IF2

^aDue to scheduling and source rise and set times some sources may have more or fewer scans resulting in total integration times varying.

^bDue to radio frequency interference the 18.8 GHz observations were discarded.

3.2. ATCA data processing

ATCA data were calibrated and imaged using MIRIAD (Sault, Teuben, & Wright 1995) software package. The PGFLAG and BLFLAG flagging routines within MIRIAD were used for automated and manual flagging, respectively, in conjunction with the traditional MIRIAD calibration tasks to perform the initial data reduction. A frequency dependent calibration solution was determined using the NFBIN option given the wide bandwidth of the Compact Array Broadband Backend (CABB) system (Wilson et al. 2011). Flux and bandpass calibration were performed on PKS 1934-638 whilst gain calibration was transferred from the phase calibrators. The calibration solutions for each individual day of observations were applied to each source. The calibrated observations different configurations at 5.5 and 9.5 GHz were then merged to provide the best u-v sampling and highest sensitivity.

3.2.1. Total flux density measurements

To measure the total flux density across all frequency bands, each source was imaged using their complete bandwidth ( $\Delta\nu = 2.048$ GHz, minus the 100 edge channels automatically flagged by ATLOD). Natural weighting (Briggs robust parameter value of 2) was used to provide the maximum signal-to-noise for total flux measurements. MFCLEAN (Sault & Wieringa 1994) was used to deconvolve the multi-frequency synthesised dirty map. Then RESTOR and LINMOS were used to restore the full bandwidth images and perform the primary beam correction whilst accounting for the spectral index of the clean components.

Figure 1. Left: The preferred model SED of GLEAM J012121-340345 with observed data points. The overlaid black line indicates the full model whilst the dotted blue line indicates the first PL component and purple dashed line indicates the second PL component which is free-free absorbed. The highlighted regions represent the 1- $\sigma$ uncertainties sampled by EMCEE. Right: The DES g-band optical image of GLEAM J012121-340345 showing the stellar extent and morphology overlaid with contours from RACS-mid at 1.37 GHz in blue and ATCA 9.5 GHz in green. Radio contours for both frequencies start at the 4 $\sigma$ level and increase by factors of $\sqrt{3}$ . The FWHM beams for RACS-mid and ATCA are given by the blue and green ellipses, respectively. The scale bar at the bottom left denotes 5 kpc.

Figure 2. Left: The preferred model SED of GLEAM J184747-602054 with observed data points. The overlaid black line indicates the full FFA_PL model. The highlighted region represents the 1- $\sigma$ uncertainties sampled by EMCEE. Right: The g-band optical image of GLEAM J184747-602054 showing the stellar extent and morphology overlaid with contours from RACS-mid at 1.37 GHz in red and ATCA 9.5 GHz in pink. Radio contours for both frequencies start at the 4 $\sigma$ level and increase by factors of $\sqrt{3}$ . The FWHM beams for RACS-mid and ATCA are given by the red and pink ellipses, respectively. The scale bar at the bottom left denotes 5 kpc.

We used an iterative procedure to perform sub-band splitting to allow for improved spectral sampling and modelling. Each CABB band was initially imaged as described above. The peak flux density was extracted and then, provided the signal-to-noise ratio (SNR) was above 8 or 10, the CABB dataset would be split into 2 or 3 sub-bands, respectively, and reprocessed. To safely use MFCLEAN, we ensured that each sub-band maintained a fractional bandwidth larger than 10 $\%$ . This resulted in a maximum of three sub-bands at 5.5 GHz, two at 9.5 GHz and no sub-band splitting at $ \gt +10$ GHz. Sub-bands were split using the LINE option in INVERT with their new central frequencies recorded individually for each source depending on their unflagged channels. Total fluxes for each image were then measured using PROFOUND task within the ProFound software package (Robotham et al. 2018).

High resolution imaging was also separately performed at 5.5 and 9.5 GHz using the entire CABB bandwidth. A Briggs robust parameter of 0.5 was chosen to achieve a balance between sensitivity to diffuse emission and angular resolution. The radio contours for the high resolution images are presented in Figs. 1, 2, and A1.

In order to account for flux scale errors between GLEAM and the other radio data, we applied a conservative 10% error in quadrature to all non-GLEAM flux density measurements. This 10% error was chosen such that it encompasses the $\approx+8\%$ flux scale error measured for the GLEAM survey (Hurley-Walker 2017) and accounts for possible source variability, although SFGs generally do not show significant variability at radio frequencies on timescales of decades at the sensitivities observed in our sample (Mooley et al. 2016). Lastly we do not attempt to match the u-v coverage or resolution of multi-frequency images used when modelling the radio SEDs due to the large range in resolutions between GLEAM observations ( $\sim+2'$ ) and other radio data ( $\sim$ 10–45”). ATCA observations are performed in a number of configurations to ensure we have sufficient short baseline u-v coverage for extended flux measurements. Measured ATCA flux densities are given in Table A2 with other radio fluxes presented in Table A1.

4. SED modelling

4.1. Modular radio continuum models

We choose to fit a series of increasingly complex models to our sources loosely following Galvin et al. (2018) with all modelling being performed in the frequency rest-frame ^c with a reference frequency of $\nu_{0}$ = 1.4 GHz. We deviate from the method of Galvin et al. (2018) such that radio continuum models are constructed using a modular approach whereby base single and two component power-law models are modified by prefix models at GLEAM frequencies ( $\nu \leq$ 300 MHz) and suffix models non-GLEAM frequencies ( $\nu \geq$ 300 MHz). The prefix and suffix models encapsulate the physical loss and absorption processes occurring within SFGs including FFA, IC losses and synchrotron losses. These models are outlined in Table 3.

We also make no attempt in this work to include a FIR dust heating component in our current SED modelling. FIR emission increases model complexity and does not help to constrain our radio SED model parameters, as, at our current highest radio frequency measurements (17 GHz) warm dust contributes $\ll+1\%$ to the total flux.

4.1.1. Base model: Power Law (PL)

Firstly we fit a simple single component power law to radio flux density measurements with the form of:

(2) \begin{equation} S_{\nu} = A\bigg(\frac{\nu}{\nu_{0}}\bigg)^{\alpha},\end{equation}

where A is a normalisation component and the spectral index, $\alpha$ , is the gradient in logarithmic space with both being treated as free parameters.

4.1.2. Base model: Synchrotron and free-free emission (SFG)

The radio continuum emission is usually modelled as the sum of two distinct power laws with one representing the flat spectrum thermal free-free emission and the second representing the steep spectrum non-thermal synchrotron emission. This model takes the form:

(3) \begin{equation} S_{\nu} = A\bigg(\frac{\nu}{\nu_{0}}\bigg)^{\alpha} + B\bigg(\frac{\nu}{\nu_{0}}\bigg)^{-0.1},\end{equation}

where the free parameters A and B are the synchrotron and free-free normalisation components, respectively. The free parameter $\alpha$ is the synchrotron spectral index which depends on the cosmic ray injection history and is known to vary (Niklas et al. 1997). The spectral index of the free-free emission is well approximated by $-0.1$ over the range of interest (Condon 1992).

Table 3. Modular radio continuum models.

Note: Column (1-3): The prefix, base and suffix model names. Column (4): The full model label. Column (5): The function representing the modelled radio spectrum. Column (6): Free parameters. Column (7): Galvin et al. (2018) equivalent model.

4.1.3. Prefix: Free-free Absorption (FFA_)

Synchrotron emission can be attenuated by FFA processes when in a coextensive environment with free-free emission producing spectral curvature at primarily low frequencies. The attenuation is influenced by the density, flux density and spatial distribution of the ionised free-free emission in comparison to the synchrotron emission. The free-free optical depth can be approximated as $\tau_1$ = $(\nu/\nu_{t,1})^{-2.1}$ with $\nu_{t,1}$ being the turn-over frequency where the optical depth reaches unity (see Section 5.1.1). This model modifies the SFG base model resulting in the full model FFA_SFG:

(4) \begin{equation} S_{\nu} = (1-e^{-\tau_1})\bigg(B + A\bigg(\frac{\nu}{\nu_{t,1}}\bigg)^{0.1+\alpha}\bigg)\bigg(\frac{\nu}{\nu_{t,1}}\bigg)^{2},\end{equation}

following Condon (1992) and Clemens et al. (2010). $\nu_{t,1}$ is limited to GLEAM frequencies (i.e. $\leq$ 300 MHz). A, B, and $\alpha$ are the synchrotron and free-free normalisation components and synchrotron spectral index, respectively, which are fitted simultaneously with $\nu_{t,1}$ . Model degeneracy is minimised by replacing the $\nu_{0}$ term with the turnover frequency parameter for each component following Galvin et al. (2018). Setting $B = 0$ gives the full model FFA_PL.

4.1.4. Suffix: Free-free Absorption (_FFA)

The models only including the FFA_ prefix assume a single volume of thermal free-free plasma mixed with relativistic electrons which produce synchrotron emission. This model was derived from observations by Condon & Yin (1990) of the irregular clumpy galaxy Markarian 325; however, Clemens et al. (2010), Galvin et al. (2018), and Dey et al. (2022); Dey et al. (2024) present a set of LIRGS/ULIRGS which show higher frequency ‘kinks’ in their radio SEDs which could also be attributed to FFA. Their interpretations suggest that when multiple SF regions of different geometric orientations or composition are integrated over a large synthesised beam the observed radio continuum can be complex. Thus following Galvin et al. (2018) we introduce a model which includes a single relativistic electron population which produces synchrotron emission that is inhomogenously mixed with two distinct star-forming regions of two distinct optical depths of which only $\tau_2$ becomes optically thick within the observed frequency range. This situation is described by the full model SFG_FFA:

(5) \begin{align} S_{\nu} &= A\bigg(\frac{\nu}{\nu_{0}}\bigg)^{\alpha} + B\bigg(\frac{\nu}{\nu_{0}}\bigg)^{-0.1}\nonumber\\[5pt] &+ (1-e^{-\tau_{2}})\bigg(D+C\bigg(\frac{\nu}{\nu_{t,2}}\bigg)^{0.1+\alpha}\bigg)\bigg(\frac{\nu}{\nu_{t,2}}\bigg)^{2},\end{align}

where the second component of the SFG_FFA model represents the second distinct star-forming region which becomes optically thick below a certain frequency. $\nu_0$ is the reference frequency of 1.4 GHz and $\tau_2$ describes the optical depth of the second component parameterised by the turnover frequency $\nu_{t,2}$ which is limited to non-GLEAM frequencies ( $\geq$ 300 MHz). A and C are the normalisation parameters for the synchrotron emission and B and D govern the free-free emission component. $\alpha$ is the spectral index of a single synchrotron emission population. Setting B = 0 and D = 0 gives the full model PL_FFA.

We can replace the first line in this model with Equation (4) to account for the situation where there is also a LFTO caused by FFA at low frequencies resulting in the FFA_SFG_FFA and (when B = 0 and D = 0) FFA_PL_FFA full models.

We then allow for two distinct electron populations by relaxing the single spectral index constraint. This is physically motivated by galaxy merger situations which can trigger new bursts of SF such that the electron distribution can be made up of a populations of newly injected and older non-thermal relativistic electrons. This is done by allowing each component in the FFA_SFG_FFA full model to have different synchrotron spectral index values $\alpha$ and $\alpha_{2}$ which are free parameters in the full model FFA_SFG_FFA2.

4.1.5. Suffix: Synchrotron and Inverse-Compton Losses (_SIC)

The spectral steepening observed at higher frequencies may not be related to FFA at different optical depths and instead be due to synchrotron or IC losses. Synchrotron losses take place as CREs age and lose energy as they propagate in a galaxy’s large scale magnetic field whilst IC losses are generally dependent on FIR or CMB photon scattering. Synchrotron and IC losses under constant electron injection from massive SF both act to gradually steepen the synchrotron spectral index by $\Delta\alpha = -0.5$ around a ‘break’ frequency $\nu_b$ . This results in the full model SFG_SIC having the form:

(6) \begin{equation} S_{\nu} = A\bigg(\frac{\nu}{\nu_{0}}\bigg)^{\alpha}\bigg(\frac{1}{1+(\nu/\nu_b)^{\Delta\alpha}}\bigg) + B\bigg(\frac{\nu}{\nu_{0}}\bigg)^{-0.1},\end{equation}

where the free parameters A and B are the synchrotron and free-free normalisation components, respectively. The free parameter $\alpha$ is the low-frequency synchrotron spectral index which gradually steepens by $\Delta\alpha$ = -0.5 around the ‘break’ frequency $\nu_b$ which is limited to non-GLEAM frequencies ( $\geq$ 300 MHz). Setting B = 0 gives the full model PL_SIC.

We also model the situation where there is a LFTO caused by FFA and gradual spectral steepening caused by synchrotron and IC losses in the FFA_SFG_SIC full model which takes the form:

(7) \begin{equation} S_{\nu} = (1-e^{-\tau_1})\bigg(B + A\bigg(\frac{\nu}{\nu_{t,1}}\bigg)^{0.1+\alpha}\bigg(\frac{1}{1+(\nu/\nu_b)^{\Delta\alpha}}\bigg)\bigg)\bigg(\frac{\nu}{\nu_{t,1}}\bigg)^{2},\end{equation}

where the free parameters are the same as in Equation (6) with the addition that $\nu_{t,1}$ is the turnover frequency where the optical depth ( $\tau_1$ ) reaches unity and is limited to GLEAM frequencies ( $\leq$ 300 MHz). Setting B = 0 gives the full model FFA_PL_SIC.

We do not make any attempt to model both the Suffix models simultaneously (i.e. two synchrotron components, one with a spectral break and both with FFA) as we do not have the spectral sampling frequency or range to be able to draw accurate conclusions. It is likely, however, that both FFA and loss processes are occurring in most SFGs and that the two separate Suffix models will improve the fitting of radio SEDs which display an obvious ‘kink’ and gradual spectral steepening at higher frequencies.

4.2. Fitting and selection

4.2.1. Model fitting

Briefly we use the ‘affine invariant’ Markov chain Monte Carlo ensemble sampler (Goodman & Weare 2010) implemented as the EMCEE PYTHON package (Foreman-Mackey et al. 2013) to constrain each of the radio continuum models for each source in our sample. The log likelihood function that EMCEE attempts to minimise relies on the assumption that measurements are independent with normally distributed errors.

4.2.2. Model priors

We choose physically motivated uninformative (uniform) priors to constrain our models within a Bayesian framework. Throughout our model fitting we ensure that the normalisation parameters A, B, C, and D remain positive (A/B/C/D > 0). The spectral index parameters $\alpha$ and $\alpha_{2}$ remain in the range of $-1.8 \leq \alpha \leq -0.2$ . The LFTO ( $\nu_{t,1}$ ) and spectral ‘kink’ ( $\nu_{t,2}$ ) frequencies are limited to between 10 MHz to 300 MHz and 300 MHz to 17 GHz, respectively. The spectral ‘break’ frequency ( $\nu_{b}$ ) is also limited to between 300 MHz and 17 GHz. These priors are founded on the sound assumptions that flux densities are positive emission processes and we can only constrain turnovers within the frequency range that we have data; however, we note that some SEDs may begin to flatten before the optical depth reaches unity. The spectral index parameters are limited to allow for the range of values found for synchrotron dominated emission in literature (Condon & Yin 1990; Clemens et al. 2010; Galvin et al. 2018; An et al. 2021; Dey et al. 2022).

4.2.3. Model selection

In order to objectively test whether the introduction of additional model complexity is justified by an improved fit and is not just a symptom of overfitting we make use of an estimate of the evidence value Z. The evidence value is defined as the ratio between the integral of the posterior volume over the prior volume (Skilling 2004) and is computationally difficult to compute but can be reliably estimated using recent algorithms. DYNESTY (Speagle 2020) uses a dynamic nested sampling method (Higson et al. 2019) to obtain an estimate of the evidence value. Given the evidence values of competing models, one is able to determine whether a model is preferred over another for a given set of data. The natural logarithm of the Bayes odds ratio between evidence values $Z_{1}$ and $Z_{2}$ for models $M_{1}$ and $M_{2}$ is described by:

(8) \begin{equation} \textrm{ln}(\Delta Z) = \textrm{ln}(Z_1) - \textrm{ln}(Z_2).\end{equation}

A value of $\textrm{ln}(\Delta Z)$ < -5 provides very strong evidence for $M_{1}$ over $M_{2}$ whilst -5 < $\textrm{ ln}(\Delta Z)$ < -3 and -3 < $\textrm{ln}(\Delta Z)$ < -1.1 provide strong and positive evidence, respectively. When $\textrm{ ln}(\Delta Z)$ > -1.1 the models are indistinguishable from each other. This scale was established by Kass & Raftery (1995) and is considered the standard for model selection.

The prior parameter space searched by DYNESTY is limited to a uniform distribution within the uncertainties given by the 1st and 99th percentiles of the samples posterior distribution as determined by EMCEE. This limitation of priors is necessary as DYNESTY requires bounded priors and different types of models are explored. Because the evidence value is entirely dependent on the ‘size’ of the prior volume (Skilling 2004) setting arbitrarily large priors on normalisation components would heavily bias the evidence values against models with extra normalisation parameters. Within DYNESTY we also select ‘Single ellipsoid bounds’ and ‘random walk’ samplers (rather than uniform samplers) to improve processing time and maintain consistency with the posteriors estimated by EMCEE. The results of the Bayes odds ratio tests for all models as defined relative to the lowest value are summarised in Table 4 with the most preferred model and any other competing models highlighted, that is, $ln(\Delta Z)$ > -1.1.

Table 4. Bayes odds ratio table.

An overview of the natural log of the Bayes odds ratio from the DYNESTY fitting of each model to every source. For each source, the values presented below are the evidence values for each model divided by the most preferred model. The natural log of the ratio is presented, such that the most preferred models have values in this table equal to ln(1) $=0$ (highlighted green). Models that are indistinguishable from the most preferred model correspond to $ \lt $ ln(3) $=1.1$ (highlighted grey). Less preferred models therefore have more negative numbers. Initial SFG sample membership is denoted by LFTO or control. Control galaxies with favoured models that do not include an LFTO have their sample membership highlighted in blue. LFTO galaxies that with FFA_ prefix models are shown in red. LFTO galaxies which do not have their most favoured model with the FFA_ prefix (i.e. the favoured model does not include a LFTO) are shown in pink. An * indicates a competing model of the opposite class.

Table 5. Preferred model parameter table.

The most preferred models, as judged by their evidence value with their constrained parameter values and 1 $\sigma$ uncertainties. We use the 50th percentile of the samples posterior distribution as the nominal value, and use the 16th and 84th percentiles to provide the 1 $\sigma$ uncertainties. Parameters not included in the models are left blank.

^aBy construction the final values of $\alpha$ for this model are $\alpha + \Delta\alpha$ (i.e the steeper value at higher frequencies). In brackets we include the value of $\alpha$ at the break frequency $\nu_b$ .

5. Results

5.1. Model results

We find that all eight galaxies selected to be controls are favoured to be modelled with no LFTO (highlighted blue in Table 4) and no control galaxies have models with the prefix FFA_ as being indistinguishable suggesting that their selection criteria were appropriate. Of these eight control galaxies three are best modelled by simple PL models, four by PL_FFA models and one by a PL_SIC model. There is, however, competition between primarily the PL and PL_SIC models with many being indistinguishable from each other. The PL_SIC models often provide a better fit to the data by accounting for both flattening and steepening of the SED; however, the addition of an extra free parameter which is difficult to constrain ( $\sigma(\nu_{b}) \sim5$ GHz) results in PL models being more favoured as selected purely by their log-likelihood values. It is possible that models which include synchrotron or IC losses are more physically motivated at the expense of an additional parameter. We discuss the implications of this in the discussion. PL_FFA models do not have any competing models in the four sources they are preferred for. Indeed they do all visually have distinctly offset power law emission of a similar spectral index suggesting either distinct components as suggested in Dey et al. (2024) or perhaps GLEAM flux scale discrepancies.

Six of the 11 galaxies in the LFTO sample have models with LFTOs (FFA_ prefix, highlighted in red in Table 4) of which only four do not have competing simpler models with no curvature. Despite their initial selection as sources with LFTOs five sources have no FFA_ prefix in their most favoured models (highlighted in pink). It is often the case that the LFTO galaxies are better fit by FFA_ prefix models; however, this is for sources with smaller GLEAM flux density uncertainties or more spectral curvature. Simpler models compete well due to the large impact of adding free parameters. We find that GLEAM_J072121-690005 is the only source to have a competing full model with both a prefix and suffix components, this source has a distinct LFTO and steepening spectral index as shown in Fig. A1.

We find our different models are more similar (i.e. the Bayes odds ratios are smaller) than those in Galvin et al. (2018) and Dey et al. (2022); Dey et al. (2024). This is primarily a result of model construction such that our prior bounds are smaller with the overall frequency space being a factor of 2.5 times smaller which is then further constrained by the prefix and suffix model conditions only allowing $\nu_{t,1}$ to take values between 10–300 MHz and 300 MHz $\leq$ $\nu_{t,2}$ or $\nu_{b}$ $\leq$ 17 GHz. In the case of Galvin et al. (2018) we have fewer free parameters due to their inclusion of FIR dust and GLEAM covariance models. We are also lacking data above 17 GHz which helps constrain and distinguish models, especially those with thermal emission components (Galvin et al. 2018; Dey et al. 2022; Dey et al. 2024).

All favoured models do not include thermal emission, this is likely due to the lack of high-frequency data needed to constrain the thermal emission component. GLEAM_J184747-602054 is the only source that has a competitive SFG based model suggesting that it may be the only source in our sample that has a significant thermal fraction at 17 GHz. We do not observe spectral flattening in our other sources towards higher frequencies and generally find the modelled spectral indices for our SFG sample agree with the commonly accepted range of $-0.8 \geq \alpha \leq -0.7$ (Condon 1992). These results are summarised in Table 4 and Figs. 1, 2, and A1.

5.1.1. Emission measure

The turnover frequency due to FFA is dependent on where the optical depth reaches unity which occurs at the turnover frequency ( $\nu_{t,1}$ or $\nu_{t,2}$ ) measured during SED modelling. The emission measure (EM) is then calculated assuming that the emitting HII regions form a cylinder orientated along of line of sight with constant temperature and density (Condon 1992). The free–free opacity is then well approximated by:

(9) \begin{equation} \tau_{\nu} = 3.28\times10^{-7} \bigg(\frac{T_e}{10^4\,\textrm{K}}\bigg)\bigg(\frac{\nu}{\textrm{ GHz}}\bigg)^{-2.1}\bigg(\frac{\textrm{EM}}{\mathrm{pc\;cm^{-6}}}\bigg),\end{equation}

where $T_e$ is the electron temperature of the HII emitting region, typically taken as $10^4$ K, and EM is the emission measure at a depth s, defined as:

(10) \begin{equation} \frac{\textrm{EM}}{\mathrm{pc\;cm^{-6}}} = \int_{los} \bigg(\frac{n_e}{\mathrm{cm^{-3}}}\bigg)^2d\bigg(\frac{s}{\textrm{ pc}}\bigg).\end{equation}

Where the EM is the integral of the electron density, $n_e$ , along the line of sight, los, of a HII region of depth s. Using the turnovers constrained by our modelling, we have estimated the EMs of our sources, outlined in Table 6, using Equation (9). We label the corresponding EM of $\nu_{t,1}$ and $\nu_{t,2}$ for all models as EM $_1$ and EM $_2$ , respectively. We find that the EM values obtained for our sources are consistent with the more luminous LIRG and ULIRG samples of Clemens et al. (2010), Galvin et al. (2016), Dey et al. (2022); Dey et al. (2024) suggesting that the depth and orientation of the HII regions within these galaxies are similar despite their lower starburst activity.

Table 6. Emission measures.

An overview of the emission measures (EM) derived for each source from the model most supported by the evidence. Objects without an emission measure constrained are not listed.

5.2. Global radio properties

5.2.1. Radio star-formation rates

RACS-mid flux densities at 1.37 GHz from the 25^′′ catalogue (Duchesne et al. 2024) are used to estimate the 1.4 GHz radio-SFR as they are available for all of our SFGs at a higher resolution than NVSS and its flux scale is consistent with that of NVSS within $ \lt $ 1%. We use the relationship from Molnár et al. (2021):

(11) \begin{equation} \log\bigg(\frac{\mathrm{SFR_\textrm{1.4 GHz}}}{\mathrm{M_{\odot}\,yr^{-1}}}\bigg)= (0.823\pm0.009)\log\bigg(\frac{L_{1.4}}{\mathrm{W\,Hz^{-1}}}\bigg) - (17.5\pm0.2)\end{equation}

where $L_{1.4}$ is the 1.4 GHz radio luminosity in W Hz^-1 and is calibrated against the $q_\textrm{TIR}$ derived SFR to estimate the radio SFR. The 1.4 GHz luminosity is calculated using:

(12) \begin{equation} \bigg(\frac{\mathrm{L_\textrm{1.4 GHz}}}{\mathrm{ W\,Hz^{-1}}}\bigg)=9.52\times10^{15}\bigg(\frac{4\pi}{(1+z)^{1+\alpha}}\bigg)\,\bigg(\frac{\mathrm{D_{L}}}{\mathrm{ Mpc}}\bigg)^{2}\bigg(\frac{S_\textrm{1.4 GHz}}{\textrm{mJy}}\bigg)\end{equation}

where z is the redshift, $\alpha$ is the modelled spectral index, D $_\textrm{L}$ is the luminosity distance from the WXSC in Mpc and $S_\textrm{1.4 GHz}$ is the RACS-mid integrated flux density in mJy.

The radio SFR surface density is given by:

(13) \begin{equation} \frac{\Sigma \textrm{SFR}_\textrm{1.4 GHz}}{\mathrm{M_{\odot}\,yr^{-1}\,kpc^{-2}}}= \bigg(\frac{\mathrm{SFR_\textrm{ 1.4 GHz}}}{\mathrm{M_{\odot}\,yr^{-1}}}\bigg)\bigg(\frac{\pi k^2 A B}{\mathrm{kpc^2}}\bigg)^{-1}\end{equation}

where k is the conversion factor from angular to physical scale in kpc arcsec^-1 as calculated based on D $_L$ and A and B are the K-band light major and minor axis radii in arcseconds. All these values are presented in Table 7.

Table 7. Derived radio continuum properties.

Note: Column (1): GLEAM source ID. Column (2): GLEAM (72–231 MHz) spectral index. Column (3): GLEAM-RACS-mid (200 MHz to 1.37 GHz) spectral index. Column (4): ATCA (5.5–17 GHz) spectral index. Column (5): RACS-mid luminosity. Column (6): Total radio SFR from RACS-mid luminosity and Equation (11). Column (7): Total radio star formation rate surface density calculated based on radio SFR and K-band light axis ratios (see Equation 13).

^aRecalculated values of $\alpha^\textrm{GLEAM}_\textrm{low}$ due to a GLEAM sub-band containing negative flux.

5.2.2. Spectral indices

To determine the relationship between modelled spectral index $\alpha$ (hereafter $\alpha^{model}$ ) and the GLEAM, GLEAM to RACS-mid, and ATCA spectral indices ( $\alpha^\textrm{GLEAM}_\textrm{low}$ , $\alpha^\textrm{RACS}_\textrm{mid}$ , and $\alpha^\textrm{ATCA}_\textrm{high}$ , respectively) for our samples we perform a Spearman’s rank correlation test. We find that $\alpha^{model}$ is most strongly correlated with $\alpha^\textrm{ATCA}_\textrm{high}$ and slightly correlated with $\alpha^\textrm{RACS}_\textrm{mid}$ as based on their p-values which reject the null hypothesis of no correlation. This is likely due to the lower error budget and increased sampling weighting the modelled spectral index towards $\alpha^\textrm{ ATCA}_\textrm{high}$ in comparison to $\alpha^\textrm{RACS}_\textrm{mid}$ .

Comparing our modelled spectral index (as presented in Table 5) to $\alpha^\textrm{GLEAM}_\textrm{low}$ , $\alpha^\textrm{RACS}_\textrm{mid}$ , and $\alpha^\textrm{ATCA}_\textrm{high}$ at GLEAM, RACS-mid and ATCA frequencies in Fig. 3 we see positive correlation between the modelled spectral index and $\alpha^\textrm{RACS}_\textrm{mid}$ and $\alpha^\textrm{ATCA}_\textrm{high}$ . The modelled spectral index tends to be steeper than $\alpha^\textrm{RACS}_\textrm{mid}$ but flatter than $\alpha^\textrm{ATCA}_\textrm{high}$ . By construction the modelled spectral index will differ based on the most favoured model as they each account for spectral flattening or steepening differently to simple unmodified power-law models. The overall effect of this is that basic PL only models have flatter spectral indices as they do not account for low-frequency flattening or include a second emission component as in PL_FFA models.

Figure 3. Comparisons between the modelled spectral index and GLEAM, GLEAM to RACS-mid and ATCA spectral indices in panels (a), (b), and (c), respectively. The slope of the weighted linear fit and its 1 $\sigma$ uncertainty and the Spearman’s rank correlation test $\rho$ and p-values are given inside each panel. PL_SIC models have had their $\alpha^{model}$ values increased by 0.25 to be comparable due to their model construction. We do not include the outlier GLEAM J003652-333315 in our statistical analysis as it is Haro-11 the Lyman-continuum leaker with extreme IR properties.

Figure 4. Comparisons between the radio SFR versus redshift and stellar mass in panels (a) and (b) and modelled spectral index versus redshift and stellar mass in panels (c) and (d), respectively. The slope of the weighted linear fit and its 1 $\sigma$ uncertainty and the Spearman’s rank correlation test $\rho$ and p-values are given inside each panel. We do not include the outlier source GLEAM J003652-333315 in our statistical analysis. The open points are compiled from the LIRG/ULIRG samples of Galvin et al. (2018), Dey et al. (2022) and (2024) and are separated based on their most favoured radio SED model with red stars being sources which include LFTOs in their radio SEDs and blue circles being all other sources. Grey dashed lines and relationships are the fits from Dey et al. (2024) to compare to our SFG sample. The blue dashed line is from Heesen et al. (2022) and measures $\alpha^{0.15}_{1.4}$ against total galaxy mass (i.e. it probes a flatter part of the radio SED). PL_SIC models have had their $\alpha^{model}$ values increased by 0.25 to be comparable due to their model construction.

The low-frequency spectral index is not related to the modelled spectral index for the total sample due to the construction of the FFA_PL model including a LFTO. For the non-FFA_PL sample $\alpha^{model}$ is not correlated with $\alpha^\textrm{GLEAM}_\textrm{low}$ and we see that $\alpha^\textrm{GLEAM}_\textrm{ low}$ is much flatter. All but one source exhibits some spectral index flattening between $\alpha^\textrm{ATCA}_\textrm{high}$ and $\alpha^\textrm{GLEAM}_\textrm{low}$ indicating that loss and absorption processes play a key role in modifying the synchrotron spectral index. Merging systems are found to have flatter modelled spectral indices than non-merging systems.

5.2.3. Modelled spectral index correlations

We compare the modelled spectral index with a number of properties to provide insight as to the possible different physical processes occurring within our LFTO and control samples. Our sample lies within the range of values of the spectral index for LIRGS/ULIRGS with higher SFR galaxies being more easily detected at higher redshifts (Galvin et al. 2018; Dey et al. 2024). Fig. 4 also shows that there is no significant relationship (as determined by the p-value of the Spearman rank correlation test) between the modelled spectral index and redshift indicating that IC losses due to CMB photons are unlikely to be a dominant mechanism acting to steepen the spectral index as expected due to the low redshifts of our SFG sample. Galvin et al. (2018) and Dey et al. (2024) also find no significant relationship between spectral index and redshift for LIRGS/ULIRGS even out to a redshift of 0.4.

On the other hand, we do find a statistically significant relationship between stellar mass and the modelled spectral index whereby more massive galaxies have a higher radio SFR and steeper modelled spectral index. This effect is seen in An et al. (2021), Heesen et al. (2022), An et al. (2024), and Dey et al. (2024) with negative correlations between stellar mass and spectral index. The steepening of the spectral index with increasing stellar mass is a result of higher mass galaxies generally being larger (Gürkan et al. 2018) which causes the CREs to take longer to escape from the galaxy thus synchrotron cooling losses become important. This effect is highlighted by the negative correlation in Fig. 5 between $\alpha^{model}$ and the galactic radius which is also seen in Heesen et al. (2022) (however, their results probe a flatter part of the radio SED hence the vertical offset).

The modelled spectral index decreases with increasing radio SFR as shown in Fig. 6; however, the correlation is not statistically significant similar to the result of Dey et al. (2024) for the IR derived SFR. This modelled spectral index-SFR relationship is secondary with galaxy size and synchrotron cooling primarily driving the correlation between these variables as shown by the Spearman rank correlation coefficient for the $\alpha^{model}$ -Radius correlation being -0.51 $(p = 0.06)$ compared to -0.37 $(p=0.2)$ for the $\alpha^{model}$ -SFR $_{1.4 GHz}$ . We do not find a statistically significant relationship between radio SFR surface density and modelled spectral index similar to Heesen et al. (2022). These results are contrary to the expectation (and results in Tabatabaei et al. 2017) that galaxies with higher SFRs and younger more energetic CREs have a flatter spectral index. Overall we do not see any separation in radio derived properties between LFTO and non-LFTO (or FFA_PL and non-FFA_PL) samples or those in Galvin et al. (2018) or Dey et al. (2022); Dey et al. (2024).

Figure 5. The modelled spectral index in comparison to the K-band light major axis radius. The slope of the weighted linear fit and its 1 $\sigma$ uncertainty and the Spearman’s rank correlation test $\rho$ and p-values are given. We do not include the outlier source GLEAM J003652-333315 in our statistical analysis. The blue dashed line is from Heesen et al. (2022) and measures $\alpha^{0.15}_{1.4}$ against the star formation radius (i.e. it probes a flatter part of the radio SED). PL_SIC models have had their $\alpha^{model}$ values increased by 0.25 to be comparable due to their model construction.

Figure 6. Comparisons between the modelled spectral index with radio SFR and radio star formation rate surface density in panels (a) and (b), respectively. The slope of the weighted linear fit and its 1 $\sigma$ uncertainty and the Spearman’s rank correlation test $\rho$ and p-values are given inside each panel. The blue dashed line is from Heesen et al. (2022) and measures $\alpha^{0.15}_{1.4}$ against TIR SFR (i.e. it probes a flatter part of the radio SED). PL_SIC models have had their $\alpha^{model}$ values increased by 0.25 to be comparable due to their model construction. We do not include the outlier source GLEAM J003652-333315 in our statistical analysis.

5.3. Global mid-infrared properties

The WISE colour-colour diagram (see Fig. 7) provides a diagnostic for the dominant mid-IR emission mechanism and activity of nearby galaxies. Our sample primarily lies in the region occupied by star-forming disc galaxies at the top end of the SF sequence identified in Jarrett et al. (2019) which is expected based on their morphology and measured SFRs. There is no significant separation between the LFTO/control or FFA_PL/non-FFA_PL samples in this parameter space. Two control galaxies GLEAM J040226-180247 and GLEAM J002238-240737 have warmer mid-IR colours but show no evidence for significant AGN emission in their mid-IR, radio or optical spectra. GLEAM J074515-712426 lies below the SF sequence likely due to contamination by a saturated foreground star which causes a deficit of W1-W2 emission throughout the field. Lastly GLEAM J003652-333315 is found to be an extremely mid-IR bright dust-obscured source also known as Haro 11. Haro 11 has been shown to be a starbursting blue compact galaxy with both dust obscured star-forming regions and Lyman- $\alpha$ leakage (Östlin et al. 2015; Östlin et al. 2021). GLEAM J003652-333315 will therefore have somewhat overestimated (but still comparable to those of Östlin et al. 2015) mid-IR based parameters including M $_*$ , SFR $_{mircor}$ and sSFR $_{mircor}$ due to the dust-obscuration.

Our sample lies above the SFR-M $_{*}$ main sequence in the WISE and GALEX Atlas of Local Galaxies (Leroy et al. 2019). Our sample also has high sSFRs (see Fig. 8) primarily due to their selection criteria requiring that their fluxes are measurable in C and X band ATCA radio observations. The minor separation in M $_*$ between the LFTO and control samples is observable with the LFTO sample probing masses up to an order of magnitude lower than the control sample. GLEAM J003652-333315 has SFR $_{mircor}$ and sSFR $_{mircor}$ more comparable to the more distant LIRGs/ULIRGs in Dey et al. (2022); Dey et al. (2024) due to its dusty obscured starbursts. Overall we see no distinct separation between galaxies which contain LFTOs and those which do not in both our sample and those in Dey et al. (2022); Dey et al. (2024) (selected based on their most favoured model containing a low- $\nu$ FFA component).

5.4. Inclination and mergers

We investigate whether there are any relationships between the morphology or inclination and the radio SED features observed. Inclinations for these galaxies are estimated using the ratio between major and minor K-band light axes where axes are compiled from NED with sources given in Table 1. Inclinations are compared to the GLEAM and modelled spectral index with merging systems flagged in Fig. 9. Overall we do not see any significant correlations between inclination and spectral index (either GLEAM or modelled) or whether edge on SFGs are more likely to contain LFTOs agreeing with the findings of Hummel (1991), Chyży et al. (2018), and Heesen et al. (2022). The two galaxies best modelled by PL_SIC models (square points) including synchrotron and IC losses are both some of the highest inclination sources in this SFG sample. We also examine whether inclination is correlated with the global SFR surface density and find no significant relationship. There is no separation between our samples based on inclination.

5.5. Global radio and infrared properties

We see that our sample shows a statistically significant linear relationship between radio-SFR and SFR $_{mircor}$ in Fig. 10 providing further evidence for the relationship between the IR and radio emission mechanisms. The SFR $_{mircor}$ -excess outlier source GLEAM J003652-333315 is the dust obscured blue compact galaxy Haro 11 (Östlin et al. 2015). The radio-SFR for this source is potentially underestimated as its compact size allows for the escape of CREs into the galactic halo on timescales shorter than the synchrotron loss timescale.

The $q_{FIR}$ parameter, which is the logarithmic ratio between the far infrared flux and 1.4 GHz flux density of an object, is a parameterisation of the FRC, where $q_{FIR}$ is defined as:

(14) \begin{equation} q_{FIR} = \log\bigg(\frac{\textrm{FIR}}{\mathrm{3.75\times10^{12}\;W\,m^{-2}}}\bigg) - \log\bigg(\frac{S_{1.4}}{\mathrm{ W\,m^{-2}\,Hz^{-1}}}\bigg).\end{equation}

Where S $_{1.4}$ is the flux density at $\nu$ = 1.4 GHz, and FIR is defined as:

(15) \begin{equation} \textrm{FIR} = 1.26\times10^{-14} (2.58S_{60} + S_{100})\,\textrm{W m}^{-2},\end{equation}

where S $_{60}$ and S $_{100}$ are the 60 and 100 $\mu$ m band flux densities from IRAS in Jy (Helou, Soifer, & Rowan-Robinson 1985). The mean $q_{FIR}$ value between 60 $\mu$ m and 1.4 GHz is typically taken as 2.34 for SFGs (Yun et al. 2001). $q_{FIR}$ is known to decrease with redshift (Magnelli et al. 2015; Delhaize et al. 2017) ^d and luminosity (Molnár et al. 2021). Deviations from the typical $q_{FIR}$ value can be a critical diagnostic of the physical processes driving the radio and IR emission. Radio excess objects ( $q_{FIR}$ < 1.6) are usually associated with AGN emission whilst IR excess objects ( $q_{FIR}$ > 3) may be dust obscured AGN or compact starbursts. With much of the scatter being influenced by the variation in dust temperature, extinction, and the different timescales associated with different SFR indicators.

Table 8. SFG sample derived IR properties.

Note: Column (1): GLEAM source ID. Column (2): WISE band W1 subtracted from WISE band W2 magnitude. Column (3): WISE band W2 subtracted from WISE band W3 magnitude. Column (4): Mid-IR + UV corrected SFR from Cluver et al. (2024). Column (5): Mid-IR + UV corrected specific SFR from Cluver et al. (2024). Column (6): Bolometric 60 $\mu$ m luminosity from Moshir et al. (1990). Column (7): q $_{FIR}$ as calculated following Yun et al. (2001) using RACS-mid 1.4 GHz radio flux densities.

^aThis source is undetected in IRAS source catalogues.

Our SFG sample lies within the large scatter of previous $q_{FIR}$ observations shown in Fig. 11 (Yun et al. 2001; Magnelli et al. 2015; Delhaize et al. 2017; Galvin et al. 2018; Dey et al. 2024) with no FIR or radio excess sources. We extend the samples of SFGs with spectral curvature down 2 orders of magnitude in L $_{60\,\mu {\rm m}}$ to non-LIRG sources. When comparing $q_{FIR}$ we find the clearest separation between the control (blue) and LFTO (pink and red) galaxy samples with the LFTO sample having a significantly higher mean $q_{FIR}$ value than the control galaxies (see Table 9). This indicates either an excess in IR emission in LFTO galaxies or a radio deficit at 1.4 GHz compared to the control sample. This separation is not seen in the sample of LIRGS/ULIRGS from Galvin et al. (2018) and Dey et al. (2024). There is a significant negative correlation between $q_{FIR}$ and FIR luminosity or stellar mass which agrees with the findings in Molnár et al. (2021). We also see a significant decrease in $q_{FIR}$ with redshift even over the small redshift range observed, however, due to our short lookback times this is likely a result of scatter and is not seen in previous research (Magnelli et al. 2015; Delhaize et al. 2017; Galvin et al. 2018; Dey et al. 2024). These results are summarised in Fig. 11.

Figure 7. WISE colour-colour diagram for our SFG sample. Magnitudes are in the Vega system with calibration described in Jarrett et al. (2011). Regions roughly delineate source types into the labelled categories with AGN and extrema including luminous dust-obscured starbursts (GLEAM J003652-333315). The grey dotted line indicates the ‘star formation sequence’ identified by the 100 largest galaxies in the WXSC (Jarrett et al. 2019).

Lastly we compare the modelled spectral index to a number of IR derived physical properties in Fig. 12. In panels (a) and (c) we see that there is no significant correlation between the modelled spectral index and the mid-IR + UV derived SFR or sSFR. We observe no significant separation between the LFTO and control or FFA_PL and non-FFA_PL samples but can clearly see the flatter spectral indices in merging systems as well as their elevated sSFRs. This suggests mergers can trigger new starburst activity which raises the sSFR and injects new CREs with a flatter spectral index (Murphy et al. 2009; Murphy 2013; Donevski 2015). In panels (b) and (d) we examine the relationship between the GLEAM and modelled spectral indices and $q_{FIR}$ and find that while there is a positive correlation it is not statistically significant above the 5% level. Panel (b) shows the most significant separation between LFTO and control sources whereby LFTO sources occupy the top right region of the plot having flatter GLEAM spectral indices and higher $q_{FIR}$ values. This provides evidence that the spectral flattening at low frequencies is more strongly tied to the ISM properties of the galaxy than the high-frequency spectral index ( $\alpha^{model}$ ).

We present p-values for Welch’s T-tests between the LFTO and control, FFA_PL and non-FFA_PL, and merger and non-merger samples for a number of parameters in Table 9. We verify that, by selection, the LFTO and control samples have $\alpha_\textrm{L}$ mean values that are different at a 1% significance level. This separation also is observed in the FFA_PL and non-FFA_PL samples (see Fig. 3) as SFGs with the flattest $\alpha_\textrm{L}$ have most preferred models which include the FFA. We find a 5% significant separation between LFTO and control galaxies in M $_*$ and SFR $_{1.4GHz}$ with control galaxies having higher values and low variance. This is likely a selection and small sample effect as we see that LIRGs/ULIRGs with high stellar mass and SFRs often contain spectral curvature (Galvin et al. 2018; Dey et al. 2022; Dey et al. 2024). We find a separation between LFTO and control galaxies in $q_{FIR}$ which does not carry through to the FFA_PL and non-FFA_PL samples. This suggests that a clearly modelled LFTO doesn’t necessarily indicate an elevated $q_{FIR}$ value but SFGs with spectral flattening generally have a higher $q_{FIR}$ as shown in Fig. 12. We find a separation between merging and non-merging SFGs in $\alpha^\textrm{ATCA}_\textrm{high}$ and $\alpha^{model}$ (at $P \lt 0.05$ ) with merging systems having flatter spectral indices as seen in Fig. 3. It has been shown that mergers can trigger star formation flattening the spectral index (Donevski & Prodanović 2015) and increasing the sSFR (Murugeshan et al. 2021) which shows a 5% significant separation in sSFR $_{mircor}$ .

Figure 8. (a): WISE mid-IR+UV corrected SFR (Cluver et al. 2024) versus stellar mass for our SFG sample. (b): WISE mid-IR+UV corrected specific SFR (Cluver et al. 2024) versus stellar mass for our SFG sample. The orange background sample and black dashed SFG main sequence best fit is from the WISE and GALEX Atlas of Local Galaxies (Leroy et al. 2019). The open points are compiled from the LIRG/ULIRG samples of Dey et al. (2022) and (2024) and are separated based on their most favoured radio SED model with red stars being sources which include LFTOs in their radio SEDs and blue circles being all other sources.

Figure 9. Comparisons between the inclination and GLEAM spectral index, modelled spectral index and star formation rate surface density in panels (a), (b), and (c), respectively. Edge-on sources have cos(i) $\sim$ 0 whilst face-on sources have cos(i) $\sim$ 1. PL_SIC models have had their $\alpha^{model}$ values increased by 0.25 to be comparable due to their model construction.

6. Discussion

6.1. SED modelling

Accurately modelling the physical processes which occur in SFGs and how they manifest in the radio continuum emission observed is a complex problem caused primarily by the lack of accurate radio data across a wide frequency range. We have observed that over the last few decades the preferred radio SED models have gone from simple single power laws to models encompassing thermal and non-thermal emission (Pacholczyk 1970; Condon 1992) with loss or absorption processes occasionally being invoked in sources which are not well modelled by simple power laws (Clemens et al. 2010; Galvin et al. 2018; Dey et al. 2022; Dey et al. 2024). The limitations in this research follow a similar vein such that it is difficult to conclusively confirm whether the features observed in our radio SEDs are the result of these loss and absorption processes or due to lack of high quality data. For example, of the initial sample of 427 SFGs in the GLEAM-6dFGS catalogue, 54 have negative flux values in at least one subband, whilst the remainder may have photometric errors which can introduce a false sense of curvature in the radio SED. Our sample selection aimed to mitigate these issues by selecting brighter sources with distinct LFTOs such that the photometry is more reliable; however, it is likely that more accurate observations with GLEAM-X and the SKA will reveal a larger population of SFGs with LFTOs highlighting the importance of including physical loss and absorption mechanisms during radio SED modelling. Table 4 shows that synchrotron only, PL based models are most preferred for all of our sources (except J1847) despite the consensus that radio emission in SFGs is due to both non-thermal synchrotron and thermal free-free emission (Condon 1992). We attribute this to lack of high-frequency radio data and limitations of the modelling framework as we discuss below.

Figure 10. The 1.4 GHz radio-SFR compared to the mid-IR+FUV corrected SFR. The slope of the weighted linear fit with its 1 $\sigma$ uncertainty and the Spearman’s rank correlation test $\rho$ and p-values are presented.

Thermal free-free emission generally does not contribute significantly to radio emission at low frequencies contributing between 1–10% at 1.4 GHz (Condon 1992); however, more recent work has been leaning towards 2–4% (Galvin et al. 2018; Dey et al. 2022; Dey et al. 2024). The flat spectral index of thermal emission is also typically thought to become dominant above $\sim$ 30 GHz (e.g. in M82; Condon 1992) but again with recent work we find that in other starburst galaxies this is likely an overestimation with Galvin et al. (2018), Dey et al. (2022); Dey et al. (2024) often finding values of TF $_{40GHz}$ < 0.5. There is a large scatter in thermal fractions which is linked to the star formation timescales with thermal emission being tied to instantaneous star formation with timescales of $\sim$ 10 Myr (Condon 1992) whilst non-thermal emission is typically delayed and has a timescale of $ \gt $ 30 Myr (Condon 1992; Vollmer et al. 2020; Heesen et al. 2024). Thus starbursts with increasing or decreasing SFRs will result in different thermal fractions (and radio SFR estimates; Cook et al. 2024) emphasising the importance of knowing the SFR histories of SFGs. We will explore the impact of SFR histories and their connection to thermal fractions, radio SED curvature and the FRC/TRC in a future paper in this series.

At our current highest frequency point at 17 GHz the flattening we would expect to see at higher frequencies is typically not observed meaning that the thermal fractions are likely low. This, coupled with the fact that the number of free parameters in our models are currently low without the inclusion of GLEAM flux covariance modelling and FIR emission, means that the addition of model complexity by including a thermal component has a huge impact on the selected model. For example, the difference between a model with 2 and 3 free parameters is far greater than the difference between SED models with 10 and 11 free parameters. Thus because we do not observe this spectral flattening and including unconstrained thermal emission in our models increases model complexity a relatively large amount we do not favour SFG based models despite their more physically realistic nature. The increased spectral sampling afforded by ALMA 40 and 115 GHz ^e and FIR observations compounded with the increased model complexity when including covariance and FIR emission models would likely result in SFG based models being favoured in a future paper in this series. This may also cause the measured model spectral indices to be steeper and more comparable to the work of Galvin et al. (2018) especially for sources with a higher thermal fraction indicating a higher contribution of synchrotron or IC based loss processes.

Our modelled synchrotron spectral index values lie within the range of values found in Galvin et al. (2018) and Dey et al. (2022); Dey et al. (2024) with them being found to be, on average, steeper than the canonical value assumed for SFGs $\alpha$ = -0.8 (Condon 1992) especially when considering models which account for spectral steepening (PL_FFA and PL_SIC). PL_SIC models are competitive in 12 of the 19 SFGs in our sample indicating that synchrotron losses are potentially an important physical loss mechanism. This mechanism should be considered during radio SED modelling above 1.4 GHz, especially for large, massive galaxies. One outcome of this is that higher frequency flux estimates derived from two point spectral indices at low frequencies will likely be overestimated and the spectral index used for radio K-corrections may have a mass dependence. Interestingly, however, PL_SIC models which take into account spectral steepening at high frequencies are less likely to be completely preferred than PL_FFA models which have an extra free component. This is likely due to the inability to constrain the break frequency when steepening is so gradual. Despite this we speculate along the same lines as Dey et al. (2022); Dey et al. (2024) that it is likely to be synchrotron losses causing this steepness rather than IC losses or FFA.

Figure 11. $q_{FIR}$ compared to the IRAS 60 $\mu$ m luminosity with the relationship from Yun et al. (2001) shown in panel (a). Panels (b) and (c) show the comparison between $q_{FIR}$ and the stellar mass and redshift, respectively. The slope of the weighted linear fit with its 1 $\sigma$ uncertainty and the Spearman’s rank correlation test $\rho$ and p-values are given inside each panel. PL_SIC models have had their $\alpha^{model}$ values increased by 0.25 to be comparable due to their model construction. We do not include the outlier source GLEAM J003652-333315 in our statistical analysis.

Over half of our sources are best fit by models including FFA whether at low or high frequencies which supports the growing evidence for more complex radio SED model requirements. The inclusion of FFA accounts for both curvature at low frequencies or kinks at higher frequencies which are often evident in well-sampled radio SEDs. In Section 5.5, see no distinct separation in global properties (except $q_{FIR}$ ) for SFGs which contain either low- or high-frequency FFA suggesting that the regions in which FFA is occurring are being averaged out and require a resolved analysis. It was also found that a number of our SFG galaxy sample that were initially selected as having LFTOs are best modelled without the inclusion of the FFA_ prefix model meant to model this spectral curvature. This is attributed to the large errors on the GLEAM flux densities for some sources which results in simpler models being preferred over the increased complexity added by including a relatively unconstrained FFA turnover.

Ionisation losses are also theorised to be a significant loss process for low energy non-thermal emission in SFGs (Lacki & Thompson 2010; Longair 2011; Basu et al. 2015; Roth et al. 2023; Roth et al. 2024) which under steady continuous injection of electrons act to gradually flatten the low-frequency spectral index by $\Delta \alpha$ = +0.5 towards $\alpha$ = -0.1 (Basu et al. 2015). However, for single injection events the resulting spectrum is dependent on the ratio of timescales between injection and ionisation loss (Basu et al. 2015). We chose not to include ionisation losses in this work as the shape of this model and its rate of curvature is currently unconstrained, however, recent work suggests ionisation losses do play a key role in flattening the spectral index of galaxies, particularly those with high SFRs (Roth et al. 2023; Roth et al. 2024). As this is the case we cannot currently disentangle the contributions of ionisation losses and FFA to the radio SEDs of our sources; however, only absorption processes allow for a decrease in the measured flux as we move to lower frequencies (i.e. $\alpha$ > 0) so it can still be safely concluded that FFA is occurring in some of these SFGs.

We also do not include a number of other loss processes within our models as they do not significantly impact the observed spectral index these include adiabatic losses, diffusion losses, and bremsstrahlung losses. These three loss mechanisms play a role in the observed radio emission of SFGs over five decades in SFR but due to their relative independence with electron energy in the non-thermal radio regime they do not act to modify the spectral index at different frequencies and result in a constant power-law spectral index with $\alpha$ = -0.6 (Lacki & Thompson 2010; Roth et al. 2023; Roth et al. 2024).

Lastly spectral curvature in the radio SED can also be caused by synchrotron self absorption or the Tsytovitch-Razin effect when the refractive index of the medium is less than unity (Israel & Mahoney 1990). Both these effects, however, occur below the observed frequency range of our GLEAM observations when considering the peak flux and size of our sources or would require unreasonably high magnetic field strengths (Kellermann & Pauliny-Toth 1969) so they can be safely ruled out as the cause of LFTOs in our sample following (Dey et al. 2024).

6.2. Global properties

We look to connect the features of the radio SED to the global properties and allow us to predict the physical processes occurring within SFGs for unresolved galaxies in future large scale radio surveys. With the improved frequency sampling and sensitivity afforded by future surveys we will begin to be able to disentangle the dominant cooling and loss mechanisms and infer the ISM properties of SFGs. We find that the parameter most strongly correlated with the modelled spectral index is the stellar mass whereby more massive galaxies have steeper spectral indices, and this however is a secondary effect as it is in fact galactic size which causes this correlation with more massive galaxies tending to be larger and more star-forming (Gürkan et al. 2018). Due to the diffusion scale length and lifetime of CREs larger galaxies are able to retain electrons as they undergo synchrotron losses with higher energy CREs losing energy faster than low energy CREs. This causes the spectral steepness we observe particularly above 1.4 GHz. This effect is more pronounced in our SFG sample than those of Heesen et al. (2022) and Dey et al. (2024) as shown by our steeper relationship between $\alpha^{model}$ and $M_{*}$ (see Fig. 4) with SFGs that have preferred models that account for the spectral steepening (_SIC and _FFA) being generally higher mass. However it is important to note that the spectral index measured in Heesen et al. (2022) is a two point spectral index at a lower frequency so it does not account for spectral curvature and will be intrinsically flatter than our modelled spectral index. Overall this suggests that synchrotron losses are the primary driver for spectral steepness above 1.4 GHz for larger nearby SFGs.

Table 9. T-Test p-values between samples.

Welch’s T-test p-values with columns indicating the two samples compared after removal of outlier source GLEAM J003652-333315. Samples which reject the null hypothesis that their means are the same at the 1% and 5% level are highlighted in green and grey, respectively. We do not compare inclinations between merging and non-merging systems as measuring the inclination of merging systems with axes light ratios is unreliable.

Typically the impact of IC losses is unable to be disentangled from synchrotron losses as they result in a similar spectral shape and their loss timescales have similar CRE energy dependence. IC losses due to CMB photons however primarily occur at high redshifts in low density galaxies due to increased CMB photon density (Lacki & Thompson 2010; Klein, Lisenfeld, & Verley 2018). We however probe low redshift galaxies and find no statistically significant correlation to the spectral index such that we can likely rule out the impact of CMB based IC losses for our sample. Dey et al. (2024) do not find a correlation with redshift out to $z\sim 0.4$ due to small sample size and that they also do not probe out to sufficient redshifts for the CMB photon density to become significant. Magnelli et al. (2015) and An et al. (2024) also do not find evidence for IC losses finding flat constant (two point) spectral indices out to $z\sim 2$ , however the spectral indices probed in these papers are below 1.4 GHz where IC losses have long timescales and so the lack of evidence for IC losses are unsurprising. Another possibility is that IC losses are occurring due to high FIR photon energy densities in optically thick starburst galaxies; Lacki & Thompson (2010) explore this but find these galaxies begin to violate the FRC, whilst Basu et al. (2015) and Roth et al. (2024) show that IC losses are more important for low density systems. A potential test that could be performed to see whether IC losses in starbursts play a role in steepening the radio SED is by testing whether we see enhanced X-ray emission in starbursts which has been upscattered, however disentangling this emission from other sources of X-ray emission is currently not feasible. IC losses will however play a role in steepening the observed radio SED at high frequencies so will need to be considered during high redshift radio surveys particularly above rest frame 1.4 GHz, and, radio K-corrections which rely on the radio spectral index above 1.4 GHz will likely need to have a frequency dependence which acts to steepen the spectral index at higher frequencies.

One would expect however that the spectral index would flatten with increasing SFR as there is a constant injection of young electrons, or that lower SFR galaxies would have a higher proportion of aged CREs. The relationship observed between modelled spectral index and 1.4 GHz radio SFR (which is a delayed tracer) of star formation, however, presents the opposite effect (see Fig. 6). This result is consistent with Heesen et al. (2022), Galvin (2019), and Dey et al. (2022); Dey et al. (2024) which all rely on more instantaneous IR SFR tracers. However, as suggested before, this trend is likely a secondary effect to galaxy size or mass with the correlation of spectral index being stronger and more significant with both of these parameters. This observed global correlation of steepening spectral index with increasing SFR also breaks down within galaxies as we see the higher SFR or SFR surface density regions are those with the flattest spectral index (Tabatabaei et al. 2017; Heesen et al. 2022). This global averaging also likely explains the lack of correlation between spectral index and total SFR surface density for our sample, similarly to Figure 7 in Heesen et al. (2022) which shows that size is the driving factor. We do however find a slightly stronger correlation between sSFR and spectral index with higher sSFR (related to source compactness Elbaz et al. 2011) sources having flatter spectral index, consistent with Murphy (2013) and An et al. (2021). This again is consistent with larger more massive galaxies having a steeper spectral index due to retaining ageing CREs for longer timescales such that they can undergo synchrotron losses.

Figure 12. Comparisons between the modelled spectral index and the IR SFR, sSFR and $q_{FIR}$ in panels (a), (c), and (d), respectively. Panel (b) compares the GLEAM spectral index to $q_{FIR}$ . The slope of the weighted linear fit and its 1 $\sigma$ uncertainty and the Spearman’s rank correlation test $\rho$ and p-values are given inside each panel. PL_SIC models have had their $\alpha^{model}$ values increased by 0.25 to be comparable due to their model construction. We do not include the outlier source GLEAM J003652-333315 in our statistical analysis.

The only parameter with a large separation between the LFTO and control samples is $q_{FIR}$ (see Fig. 11 panel (b)) which implies there may be differences in their SFR timescales, emission sources or ISM properties. Curvature in the low-frequency spectral index is also related to the ISM properties of galaxies with both FFA and ionisation losses being dependent on ionised and neutral gas densities, respectively (Lacki & Thompson 2010). This separation then suggests that these LFTO sources that have higher $q_{FIR}$ and flattened spectral indices have different ISM or star formation properties than the control sources. We can likely rule out AGN emission as a reason for this separation as we have no evidence for AGN emission in their optical spectra or significant $q_{FIR}$ deficits. The spectral flattening and turnovers in the LFTO sample imply a higher gas density of ionised gas which is related to HMS formation and/or neutral gas which is the fuel for star formation in both molecular clouds and HI reservoirs. Thus we would expect these galaxies to have high instantaneous SFRs (as the HMS lifetime is short) and/or be gas rich. We unfortunately do not have resolved HI or H $\alpha$ observations for most of these sources so we can not currently determine whether they have enhanced neutral or ionised gas densities. ^f

We do have evidence of these galaxies having high instantaneous SFRs, in their high mid-IR SFRs and their elevated $q_{FIR}$ values. Elevated $q_{FIR}$ values are a result of an excess of FIR compared to radio emission which would be the case in young starbursts or galaxies with increasing SFRs (Galvin 2019; Cook et al. 2024) as the IR emission is produced approximately 10 Myr before the synchrotron dominated 1.4 GHz radio emission. These differences in SFR indicator timescales and contamination by AGN emission are partly responsible for scatter in the $q_{FIR}$ relationship (Galvin 2019; Molnár et al. 2021) and would be somewhat mitigated by considering the short timescale thermal emission at 40 GHz (however, this introduces errors with radio SED decomposition or H $\alpha$ scaling relationships). Overall this paints a picture of recent compact starbursts in regions of high gas density being the cause of the spectral curvature and flattening in the LFTO sample. Thus comparisons between resolved spectral index, FRC and gas density maps as well as measuring SFR histories are vital to understanding complex radio SED features and will be explored in future papers in this series.

6.3. Inclination and mergers

In Section 5.4, we find no relationship between inclination and either modelled or GLEAM spectral index for our SFG sample. This suggests that the spectral index is not dependent on the viewing angle of galaxies and that spectral curvature or flattening is not an effect of looking through the galactic disc. Hummel (1991) and Heesen et al. (2022) also find no connection between the low-frequency spectral index and inclination. For FFA to occur a sufficient electron density along a line of sight is required, dense ionised regions typically occur within individual star-forming molecular clouds which exist throughout the galactic disc. It is then the superposition of these individual regions within the large synthesised GLEAM beam that causes us to observe this spectral curvature. Interestingly if we remove the merging systems from Fig. 9 panel (b) we see that the modelled spectral index (which is most strongly related to $\alpha^\textrm{ATCA}_\textrm{high}$ ) is steeper for edge on galaxies. This may be an effect of ATCA sensitivity or CRE diffusion into the halo and would be the case if off-planar radio emission from a radio halo due to CRE transport was preferentially detected at low frequencies (Vollmer et al. 2020) due to their longer lifetimes (Heesen et al. 2024).

As a number of sources in our SFG sample are merging we investigate the effects of mergers on the radio SED and other global properties. Mergers have been shown to enhance the sSFR of galaxies as turbulence and gas inflows can trigger star formation (Ellison, Catinella, & Cortese 2018). These triggered starbursts inject young CREs with flatter spectral indices which have not had the time to undergo synchrotron cooling. We find that the merging systems in our sample have systematically higher sSFRs and flatter spectral indices both modelled and at higher frequencies indicating that they may have a younger population of CREs or enhanced thermal fractions. Four of these merging systems are members of the LFTO sample supporting the evidence that merger driven starbursts can result in FFA and cause spectral flattening or curvature at low frequencies.

7. Summary and conclusions

We created a sample of 19 GLEAM selected SFGs of which 11 display LFTOs at GLEAM frequencies and eight do not. These SFGs are observed with ATCA between 5.5 and 17 GHz. Their radio SEDs are then modelled between 71 MHz and 17 GHz using a modular radio continuum models within a Bayesian framework which encompass combinations of thermal and non-thermal emission processes as well as FFA, synchrotron and IC losses. We find that:

1. 1. PL based models are preferred for our entire sample of SFGs over models containing thermal free-free emission due to lack of high-frequency radio data displaying spectral flattening and the relatively large increase in model complexity being unfavoured.

2. 2. Of the originally selected 11 LFTO galaxies five favour SED models with no FFA based LFTO which is attributed to large uncertainties in GLEAM flux densities resulting in the inability to accurately constrain the turnover frequency.

3. 3. Radio SED models with spectral ‘kinks’ or gradual steepening are competitive with simple PL only models.

4. 4. The radio emission is most strongly related to the stellar mass (correlated with physical size) of these galaxies such that more massive or larger galaxies generally have higher SFRs and retain synchrotron emitting CREs for longer steepening the spectral index and lowering $q_{FIR}$ .

5. 5. As both GLEAM spectral index and $q_{FIR}$ depend on SFR timescale with short lived HMSs being responsible for increasing the ionised gas density and FIR emission it is likely that LFTOs are transient and dependent on recent starburst activity.

6. 6. The merger systems in our sample have elevated sSFRs and flatter spectral indices than our non-merging SFGs indicating that they are undergoing recently triggered starburst activity. Four (of five) are members of the LFTO sample with three of these sources having elevated $q_{FIR}$ values suggesting that LFTOs may be caused by recently triggered compact starbursts.

7. 7. We find no relationship between the inclination and whether a galaxy contains an LFTO suggesting FFA occurs within specific star-forming regions rather than because we are looking through the galactic disc. Galaxies with favoured models that include synchrotron and IC losses tend to be the most inclined suggesting that CRE diffusion into the galactic halo is partly responsible for the steep spectral index.

Overall we show that for our sample single power-law radio SED models are often appropriate to model radio continuum emission from SFGs especially with limited sampling below 300 MHz and above 10 GHz. However, with improved sampling, especially at low and high frequencies more complex models including loss and absorption processes are often needed to explain the spectral features observed. We also find that no specific global properties act as a predictor of LFTOs for our sample suggesting that resolved observations are required to measure the physical properties in specific star-forming regions in which FFA may be occurring. The relationship between the radio SED and resolved emission properties and star-formation history will be explored in future papers in this series.