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

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 data set obtained on 2019 December 19 as "Night1." We also used UHF-band spectra from the data set obtained on 2020 July 13 and presented in Combes et al. (2021). In addition to these previously published data sets, we observed PKS 1830 a second time in the 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⁻¹ at the center of the band. For the UHF-band observations, 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⁻¹ 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. Because PKS 1830 is a bona fide VLA gain calibrator at this spatial resolution, a separate gain calibrator was not observed. The total on-source times on PKS 1830 are: 40 minutes (L-band Night1), 90 minutes (UHF-band), and 90 minutes (L-band Night2).

All MALS data sets 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 data sets 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 the L Night1 and UHF bands with robust = 0 weighting are 11.245 ± 0.001 Jy at 1270 MHz and 11.40 ± 0.01 Jy at 832 MHz, respectively. The radio emission is unresolved in these images, with spatial resolutions of 129 × 81 (position angle = −763) and 174 × 131 (position angle = +690), respectively (Gupta et al. 2021; Combes et al. 2021). The quoted uncertainties 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 data set 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 data sets 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 Combes et al. 2021; Gupta et al. 2021, for details of SPW boundaries). The pixel sizes for L- and UHF-band image cubes are 20 and 30, 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 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 data set 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 minutes. 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 data set following the procedures described above and use the resultant spectra to establish the genuineness of the results obtained for PKS 1830.

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. ¹⁶ At z = 0.89, for example, the observations cover 31α recombination lines in the L band, spanning principal quantum numbers n = 128–158, and 36 α recombination lines in the 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⁻¹ 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-order (≤5) polynomial—with an order determined by minimizing the reduced χ²—could be fit over a well-behaved bandpass. RRLs that fell within ±1000 km s⁻¹ 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 H i 21 cm and OH 18 cm line features (at relevant redshifts; see Combes et al. 2021; Gupta et al. 2021), and we discarded the full spectral line chunk when flagged channels fell within ∣v∣ ∼ 500 km s⁻¹ 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 broadband RFI and for which 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 Section 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⁻¹ to +55 km s⁻¹.

Figure 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.

Zoom In Zoom Out Reset image size

We next interpolated each spectral chunk to a common velocity grid with channel widths of 1 km s⁻¹, 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). ¹⁷ We next smoothed the channel resolution to 8 km s⁻¹ 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.

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 toward lower frequency, thereby providing information 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⁻¹ and FWHM_H163α = 155 ± 24 km s⁻¹ are consistent within errors. For the SW component, widths of FWHM_H144α = 62.3 ± 9.8 km s⁻¹ and FWHM_H163α = 80 ± 28 km s⁻¹ are also consistent within errors. This indicates that the lines are Doppler-broadened. Gaussian distributions do fit our observed line profiles reasonably well (see Figures 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⁻¹, and for the SW component it is 64 km s⁻¹. We then use the width at the lowest frequency to put a firm upper limit on the gas density, assuming pressure broadening dominates. It should be noted that 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⁻³, a = −7.386, γ = 4.439, and n is the principal quantum level. We place an upper limit of n_e ≲ 7900 cm⁻³ for the NE component and n_e ≲ 3700 cm⁻³ for the SW component.

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

$$\begin{eqnarray}&&{S}_{{\mathsf{n}},\nu }\approx -{\tau }_{{\mathsf{n}},\nu }^{* }\left({b}_{{\mathsf{n}}}\beta \right){S}_{\mathrm{bkg},\nu },\end{eqnarray} \tag{ 1 }$$

where ${\tau }_{{\mathsf{n}},\nu }^{* }$ 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. Here, 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 the path length. 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) at high principal quantum numbers takes the form

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \int {\tau }_{{\mathsf{n}}}\,{\mathsf{d}}v & \approx & 6.13\times {10}^{-4}\,\mathrm{km}\,\ {{\rm{s}}}^{-1}\left({b}_{{\mathsf{n}}}\beta \right)\\ & & \times \,\left(\displaystyle \frac{\mathrm{EM}}{{10}^{4}\,{\mathrm{cm}}^{-}6\,\mathrm{pc}}\right){\left(\displaystyle \frac{{T}_{e}}{{10}^{4}\,{\rm{K}}}\right)}^{-2.5}{\left(\displaystyle \frac{\nu }{\mathrm{GHz}}\right)}^{-1},\end{array}\end{eqnarray} \tag{ 2 }$$

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 path length ℓ. Because the departure coefficients at each principal quantum number are strongly dependent upon density, they are key to breaking the degeneracy between density and path length in the emission measure.

In Figure 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⁻¹ 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⁻¹ and ∫τ_H163α = −0.027 ± 0.002 km s⁻¹. 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.

Zoom In Zoom Out Reset image size

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, 2017b). 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 (Equation (2)), assuming a Gaussian distribution function with dispersion equal to the measurement uncertainty (see Table 1). The model is taken to be Equation (2), in which the values of ν, n_e , T_e , and corresponding b_n β are input, and the path length 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, $\mathrm{log}({n}_{e})\,=\,$[−1, 6] cm⁻³ for the Total component, and we used the upper limits derived from the line width for the NE and SW components, $\mathrm{log}({n}_{e})\,=\,$[−1, 3.4] and $\mathrm{log}({n}_{e})\,=\,$[−1, 3.6], respectively, in units of cm⁻³; for path length, $\mathrm{log}({\ell })=$[−9, 6] pc; and for temperature, we select a strict range of $\mathrm{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⁵ 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 Figure 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⁻³ and 68.3% credible interval of 130 cm⁻³ ≤ n_e ≤ 2000 cm⁻³, 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 path length for the Total component is ℓ = 0.025 pc and the 68.3% credible interval is 7.9 × 10⁻³ 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 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 scii regions (for further information, see Section 5.1), etc.

Zoom In Zoom Out Reset image size

In Figure 9, we show the marginal posterior distributions for all three components and each of the parameters, electron density (n_e ), electron temperature (T_e ), path length (ℓ), 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 Figure 7, in order to demonstrate the variation in the predicted integrated optical depth for different models.

Zoom In Zoom Out Reset image size

Table 2. Constraints on Physical Conditions. Temperature Is Not Well-constrained within the Allowed Range $\mathrm{log}({T}_{e})=[3,4.3]$ K

	vcen	$\mathrm{log}({n}_{e})$	$\mathrm{log}({\ell })$	$\mathrm{log}({EM})$
	(km s−1)	(cm−3)	(pc)	(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}$

Download table as:  ASCIITypeset image

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 they 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

$$\begin{eqnarray}\begin{array}{rcl}{M}_{\mathrm{ion}} & \approx & 1.36{m}_{{\rm{H}}}\cdot {n}_{e}\cdot (\pi {r}_{c}^{2}){\ell }\\ {M}_{\mathrm{ion}} & \approx & 0.11\,{M}_{\odot }\left(\displaystyle \frac{{n}_{e}}{{\mathrm{cm}}^{-3}}\right){\left(\displaystyle \frac{{r}_{c}}{\mathrm{pc}}\right)}^{2}\left(\displaystyle \frac{{\ell }}{\mathrm{pc}}\right),\end{array}\end{eqnarray} \tag{ 3 }$$

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 01 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 005 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}_{\mathrm{ion}}\approx {10}^{{6.0}_{-0.7}^{+0.6}}$ M_⊙ and for the SW component, ${M}_{\mathrm{ion}}\approx {10}^{{5.5}_{-0.4}^{+0.8}}$ 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}_{\mathrm{ion}}\approx {10}^{{6.4}_{-0.5}^{+0.7}}$ M_⊙.

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

$$\begin{eqnarray}&&{{\rm{\Sigma }}}_{\mathrm{ion}}\approx 0.033\,{M}_{\odot }\,{\mathrm{pc}}^{-2}\left(\displaystyle \frac{{n}_{e}}{{\mathrm{cm}}^{-3}}\right)\left(\displaystyle \frac{{\ell }}{\mathrm{pc}}\right).\end{eqnarray} \tag{ 4 }$$

For the NE, SW, and Total components, we calculate Σ_ion to be ${10}^{{0.3}_{-0.7}^{+0.6}}$ M_⊙ pc⁻², ${10}^{-{0.8}_{-0.4}^{+0.8}}$ M_⊙ pc⁻², and ${10}^{{0.0}_{-0.5}^{+0.7}}$ M_⊙ pc⁻², 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 $\mathrm{log}({M}_{\mathrm{ion},{\rm{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, given that we do not trace the bulk of the ionized gas mass of the galaxy contained in the Warm Ionized Medium (WIM; see Tielens 2005).

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}\begin{array}{rcl}\displaystyle \frac{{Q}_{o}}{\mathrm{area}} & \approx & \mathrm{EM}\cdot {\alpha }_{B}\\ \displaystyle \frac{{Q}_{o}}{\mathrm{area}} & \approx & 7.6\times {10}^{45}\,\mathrm{photons}\ \,\ {{\rm{s}}}^{-1}\,{\mathrm{pc}}^{-2}\displaystyle \frac{{EM}}{{10}^{3}\,{\mathrm{cm}}^{-6}\,\mathrm{pc}},\end{array}\end{eqnarray} \tag{ 5 }$$

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⁻¹ pc⁻² of $\mathrm{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 (Kado-Fong et al. 2020, and references therein).

While the gas mass estimates (Section 4.3) are likely lower limits, the ionized photon fluxes are closer to realistic values, because they are dominated by the large emission measures that we are sensitive to.

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.8 cm⁻⁶ pc, within a typical error of 0.8 dex, and therefore we adopt an ionizing photon flux of Q_o /area ≈ 10^46.9±0.8 photons s⁻¹ pc⁻². 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⁻¹. 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⁻¹ (Q_o /1.4 × 10⁵³ photons s⁻¹). Using this relation, an estimated SFR for the z = 0.89 galaxy is SFR ≈ 10^1.7±0.8 M_⊙ yr⁻¹ and the SFR per unit area is Σ_SFR ≈ 10^−0.2±0.8 M_⊙ yr⁻¹ kpc⁻². A galaxy on the main sequence at z = 0.89 with SFR ∼ 50 M_⊙ yr⁻¹ has a typical stellar mass of ∼10¹¹ M_⊙ (e.g., Schreiber et al. 2015). Current lensing models estimate the total mass within the Einstein ring as M_E ≈ 4 × 10¹¹ M_⊙ (Muller et al. 2020), and therefore a stellar mass of ∼8 × 10¹⁰ M_⊙ given a typical mass-to-stellar-light ratio of ∼5 (Treu & Koopmans 2004). This is likely a main-sequence galaxy.

In Section 4.3, we estimated the ionized gas mass per unit area of the NE component to be Σ_ion ∼ 2.1 M_⊙ pc⁻². For comparison, the H i column density is estimated to be N_{H I} ≈ 5 × 10²¹ cm⁻², 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⁻². We use the OH column density of N_OH ≈ 1.5 × 10¹⁵ cm⁻² (Gupta et al. 2021) and assume an abundance of 10⁻⁷ (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⁻². This H₂ column density is higher than the ${N}_{{{\rm{H}}}_{2}}\sim {10}^{21}$ cm⁻² derived from H₂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 ${{\rm{\Sigma }}}_{{{\rm{H}}}_{2}}\approx 240$ M_⊙ pc⁻².

With these estimates, the neutral gas mass per unit area is ${{\rm{\Sigma }}}_{{\rm{H}}{\rm\small{I}}+{{\rm{H}}}_{2}}\sim 290$ M_⊙ pc⁻². 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) has long established a direct correlation between the surface densities of the neutral gas mass ${{\rm{\Sigma }}}_{{\rm{H}}{\rm\small{I}}+{{\rm{H}}}_{2}}$ 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⁻¹ kpc⁻² has a neutral gas mass per unit area of ${{\rm{\Sigma }}}_{{\rm{H}}{\rm\small{I}}+{{\rm{H}}}_{2}}\sim 270$ M_⊙ pc⁻² (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⁻³, is typical of young, compact H scii regions. If we assume the H scii 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⁻⁴. The total number of H scii regions, N_{H II} , is calculated from the surface density estimates, $\pi {r}_{C}^{2}={N}_{{\rm{H}}\,{\rm\small{II}}}/{f}_{c}\cdot (\pi \,{r}_{{\rm{H}}\,{\rm\small{II}}}^{2})$, which computes to N_{H II} ∼ 116, where we recall from Section 4.3 that r_C = 786 pc for the SW component. Equivalently, assuming all H scii regions have n_e ≈ 1600 cm⁻³ and r_{H II} ≈ 2 pc, the ionized gas mass per H scii region is estimated to be M_{H II} = 1860 M_⊙, and that 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⁵ M_⊙)/(1860 M_⊙) ∼ 113. However, this calculation warrants caution because the number of H scii 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 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 scii regions thermalize with the ISM (Gutiérrez & Beckman 2010), i.e., take on a pressure balance with others phases, then the volume-averaged 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 − Σ_SFR relation that serves as an important test bed for pressure-regulated, feedback-modulated star formation (e.g., Kim et al. 2013; Ostriker & Kim 2022). Using doublet ratios of [S ii] and [O ii] in the optical (Kewley et al. 2019), Shimakawa et al. (2015) and Jiang et al. (2019) find a relation of ${{\rm{\Sigma }}}_{\mathrm{SFR}}\propto {n}_{e}^{1.7\pm 0.3}$ in a sample of z ∼ 1–3 starburst galaxies. In the far-IR, the ratio of [N ii] fine structure lines from a sample of nearby normal galaxies and (ultra)luminous infrared galaxies ((U)LIRGs; L_IR ≥ 10¹¹ L_⊙) establish ${{\rm{\Sigma }}}_{\mathrm{SFR}}\propto {n}_{e}^{1.5}$ (Herrera-Camus et al. 2016).

For Σ_SFR ∼ 0.6 M_⊙ yr⁻¹ kpc⁻², the electron density predicted by the Shimakawa et al. (2015) relation is ∼110 cm⁻³ and that by the Herrera-Camus et al. (2016) relation is ∼200 cm⁻³. 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 toward smaller galactic radii, r, as $\unicode{x03008}{n}_{e}\unicode{x03009}={\unicode{x03008}{n}_{e}\unicode{x03009}}_{o}\exp (-r/{R}_{g})$, where R_g is the scale length of the disk at which star formation and density drops off and ${\unicode{x03008}{n}_{e}\unicode{x03009}}_{o}$ is the innermost density. With R_g = 5.3 kpc and $\unicode{x03008}{n}_{e}\unicode{x03009}=100$ cm⁻³ of the NE component, at r = 2.4 kpc the expected density is n_e ∼ 160 cm⁻³. 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⁻³ does fall within the 68% credible interval of the SW component. The smaller path length 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.

The radio continuum emission from PKS 1830 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, 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 became optically thick, it 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}\,\mathrm{EM}\,{T}_{e}^{-1.323}{(\nu /\mathrm{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 Figure 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 ${S}_{\nu }\,={S}_{o}{\left(\tfrac{\nu }{{\nu }_{o}}\right)}^{\alpha }\exp (-{\tau }_{\nu })$ with ${\tau }_{\nu }={\tau }_{o}{\left(\tfrac{\nu }{{\nu }_{o}}\right)}^{-2.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⁻⁵.

Zoom In Zoom Out Reset image size

In Figure 10, we also plot how ionized gas that has physical properties constrained by the RRL models—using the same selection of models presented in Figure 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 those observed in PKS 1830. For the z = 0.89 galaxy to cause the free–free absorption, volume-averaged path lengths a factor of 5 larger are needed. A smaller filling factor and larger path length intercepting the radio emission would not be 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 EMs 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.