Discovery of Hydrogen Radio Recombination Lines at z=0.89 towards PKS  1830-211

Kimberly L. Emig

0000-0001-6527-6954]Kimberly L. Emig
Jansky Fellow of the National Radio Astronomy Observatory
National Radio Astronomy Observatory, 520 Edgemont Road, Charlottesville, VA 22903, USA

0000-0001-7547-4241]Neeraj Gupta
Inter-University Centre for Astronomy and Astrophysics, Post Bag 4, Ganeshkhind, Pune 411 007, India

Pedro Salas
Green Bank Observatory, 155 Observatory Road, Green Bank, WV 24915, USA

Sébastien Muller
Department of Space, Earth and Environment, Chalmers University of Technology, Onsala Space Observatory, SE-43992 Onsala, Sweden

Sergei A. Balashev
Ioffe Institute, Politekhnicheskaya 26, 194021 Saint Petersburg, Russia
HSE University, Saint Petersburg, 190121, Russia

Françoise Combes
Observatoire de Paris, Collège de France, PSL University, Sorbonne University, CNRS, LERMA, Paris, France

Emmanuel Momjian
National Radio Astronomy Observatory, 1003 Lopezville Road, Socorro, NM 87801, USA

Yiqing Song
Department of Astronomy, University of Virginia, 530 McCormick Road, Charlottesville, VA 22903, USA
National Radio Astronomy Observatory, 520 Edgemont Road, Charlottesville, VA 22903, USA

Preshanth Jagannathan
National Radio Astronomy Observatory, 1003 Lopezville Road, Socorro, NM 87801, USA

0000-0001-9174-1186]Partha P. Deka
Inter-University Centre for Astronomy and Astrophysics, Post Bag 4, Ganeshkhind, Pune 411 007, India

Gyula I. G. Józsa
Max-Planck Institut für Radioastronomie, Auf dem Hügel 69, 53121 Bonn, Germany
Department of Physics and Electronics, Rhodes University, PO Box 94, Makhanda, 6140, South Africa

Hans-Rainer Klöckner
Max-Planck Institut für Radioastronomie, Auf dem Hügel 69, 53121 Bonn, Germany

Abhisek Mohapatra
Inter-University Centre for Astronomy and Astrophysics, Post Bag 4, Ganeshkhind, Pune 411 007, India

Pasquier Noterdaeme
Institut d’Astrophysique de Paris, Sorbonne Université and CNRS, 98bis boulevard Arago, F-75014 Paris, France
Franco-Chilean Laboratory for Astronomy, IRL 3386, CNRS and U. de Chile, Casilla 36-D, Santiago, Chile

Patrick Petitjean
Institut d’Astrophysique de Paris, Sorbonne Université and CNRS, 98bis boulevard Arago, F-75014 Paris, France

Raghunathan Srianand
Inter-University Centre for Astronomy and Astrophysics, Post Bag 4, Ganeshkhind, Pune 411 007, India

Jonah D. Wagenveld
Max-Planck Institut für Radioastronomie, Auf dem Hügel 69, 53121 Bonn, Germany

We report the detection of stimulated hydrogen radio recombination line (RRL) emission from ionized gas in a $z=0.89$ galaxy using 580–1670 MHz observations from the MeerKAT Absorption Line Survey (MALS). The RRL emission originates in a galaxy that intercepts and strongly lenses the radio blazar $-$211 ($z=2.5$). This is the second detection of RRLs outside of the local universe and the first clearly associated with hydrogen. We detect effective H144$\alpha$ (and H163$\alpha$) transitions at observed frequencies of 1156 (798) MHz by stacking 17 (27) RRLs with 21$\sigma$ (14$\sigma$) significance. The RRL emission contains two main velocity components and is coincident in velocity with  21 cm and OH 18 cm absorption. We use the RRL spectral line energy distribution and a Bayesian analysis to constrain the density ($n_e$) and the volume-averaged pathlength ($\ell$) of the ionized gas. We determine $\log( n_e ) = 2.0_{-0.7}^{+1.0}$  and $\log( \ell ) = -0.7_{-1.1}^{+1.1}$ pc towards the north east (NE) lensed image, likely tracing the diffuse thermal phase of the ionized ISM in a thin disk. Towards the south west (SW) lensed image, we determine $\log( n_e ) = 3.2_{-1.0}^{+0.4}$  and $\log( \ell ) = -2.7_{-0.2}^{+1.8}$ pc, tracing gas that is more reminiscent of  regions.
We estimate a star formation (surface density) rate of $\Sigma_{\rm SFR} \sim 0.6$  kpc$^{-2}$ or SFR $\sim 50$ , consistent with a star-forming main sequence galaxy of  $M_{\star} \sim 10^{11}$ .
The discovery presented here opens up the possibility of studying ionized gas at high redshifts using RRL observations from current and future (e.g., SKA and ngVLA) radio facilities.

§ INTRODUCTION

Radio recombination lines (RRLs) result from the radiative de-excitation of electrons at high excitation levels of atoms. RRLs with frequencies $\nu_{\mathrm{rest}} \lesssim $10 GHz in extragalactic sources probe gas with relatively low density that can be stimulated by radio continuum <cit.>, upon which the line emission is significantly enhanced compared with the local thermodynamic equilibrium (LTE) Boltzmann distribution. The physical conditions of the gas strongly influence which principal quantum numbers (thus frequencies) show enhanced emission. Accordingly, the relative intensities of RRLs, the so-called spectral line energy distribution (SLED), carry information on the temperature, density, and pathlength of the diffuse gas component <cit.>. Stimulated emission has the added benefit that its intensity is proportional to the background continuum intensity at a given frequency, $S_{\mathrm{RRL}} \propto S_{\mathrm{bkg \, cont}}$. Therefore, it is conceivable to observe these RRLs at cosmological distances wherever bright radio sources are present [Shaver, 1978]. In contrast, RRLs at high radio frequencies ($\gtrsim 10$ GHz) typically trace spontaneous recombination emission (in LTE), making them a great direct measure of ionizing photons and therefore star formation rates. However, the flux of spontaneous transitions falls off with distance ($D$) and frequency ($\nu$) as $S \propto D^{-2} \nu^{-1}$, limiting the observation of this faint emission to galaxies in the nearby universe.

The very first extragalactic detections of RRLs, in M82 and NGC253, found contributions from stimulation [Shaver et al., 1977, Seaquist & Bell, 1977, Shaver et al., 1978]. In the case of M82, [Bell & Seaquist, 1978] modeled the SLED of hydrogen RRLs and showed that the increasing intensity of the RRL emission at $\nu < 10$ GHz can only result due to stimulation by the large-scale synchrotron continuum of the galaxy. They determined that the hydrogen RRLs originate in ionized gas with $n_e \approx 150$  and pathlength $\ell \approx 110$ pc. Soon after, [Churchwell & Shaver, 1979] used the Arecibo 300 m telescope to search 21 galaxies and active galactic nuclei (AGN) for RRL emission with the 1.4 GHz receiver and three AGN with the 430 MHz receiver, with the set-up covering just a single RRL transition. They did not detect emission with line-to-continuum ratios of $\tau_{\mathrm{RRL}} > 10^{-3}$. To similar sensitivities, Bell et al., 1984 used the Effelsberg 100 m dish at 4.8 GHz to search ten galaxies without clear success. Bell & Seaquist, 1980 discovered the H83$\alpha$ and H99$\alpha$ lines at 10.5 GHz and 6.2 GHz, respectively, in the GHz peaked-spectrum source OQ 208 at $z=0.0763$, showing that this RRL emission could also only arise due to stimulation. These studies used narrow bandwidth receivers and were only sensitive to one RRL spectral line per observation.

To date, 8 of the 23 external galaxies with detected RRL emission show evidence for stimulated emission by non-thermal emission <cit.>[We also refer the reader for an update collection of extragalactic RRL detections at <www.tinyurl.com/DatabaseForExtragalacticRRLs>]. These were all local ($D < 350$ Mpc) and mostly star forming galaxies, until recently, [Emig et al., 2019] used the Low Frequency Array to detect stimulated RRLs at $z = 1.124$ with a rest-frame frequency of 284 MHz. They argue that the RRL emission most likely arises from carbon in an intervening galaxy along the line-of-sight to 3C 190. However, they could not clearly discern whether the emission was from carbon or hydrogen since the lines are separated by 150  and have comparable intensities at those frequencies in the Milky Way <cit.>. The improved sensitivity of high-resolution interferometers and the large fractional bandwidths that enable deeper searches through line stacking are making extragalactic RRL detections now feasible. Furthermore, given that the frequency separation between each RRL is unique, the development of new cross-correlation techniques are enabling blind searches of RRLs across redshift space [Emig et al., 2020].

The principal quantum numbers, $\mathsf{n}$, of Hydrogen RRLs at which the stimulated-only emission of gas with a given density peaks. The thick black region marks the maximum $\mathsf{n}$ in the peak of emission, for ionized gas temperatures $ 5000~\mathrm{K} \leq T_e \leq 12\,000$ K, and the thin black region indicates the $\mathsf{n}$ for which the peak intensity is more than one half of the maximum. The shaded blue region shows the parameter space covered by MALS for RRL redshifts of $0 \leq z \leq 3$, while the red hatched region shows the coverage for a $z=0.89$ observation. To calculate the peak optical depths we do not take into account an external radiation field (which has a minimal effect), but do take into account collisional broadening of the line profile occurring in dense gas.

Current large-bandwidth spectral line surveys, such as the MeerKAT Absorption Line Survey <cit.>, the First Large Absorption line Survey <cit.>, the Search for HI Absorption with APERTIF <cit.> – and in the future with the Square Kilometer Array <cit.> – can observe tens of RRLs simultaneously, opening the way for ionized gas studies with optimum sensitivity to gas with electron densities of $1 ~\mathrm{cm^{-3}} \lesssim n_e \lesssim 10^4$  (e.g., see Fig. <ref>). These surveys can explore RRLs as a tracer of (diffuse) ionized gas in external galaxies for the first time in a large systematic way and address at what level stimulated RRLs are present in galaxies. Ultimately, these RRL observations will bring new insight into the evolution of the (ionized) interstellar medium of galaxies and the environment of AGN.

In particular, MALS is carrying out the most sensitive search to date ($\sigma \sim 0.6$  per 6  channel) of  21 cm and OH 18 cm absorption lines at $0 \lesssim z \lesssim 2$ [Gupta et al., 2016]. MALS is observing $\sim$500 pointings centered on the brightest ($S_{\mathrm{1\, GHz}}>0.2$ Jy) radio sources at declination $\delta \lesssim +30^{\circ}$ <cit.>, using the MeerKAT [Jonas & MeerKAT Team, 2016] L band, nominally covering 900–1670 MHz, and UHF band, nominally covering 580–1015 MHz. Towards each sight line, the survey is sensitive to peak RRL line to continuum ratios of $\tau_{\rm RRL} = (0.08 -1) \times 10^{-3}$ through line stacking, reaching the range of optical depths observed from ionized gas in the Milky Way <cit.>. For an electron temperature of 8000 K, these optical depths convert into emission measures $10^{2.6} \gtrsim EM / \mathrm{cm^{-6}~pc} \gtrsim 10^{4.7}$ for densities $1 ~\mathrm{cm^{-3}} \lesssim n_e \lesssim 10^4$  (see Fig. <ref>).

The first science verification observations of MALS [Gupta et al., 2021, Combes et al., 2021] focused on the bright ($S_{1\,\mathrm{GHz}} \approx 11$ Jy) $z = 2.507$ [Lidman et al., 1999] blazar, PKS 1830$-$211 (referred to hereafter as ).  has a radio spectral index of $\alpha_{1-15~\rm GHz} \approx -0.26$, classifying it as a flat spectrum radio quasar [Pramesh Rao & Subrahmanyan, 1988, Subrahmanyan et al., 1990], and at lower frequencies it is known to have a spectral turnover [Lovell et al., 1996].  is strongly gravitationally lensed [Patnaik et al., 1993, Nair et al., 1993] by a galaxy at $z=0.89$ [Wiklind & Combes, 1996]. Its morphology reveals two compact radio components, referred to as northeast (NE) and southwest (SW), approximately 1 apart and surrounded by a low surface-brightness Einstein ring <cit.>. While the NE and SW components are images of the blazar core, the ring is mainly due to the jet and a bright knot in the jet [Jin et al., 2003].  is known to be variable, by up to a factor of two in radio, on timescales of hours to years <cit.>, and these variations are seen in all three lensed components — the NE, SW, and ring.  The continuum flux is dominated by the NE and SW components at least down to 1.4 GHz [Verheijen et al., 2001, Koopmans & de Bruyn, 2005], where the Einstein ring contributes $\sim$1% to the continuum flux [Combes et al., 2021, Patnaik et al., 1993]. For an overview image of the system, we refer the reader to [Nair et al., 1993].

Two absorption line systems are present along the line-of-sight to . The lensing galaxy at $z=0.89$ has become an extragalactic-prototype absorption system, with the most molecular species detected to date (towards the SW image) <cit.>. The $z=0.89$ galaxy has been directly imaged with Hubble Space Telescope and appears to be a barred-spiral [Courbin et al., 1998, Courbin et al., 2002], but remains weak and elusive. From a well-constrained lensing model [Nair et al., 1993, Koopmans & de Bruyn, 2005, Muller et al., 2020, Combes et al., 2021], the kinematics ($v_{rot} \sim 260$ ) and orientation also suggest a nearly face-on ($i \sim 25$) barred-spiral galaxy of $\sim$10$^{11}$ . The SW image of the blazar core passes through the galaxy at a radius of $R_g \sim 2.4$ kpc and a central velocity (from spatially unresolved emission) at $v_{\rm cen} \sim 0$ , and the NE image passes through at $R_g \sim 5.3$ kpc and $v_{\rm cen} \sim - 150$ .
The SW image intercepts a spiral arm of the $z=0.89$ galaxy and traces dense ($n_{H_2}\sim 2000$ cm$^{-3}$) molecular gas [Wiklind & Combes, 1996, Courbin et al., 2002, Muller et al., 2013]. The gas along the line of sight to the NE image arises in the diffuse ISM [Muller et al., 2011] and is bright in  21 cm and main OH 18 cm absorption [Chengalur et al., 1999, Koopmans & de Bruyn, 2005, Gupta et al., 2021, Combes et al., 2021]. No time variation is visible in the cm-line spectra from the $z=0.89$ galaxy [Combes et al., 2021], which contrasts strongly with the variations detected in the mm-wave absorption spectra [Muller & Guélin, 2008, Muller et al., 2014, Schulz et al., 2015]. The second absorption system towards  at $z=0.19$ has been seen only in  21 cm absorption so far [Lovell et al., 1996, Allison et al., 2017, Gupta et al., 2021].

MALS observations of  verified system performance and led to the first detection of OH 18 cm satellite lines at $z = 0.89$, which had previously only been detected at $z \lesssim 0.25$ [Gupta et al., 2021, Combes et al., 2021]. In this article, we use these MALS observations to search for radio recombination line emission. We aim to understand whether RRLs are present and detectable and what they can tell us about photoionized gas in galaxies and AGN. In Sec. <ref>, we describe the methods used to process the data. In Sec. <ref>, we report (i) detections of hydrogen RRLs in both the L and UHF bands originating from ionized gas in the $z=0.89$ galaxy, (ii) tests we performed to verify these results, and (iii) non-detections at the additional redshifts searched. In Sec. <ref>, we constrain physical conditions of the ionized gas by modeling the stimulated RRL emission. Finally, we discuss implications of the ionized gas measured by the RRLs in Sec. <ref> and conclude in Sec. <ref>. In this article, velocities are reported using the heliocentric frame, with respect to $z=0.88582$ unless stated otherwise, and are converted from frequency using the relativistic definition. We use $\Lambda$CDM cosmology with $\Omega_m = 0.29$, $\Omega_{\Lambda} = 0.71$, and $H_o = 70$  Mpc$^{-1}$, for which 1 $\sim$ 7.860 kpc at $z=0.88582$.

§ DATA AND PROCESSING

§.§ Observations and Data Reduction

was observed with the MeerKAT Radio Telescope [Jonas & MeerKAT Team, 2016] as the first science verification target of MALS [Gupta et al., 2016].
For the work presented here, we used MALS L band spectra originally presented by [Gupta et al., 2021]. Hereafter, we refer to this L band dataset obtained on 2019 December 19 as “Night1”. We also used UHF band spectra from the dataset obtained on 2020 July 13 and presented in [Combes et al., 2021]. In addition to these previously published datasets, we observed  a second time in L band on 2020 September 18 using 59 antennas, which we will refer to as “Night2” in the article.

For both L band observations, the total bandwidth of 856 MHz was centered at 1283.987 MHz, covering 856–1712 MHz, and split into 32 768 frequency channels. This delivered a frequency resolution of 26.123 kHz, which is 6.1  at the center of the band. For UHF band, the total observable bandwidth of 544 MHz covering 544–1088 MHz was also split into 32 768 frequency channels, providing a channel resolution of 16.602 kHz, or 6.1  at the center of the band, i.e., 815.992 MHz. For all observations the correlator dump time was 8 s and the data were acquired for all four polarization products, labeled as XX, XY, YX and YY. We also observed PKS 1934–638 and/or 3C 286 for flux density scale and bandpass calibrations. Since  is a bonafide VLA gain calibrator at this spatial resolution, a separate gain calibrator was not observed. The total on-source times on  are: 40 min (L band Night1), 90 min (UHF band) and 90 min (L band Night2).

All MALS datasets have been processed using , the Automated Radio Telescope Imaging Pipeline [Gupta et al., 2021], a Python-based pipeline using tasks and tools from the Common Astronomy Software Applications <cit.>. The specific details of the observations, calibration, and imaging of L and UHF band datasets can be found in Gupta et al., 2021 and Combes et al., 2021, respectively. The Stokes-$I$ continuum flux densities of  obtained from wideband radio continuum images in L Night1 and UHF band with robust=0 weighting are $11.245 \pm 0.001$ Jy at 1270 MHz and $11.40 \pm 0.01$ Jy at 832 MHz, respectively. The radio emission is unresolved in these images with a spatial resolution of $12\farcs9\times8\farcs1$ (position angle = $-76\fdg3$) and $17\farcs4\times13\farcs1$ (position angle = $+69\fdg0$), respectively <cit.>. The quoted uncertainty of the flux densities are based on a single 2D Gaussian fit to the continuum emission. The continuum flux density of  obtained using the Night2 dataset is $S_{1.27~\mathrm{GHz}} = 11.86 \pm 0.02$ Jy. This matches with the Night1 measurement within the flux density uncertainty ($\sim$5%) expected at these low frequencies. Therefore, throughout the article we use the average flux density from the Night1 and Night2 datasets as the representative L band flux density, i.e., $S_{1.27~\mathrm{GHz}} \approx 11.5$ Jy.

The spectral line data products from  for RRL analysis are continuum-subtracted XX and YY parallel hand image cubes obtained with robust=0 weighting. We also note that for spectral line processing  splits the L and UHF bands into 15 spectral windows (hereafter SPWs) <cit.>. The pixel sizes for L and UHF band image cubes are $2\farcs0$ and $3\farcs0$, respectively. XX and YY spectra were extracted in all SPWs from a single pixel at the location of  determined from the continuum images. The residual flux density in each continuum-subtracted spectrum is on the order of 0.5% and required additional spectral-baseline subtraction.

Further details of RRL specific spectral line processing are provided in the next section. In passing, we note that we also made use of a MALS L band dataset of another radio source, PKS 1740-517 (hereafter PKS 1740), also known as J1744-5144, observed on 2020 September 20 (two days after Night2 observations of ) with an on-source time of 56 min. This observation also used PKS 1934-638 and 3C 286 for flux density and bandpass calibrations, and the unresolved radio source has a flux density of $\sim$7 Jy at 1270 MHz, comparable to . We process this dataset following the procedures described above and use the resultant spectra to establish the genuineness of the results obtained for .

§.§ RRL Spectral Processing

Spectral properties prior to stacking RRL transitions for L band spectra from Night 2. Velocities are shown with respect to $z=0.88582$ and spectra are shifted along the ordinate for display purposes. The principal quantum number of each spectrum is given on the right hand side ordinate. The shaded gray region in the left most panel indicates the line blank region. “Noise”, $\sigma$, is the rms (outside of the line blank region) per 8  channel of the XX (yellow circles) or YY (purple circles) spectrum. “Integrated Signal” $= \Delta v \Sigma_{i=0}^{N}\, S_i$, the sum of emission from the channels inside the line blank region. “Integrated S/N” $= (\Sigma_{i=0}^{N}\, S_i) / (\sqrt{N}\sigma)$, the integrated signal divided by the integrated noise; we note that corresponding to negative values of Integral Signal, the Integrated S/N is also negative. The vertical dashed line in the three right panels indicates the median value. There are no obvious and significant spurious features that could contaminate the RRL spectrum in our final stacking.

Considering the rest frequencies of RRLs, 38 and 44 of the $\alpha$ ($\Delta \mathsf{n} =1$) recombination lines (per element) fall within the MALS L and UHF band coverage, respectively.[We refer to the reader to the CRRLpy module [Salas et al., 2016] found at <https://github.com/astrofle/CRRLpy> for a list of RRL line frequencies.] At $z=0.89$, for example, the observations cover 31 $\alpha$ recombination lines in L band spanning principal quantum numbers $\mathsf{n} = 128 - 158$ and 36 $\alpha$ recombination lines in UHF band spanning $\mathsf{n} = 148 - 183$. Because line properties of RRLs are correlated over a sizeable frequency range, we stacked the spectral lines to drive down the noise and increase the signal to noise ratio, as is commonly done in Galactic and extragalactic RRL observations [Balser, 2006, Emig et al., 2020].

We began spectral processing by identifying the observed frequency of an RRL and extracting a spectrum equivalent to $v_{\mathrm{sys}} \pm 1500$  centered on the line. We selected this velocity chunk and discarded coverage outside of it in order to (i) have a sufficient number of channels to minimize errors in the estimation of the spectral-baseline continuum level [Sault, 1994], while (ii) ensuring that a low ($\leq5$) order polynomial – with an order determined by minimizing the reduced $\chi^2$ – could be fit over a well-behaved bandpass.  RRLs that fell within $\pm 1000$  from the edge of a spectral window were excluded from subsequent processing. We flagged channels with persistent radio frequency interference (RFI) and/or at the HI 21 cm and OH 18 cm line features (at relevant redshifts) <cit.>, and we discarded the full spectral line chunk when flagged channels fell within $|v| \sim 500$  from the line center in order to ensure a reliable fit to the spectral-baseline. We also flagged spectral line chunks that were clearly contaminated by broad-band RFI and reliable estimates of the baseline could not be attained.

After flagging, we next fit a low order ($\leq5$) polynomial to (line-free) channels in each spectral chunk and subtracted the fit. For the results presented in Sec. <ref>, we did not impose a line-blank region, except for the $z=0.89$ stacks. For the $z=0.89$ results, we first stacked the spectrum without a line-blank region. Based on the significant feature we obtained in that spectrum, we set the line blank as -230  to +55 .

Fig. <ref> shows the baseline-subtracted spectral chunks, as an example, from L band Night2 observations processed for $z=0.89$ RRLs. The spectral noise has a median of 3.4 mJy across the 17 RRL spectral chunks used in the stack. The noise properties and number of lines used in the final stacks can be found in Table <ref>. The noise properties are consistent across the bands and between parallel hand spectra.

We next interpolated each spectral chunk to a common velocity grid with channel widths of 1 , intentionally oversampling the spectral channels to avoid artificially smoothing-out spectral features. The spectral line chunks were then averaged together using the inverse noise squared in line-free channels as a weight <cit.>.[ We tested combining the spectral chunks with an additional weight that depended upon the (inverse) line frequency (raised to a power), but this did not significantly change the S/N properties, indicating that the line properties are similar and well-correlated across the observing bands.] We next smoothed the channel resolution to 8  using a boxcar averaging kernel, to better match the lowest resolution achieved at the low frequency end of the bands. At this stage, we had obtained a single H$n\alpha$ spectrum for each parallel hand XX and YY. Finally, we averaged the XX and YY spectra to create a Stokes-$I$ spectrum. We chose to combine the stacked XX spectrum and the stacked YY spectrum rather than creating a Stokes-$I$ spectrum for each line before combining polarizations because it (i) resulted in better-behaved, i.e., Gaussian-like noise properties and (ii) had the benefit of independently examining the differences in the signal detected in two parallel hand spectra. The latter is particularly useful in distinguishing true astrophysical signal from RFI, which is often linearly polarized.

§ RESULTS

We applied the spectral processing procedures described in the previous section to  observations at the redshifts (i) $z=0$ for Galactic emission, (ii) $z=0.19259$ for the low redshift intervening absorber, (iii) $z=0.88582$ for the high redshift intervening absorber, and (iv) $z=2.507$ for the intrinsic redshift of . We detect RRL emission from the $z=0.89$ absorber in both Nights of L band and in UHF band observations. We report non-detections and upper limits at all other redshifts and bands.

Spectral and Line Properties

$N_{\rm lines}$
$S_{\rm peak}$
$\int S_{\rm RRL} \, \mathsf{d}v$
$\int \tau \, \mathsf{d}v$

(mJy )

PKS 1830-211        0.88582       L     17        H 144 $\alpha$        0.34      -117.4 ±5.3        1.86 ±0.12     131 ±14        258 ±33        -0.045 ±0.006

7.7 ±4.2       1.54 ±0.19         62.3 ±9.8      102 ±20      -0.018 ±0.003

UHF       27        H 163 $\alpha$        0.57      -124.4 ±9.4           1.59 ±0.15        155 ±24     262 ±47      -0.046 ±0.008

7.7      0.81 ±0.25         80 ±28         70 ±25       -0.012 ±0.004

0.0           L     17        H 175 $\alpha$        0.36        ...        ...        ...        <22.7       <|0.0020|

UHF     28        H 203 $\alpha$        0.41        ...       ...       ...       <26.5       <|0.0023|

0.19259       L     17        H 166 $\alpha$        0.34       ...       ...       ...        23.6 ±7.3     - 0.002 1 ±0.000 6

UHF       33        H 189 $\alpha$        0.43       ...       ...       ...       <27.2       <|0.0024|

2.507         L     14        H 116 $\alpha$        0.40      ...       ...       ...       <25.5       <|0.0023|

UHF       19        H 132 $\alpha$     0.55      ...       ...       ...       <34.9       <|0.0031|

PKS 1740-517       0.88582       L      17        H 142 $\alpha$     0.39      ...       ...       ...       <24.8     <|0.004|

Uncertainties of the line properties are determined from the variance of each parameter as determined by the fit. $N_{\rm{lines}}$ is the effective number of recombination lines stacked in the final spectrum. “RRL” refers to the effective radio recombination line transition of the stacked spectrum. “noise" is the weighted standard deviation of line-free channels in units of mJy per 8 channel. $\rm{v_{center}}$ is the central velocity of the best fit Gaussian and uncertainty. $S_{\rm{peak}}$ is the peak amplitude of the best fit Gaussian fit. FWHM is the full-width half maximum of the best fit Gaussian. $\int S_{\rm RRL} \, \mathsf{d}v$ is the velocity-integrated flux density of the best-fit Gaussian profile, or in the case of an upper limit, the reported value is equal to an integrated signal to noise ratio of 3 assuming a line width of 60 . $\int \tau \, \mathsf{d}v$ is the velocity-integrated optical depth computed as $- \int S_{\rm RRL}/ S_{\rm c} ~ \mathsf{d}v$ where the values used for the continuum flux density, $S_{\rm c}$, are described in Sec. <ref>.

In Fig. <ref>, we show the L band detections obtained from the $z=0.89$ absorber. We overlay the parallel hand spectra, showing that they are consistent within the noise in both nights of observations and therefore not due to low-level linearly-polarized RFI. We also overlay the Stokes-$I$ spectrum obtained from each night of L band observations; they are consistent within the noise, giving further evidence that the signal is astrophysical in nature.
Lastly, we averaged L band Stokes-$I$ spectra from Nights 1 and 2, producing a spectrum with the highest signal-to-noise ratio attainable which we refer to as the total-$I$ spectrum. The reduction in spectral noise of the incrementally combined spectra follows root N statistics. The total integrated flux density of $\int S_{\rm H144\alpha} \, \mathsf{d}v = 337 \pm 16$ mJy  with a maximum integrated signal-to-noise ratio (S/N) $= (\Sigma_{i=0}^{N}\, S_i) / (\sqrt{N}\sigma)$ of 21 is computed from the total-$I$ spectrum by integrating over the $N$ channels covering the velocity range $-230$ to 55 . The effective transition of the total-$I$ spectrum presented in the bottom panel of Fig. <ref> is H144$\alpha$, at a rest-frame frequency of 2179.6 MHz and observed at 1155.8 MHz. The effective transition is determined from the noise-weighted average of the line transitions included in the stack. The line properties of the observed emission are consistent with $S_{\mathrm{H}136\alpha} < 5$ mJy upper limits to the H136$\alpha$ line at $\nu_{\mathrm{rest}} = 2585.7$ MHz ($\nu_{\mathrm{obs}} = 1371.1$ MHz) — a transition which is included in our stack — obtained by [Mohan et al., 2002].

Comparison of RRL stacked spectra in L band at $z=0.89$. The top three panels compare spectra from parallel hands and observing nights. The bottom panel shows (i) the final spectral result for the band with Gaussian profiles fit to two components and (ii) underneath, the final spectrum with the Gaussian profiles subtracted. The vertical dashed lines mark the velocities $-230$  and 55  within which we integrate the signal of the Total component.

In addition to (1) multiple nights of observations and (2) comparing XX and YY parallel hand spectra, we verified additional evidence that the signal is true recombination line emission by (3) observing an independent detection of RRL emission in UHF band (more details below), (4) performing jackknife tests, in which one line spectrum at a time is omitted from the stacked spectrum, showing that the signal is not dominated by a single outlying spectral chunk[We also refer the reader to Figure <ref>, where the noise, integrated signal, and integrated signal-to-noise properties also demonstrate no single or few outlying spectra dominate.], (5) splitting the lines in the band into two groups creating two sub-stacks and this resulted in consistent line properties, (6) verifying a signal is not reproducible by stacking RRLs at other redshifts (more details at the end of the Section and see Fig. <ref>), and (7) finding that an RRL spectrum of PKS 1740 stacked at $z=0.89$ does not systematically produce a signal. We made use of additional MALS L band observations of PKS 1740 and followed the spectral processing procedures described in Sec. <ref>. We created an average RRL spectrum at $z=0.89$ with an effective transition of H142$\alpha$ and show it in Fig. <ref>. This spectrum reached a noise, $\sigma = 0.39$ mJy, comparable to the  stack. However, no emission or significant spectral features are present in the PKS 1740 spectrum, and it is consistent with noise. This lends additional support to the physical and real nature of the emission from the $z=0.89$ absorber in .

Non-detection RRL spectra in the spectrum of  (top three) and in PKS 1740 (bottom). Spectra have been given an arbitrary offset in intensity in multiples of 4 mJy. Velocities are defined with respect to the labeled redshift.

Two velocity components dominate the L band H144$\alpha$ emission from . In the bottom panel of Fig. <ref>, we show the best fit of two Gaussian profiles and the spectrum that results when the two Gaussian profiles are subtracted. The properties of the Gaussian fits are listed in Table <ref>; the errors of each quantity are determined from the variance of each variable as determined by the fit. The H144$\alpha$ component centered on $-117.4 \pm 5.3$  arises from the NE sight-line (and thus we will refer to this velocity component as the NE component hereafter). The H144$\alpha$ component centered on $7.7 \pm 4.2$  arises from the SW sight-line (and thus we will refer to this velocity component as the SW component hereafter). In Table <ref>, we also include the integrated optical depth equal to $\int \tau \, \mathsf{d}v \approx - \int S_{\mathrm{RRL}} / S_{\rm c} ~ \mathsf{d}v$ computed by letting $S_c$ of each component equal half the total continuum flux density in the band, $S_c \approx 5.75$ Jy. We assume the continuum flux is equally split between the two core components following [Koopmans & de Bruyn, 2005] and [Combes et al., 2021] which show this to be the case and that the core components dominate the emission at least down to 1.4 GHz. Furthermore, ALMA observations measure a NE/SW flux density ratio close to one in July 2019 [Muller et al., 2021].

Fig. <ref> overlays the H144$\alpha$ spectrum with the  21 cm and OH 18 cm absorption spectra that have been rescaled by factors of $-240$ and $-30$, respectively, in flux density units and obtained from the same MALS datasets. The RRL emission appears to span a similar velocity range as these other diffuse gas tracers and likewise, is also dominated by two main velocity components.

The H144$\alpha$ spectrum detected in L band overlaid by the rescaled  and OH absorption spectra. The right-hand ordinate shows the RRL flux density normalized by the continuum.

The RRL line centroids of the NE and SW components agree within error with the OH 18 cm profiles, which have values of $-110 \pm 3$  and $6 \pm 3$ , respectively. However, the RRL centroids are significantly offset from the dominant  components at $\sim -150$  and 0 , respectively, albeit the  profile is more complex with perhaps five velocity components best-fitting the observed profile. The variation in central velocity with  could arise due to different filling factors, intrinsic distributions, and shapes of the continuum with the higher frequency tracers (OH and RRLs) finding better agreement.

The RRL line width of the NE component also agrees within error with the OH absorption width, but the widths are significantly different for the SW component, with the RRL FWHM of $63 \pm 10$  and the OH FWHM of $94.2 \pm 5$ . An estimate of the  width is close to $\sim 100$  for both components, which would be inconsistent with the RRL widths from either component. Given the smaller line width of the warmer gas traced by the RRLs in the SW, the emission may likely have a smaller filling factor, which is to say, fewer individual components contribute to the total profile. This can be expected since the SW line of sight is dominated by dense molecular gas.

While the two dominant components of H144$\alpha$ emission generally agree with the  and OH profiles, more than two velocity components are discernible in the higher S/N spectra of  and OH. For example, OH absorption has an additional component centered at $-211 \pm 3$  with a FWHM of $28 \pm 9$ ; the RRL line profile fit to the NE component encompasses some emission in this velocity range. The  absorption spectrum also shows two velocity features that are blueshifted with respect to the main $\sim -150$  peak.

In Fig. <ref>, we show the UHF band detections obtained at the $z=0.89$ absorber. We overlay the parallel hand spectra, showing that they are consistent within the noise and thus likely not a result of RFI. Integrating the Stokes-$I$ spectrum over the velocity channels from $-230$ to 55  results in $\int S_{\rm H163\alpha} \, \mathsf{d}v = 309 \pm 22$ mJy  and a maximum integrated S/N of 14. The effective transition of the final spectrum is H163$\alpha$, with a rest-frame frequency of 1504.6 MHz and observed at 797.85 MHz.

The H163$\alpha$ emission arises across a similar velocity range as the H144$\alpha$ emission, yet the peak intensities are slightly lower. The distinction of two components is less obvious in the H163$\alpha$ stack (UHF band), as compared with the H144$\alpha$ stack (L band). We fit two Gaussian components to the spectrum, fixing the central velocity of the SW component equal to that obtained from the high S/N spectrum in L band, $v_{\mathrm{cen}} = 7.7$ . The best fits are shown in Fig. <ref> and the fit parameters are listed in Table <ref>. As in the L band spectrum, we compute the integrated optical depth by assuming the UHF band continuum flux is equally split towards each core component, $S_c \approx 5.7$ Jy.

RRL stacked spectra in UHF band at $z=0.89$. The top panel compares spectra from the two parallel hands. The bottom panel is the same as in Fig. <ref> but for the UHF band spectrum.

In Fig. <ref>, we show the L and UHF band stacked spectra of  at each redshift where we obtained null results: $z=0$, $z=0.19$, and $z=2.5$. The spectral properties of the non-detections are included in Table <ref> and the integrated optical depth, $\int \tau \, \mathsf{d}v \approx - \int S_{\mathrm{RRL}} / S_{\rm c} ~ \mathsf{d}v$, is  computed by letting $S_c$ equal the total continuum flux density in the L and UHF bands, respectively. [Mohan et al., 2002] previously reported a 5$\sigma$ upper limit to the H158$\alpha$ line from the $z=0.19$ absorber of $S_{\rm H158\alpha} < 0.5$ mJy. In our L band spectrum of the $z=0.19$ stack at an effective transition of H166$\alpha$, we reach a spectral noise of 0.34 mJy, and our 3$\sigma$ upper limit to the peak line emission of $S_{\rm H166\alpha} < 1.0$ mJy is consistent with the results obtained by [Mohan et al., 2002].

Lastly, we note that carbon RRL emission becomes significantly enhanced only at frequencies $\lesssim$ 350 MHz, and thus we do not expect to detect it at the frequencies of our observations. The detected signal likely arises only from hydrogen emission given that it is coincident in velocity with the HI 21 cm and OH 18 cm absorption. A 3$\sigma$ upper limit to the carbon RRLs at $z=0.89$ in L band is $\int S_{\rm C144\alpha}\, \mathsf{d}v < 22.3$ mJy  and in UHF band is $\int S_{\rm C163\alpha}\, \mathsf{d}v < 37.5$ mJy , assuming a line width of 60 . We searched for H$\beta$ emission by stacking all available lines in both L and UHF band following the procedures described in Sec. <ref>. We report a 3$\sigma$ upper limit to the H$\beta$ emission of $\int S_{\rm H192\beta}\, \mathsf{d}v < 18.3$ mJy  and an in-band ratio of peak H192$\beta$/H144$\alpha <0.5$, where typical in-band $\beta/\alpha$ ratios are $\sim 0.2$ <cit.>.

§ PHYSICAL CONDITIONS OF IONIZED GAS IN THE Z=0.89 ABSORBER

Because  has a relatively low Galactic latitude and is behind the Inner Galaxy, $(\ell, b = 12.0^{\circ}, -5.7^{\circ})$, it has not yet been feasible to observe ionized gas tracers at UV through IR wavelengths from the $z=0.89$ galaxy  <cit.>. However, there have been some indications for the presence of ionized gas that could be attributed to the $z=0.89$ galaxy.
Firstly, the jet emission which forms an Einstein ring shows a complete turnover at $\sim300$ MHz in its spectral energy distribution (SED), which could be due to free-free absorption and would not arise from synchrotron self absorption
Secondly, jet emission in the Einstein ring (i.e., not the core emission) is strongly polarized [Pramesh Rao & Subrahmanyan, 1988, Subrahmanyan et al., 1990], indicating ionized plasma lies along the lines of sight.  However, it is debated whether the ionized gas originates in the Milky Way or the $z=0.89$ galaxy and if the blazar core components are free-free or synchrotron-self absorbed   [Jones et al., 1996, Guirado et al., 1999, Martí-Vidal et al., 2015].

Stimulated hydrogen radio recombination lines provide strong constraints on the gas volume density of electrons and volume-averaged pathlength. In the following sections we model the RRL emission to derive physical properties. We use these constraints to estimate the ionized gas mass per unit area and the ionizing photon flux.

§.§ RRL line width

The width of recombination lines as a function of frequency provides insight into the physical properties of the emitting gas. A Doppler-broadened profile has a constant width in velocity units as a function of frequency, and its Gaussian profile indicates broadening from the intrinsic gas particle motions (e.g., due to the temperature of the gas or turbulence) or from multiple emitting regions rotating at different velocities in a galaxy.
However, collisional and radiation broadening create a Lorentzian line profile and have an increasing line width towards lower frequency, thereby informing on the electron density of the gas or the incident radiation field strength, respectively.

For the NE component, widths of FWHM$_{\rm H144\alpha} = 131 \pm 14$  and FWHM$_{\rm H163\alpha} = 155 \pm 24$  are consistent within error. For the SW component, widths of FWHM$_{\rm H144\alpha} = 62.3 \pm 9.8$  and FWHM$_{\rm H163\alpha} = 80 \pm 28$  are also consistent within error. This indicates that the lines are Doppler broadened. Gaussian distributions do fit our observed line profiles reasonably well (see Figs. <ref> and <ref>).

Assuming the line widths of each component are equal at the two frequencies, the weighted average of the FWHM for the NE component is 137  and for the SW component is 64 . We then use the width at the lowest frequency to put a firm upper limit on the gas density, assuming pressure broadening dominates. Note, there is no indication to expect an extreme radiation field that would cause the line to be radiation broadened. [Salgado et al., 2017] provide the FWHM of a collisionally broadened profile, $\Delta \nu_{col} = n_e \mathsf{n}^{\gamma} \cdot 10^a / \pi $ where $\Delta \nu_{col}$ is in units of Hz, $n_e$ is in units of , $a = -7.386$, $\gamma = 4.439$ and $\mathsf{n}$ is the principal quantum level. We place an upper limit of $n_e \lesssim 7900$  for the NE component and $n_e \lesssim 3700$  for the SW component.

§.§ RRL SLED

As shown in [Shaver, 1975] and [Shaver, 1978], the flux density of an optically thin RRL of principal quantum number, $\mathsf{n}$, and frequency, $\nu$, observed in front of a significantly-brighter background radio source of flux density $S_{\rm{bkg},\nu}$ and which results from stimulated emission is given by

\begin{equation}
S_{\mathsf{n},\nu} \approx - \tau^*_{\mathsf{n},\nu} \left( b_{\mathsf{n}}\beta \right) S_{\rm{bkg},\nu}
\label{eq:Sline}
\end{equation}

where $\tau^*_{\mathsf{n},\nu}$ is the LTE RRL optical depth, $b_{\mathsf{n}}$ is the ratio of the number of hydrogen atoms in level $\mathsf{n}$ between the non-LTE and the LTE cases, and $\beta$ characterizes the effect of stimulated transitions. $b_{\mathsf{n}}$ and $\beta$, collectively referred to as departure coefficients, depend on the atomic physics of the hydrogen atom, and their product is dependent upon electron density, electron temperature, the radiation field, and pathlength. Computation of the departure coefficients has been thoroughly studied since the first observations of hydrogen RRLs <cit.>. Integrating over the line profile and expressing the LTE optical depth of the line, the integrated optical depth of an $\alpha$ transition ($\Delta \mathsf{n} = 1$) at high principal quantum numbers takes the form

\begin{eqnarray}
\label{eq:inttau}
\int \tau_{\mathsf{n}} \, \mathsf{d}v && \approx 6.13 \times 10^{-4} \,\mathrm{km\, s^{\text{-}1}} \left( b_{\mathsf{n}}\beta \right) \\ \nonumber
&&\left( \frac{EM}{10^4\, \mathrm{cm^{\text{-}6}\, pc}} \right)
\left( \frac{T_e}{10^4\, \mathrm{K}} \right)^{\text{-}2.5}
\left( \frac{\nu}{\rm{GHz}} \right)^{\text{-}1}
\end{eqnarray}

where $T_e$ is the electron temperature and $EM$ is the emission measure equal to $EM = \int n_e n_{ion} \, \mathsf{d}\ell$ for electron density $n_e$ and ion density $n_{ion}$ integrated over the pathlength $\ell$. Because the departure coefficients at each principal quantum number are strongly dependent upon density, they are key to breaking the degeneracy between density and pathlength in the emission measure.

The integrated optical depth as a function of frequency for the NE, SW, and Total components. We overlay predicted line intensities from RRL models, from 100 model combinations that fall within the reported uncertainties and chosen at random.  We assume a redshift $z=0.89$ for the conversion between observed and rest frequencies.

In Fig. <ref>, we plot the RRL SLED using the integrated optical depths corresponding to the NE and SW Gaussian components from from Table <ref>. We also plot (and use for analysis in this section), the integrated optical depth of the Total component computed from the velocity-integrated emission over -230 to 55  and setting the continuum flux density, $S_c$, equal to the flux density in each respective band, i.e. $\int \tau_{\rm{H}144\alpha} = - 0.029 \pm 0.001$  and $\int \tau_{\rm{H}163\alpha} = - 0.027 \pm 0.002 $ . While it would be ideal to model the two velocity components individually, the large errors of the Gaussian fit parameters warrant caution. The results from the Total component represent an average of the two lines-of-sight.

To model the recombination line emission, we calculated the departure coefficients for a range of electron densities and electron temperatures using the code and framework described in [Salgado et al., 2017, Salgado et al., 2017]. When computing $b_{n} \beta $ we only consider the effect of the cosmic microwave background on the level populations. The grid in density was sampled at an interval of 0.5 dex and the temperature in multiples of 100. We then interpolate over the grid values using a cubic bivariate spline for the following analysis.

We used Bayesian analysis to constrain the posterior distribution of the parameters $n_e$, $T_e$, and $\ell$, and we assume the gas is fully ionized with $n_{ion} = n_e$. The likelihood function describes the comparison of the observed integrated optical depth with the model (Eq. <ref>), assuming a Gaussian distribution function with dispersion equal to the measurement uncertainty (see Table <ref>). The model is taken to be Eq. <ref> in which the values of $\nu$, $n_e$, $T_e$, and corresponding $b_{n} \beta$ are input, and the pathlength is left as a free parameter. We used flat priors for the parameters expressed in logarithmic scale. We assumed reasonable ranges of the parameters: for density, $\log(n_e) =[-1, 6]$  for the Total component, and we used the upper limits derived from the line width for the NE and SW components, $\log(n_e) =[-1, 3.4]$ and $\log(n_e) =[-1, 3.6]$, respectively, in units of ; for pathlength, $\log(\ell) =[-9, 6]$ pc; and for temperature, we select a strict range of $\log( T_e ) =[3, 4.3]$ K, given that theory and observations substantiate temperatures of photoionized gas to have typical values of $T_e \approx 6000 - 10\,000$ K <cit.>. We also used the prior that the departure coefficients must be negative and thus the line appears in emission. To sample the posterior distribution we used an affine-invariant sampler within the package  [Foreman-Mackey et al., 2013]. We used 20 walkers, $10^5$ iterations, and verified convergence using an autocorrelation time analysis.

Corner plot showing the Total component constraints for the electron densities $n_e$ in units of , electron temperatures $T_e$ in units of K, and pathlength $\ell$ in units of pc. Contours are drawn at 68% and 95% intervals. The histograms show the marginalized distributions and the shaded region marks the 68% credible intervals.

The 2D and 1D marginal posterior distributions of the parameters of the Total component are plotted in Fig. <ref>. They show that density is constrained by the models to within 0.6 dex with a maximum a posteriori value of $n_e = 500$  and 68.3% credible interval of $130~\mathrm{cm^{-3}} \leq n_e \leq 2000$ , and that electron temperatures are not well constrained. Typical temperatures of photoionized gas have values of $T_e \approx 6000 - 10\,000$ K, with variations that are largely metallicity dependent <cit.>. The RRL modeling shows that typical temperatures are slightly more likely than cooler temperatures. The maximum a posteriori value of the volume-averaged pathlength for the Total component is $\ell = 0.025$ pc and the 68.3% credible interval is $7.9 \times 10^{-3}~\mathrm{pc} \leq \ell \leq 0.13$ pc. It is reasonable to infer that this ionized gas is not distributed in a very thin sheet of  $\ell \sim 0.025$ pc, but instead has a volume filling factor less than unity, and the emission arises from multiple, discrete clouds within the galaxy — either as ionized clouds, interface layers of molecular clouds,  regions (for further information see Sec. <ref>), etc.

Posterior probability distributions for the parameters constrained through the RRL observations: electron density $n_e$, electron temperature $T_e$, pathlength $\ell$, and emission measure $EM$, derived for the NE, SW, and Total components. The shaded regions represent the 68.3% credible intervals. In the left-most panel, the gray dashed (dotted) line marks the $n_e < 7900$  ($<$3700 ) constraint for the NE (SW) component obtained from the line width.

Constraints on physical conditions. Temperature is not well constrained within the allowed range $\log(T_e) = [3,4.3]$ K.

$\log( n_e ) $
$\log( \ell ) $
$\log( EM ) $

(cm$^{-6}$ pc)

NE       -120         2.0_-0.7^+1.0           -0.7_-1.1^+1.1      4.1_-1.2^+0.8

SW       8            3.2_-1.0^+0.4          -2.7_-0.2^+1.8       3.9_-1.0^+0.7

Total     (-230–55)      2.7_-0.6^+0.6      -1.6_-0.5^+0.7       4.9_-1.7^+0.4

In Fig. <ref>, we show the marginal posterior distributions for all three components and each of the parameters, electron density ($n_e$), electron temperature ($T_e$), pathlength ($\ell$), and emission measure ($EM$). The parameter constraints are listed in Table <ref>, with the maximum a posteriori values and 68.3% credible intervals. We do not list the temperature in Table <ref> because it is not well constrained by these measurements. We take 100 RRL models at random with properties that fall within the constraints of the 68.3% credible intervals and plot the SLEDs of these models in Fig. <ref>, in order to demonstrate the variation in the predicted integrated optical depth for different models.

§.§ Mass of ionized gas

Using the most likely values of the physical quantities in Table <ref>, we estimate the mass and mass per unit area of ionized gas in each component. These mass estimates are meant for qualitative comparisons as order of magnitude indications and do warrant caution. We assume the volume of gas is effectively described by a cylinder, as a circular region on the plane of the sky (core size) and with a distance into the plane equal to the path length. We also assume that the radio continuum emission is dominated by the NE and SW components, as these have been shown to account for $\sim99$% of the emission at 1.4 GHz [Koopmans & de Bruyn, 2005, Combes et al., 2021]. To calculate the mass of ionized gas, we expect

\begin{eqnarray}
M_{\mathrm{ion}} && \approx 1.36 m_{\mathrm{H}} \cdot n_e  \cdot (\pi r_{c}^2) \ell \nonumber \\
M_{\mathrm{ion}} && \approx 0.11 \, \mathrm{M}_{\odot}
\left( \frac{n_e}{ \mathrm{cm}^{-3} } \right)
\left( \frac{r_c}{ \mathrm{pc} } \right)^2
\left( \frac{\ell}{ \mathrm{pc} } \right).
\label{eq:m_ion}
\end{eqnarray}

where $m_{\rm H}$ is the hydrogen atom mass and $r_c$ is the radius of the radio continuum core. [Guirado et al., 1999] constrained the source size of the SW component to follow $\propto \nu^{-2.0}$ resulting in a size of 0.1 at 1 GHz, which corresponds to the value we set of $r_c = 786$ pc at $z=0.89$. For the NE component we adopt the separation of 0.05 between the two brightest emission peaks as the size, $ r_c = 393$ pc. We do not include uncertainties for $r_c$ when computing the mass; however, we note that their errors are small in comparison to the large range uncertainty of the credible intervals. Using the posterior distributions of $n_e$ and $\ell$, we estimate a total mass for the NE component of $M_{\mathrm{ion}} \approx 10^{6.0_{-0.7}^{+0.6}}$  and for the SW component, $M_{\mathrm{ion}} \approx 10^{5.5_{-0.4}^{+0.8}}$ . If we assume that the area of the Total component is a summation of the areas intercepted by the NE and SW lines-of-sight, the estimated ionized mass is $M_{\mathrm{ion}} \approx 10^{6.4_{-0.5}^{+0.7}}$ .

It is informative to also calculate the gas mass per unit area

\begin{equation}
\Sigma_{\mathrm{ion}} \approx 0.033 \, \mathrm{M}_{\odot} \,  \mathrm{pc}^{-2}
\left( \frac{n_e}{ \mathrm{cm}^{-3} } \right)
\left( \frac{\ell}{ \mathrm{pc} } \right).
\end{equation}

For the NE, SW, and Total components, we calculate $\Sigma_{\mathrm{ion}}$ to be $10^{0.3_{-0.7}^{+0.6}}$ , 10$^{-0.8_{-0.4}^{+0.8}}$ , and $10^{0.0_{-0.5}^{+0.7}}$ , respectively.

Using the surface densities of the NE, SW and Total components, we calculate an estimate for the total ionized gas mass of the galaxy within the Einstein ring of $R_{\rm g} \sim 5.3$ kpc in units of  as $\log( M_{\rm ion,g} ) \approx 8.2_{-0.7}^{+0.6}$,  $7.1_{-0.4}^{+0.8}$, and $7.9_{-0.5}^{+0.7}$, respectively. While the mass estimated from the SW is more equivalent to a true lower limit of the total ionized gas mass, the estimates from the NE and Total components are almost certainly lower limits as well since we do not trace the bulk of the ionized gas mass of the galaxy contained in the Warm Ionized Medium <cit.>.

§.§ Ionizing photon flux

We use the ionized gas emission measure to infer the ionizing photon flux. The ionizing photon rate, $Q_o$, per unit area is,

\begin{eqnarray}
\frac{ Q_o }{\mathrm{area}} && \approx EM \cdot \alpha_B \nonumber \\
\frac{ Q_o }{\mathrm{area}} && \approx 7.6 \times 10^{45} \,\mathrm{photons \, s^{-1} \, pc^{-2}}  \frac{EM}{10^3 \, \mathrm{cm^{-6} \, pc} }
\end{eqnarray}

where $\alpha_B$ is the case B recombination coefficient. Using the posterior distributions of the emission measure, we calculate ionizing photon fluxes in units of  of $\log{(Q_o / \mathrm{area})} = 47.0_{-1.2}^{+0.8}$, $46.8_{-1.0}^{+0.7}$, and $47.8_{-1.7}^{+0.4}$ for the NE, SW, and Total components, respectively. These values are about an order of magnitude or more higher than the ionizing photon flux in the disk of the Milky Way <cit.>.

While the gas mass estimates (Sec. <ref>) are likely lower limits, the ionized photon fluxes are closer to realistic values, since they are dominated by the large emission measures that we are sensitive to.

§ DISCUSSION

§.§ Star Formation Rate and ISM Properties

The posterior distributions for the emission measure and thus directly the ionizing photon flux are constrained with fairly similar 68.3% credible intervals. To estimate a star formation rate (SFR) of the galaxy, we adopt the log-average $EM$ within the credible intervals of the three components, $EM \approx 10^{4.0 \pm 0.8}$ , within a typical error of 0.8 dex, and therefore we adopt an ionizing photon flux is $Q_o / {\rm area} \approx 10^{46.9 \pm 0.8}$ .
We estimate the total ionizing photon rate for the galaxy within $R_g = 5.3$ kpc as $Q_o = \pi R_g^2 \cdot 10^{46.9 \pm 0.8} \approx 10^{54.8 \pm 0.8}$ photons s$^{-1}$. A Starburst99 model [Leitherer et al., 1999] of continuous star formation establishes the relation with the ionizing photon rate (which levels off after $\sim$50 Myr) of $\mathrm{SFR} = 1~\mathrm{M_{\odot} \, yr^{-1}} \, (Q_o / 1.4 \times 10^{53} \, \mathrm{photons \, s^{-1}})$. Using this relation, an estimated SFR for the $z=0.89$ galaxy is $\mathrm{SFR} \approx 10^{1.7 \pm 0.8}$  and the SFR per unit area is $\Sigma_{\mathrm{SFR}} \approx 10^{-0.2 \pm 0.8}$  kpc$^{-2}$. A galaxy on the main sequence at $z=0.89$ with SFR $\sim$50  has a typical stellar mass of $\sim10^{11}$  <cit.>. Current lensing models estimate the total mass within the Einstein ring as $M_E \approx 4 \times 10^{11}$  [Muller et al., 2020] and therefore a stellar mass of $\sim 8 \times 10^{10}$  given a typical mass to stellar light ratio of $\sim5$ [Treu & Koopmans, 2004]. This is likely a main-sequence galaxy.

In Sec. <ref> we estimated the ionized gas mass per unit area of the NE component to be $\Sigma_{\rm ion} \sim 2.1$ . For comparison, the  column density is estimated to be $N_{\mathrm{H\,I}} \approx 5 \times 10^{21}$ cm$^{-2}$ assuming half of the continuum flux comes from the NE component [Chengalur et al., 1999, Combes et al., 2021]. Assuming the average particle mass is $\approx 1.36m_{\rm H}$, the atomic gas mass per unit area is $\Sigma_{\mathrm{H\,I}} \approx 50$ . We use the OH column density of $N_{\rm OH} \approx 1.5 \times 10^{15}$ cm$^{-2}$ [Gupta et al., 2021] and assume an abundance of $10^{-7}$ [Balashev et al., 2021] given that this line-of-sight has properties of a diffuse cloud [Muller et al., 2011] in order to estimate a molecular gas column density of $N_{\rm H_2} \sim 1.5 \times 10^{22}$ cm$^{-2}$. This H$_2$ column density is higher than the $N_{\rm H_2} \sim 10^{21}$ cm$^{-2}$ derived from H$_2$O absorption, with a smaller continuum cross section at higher frequencies [Muller et al., 2014]. Assuming the average particle mass is $\approx 2m_{\rm H}$, the molecular gas mass per unit area is $\Sigma_{\rm H_2} \approx 240$ .

With these estimates, the neutral gas mass per unit area is $\Sigma_{\rm HI + H_2} \sim 290$ .
The ionized gas mass detected in RRLs of the NE component is a small fraction ($\sim$0.7%) of the total gas mass.
The Kennicutt-Schmidt law [Kennicutt, 1998] long establishes a direct correlation between the surface densities of the neutral gas mass $\Sigma_{\rm HI + H_2}$ and the star formation rate $\Sigma_{\rm SFR}$ in galaxies on spatial scales 300–500 pc or more [Schruba et al., 2010, Kruijssen & Longmore, 2014]. It suggests that a region in a galaxy with $\Sigma_{\rm SFR} \sim 0.6$  kpc$^{-2}$ has a neutral gas mass per unit area of $\Sigma_{\rm HI + H_2} \sim 270$  <cit.>. Our measured neutral gas mass agrees to within 10% of the expected value, although our measured SFR surface density has a large uncertainty.

Although the estimated densities of the NE, SW and Total components are consistent within the errors, the most likely density of the SW component, $n_e \approx 1600$ , is typical of young, compact  regions. If we assume the  regions are $r_{\rm{H\,II}} \approx 2$ pc in size, then the covering fraction through this cross section of the galaxy is $f_c \approx \ell_{\rm SW}/ (4/3 \cdot r_{\rm{H\,II}}) \approx 7.5 \times 10^{-4}$. The total number of  regions, $N_{\rm{H\,II}}$, is calculated from the surface density estimates, $\pi r_C^2 = N_{\rm{H\,II}} / f_c \cdot (\pi\, r_{\rm{H\,II}}^2) $, which computes to $N_{\rm{H\,II}} \sim 116$, where we recall from Sec. <ref> that $r_C = 786$ pc for the SW component. Equivalently, assuming all  regions have $n_e \approx 1600$  and $r_{\rm{H\,II}} \approx 2$ pc, the ionized gas mass per  region is estimated to be $M_{\rm H\,II} = 1860$ , and from the total ionized mass of the SW component, $M_{\rm SW} = N_{\rm H\,II} M_{\rm H\,II}$, also results in $N_{\rm H\,II} \approx (2.1 \times 10^5 \, \rm{M}_{\odot}) / (1860 \,  \rm{M}_{\odot}) \sim 113$. However, this calculation warrants caution because the number of  regions changes dramatically depending on the assumed size. We set $r_{\rm H\,II} \approx 2$ pc because it is the typical size of young, massive star clusters [Ryon et al., 2017] that are expected to dominate star formation in galaxies of this epoch, i.e., with high gas surface densities and star formation rates [Kruijssen, 2012].

§.§ The $n_e - \Sigma_{\rm SFR}$ relation

The star formation rate per unit area is shown to be correlated with the electron density of ionized gas in the region [Shimakawa et al., 2015, Herrera-Camus et al., 2016, Jiang et al., 2019]. Assuming  regions thermalize with the ISM [Gutiérrez & Beckman, 2010], i.e., take on a pressure balance with others phases, then the volume-average density of the ionized gas (along with fairly consistent ionized gas temperatures) indicates the thermal pressure of the medium [Jiang et al., 2019, Barnes et al., 2021]. This results in a $P-\Sigma_{\mathrm{SFR}}$ relation that serves as an important test-bed for pressure-regulated, feedback-modulated star formation <cit.>. Using doublet ratios of [SII] and [OII] in the optical [Kewley et al., 2019], [Shimakawa et al., 2015] and [Jiang et al., 2019] find a relation of $\Sigma_{\mathrm{SFR}} \propto n_e^{1.7 \pm 0.3}$ in a sample of $z \sim 1-3$ starburst galaxies. In the far IR, the ratio of [NII] fine structure lines from a sample of nearby normal galaxies and (ultra) luminous infrared galaxies ((U)LIRGs; $L_{\mathrm{IR}} \ge 10^{11}$ ) establish $\Sigma_{\mathrm{SFR}} \propto n_e^{1.5}$ [Herrera-Camus et al., 2016].

For $\Sigma_{\mathrm{SFR}} \sim 0.6$  kpc$^{-2}$, the electron density predicted by the [Shimakawa et al., 2015] relation is $\sim 110$  and by the [Herrera-Camus et al., 2016] relation is $\sim 200$ . These density estimates agree well with the NE component and lie on the cusp of the 68% credible intervals of the SW and Total components.
If the RRL emission is tracing the thermal properties of the diffuse medium in the galaxy's disk, higher pressures and densities are typically found at smaller galactic radii, i.e., the SW component, than at larger galactic radii, i.e., the NE component.
For example, [Gutiérrez & Beckman, 2010] measured the electron density to increase towards smaller galactic radii, $r$, as $\left< n_e \right> = \left< n_e \right>_o \exp{(-r / R_g)}$, where $R_g$ is the scale length of the disk at which star-formation and density drops off and $\left< n_e \right>_o$ is inner most density.  With $R_g = 5.3$ kpc and $\left< n_e \right> = 100$  of the NE component, at $r=2.4$ kpc the expected density is $n_e \sim 160$ . While the contrast in density at the two radii is not as extreme as the best fit densities of the components indicate, the general trend is consistent and a density of 160  does fall within the 68% credible interval of the SW component. The smaller pathlength in comparison to the NE component might then indicate a smaller covering fraction of the overall larger cross section of this line-of-sight, i.e. only within or close to the spiral arm.

§.§ Radio Continuum SED

The radio continuum emission from  is complex. The lensing is achromatic and contains a small-scale core-jet structure with regions of different spectral indices and opacities. In addition, the source is variable on hourly to yearly timescales [Pramesh Rao & Subrahmanyan, 1988, van Ommen et al., 1995, Lovell et al., 1996, Lovell et al., 1998, Garrett et al., 1997, Jin et al., 2003, Martí-Vidal et al., 2013, Allison et al., 2017], which makes it difficult to compare observations and model the radio continuum emission [Muller et al., 2020].

Ionized gas that is detectable in RRLs and has a large covering fraction would also emit free-free emission and when it becomes optically thick, would absorb any background radio continuum. A reliable model and knowledge of the spatially-resolved radio SED could independently constrain the physical conditions of the gas through free-free absorption. The frequency at which radio emission becomes optically thick to free-free absorption is defined by $\tau_\nu = 6.67 \times 10^{-2} \, EM \, T_e^{-1.323} (\nu / \mathrm{GHz})^{-2.118}$ <cit.>.

Compilation of continuum flux density measurements from the bright and highly variable (factors of $\sim$1.5 on weeks and years timescales) . The gray shaded region encompasses the 1$\sigma$ confidence region of a best-fitting power-law that is attenuated by free-free absorption. We also include the continuum SEDS (with colors corresponding to the components on the left hand plot) from a fixed power-law and that is attenuated by ionized gas with properties constrained by the RRL models. We assume a redshift $z=0.89$ for the conversion between observed and rest frequencies.

We collected radio continuum measurements of  at $\nu \lesssim 22$ GHz <cit.> and plot these in Fig. <ref>. We fit the SED of  as a power-law index with an external screen of free-free absorption, using the functional form $S_{\nu} = S_o (\frac{\nu}{\nu_o})^{\alpha} \exp (- \tau_{\nu} ) $ with $\tau_{\nu} = \tau_o (\frac{\nu}{\nu_o})^{-2.118}$. Setting $\nu_o = 40$ GHz with respect to observed frequencies, the best fit parameters with 1$\sigma$ uncertainties are $S_o = 4.8 \pm 0.4$ Jy, $\alpha = -0.24 \pm 0.03$ and $\tau_o = (1.1 \pm 0.2) \times 10^{-5}$.

In Fig. <ref>, we also plot how ionized gas that has physical properties constrained by the RRL models — using the same selection of models presented in Fig. <ref> — would attenuate a power-law continuum SED (normalized by $S_o$ and $\alpha$ of our fit). The emission from the RRLs we detect at $z=0.89$ would result in free-free absorption at lower frequencies than is observed in .
For the $z=0.89$ galaxy to cause the free-free absorption, volume-average pathlengths a factor of 5 larger are needed. A smaller filling factor and larger pathlength intercepting the radio emission is not unreasonable.

It would still be possible for ionized gas in the environment of the blazar at $z=2.5$ to create the absorption in the radio SED. Even though we do not detect hydrogen RRL emission at $z=2.5$, ionized gas with different properties, for example higher densities, could be present and still be consistent with the RRL constraints.
We also note that ionized gas in the Milky Way with $EM$s that result in a turnover at the observed frequencies (especially, on $<1$ scales) would be observable in hydrogen RRL emission at $z=0$, and we do not detect any RRL emission from the Milky Way.

§ CONCLUSIONS

We used MALS observations to detect RRL emission in the spectrum of the radio blazar . The RRL emission is observed at $z=0.89$ from a galaxy that lies along the line of sight and strongly lenses . This is the second detection of RRLs outside of the local universe (i.e., at $z \geq 0.076$) and the first clearly associated with hydrogen <cit.>. We detect H144$\alpha$ by stacking 17 RRLs covered by the L band (856–1712 MHz) with a S/N of 21 (see Fig. <ref>), and we detect H163$\alpha$ by stacking 27 lines in the UHF band (544–1088 MHz) with a S/N of 14 (see Fig. <ref>). Emission from the H144$\alpha$ line is consistent over two separate observations, when comparing the spectra in parallel hand polarizations, and is robust against additional spectral stacking verification methods. Like the  21 cm and OH 18 cm absorption spectra (see Fig. <ref>), the H144$\alpha$ and H163$\alpha$ emission profiles span $\sim$250  in velocity and are dominated by two velocity components associated with two physically distinct regions of the galaxy, the NE and SW lines of sight. We do not detect RRL emission in either band intrinsic to  ($z=2.5$), from the $z=0.19$ absorption system along this line of sight, or from the Milky Way (see Fig. <ref>).

Hydrogen RRL emission typically arises from fully ionized gas and only stimulated emission is observable outside of the local universe. The maser-like properties of stimulated emission enable the RRL SLED to constrain the density and pathlength of the ionized gas (see Table <ref> and Fig. <ref>). Considering the total integrated line intensity, referred to as the Total component, we used a Bayesian analysis to constrain the electron density of the gas $\log( n_e ) = 2.6 \pm 0.6$  and a volume-averaged pathlength of $\log(\ell) = -1.6_{-0.5}^{+0.7}$ pc, which likely has a non-unity filling factor. Analyzed separately, the NE line-of-sight appears to harbor less dense gas with $\log( n_e ) = 2.0_{-0.7}^{+1.9}$  and $\log( \ell ) = -0.7 \pm 1.1 $ pc, and the SW line-of-sight appears to intercept dense gas that is more typical of  regions with $\log( n_e ) = 3.2_{-1.0}^{+0.4}$  and $\log( \ell ) = -2.7_{-0.2}^{+1.8}$ pc. These scenarios are consistent with the NE line-of-sight passing through diffuse clouds at a larger galactic radius, and the SW component directly intercepting a spiral arm, as has previously been determined.

The RRL components measure an ionizing photon flux of $Q_o / \mathrm{area} \approx 10^{46 \pm 0.8}$  and star formation rate surface density of $\Sigma_{\rm SFR} \sim 10^{-0.2 \pm 0.8}$  kpc$^{-2}$. Taken over the $z=0.89$ galaxy within $R_g \sim 5.3$ kpc, the ionizing photon rate of $Q_o \sim 10^{54.8}$ photons s$^{-1}$ yields an average star formation rate of SFR $\sim 50$ . Despite the plethora of molecular species observed, the ionized gas content and SFR have not been previously measured for this source, largely due to the highly reddened nature of  at optical and NIR wavelengths. In comparing the SFR and the galaxy's mass (from lensing), the $z=0.89$ system is likely on the main sequence.

The ionized gas mass per unit area of the diffuse NE component as measured by the RRL emission is $\Sigma_{\rm ion} \approx 2.1$ , in comparison with gas masses of $\Sigma_{\mathrm{H\,I}} \approx 50$  and $\Sigma_{\rm H_2} \approx 240$  (via OH) estimated from only the MALS observations. Given our estimated SFR, the +H$_2$ gas mass surface density is close to the gas content predicted by the Kennicutt-Schmidt law. Our measured electron densities also match reasonably well with the $n_e - \Sigma_{\mathrm{SFR}}$ relation determined from optical and FIR line ratios.

is the first source investigated with MALS, and the detection of RRLs in the source is promising for the remaining $\sim$500 targets of the survey. With the first hydrogen RRL detection that breaks the redshift-barrier, we show that this tracer can be an important tool for investigating (a) the electron density (thermal pressure) of ionized gas in the ISM of galaxies (and the $n_e - \Sigma_{\rm SFR}$ relation) and (b) the SED of AGN, thus eventually AGN evolution. We have also demonstrated the unique science that can be achieved through  21 cm, OH 18 cm, and RRL measurements that are simultaneously observed in the MALS survey. The ionized gas properties in the $z=0.89$ galaxy will be substantially improved through RRL observations at higher and lower radio frequencies and at higher ($<$1) spatial resolutions which can separate the two main (velocity) components of emission. The new science afforded by high-redshift RRL studies is accessible with on-going wide-bandwidth spectral line surveys and will be explored in unprecedented capacities with future facilities such as the next generation Very Large Array <cit.> and the SKA [Carilli, 2015].

The authors acknowledge and appreciate the efforts and input of the anonymous reviewer of this article. The authors thank Peter Shaver for comments on the article and for the inspiration and motivation to carry through with this research.

SAB was supported by RSF grant 18-12-00301. The MeerKAT telescope is operated by the South African Radio Astronomy Observatory, which is a facility of the National Research Foundation, an agency of the Department of Science and Innovation.

The MeerKAT data were processed using the MALS computing facility at IUCAA (https://mals.iucaa.in/releases)

The National Radio Astronomy Observatory is a facility of the National Science Foundation operated under cooperative agreement by Associated Universities, Inc. This research has made use of the NASA/IPAC Extragalactic Database (NED), which is funded by the National Aeronautics and Space Administration and operated by the California Institute of Technology.

ARTIP [Gupta et al., 2021],
Astropy [The Astropy Collaboration, 2018, The Astropy Collaboration et al., 2022],
CASA [McMullin et al., 2007, The CASA Team et al., 2022],
ChainConsumer [Hinton, 2016],
CRRLpy [Salas et al., 2016],
emcee [Foreman-Mackey et al., 2013],
Matplotlib [Hunter, 2007],
and NumPy [Harris et al., 2020]

[Allison et al., 2017]
Allison, J. R., Moss, V. A., Macquart, J. P., et al. 2017, MNRAS, 465, 4450,

[Allison et al., 2022]
Allison, J. R., Sadler, E. M., Amaral, A. D., et al. 2022, PASA, 39, 10,

[Anantharamaiah, 1985]
Anantharamaiah, K. R. 1985, JApA, 6, 177

[Balashev et al., 2021]
Balashev, S. A., Gupta, N., & Kosenko, D. N. 2021, MNRAS, 504, 3797,

[Balser, 2006]
Balser, D. S. 2006, AJ, 132, 2326, 10.1086/508515

[Barnes et al., 2021]
Barnes, A. T., Glover, S. C. O., Kreckel, K., et al. 2021, MNRAS, 508, 5362,

[Bell & Seaquist, 1978]
Bell, M. B., & Seaquist, E. R. 1978, ApJ, 223, 378, 10.1086/156272

[Bell & Seaquist, 1980]
—. 1980, ApJ, 238, 818, 10.1017/CBO9781107415324.004

[Bell et al., 1984]
Bell, M. B., Seaquist, E. R., Mebold, U., Reif, K., & Shaver, P. A. 1984,
A&A, 130, 1.

[Blyth et al., 2015]
Blyth, S., van der Hulst, T. M., Verheijen, M. A. W., et al. 2015, in
Proceedings of Advancing Astrophysics with the Square Kilometre
Array, eprint: arXiv:1501.01295, 128.

[Carilli, 2015]
Carilli, C. 2015, in AASKA, Giardini Naxos, Italy, 171.

[Chengalur et al., 1999]
Chengalur, J. N., de Bruyn, A. G., & Narasimha, D. 1999, A&A, 343, L79.
<https://ui.adsabs.harvard.edu/abs/1999A   A...343L..79C/abstract>

[Churchwell & Shaver, 1979]
Churchwell, E., & Shaver, P. A. 1979, A&A, 77, 316.
<https://ui.adsabs.harvard.edu/abs/1979A   A....77..316C/abstract>

[Combes et al., 2021]
Combes, F., Gupta, N., Muller, S., et al. 2021, A&A, 648, 116,

[Courbin et al., 1998]
Courbin, F., Lidman, C., Frye, B. L., et al. 1998, ApJ, 499, L119,

[Courbin et al., 2002]
Courbin, F., Meylan, G., Kneib, J. P., & Lidman, C. 2002, ApJ, 575, 95,

[Djorgovski et al., 1992]
Djorgovski, S., Meylan, G., Klemola, A., et al. 1992, MNRAS, 257, 240,

[Emig, 2021]
Emig, K. L. 2021, PhD thesis, Universiteit Leiden.

[Emig et al., 2020]
Emig, K. L., Salas, P., de Gasperin, F., et al. 2020, A&A, 634, 138,

[Emig et al., 2019]
—. 2019, A&A, 622, 7, 10.1051/0004-6361/201834052

[Emig et al., 2022]
Emig, K. L., White, G. J., Salas, P., et al. 2022, A&A, 664, 88,

[Foreman-Mackey et al., 2013]
Foreman-Mackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, PASP, 125, 306,

[Garrett et al., 1997]
Garrett, M. A., Nair, S., Porcas, R. W., & Patnaik, A. R. 1997, VA, 41, 281,

[Guirado et al., 1999]
Guirado, J. C., Jones, D. L., Lara, L., et al. 1999, A&A, 346, 392.
<https://ui.adsabs.harvard.edu/abs/1999A   A...346..392G/abstract>

[Gupta et al., 2016]
Gupta, N., Srianand, R., Baan, W., et al. 2016, in MKS Conf, 14.

[Gupta et al., 2021]
Gupta, N., Jagannathan, P., Srianand, R., et al. 2021, ApJ, 907, 11,

[Gupta et al., 2022]
Gupta, N., Shukla, G., Srianand, R., et al. 2022, ApJ, 929, 108,

[Gutiérrez & Beckman, 2010]
Gutiérrez, L., & Beckman, J. E. 2010, ApJ, 710, L44,

[Harris et al., 2020]
Harris, C. R., Millman, K. J., van der Walt, S. J., et al. 2020, Natur, 585,
357, 10.1038/s41586-020-2649-2

[Henkel et al., 2008]
Henkel, C., Braatz, J. A., Menten, K. M., & Ott, J. 2008, A&A, 485, 451,

[Herrera-Camus et al., 2016]
Herrera-Camus, R., Bolatto, A., Smith, J. D., et al. 2016, ApJ, 826, 175,

[Hinton, 2016]
Hinton, S. R. 2016, JOSS, 1, 00045, 10.21105/joss.00045

[Hummer & Storey, 1987]
Hummer, D. G., & Storey, P. J. 1987, MNRAS, 224, 801

[Hunter, 2007]
Hunter, J. D. 2007, Computing in Science & Engineering, 9, 90

[Intema et al., 2017]
Intema, H. T., Jagannathan, P., Mooley, K. P., & Frail, D. A. 2017, A&A, 598,
A78, 10.1051/0004-6361/201628536

[Jauncey et al., 1991]
Jauncey, D. L., Reynolds, J. E., Tzioumis, A. K., et al. 1991, Natur, 352,
132, 10.1038/352132a0

[Jiang et al., 2019]
Jiang, T., Malhotra, S., Yang, H., & Rhoads, J. E. 2019, ApJ, 872, 146,

[Jin et al., 2003]
Jin, C., Garrett, M. A., Nair, S., et al. 2003, MNRAS, 340, 1309,

[Jonas & MeerKAT Team, 2016]
Jonas, J., & MeerKAT Team. 2016, in MeerKAT Science: On the Pathway
to the SKA.

[Jones et al., 1996]
Jones, D. L., Preston, R. A., Murphy, D. W., et al. 1996, ApJ, 470, L23,

[Kado-Fong et al., 2020]
Kado-Fong, E., Kim, J.-G., Ostriker, E. C., & Kim, C.-G. 2020, ApJ, 897, 143,

[Kennicutt, 1998]
Kennicutt, Jr., R. C. 1998, ApJ, 498, 541, 10.1086/305588

[Kennicutt & De Los Reyes, 2021]
Kennicutt, Jr., R. C., & De Los Reyes, M. A. C. 2021, ApJ, 908, 61,

[Kewley et al., 2019]
Kewley, L. J., Nicholls, D. C., & Sutherland, R. S. 2019, ARA&A, 57, 511,

[Kim et al., 2013]
Kim, C.-G., Ostriker, E. C., & Kim, W.-T. 2013, ApJ, 776, 1,

[Koopmans & de Bruyn, 2005]
Koopmans, L. V. E., & de Bruyn, A. G. 2005, MNRAS, 360, L6,

[Kruijssen, 2012]
Kruijssen, J. M. D. 2012, MNRAS, 426, 3008,

[Kruijssen & Longmore, 2014]
Kruijssen, J. M. D., & Longmore, S. N. 2014, MNRAS, 439, 3239,

[Leitherer et al., 1999]
Leitherer, C., Schaerer, D., Goldader, J. D., et al. 1999, ApJS, 123, 3,

[Lidman et al., 1999]
Lidman, C., Courbin, F., Meylan, G., et al. 1999, ApJ, 514, L57,

[Lovell et al., 1998]
Lovell, J. E. J., Jauncey, D. L., Reynolds, J. E., et al. 1998, ApJ, 508,
L51, 10.1086/311723

[Lovell et al., 1996]
Lovell, J. E. J., Reynolds, J. E., Jauncey, D. L., et al. 1996, ApJ, 472, L5,

[Marti-Vidal & Muller, 2019]
Marti-Vidal, I., & Muller, S. 2019, A&A, 621, A18,

[Martí-Vidal et al., 2015]
Martí-Vidal, I., Muller, S., Vlemmings, W., Horellou, C., & Aalto, S.
2015, Sci, 348, 311, 10.1126/science.aaa1784

[Martí-Vidal et al., 2013]
Martí-Vidal, I., Muller, S., Combes, F., et al. 2013, A&A, 558, A123,

[McMullin et al., 2007]
McMullin, J. P., Waters, B., Schiebel, D., Young, W., & Golap, K. 2007, in
Astronomical Data Analysis Software and Systems XVI ASP
Conference Series, Vol. 376 (Astronomical Society of the Pacific), 127

[Mohan et al., 2002]
Mohan, N. R., Anantharamaiah, K. R., & Goss, W. M. 2002, in The Universe at
Low Radio Frequencies, Proceedings of IAU Symposium, ed.
A. Pramesh Rao, G. Swarup, & Gopal-Krishna, Vol. 199 (Pune, India:
Cambridge University Press), 116.

[Morganti & Oosterloo, 2018]
Morganti, R., & Oosterloo, T. 2018, AARv, 26, 4,

[Muller & Guélin, 2008]
Muller, S., & Guélin, M. 2008, A&A, 491, 739,

[Muller et al., 2020]
Muller, S., Jaswanth, S., Horellou, C., & Martí-Vidal, I. 2020, A&A, 641,
L2, 10.1051/0004-6361/202038978

[Muller et al., 2021]
Muller, S., Ubachs, W., Menten, K. M., Henkel, C., & Kanekar, N. 2021, A&A,
652, A5, 10.1051/0004-6361/202140531

[Muller et al., 2011]
Muller, S., Beelen, A., Guélin, M., et al. 2011, A&A, 535, A103,

[Muller et al., 2013]
Muller, S., Beelen, A., Black, J. H., et al. 2013, A&A, 551, A109,

[Muller et al., 2014]
Muller, S., Combes, F., Guélin, M., et al. 2014, A&A, 566, A112,

[Murphy et al., 2018]
Murphy, E. J., Bolatto, A., Chatterjee, S., et al. 2018, in ASP
Conference Series, Vol. 517, 3.

[Nair et al., 1993]
Nair, S., Narasimha, D., & Rao, A. P. 1993, ApJ, 407, 46,

[Oonk et al., 2017]
Oonk, J. B. R., van Weeren, R. J., Salas, P., et al. 2017, MNRAS, 465, 1066,

[Ostriker & Kim, 2022]
Ostriker, E. C., & Kim, C.-G. 2022, ApJ.

[Patnaik et al., 1993]
Patnaik, A. R., Muxlow, T. W. B., & Jauncey, D. L. 1993, in LIACo, Vol. 31,

[Pramesh Rao & Subrahmanyan, 1988]
Pramesh Rao, A., & Subrahmanyan, R. 1988, MNRAS, 231, 229,

[Prozesky & Smits, 2018]
Prozesky, A., & Smits, D. P. 2018, MNRAS, 478, 2766,

[Roshi & Anantharamaiah, 2000]
Roshi, D. A., & Anantharamaiah, K. R. 2000, ApJ, 535, 231,

[Ryon et al., 2017]
Ryon, J. E., Gallagher, J. S., Smith, L. J., et al. 2017, ApJ, 841, 92

[Salas et al., 2016]
Salas, P., Morabito, L., Salgado, F., Oonk, J. B. R., & Tielens, A. G. G. M.
2016, CRRLpy: First Pre-release.

[Salas et al., 2017]
Salas, P., Oonk, J. B. R., van Weeren, R. J., et al. 2017, MNRAS, 467, 2274,

[Salgado et al., 2017]
Salgado, F., Morabito, L. K., Oonk, J. B. R., et al. 2017a, ApJ,
837, 141, 10.3847/1538-4357/aa5d9e

[Salgado et al., 2017]
—. 2017b, ApJ, 837, 142, 10.3847/1538-4357/aa5d9a

[Sault, 1994]
Sault, R. J. 1994, A&AS, 107, 55.

[Schreiber et al., 2015]
Schreiber, C., Pannella, M., Elbaz, D., et al. 2015, A&A, 575, A74,

[Schruba et al., 2010]
Schruba, A., Leroy, A. K., Walter, F., Sandstrom, K., & Rosolowsky, E. 2010,
ApJ, 722, 1699, 10.1088/0004-637X/722/2/1699

[Schulz et al., 2015]
Schulz, A., Henkel, C., Menten, K. M., et al. 2015, A&A, 574, A108,

[Seaquist & Bell, 1977]
Seaquist, E. R., & Bell, M. B. 1977, A&A, 60, L1.

[Shaver, 1975]
Shaver, P. A. 1975, Pra, 5, 1, 10.1007/BF02875147

[Shaver, 1978]
—. 1978, A&A, 68, 97.

[Shaver et al., 1977]
Shaver, P. A., Churchwell, E., & Rots, A. H. 1977, A&A, 55, 435.

[Shaver et al., 1978]
Shaver, P. A., Churchwell, E., & Walmsley, C. M. 1978, A&A, 64, 1

[Shaver et al., 1983]
Shaver, P. A., McGeem, R. X., Newton, M. L., et al. 1983, MNRAS, 204, 53,

[Shimakawa et al., 2015]
Shimakawa, R., Kodama, T., Steidel, C. C., et al. 2015, MNRAS, 451, 1284,

[Subrahmanyan et al., 1990]
Subrahmanyan, R., Narasimha, D., Pramesh Rao, A., & Swarup, G. 1990, MNRAS,
246, 263.

[Tercero et al., 2020]
Tercero, B., Cernicharo, J., Cuadrado, S., de Vicente, P., & Guélin, M.
2020, A&A, 636, L7, 10.1051/0004-6361/202037837

[The Astropy Collaboration>, 2018]
The Astropy Collaboration. 2018, AJ, 156, 123

[The Astropy Collaboration> et al., 2022]
The Astropy Collaboration, Price-Whelan, A. M., Lim, P. L., et al. 2022,
ApJ, 935, 167, 10.3847/1538-4357/ac7c74

[The CASA Team> et al., 2022]
The CASA Team, Bean, B., Bhatnagar, S., et al. 2022, PASP, eprint
arXiv:2210.02276, 20, 10.1088/1538-3873/ac9642

[Tielens, 2005]
Tielens, A. G. G. M. 2005, The Physics and Chemistry of the Interstellar
Medium (Cambridge University Press).

[Treu & Koopmans, 2004]
Treu, T., & Koopmans, L. V. E. 2004, ApJ, 611, 739, 10.1086/422245

[van Ommen et al., 1995]
van Ommen, T. D., Jones, D. L., Preston, R. A., & Jauncey, D. L. 1995, ApJ,
444, 561, 10.1086/175630

[Verheijen et al., 2001]
Verheijen, M. A. W., Carilli, C. L., & Yun, M. S. 2001, in ASPC, Vol. 240,
San Francisco, 69.

[Wenger et al., 2019]
Wenger, T. V., Balser, D. S., Anderson, L. D., & Bania, T. M. 2019, ApJ, 887,

[Wiklind & Combes, 1996]
Wiklind, T., & Combes, F. 1996, Natur, 379, 139, 10.1038/379139a0