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

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. The RRL emission originates in a galaxy that intercepts and strongly lenses the radio blazar PKS 1830−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α (and H163α) transitions at observed frequencies of 1156 (798) MHz by stacking 17 (27) RRLs with 21σ (14σ) significance. The RRL emission contains two main velocity components and is coincident in velocity with H i 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 path length (ℓ) of the ionized gas. We determine log(ne)=2.0−0.7+1.0 cm−3 and log(ℓ)=−0.7−1.1+1.1 pc toward the northeast (NE) lensed image, likely tracing the diffuse thermal phase of the ionized ISM in a thin disk. Toward the southwest (SW) lensed image, we determine log(ne)=3.2−1.0+0.4 cm−3 and log(ℓ)=−2.7−0.2+1.8 pc, tracing gas that is more reminiscent of H scii regions. We estimate a star formation (surface density) rate of ΣSFR ∼ 0.6 M ⊙ yr−1 kpc−2 or SFR ∼ 50 M ⊙ yr−1, consistent with a star-forming main-sequence galaxy of M ⋆ ∼ 1011 M ⊙. 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 ν rest 10 GHz in extragalactic sources probe gas with relatively low density that can be stimulated by radio continuum (for an overview, see Emig 2021), 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 (e.g., Shaver 1975;Salgado et al. 2017a;Oonk et al. 2017).Stimulated emission has the added benefit that its intensity is proportional to the background continuum intensity at a given frequency, S RRL ∝ S 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 ( 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 (ν) as S ∝ D −2 ν −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 ν < 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 ≈ 150 cm −3 and pathlength ≈ 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 τ 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α and H99α 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 RRLs at which the stimulated-only emission of gas with a given density peaks.The thick black region marks the maximum n in the peak of emission, for ionized gas temperatures 5000 K ≤ Te ≤ 12 000 K, and the thin black region indicates the 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 ≤ z ≤ 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.
non-thermal emission (for an overview, see Emig 2021) 1 .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-ofsight to 3C 190.However, they could not clearly discern whether the emission was from carbon or hydrogen since the lines are separated by 150 km s −1 and have comparable intensities at those frequencies in the Milky Way (e.g., Anantharamaiah 1985).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).
Current large-bandwidth spectral line surveys, such as the MeerKAT Absorption Line Survey (MALS; Gupta et al. 2016), the First Large Absorption line Survey (FLASH; Allison et al. 2022), the Search for HI Ab-sorption with APERTIF (SHARP; e.g., see Morganti & Oosterloo 2018) -and in the future with the Square Kilometer Array (SKA; e.g., Blyth et al. 2015) -can observe tens of RRLs simultaneously, opening the way for ionized gas studies with optimum sensitivity to gas with electron densities of 1 cm −3 n e 10 4 cm −3 (e.g., see Fig. 1).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 (σ ∼ 0.6 mJy beam −1 per 6 km s −1 channel) of H i 21 cm and OH 18 cm absorption lines at 0 z 2 (Gupta et al. 2016).MALS is observing ∼500 pointings centered on the brightest (S 1 GHz > 0.2 Jy) radio sources at declination δ +30 • (see Gupta et al. 2022, for the survey footprint), 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 τ RRL = (0.08 − 1) × 10 −3 through line stacking, reaching the range of optical depths observed from ionized gas in the Milky Way (e.g., Roshi & Anantharamaiah 2000).For an electron temperature of 8000 K, these optical depths convert into emission measures 10 2.6 EM/cm −6 pc 10 4.7 for densities 1 cm −3 n e 10 4 cm −3 (see Fig. 1).
The first science verification observations of MALS (Gupta et al. 2021;Combes et al. 2021) focused on the bright (S 1 GHz ≈ 11 Jy) z = 2.507 (Lidman et al. 1999) blazar, PKS 1830−211 (referred to hereafter as PKS 1830).PKS 1830 has a radio spectral index of α 1−15 GHz ≈ −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).PKS 1830 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 (e.g., Jauncey et al. 1991).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).PKS 1830 is known to be variable, by up to a factor of two in radio, on timescales of hours to years (e.g., Pramesh Rao & Subrahmanyan 1988;Martí-Vidal et al. 2013;Allison et al. 2017;Marti-Vidal & Muller 2019), 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 ∼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 PKS 1830.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) (e.g., Wiklind & Combes 1996;Muller et al. 2011;Tercero et al. 2020).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 wellconstrained lensing model (Nair et al. 1993;Koopmans & de Bruyn 2005;Muller et al. 2020;Combes et al. 2021), the kinematics (v rot ∼ 260 km s −1 ) and orientation also suggest a nearly face-on (i ∼ 25 • ) barred-spiral galaxy of ∼10 11 M .The SW image of the blazar core passes through the galaxy at a radius of R g ∼ 2.4 kpc and a central velocity (from spatially unresolved emission) at v cen ∼ 0 km s −1 , and the NE image passes through at R g ∼ 5.3 kpc and v cen ∼ −150 km s −1 .The SW image intercepts a spiral arm of the z = 0.89 galaxy and traces dense (n H2 ∼ 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 H i 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 PKS 1830 at z = 0.19 has been seen only in H i 21 cm absorption so far (Lovell et al. 1996;Allison et al. 2017;Gupta et al. 2021).
MALS observations of PKS 1830 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 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. 2, we describe the methods used to process the data.In Sec. 3, 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. 4, 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. 5 and conclude in Sec. 6.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 ΛCDM cosmology with Ω m = 0.29, Ω Λ = 0.71, and H o = 70 km s −1 Mpc −1 , for which 1 ∼ 7.860 kpc at z = 0.88582.

Observations and Data Reduction
PKS 1830 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 PKS 1830 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 km s −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 km s −1 at the center of the band, i.e., 815.992MHz.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 PKS 1830 is a bonafide VLA gain calibrator at this spatial resolution, a separate gain calibrator was not observed.The total on-source times on PKS 1830 are: 40 min (L band Night1), 90 min (UHF band) and 90 min (L band Night2).
All MALS datasets have been processed using ARTIP, the Automated Radio Telescope Imaging Pipeline (Gupta et al. 2021), a Python-based pipeline using tasks and tools from the Common Astronomy Software Applications (CASA;McMullin et al. 2007;The CASA Team et al. 2022).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 PKS 1830 obtained from wideband radio continuum images in L Night1 and UHF band with robust=0 weighting are 11.245±0.001Jy at 1270 MHz and 11.40±0.01Jy at 832 MHz, respectively.The radio emission is unresolved in these images with a spatial resolution of 12. 9 × 8. 1 (position angle = −76. • 3) and 17. 4 × 13. 1 (position angle = +69.• 0), respectively (Gupta et al. 2021;Combes et al. 2021).The quoted uncertainty of the flux densities are based on a single 2D Gaussian fit to the continuum emission.The continuum flux density of PKS 1830 obtained using the Night2 dataset is S 1.27 GHz = 11.86 ± 0.02 Jy.This matches with the Night1 measurement within the flux density uncertainty (∼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 GHz ≈ 11.5 Jy.
The spectral line data products from ARTIP 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 ARTIP splits the L and UHF bands into 15 spectral windows (hereafter SPWs) (see Gupta et al. 2021;Combes et al. 2021, for details of SPW boundaries).The pixel sizes for L and UHF band image cubes are 2. 0 and 3. 0, respectively.XX and YY spectra were extracted in all SPWs from a single pixel at the location of PKS 1830 determined from the continuum images.The residual flux density in each continuum-subtracted spectrum is on the order of 0.5% and required additional spectralbaseline 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 PKS 1830) 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 ∼7 Jy at 1270 MHz, comparable to PKS 1830.We process this dataset following the procedures described above and use the resultant spectra to establish the genuineness of the results obtained for PKS 1830.

RRL Spectral Processing
Considering the rest frequencies of RRLs, 38 and 44 of the α (∆n = 1) recombination lines (per element) fall within the MALS L and UHF band coverage, respectively.2At z = 0.89, for example, the observations cover 31 α recombination lines in L band spanning principal quantum numbers n = 128 − 158 and 36 α recombination lines in UHF band spanning 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 sys ± 1500 km s −1 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 (≤ 5) order polynomial -with an order determined by minimizing the reduced χ 2 -could be fit over a well-behaved bandpass.RRLs that fell within ±1000 km s −1 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) (see Gupta et al. 2021;Combes et al. 2021), and we discarded the full spectral line chunk when flagged channels fell within |v| ∼ 500 km s −1 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 (≤ 5) polynomial to (line-free) channels in each spectral chunk and subtracted the fit.For the results presented in Sec. 3, 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 km s −1 to +55 km s −1 .
Fig. 2 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 1.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 km s −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 (e.g., Emig et al. 2020). 3We next smoothed the channel resolution to 8 km s −1 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 Hnα 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 PKS 1830 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 PKS 1830.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.
In Fig. 3, 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 spec-  Note-Uncertainties of the line properties are determined from the variance of each parameter as determined by the fit.N 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 km s −1 channel.vcenter is the central velocity of the best fit Gaussian and uncertainty.S peak is the peak amplitude of the best fit Gaussian fit.FWHM is the full-width half maximum of the best fit Gaussian.SRRL dv 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 km s −1 .τ dv is the velocity-integrated optical depth computed as − SRRL/Sc dv where the values used for the continuum flux density, Sc, are described in Sec. 3.
tra 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 S H144α dv = 337±16 mJy km s −1 with a maximum integrated signal-to-noise ratio (S/N) = (Σ N i=0 S i )/( √ N σ) of 21 is computed from the total-I spectrum by integrating over the N channels covering the velocity range −230 to 55 km s −1 .The effective transition of the total-I spectrum presented in the bottom panel of Fig. 3 is H144α, 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 H136α < 5 mJy upper limits to the H136α line at ν rest = 2585.7 MHz (ν obs = 1371.1 MHz) -a transition which is included in our stack -obtained by Mohan et al. (2002).
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 4 , (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. 4), 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. 2. We created an average RRL spectrum at z = 0.89 with an effective transition of H142α and show it in Fig. 4.This spectrum reached a noise, σ = 0.39 mJy, comparable to the PKS 1830 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 PKS 1830.
Two velocity components dominate the L band H144α emission from PKS 1830.In the bottom panel of Fig. 3, we show the best fit of two Gaussian profiles and the 4 We also refer the reader to Figure 2, where the noise, integrated signal, and integrated signal-to-noise properties also demonstrate no single or few outlying spectra dominate.spectrum that results when the two Gaussian profiles are subtracted.The properties of the Gaussian fits are listed in Table 1; the errors of each quantity are determined from the variance of each variable as determined by the fit.The H144α component centered on −117.4 ± 5.3 km s −1 arises from the NE sight-line (and thus we will refer to this velocity component as the NE component hereafter).The H144α component centered on 7.7 ± 4.2 km s −1 arises from the SW sightline (and thus we will refer to this velocity component as the SW component hereafter).In Table 1, we also include the integrated optical depth equal to τ dv ≈ − S RRL /S c dv computed by letting S c of each component equal half the total continuum flux density in the band, S c ≈ 5.75 Jy.We assume the continuum flux is equally split between the two core components following Koopmans & de Bruyn (2005)  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. 5 overlays the H144α spectrum with the H i 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 RRL line centroids of the NE and SW components agree within error with the OH 18 cm profiles, which have values of −110 ± 3 km s −1 and 6 ± 3 km s −1 , respectively.However, the RRL centroids are significantly offset from the dominant H i components at ∼ −150 km s −1 and 0 km s −1 , respectively, albeit the H i profile is more complex with perhaps five velocity components best-fitting the observed profile.The variation in central velocity with H i 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 ± 10 km s −1 and the OH FWHM of 94.2 ± 5 km s −1 .An estimate of the H i width is close to ∼ 100 km s −1 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α emission generally agree with the H i and OH profiles, more than two velocity components are discernible in the higher S/N spectra of H i and OH.For example, OH absorption has an additional component centered at −211 ± 3 km s −1 with a FWHM of 28 ± 9 km s −1 ; the RRL line profile fit to the NE component encompasses some emission in this velocity range.The H i absorption spectrum also shows two velocity features that are blueshifted with respect to the main ∼ −150 km s −1 peak.
In Fig. 6, 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 km s −1 results in S H163α dv = 309 ± 22 mJy km s −1 and a maximum integrated S/N of 14.The effective transition of the final spectrum is H163α, with a rest-frame frequency of 1504.6 MHz and observed at 797.85 MHz.
The H163α emission arises across a similar velocity range as the H144α emission, yet the peak intensities are slightly lower.The distinction of two components is less obvious in the H163α stack (UHF band), as compared with the H144α 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 cen = 7.7 km s −1 .The best fits are shown in Fig. 6 and the fit parameters are listed in Table 1.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 ≈ 5.7 Jy.
In Fig. 4, we show the L and UHF band stacked spectra of PKS 1830 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 1 and the integrated optical depth, τ dv ≈ − S RRL /S c dv, 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σ upper limit to the H158α line from the z = 0.19 absorber of S H158α < 0.5 mJy.In our L band spectrum of the z = 0.19 stack at an effective transition of H166α, we reach a spectral noise of 0.34 mJy, and our 3σ upper limit to the peak line emission of S H166α < 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 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σ upper limit to the carbon RRLs at z = 0.89 in L band is S C144α dv < 22.3 mJy km s −1 and in UHF band is S C163α dv < 37.5 mJy km s −1 , assuming a line width of 60 km s −1 .We searched for Hβ emission by stacking all available lines in both L and UHF band following the procedures described in Sec. 2. We report a 3σ upper limit to the Hβ emission of S H192β dv < 18.3 mJy km s −1 and an in-band ratio of peak H192β/H144α < 0.5, where typical in-band β/α ratios are ∼ 0.2 (e.g., Salas et al. 2017).

PHYSICAL CONDITIONS OF IONIZED GAS IN THE Z=0.89 ABSORBER
Because PKS 1830 has a relatively low Galactic latitude and is behind the Inner Galaxy, ( , b = 12.0 • , −5.7 • ), it has not yet been feasible to observe ionized gas tracers at UV through IR wavelengths from the z = 0.89 galaxy (e.g., Djorgovski et al. 1992;Courbin et al. 1998).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 ∼ 300 MHz in its spectral energy distribution (SED), which could be due to free-free absorption and would not arise from synchrotron self absorption (Pramesh Rao & Subrahmanyan 1988;Jauncey et al. 1991;Jones et al. 1996;Lovell et al. 1996, and for more details, see Sec 5.3).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 H144α = 131±14 km s −1 and FWHM H163α = 155±24 km s −1 are consistent within error.For the SW component, widths of FWHM H144α = 62.3±9.8 km s −1 and FWHM H163α = 80 ± 28 km s −1 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. 3 and 6).
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 km s −1 and for the SW component is 64 km s −1 .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. (2017b) provide the FWHM of a collisionally broadened profile, ∆ν col = n e n γ • 10 a /π where ∆ν col is in units of Hz, n e is in units of cm −3 , a = −7.386,γ = 4.439 and n is the principal quantum level.We place an upper limit of n e 7900 cm −3 for the NE component and n e 3700 cm −3 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, n, and frequency, ν, observed in front of a significantly-brighter background radio source of flux density S bkg,ν and which results from stimulated emission is given by where τ * n,ν is the LTE RRL optical depth, b n is the ratio of the number of hydrogen atoms in level n between the non-LTE and the LTE cases, and β characterizes the effect of stimulated transitions.b n and β, 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 (e.g., Shaver 1975;Hummer & Storey 1987;Salgado et al. 2017a;Prozesky & Smits 2018).Integrating over the line profile and expressing the LTE optical depth of the line, the integrated optical depth of an α transition (∆n = 1) 10 0 10 1 Rest Frequency (GHz) 10 2 dv (km s 1 ) NE SW Total 10 0 Observed Frequency (GHz) Figure 7.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.
at high principal quantum numbers takes the form where T e is the electron temperature and EM is the emission measure equal to EM = n e n ion d for electron density n e and ion density n ion integrated over the pathlength .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.
In Fig. 7, we plot the RRL SLED using the integrated optical depths corresponding to the NE and SW Gaussian components from from Table 1.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 km s −1 and setting the continuum flux density, S c , equal to the flux density in each respective band, i.e.
τ H144α = −0.029± 0.001 km s −1 and τ H163α = −0.027±0.002km s −1 .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. (2017a,b).When computing b n β 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 , 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.2), assuming a Gaussian distribution function with dispersion equal to the measurement uncertainty (see Table 1).The model is taken to be Eq. 2 in which the values of ν, n e , T e , and corresponding b n β 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] cm −3 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 cm −3 ; for pathlength, log( ) = [−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 ≈ 6000 − 10 000 K (e.g., Wenger et al. 2019;Tielens 2005).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 emcee (Foreman-Mackey et al. 2013).We used 20 walkers, 10 5 iterations, and verified convergence using an autocorrelation time analysis.
The 2D and 1D marginal posterior distributions of the parameters of the Total component are plotted in Fig. 8.They show that density is constrained by the models to within 0.6 dex with a maximum a posteriori value of n e = 500 cm −3 and 68.3% credible interval of 130 cm −3 ≤ n e ≤ 2000 cm −3 , and that electron temperatures are not well constrained.Typical temperatures of photoionized gas have values of T e ≈ 6000 − 10 000 K, with variations that are largely metallicity dependent (e.g., Shaver et al. 1983;Wenger et al. 2019).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 = 0.025 pc and the 68.3% credible interval is 7.9 × 10 −3 pc ≤ ≤ 0.13 pc.It is reasonable to infer that this ionized gas is not distributed in a very thin sheet of ∼ 0.025 pc, but instead has a volume Table 2. Constraints on physical conditions.Temperature is not well constrained within the allowed range log(Te) = [3, 4.3] K.
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, H ii regions (for further information see Sec. 5.1), etc.In Fig. 9, we show the marginal posterior distributions for all three components and each of the parameters, electron density (n e ), electron temperature (T e ), pathlength ( ), and emission measure (EM ).The parameter constraints are listed in Table 2, with the maximum a posteriori values and 68.3% credible intervals.We do not list the temperature in Table 2 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. 7, 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 2, 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 ∼ 99% of the emission at 1.4 GHz (Koopmans & de Bruyn 2005;Combes et al. 2021).To calculate the mass of ionized gas, we expect where m 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 ∝ ν −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 , we estimate a total mass for the NE component of M ion ≈ 10 6.0 +0.6 −0.7 M and for the SW component, M ion ≈ 10 5.5 +0.8 −0.4 M .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 ion ≈ 10 6.4 +0.7 −0.5 M .It is informative to also calculate the gas mass per unit area For the NE, SW, and Total components, we calculate Σ ion to be 10 0.3 +0.6 −0.7 M pc −2 , 10 −0.8 +0.8 −0.4 M pc −2 , and 10 0.0 +0.7 −0.5 M pc −2 , 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 g ∼ 5.3 kpc in units of M as log(M ion,g ) ≈ 8.2 +0.6 −0.7 , 7.1 +0.8 −0.4 , and 7.9 +0.7 −0.5 , 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 (WIM; see Tielens 2005).

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, where α B is the case B recombination coefficient.Using the posterior distributions of the emission measure, we calculate ionizing photon fluxes in units of photons s −1 pc −2 of log (Q o /area) = 47.0 +0.8 −1.2 , 46.8 +0.7 −1.0 , and 47.8 +0.4  −1.7 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 (Kado-Fong et al. 2020, and references therein).
While the gas mass estimates (Sec.4.3) 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.

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 ≈ 10 4.0±0.8cm −6 pc, within a typical error of 0.8 dex, and therefore we adopt an ionizing photon flux is Q o /area ≈ 10 46.9±0.8 photons s −1 pc −2 .We estimate the total ionizing photon rate for the galaxy within R g = 5.3 kpc as Q o = πR 2 g • 10 46.9±0.8 ≈ 10 54.8±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 ∼50 Myr) of SFR = 1 M yr −1 (Q o /1.4 × 10 53 photons s −1 ).Using this relation, an estimated SFR for the z = 0.89 galaxy is SFR ≈ 10 1.7±0.8M yr −1 and the SFR per unit area is Σ SFR ≈ 10 −0.2±0.8M yr −1 kpc −2 .A galaxy on the main sequence at z = 0.89 with SFR ∼50 M yr −1 has a typical stellar mass of ∼ 10 11 M (e.g., Schreiber et al. 2015).Current lensing models estimate the total mass within the Einstein ring as M E ≈ 4 × 10 11 M (Muller et al. 2020) and therefore a stellar mass of ∼ 8×10 10 M given a typical mass to stellar light ratio of ∼ 5 (Treu & Koopmans 2004).This is likely a main-sequence galaxy.
In Sec.4.3 we estimated the ionized gas mass per unit area of the NE component to be Σ ion ∼ 2.1 M pc −2 .For comparison, the H i column density is estimated to be N H I ≈ 5×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 ≈ 1.36m H , the atomic gas mass per unit area is Σ H I ≈ 50 M pc −2 .We use the OH column density of N OH ≈ 1.5 × 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 H2 ∼ 1.5 × 10 22 cm −2 .This H 2 column density is higher than the N H2 ∼ 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 ≈ 2m H , the molecular gas mass per unit area is Σ H2 ≈ 240 M pc −2 .
With these estimates, the neutral gas mass per unit area is Σ HI+H2 ∼ 290 M pc −2 .The ionized gas mass detected in RRLs of the NE component is a small fraction (∼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 Σ HI+H2 and the star formation rate Σ 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 Σ SFR ∼ 0.6 M yr −1 kpc −2 has a neutral gas mass per unit area of Σ HI+H2 ∼ 270 M pc −2 (e.g., Kennicutt & De Los Reyes 2021).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 ≈ 1600 cm −3 , is typical of young, compact H ii regions.If we assume the H ii regions are r H II ≈ 2 pc in size, then the covering fraction through this cross section of the galaxy is f c ≈ SW /(4/3 • r H II ) ≈ 7.5 × 10 −4 .The total number of H ii regions, N H II , is calculated from the surface density estimates, πr 2 C = N H II /f c • (π r 2 H II ), which computes to N H II ∼ 116, where we recall from Sec. 4.3 that r C = 786 pc for the SW component.Equivalently, assuming all H ii regions have n e ≈ 1600 cm −3 and r H II ≈ 2 pc, the ionized gas mass per H ii region is estimated to be M H II = 1860 M , and from the total ionized mass of the SW component, M SW = N H II M H II , also results in N H II ≈ (2.1 × 10 5 M )/(1860 M ) ∼ 113.However, this calculation warrants caution because the number of H ii regions changes dramatically depending on the assumed size.We set r H II ≈ 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 − Σ 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 H ii 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) indi-cates the thermal pressure of the medium (Jiang et al. 2019;Barnes et al. 2021).This results in a P − Σ SFR relation that serves as an important test-bed for pressureregulated, feedback-modulated star formation (e.g., Kim et al. 2013;Ostriker & Kim 2022).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 Σ SFR ∝ n 1.7±0.3 e in a sample of z ∼ 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 IR ≥ 10 11 L ) establish Σ SFR ∝ n 1.5 e (Herrera-Camus et al. 2016).
For Σ SFR ∼ 0.6 M yr −1 kpc −2 , the electron density predicted by the Shimakawa et al. (2015) relation is ∼ 110 cm −3 and by the Herrera-Camus et al. ( 2016) relation is ∼ 200 cm −3 .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 n e = n e o exp (−r/R g ), where R g is the scale length of the disk at which star-formation and density drops off and n e o is inner most density.With R g = 5.3 kpc and n e = 100 cm −3 of the NE component, at r = 2.4 kpc the expected density is n e ∼ 160 cm −3 .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 cm −3 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 PKS 1830 is complex.The lensing is achromatic and contains a smallscale 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. 1996Lovell 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).10 0 10 1 Observed Frequency (GHz) Figure 10.Compilation of continuum flux density measurements from the bright and highly variable (factors of ∼1.5 on weeks and years timescales) PKS 1830.The gray shaded region encompasses the 1σ 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.
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 τ ν = 6.67 × 10 −2 EM T −1.323 e (ν/GHz) −2.118 (e.g., Emig et al. 2022).We collected radio continuum measurements of PKS 1830 at ν 22 GHz (using the NASA/IPAC Extragalactic Database (NED) and Pramesh Rao & Subrahmanyan 1988;Henkel et al. 2008;Intema et al. 2017) and plot these in Fig. 10.We fit the SED of PKS 1830 as a power-law index with an external screen of free-free absorption, using the functional form .118 .Setting ν o = 40 GHz with respect to observed frequencies, the best fit parameters with 1σ uncertainties are S o = 4.8 ± 0.4 Jy, α = −0.24± 0.03 and τ o = (1.1 ± 0.2) × 10 −5 .
In Fig. 10, we also plot how ionized gas that has physical properties constrained by the RRL models -using the same selection of models presented in Fig. 7 would attenuate a power-law continuum SED (normalized by S o and α 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 PKS 1830.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 PKS 1830.The RRL emission is observed at z = 0.89 from a galaxy that lies along the line of sight and strongly lenses PKS 1830.This is the second detection of RRLs outside of the local universe (i.e., at z ≥ 0.076) and the first clearly associated with hydrogen (e.g., Emig et al. 2019).We detect H144α by stacking 17 RRLs covered by the L band (856-1712 MHz) with a S/N of 21 (see Fig. 3), and we detect H163α by stacking 27 lines in the UHF band (544-1088 MHz) with a S/N of 14 (see Fig. 6).Emission from the H144α 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 H i 21 cm and OH 18 cm absorption spectra (see Fig. 5), the H144α and H163α emission profiles span ∼250 km s −1 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 PKS 1830 (z = 2.5), from the z = 0.19 absorption system along this line of sight, or from the Milky Way (see Fig. 4).
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 2 and Fig. 9).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 ± 0.6 cm −3 and a volume-averaged pathlength of log( ) = −1.6 +0.7 −0.5 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 +1.9 −0.7 cm −3 and log( ) = −0.7 ± 1.1 pc, and the SW line-of-sight appears to intercept dense gas that is more typical of H ii regions with log(n e ) = 3.2 +0.4  −1.0 cm −3 and log( ) = −2.7 +1.8 −0.2 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 /area ≈ 10 46±0.8 photons s −1 pc −2 and star formation rate surface density of Σ SFR ∼ 10 −0.2±0.8M yr −1 kpc −2 .Taken over the z = 0.89 galaxy within R g ∼ 5.3 kpc, the ionizing photon rate of Q o ∼ 10 54.8 photons s −1 yields an average star formation rate of SFR ∼ 50 M yr −1 .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 PKS 1830 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 Σ ion ≈ 2.1 M pc −2 , in comparison with gas masses of Σ H I ≈ 50 M pc −2 and Σ H2 ≈ 240 M pc −2 (via OH) estimated from only the MALS observations.Given our estimated SFR, the H i+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 − Σ SFR relation determined from optical and FIR line ratios.
PKS 1830 is the first source investigated with MALS, and the detection of RRLs in the source is promising for the remaining ∼500 targets of the survey.With the first hydrogen RRL detection that breaks the redshiftbarrier, 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 −Σ SFR relation) and (b) the SED of AGN, thus eventually AGN evolution.We have also demonstrated the unique science that can be achieved through H i 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 unprece-dented capacities with future facilities such as the next generation Very Large Array (ngVLA; Murphy et al. 2018) and the SKA (Carilli 2015).

Figure 1 .
Figure1.The principal quantum numbers, n, of Hydrogen RRLs at which the stimulated-only emission of gas with a given density peaks.The thick black region marks the maximum n in the peak of emission, for ionized gas temperatures 5000 K ≤ Te ≤ 12 000 K, and the thin black region indicates the 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 ≤ z ≤ 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.

Figure 2 .
Figure2.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", σ, is the rms (outside of the line blank region) per 8 km s −1 channel of the XX (yellow circles) or YY (purple circles) spectrum."Integrated Signal" = ∆vΣ N i=0 Si, the sum of emission from the channels inside the line blank region."Integrated S/N" = (Σ N i=0 Si)/( √ N σ), 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.

Figure 3 .
Figure 3.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 km s −1 and 55 km s −1 within which we integrate the signal of the Total component.

Figure 4 .Figure 5 .
Figure 4. Non-detection RRL spectra in the spectrum of PKS 1830 (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.

Figure 6 .
Figure 6.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. 3 but for the UHF band spectrum.

Figure 8 .
Figure8.Corner plot showing the Total component constraints for the electron densities ne in units of cm −3 , electron temperatures Te in units of K, and pathlength 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.

Figure 9 .
Figure9.Posterior probability distributions for the parameters constrained through the RRL observations: electron density ne, electron temperature Te, pathlength , 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 ne < 7900 cm −3 (<3700 cm −3 ) constraint for the NE (SW) component obtained from the line width.

Table 1 .
Spectral and Line Properties