© The Authors 2023

1. Introduction

Astronomical objects often appear to be of small angular size. Therefore, the study and characterization of these objects and of the inherent physical mechanisms (i.e., the whole subject of astrophysics) rely on our ability to access information corresponding to these small scales. This can be done in two ways: either with direct high-angular-resolution observations, for example by means of very long-baseline interferometry (VLBI), or, for varying sources, by studying their time variations and their timescales. Using the causality argument, the typical timescale of signal variations should indeed reflect the size of the corresponding regions.

The time variability technique has been applied to various sources and contexts, including compact transient sources like transiting exoplanets, supernovae, pulsars, and black holes, to name but a few. For example, the reverberation mapping technique within the broad-line region in active galactic nuclei (AGN) provides size estimates from the line variability, which, combined with spectroscopy measurements of the kinematics, yields a mass estimate for the central black hole (e.g., Bahcall et al. 1972; Blandford & McKee 1982). Another example is provided by the remarkable periodicity in the optical light curve of the gravitationally lensed quasar Q J 0158−4325, which with the help of microlensing led (Millon et al. 2022) to propose a model based on a pair of supermassive black holes in the coalescence stage. Yet another application is the use of temporal variability of H I, optical, or ultraviolet absorption lines to study the small-scale structures down to a few AU in the interstellar medium (ISM; see e.g., a review by Stanimirović & Zweibel 2018).

The direction toward the lensed quasar PKS 1830−211, the subject of this work, is particularly favorable for studying time variability. PKS 1830−211 is a blazar at z = 2.5 (Lidman et al. 1999) with strong gravitational lensing due to an intervening nearly face-on spiral galaxy at z = 0.89 (Wiklind & Combes 1996; Winn et al. 2002). The blazar shows significant activity on short timescales (down to hours and days) when observed from radio waves to X- and γ-rays (e.g., Abdo et al. 2015; Martí-Vidal & Muller 2019; Martí-Vidal et al. 2020). Two factors boost the apparent variability of PKS 1830−211 from our viewpoint on Earth: (i) the orientation of the blazar’s relativistic jet close to the line of sight, which makes it prone to apparent superluminal motions, and (ii) the lensing and its associated source magnification. As a result of the lensing geometry, the quasar appears as two main bright and compact images (to the northeast and the southwest, hereafter the NE and SW images), separated by ∼1″ and embedded in an Einstein ring with a steep spectral index (mostly seen at centimetre (cm) wavelengths; e.g., Jauncey et al. 1991), and a weak third lensed image at approximately the position of the lens galaxy (Muller et al. 2020a). There is a geometrical time delay of ∼27 days between the NE (leading) and SW images (see Sect. 3.3).

In addition to the peculiarities of PKS 1830−211 mentioned above, the sightlines to the two main images intercept a molecular component in the disk of the z = 0.89 lens galaxy, causing a remarkable molecular absorption-line system (Wiklind & Combes 1998; Muller et al. 2014). A large number (> 60) of molecules have been detected along the SW line of sight (e.g., Muller et al. 2011, 2014; Tercero et al. 2020), which have been used not only to investigate the chemical and isotopic composition of the gas in the absorber but also as sensitive cosmological probes; for example, of the cosmic microwave background temperature (Muller et al. 2013) and invariance of fundamental constants (e.g., Bagdonaite et al. 2013; Kanekar et al. 2015; Muller et al. 2021). However, these studies can be affected by systematic uncertainties due to the peculiarity of PKS 1830−211, with time variations being one of the major effects (e.g., Muller et al. 2021).

VLBI observations have revealed micro-changes in the morphology of the quasar on the timescale of months (Garrett et al. 1997). Furthermore, the apparent separation between the two lensed images has been measured to vary by up to 200 μas within eight months (Jin et al. 2003), which was interpreted by Nair et al. (2005) as evidence for recurrent plasmon ejections along a helical jet, with a jet precession period of one year. In addition, there is evidence for micro-lensing events in the X-ray and γ-ray light curves (e.g., Oshima et al. 2001; Neronov et al. 2015), implying that the hard emission is produced in the vicinity of the supermassive black hole.

These changes in the morphology of the quasar lead to variation in the illumination of the foreground material in the lens–absorber, and significant variation in the absorption profiles of molecular lines has been reported between spectra taken months to years apart at millimeter (mm) wavelengths (Muller & Guélin 2008; Muller et al. 2014, 2021; Schulz et al. 2015). In contrast, there seems to be no line variability at cm wavelengths: Allison et al. (2017) and Combes et al. (2021) reported no significant variation in the redshifted H I 21-cm and OH lines between observations made at up to 20 year time intervals. This could be explained if the H I gas is distributed on larger scales than that probed by the changes in the source structure at these wavelengths. Regarding variations of absorption lines, the situation is indeed more favorable for molecular lines at high frequencies (i.e., in the mm/submm part of the electromagnetic spectrum): the background continuum emission is smaller (e.g., Guirado et al. 1999) and originates closer to the base of the jet; in addition, the molecular screen is expected to be more clumpy and to have a smaller scale height than the atomic screen.

Overall, PKS 1830−211, with its sustained blazar activity, its low-inclination jet with superluminal motions, its lens magnification, potential micro-lensing events, and its intervening lens–absorber with absorption profile variability, is thus a particularly interesting source for time monitoring. In this work, we present the results of a monitoring of the continuum activity of PKS 1830−211 and of H₂O and CH submm absorption lines with the Atacama Large Millimeter Array (ALMA) over a time span of six months in 2016. The monitoring is complemented by ALMA archival data in the 2012−2022 period. Observations are presented in Sect. 2. Section 3 covers the quasar continuum variations. Results from spectral line absorptions are given in Sect. 4. A discussion in Sect. 5 brings together the outcomes of continuum and spectral-line variations. A summary is given in Sect. 6.

2. Observations

2.1. 2016 monitoring campaign

Observations were made during ALMA Cycle 3, between March and September 2016. The proposed monitoring strategy comprised one visit every 10 days, over a time span of six months. To accommodate dynamic scheduling constraints, a time window of ±7 days was set for each visit. A journal of the observations is given in Table 1. The proposed cadence was achieved – despite the scheduling challenge –, with a total of 17 successful visits. Among them, two visits were carried out on the same day (on 26 March 2016), providing a high signal-to-noise-ratio reference spectrum when combined, and the last visit (on September 8 2016) was carried out after a gap of 42 days with respect to the penultimate visit. Otherwise, the time intervals between visits range between 6 and 21 days.

During the monitoring, the array configuration changed, but the longest baseline was always longer than 450 m, such that the synthesized beam was always smaller than 0.5″ and the two main lensed images of the quasar are therefore always well resolved. The on-source integration time for each visit was slightly less than 10 min, set by our original goal, which is to be able to track changes in the spectral line profiles down to a few percent during each individual visit.

Observations were made in dual polarization mode, with two orthogonal polarizers (XX and YY) in the antenna alt-azimuthal frame. The correlator setup was chosen to simultaneously cover the ground-state transitions of H₂O (557 GHz, rest frequency) and CH (two Λ-doublets at ∼533 and 537 GHz, rest frequency) redshifted into Band 7. To cover these lines, three spectral windows were centered at ∼282.4, 284.6, and 295.4 GHz sky frequency, with a bandwidth of 938 MHz and a spectral resolution of 0.56 MHz, resulting in a velocity resolution of 0.6 km s⁻¹ after Hanning smoothing. A fourth spectral window, not including any detectable spectral line, was centered at 296.8 GHz with a coarser spectral resolution of 31.25 MHz in order to provide extra continuum data.

The calibration was done with the CASA (version 4.7.2) package (CASA Team 2022). The bandpass calibration was determined from the bright quasar J 1924−2914. The flux calibration was carried out using Titan (Butler-JPL-Horizons 2012 model), Pallas, J 1924−2914, or J 1733−130, as indicated in Table 1. We expect an absolute flux accuracy of the order of 10%−20% according to the ALMA Technical Handbook (Cortes et al. 2022). The gain calibrator was J 1832−2039 in all visits. After the standard calibration, a step of phase-only self-calibration on PKS 1830−211 was done with a short solution interval of 10 s.

The calibrated visibilities were finally fitted using the Python task UVMultiFit (Martí-Vidal et al. 2014) with a source model consisting of two point-like components, where the flux density of the NE image and the flux-density ratio, ℜ = f_NE/f_SW, were set as free parameters. Different fit runs were carried out using Stokes I data and data from the XX and YY polarizers separately in order to extract differential polarimetric information (see Sect. 3.4). The additional weak features of PKS 1830−211 reported from ALMA observations by Muller et al. (2020a), namely the third image (about 150 times weaker than the NE and SW images) and the small extensions from the NE and SW images (corresponding to the start of the Einstein ring and with a steep spectral index), have a negligible impact on the fitting process and were therefore not considered.

2.2. Complementary data

In addition to the 2016 monitoring data, we also used observations of PKS 1830−211 available in the ALMA archive to extract continuum information (flux density and flux density ratio) between 2012 and 2022, as well as previous ALMA observations of the same fundamental transitions of CH and H₂O as mentioned above, as well as of the rare isotopologs HO and HO (see Table 2). The data reduction was done with the same method as described above. All the spectral line data have a velocity resolution of ∼1.1 km s⁻¹ or better.

3. Results from continuum data

Table 1 lists the values of the flux densities, flux density ratios, polarization ratios, and spectral index differences of the two NE and SW lensed images during the 2016 monitoring and we present their time variations in Fig. 1. We used four different flux-density calibrators (Titan, Pallas, J 1924−292, and J 1733−130). The data points follow a relatively smooth light curve, which is free of significant outliers, giving confidence that potential systematic effects in the flux calibration are limited. In addition, while the absolute flux accuracy is of ∼10%−20%, we emphasize that the accuracy of relative quantities, such as the flux-density ratio of both lensed images and the polarization ratio, is much higher (e.g., at a per mil level for the flux density ratio), and can therefore offer subtle diagnostics of the quasar activity, as discussed in the following subsections.

Fig. 1. Time evolution of the continuum properties of PKS 1830−211 during our monitoring in 2016 (see Table 1): (a) Flux density of the NE and SW images (error bars of 20% are indicated for the NE image data points only; the color code indicates the source used for flux density calibration). (b) NE/SW flux density ratio. (c) Difference in flux density between the XX and YY polarizers normalized by the Stokes I flux density. (d) Polarization ratio ℛpol (values are encoded as a function of the PA of the observations). (e) Difference in spectral index. The expected time delay of ∼27 days between the two lensed images of the quasar is indicated in panel a.

3.1. Submillimeter light curve

Over our six-month monitoring of PKS 1830−211, we observe a remarkable, nearly monotonic, and relatively featureless decrease in the flux density of the NE image, from ∼1.2 Jy to ∼0.2 Jy, which is a drop of more than 80% of its flux density at the start of the monitoring. The evolution of the NE and SW images is highly correlated and on a longer timescale than the time delay between the two images, indicating that this spectacular drop is intrinsic to the quasar and not due to a differential lensing effect. On average, the flux density decreased by ∼5 mJy day⁻¹ over the six-month monitoring. This rate is small compared to previously reported flux variations of up to ∼10% on an hourly timescale (e.g., Martí-Vidal et al. 2016; Martí-Vidal & Muller 2019), but it is remarkable that it was sustained over such a relatively long period.

3.2. Flux-density ratio

In contrast to the large decrease in flux density during the 2016 monitoring, the instantaneous flux-density ratios ℜ = f_NE/f_SW show only mild (∼10%) variations around an average value of ∼1.04, with extreme values of 0.96 and 1.15 (Fig. 1b). The flux-density ratio measurements, being relative by nature, are free from instrumental effects. Given the high signal-to-noise ratio of the data, their accuracy is of the order of one per mil. As discussed by Martí-Vidal et al. (2013), the variations of the flux-density ratio over timescales of shorter than a few days to a few weeks are likely a direct signature of the intrinsic activity of the quasar modulated by the time delay between the two images, because we would expect timescales of longer than a month for milli-lensing.

The flux-density ratio is best measured at mm/submm wavelengths, where the contribution from the Einstein ring is lower due to its steep spectral index (Jauncey et al. 1991, see also Muller et al. 2020a). However, prior to ALMA it was difficult to obtain mm observations with sufficient angular resolution to separate the two lensed images (partly due to the southern declination of PKS 1830−211 and to the small number of antennas of the previous generation of mm interferometers). Nevertheless, the flux-density ratio could be estimated from the saturation level of saturated lines, such as HCO⁺ or HCN J = 2 − 1, which block the SW image almost completely (e.g., Wiklind & Combes 2001; Muller & Guélin 2008). Between 1995 and 2007, the flux-density ratio was found to vary between 1 and 3, but with only rare excursions to the extreme values. The average value of these old measurements was 1.7 ± 0.3 (Muller & Guélin 2008).

The flux-density ratios during the 2016 monitoring are among the lowest of all ALMA measurements since 2012, as shown in Fig. 2a. In 2012, the first ALMA Cycle 0 observations were taken serendipitously at the time of a γ-ray flare. The submm flux-density ratios showed remarkable temporal and chromatic variations between 1.2 and 1.5 within two months, which were interpreted by Martí-Vidal et al. (2013) with a simple model of plasmon injection and opacity effects in the jet of the quasar. Subsequent ALMA observations in 2014−2015 show slightly higher flux-density ratios between 1.4 and 2.2, contrasting with the measurements in 2016−2020, which are closer to ∼1, in particular those taken close in time to the record-breaking γ- and radio flare in 2019. The most recent measurements at hand (2022) show values in the range 1.3−2.0.

Fig. 2. History of ALMA measurements toward PKS 1830−211. (a) Flux-density ratios of PKS 1830−211 NE and SW images (with the color code corresponding to ALMA bands, as indicated at the top of the figure). (b) Flux densities of the NE image of PKS 1830−211, scaled to an equivalent frequency of 300 GHz such as S​300 GHz = S​ν × (ν/300 GHz)α with a spectral index α = −0.7 (error bars indicate a nominal flux uncertainty of 20%) overlaid on top of the Fermi-LAT γ-ray daily-monitoring light curve (in green).

With the current collection of flux-density ratio measurements, it is difficult to obtain a clear and definitive determination of the true differential magnification, ℛ_q, between the two lensed images. Taken at face value, the average value of all the ALMA flux-density ratios measured between 2012 and 2022 is 1.38 ± 0.35 (with 1σ dispersion) and the median is 1.29. However, the evolution of the flux-density ratios in Fig. 2a shows a modulation on timescales of the order of several years, which could be the signature of “milli-lensing” (e.g., as discussed by Martí-Vidal et al. 2013). On shorter timescales of a few days to weeks, micro-lensing could also affect the picture (e.g., Wambsganss 1992). This could be tested with a long and regular time monitoring, including multifrequency observations in order to probe the chromaticity of the flux-density ratio evolution.

3.3. The time delay and the differential magnification factor

The determination of the geometrical time delay between lensed images is interesting for measuring the Hubble constant, H₀, via the time-delay cosmography method (Refsdal 1964). Now that the third lensed image of the quasar has been unambiguously detected by Muller et al. (2020a) – thereby pinpointing the location of the lensing galaxy –, the major remaining sources of uncertainty for the lens modeling of the PKS 1830−211 system that can be constrained observationally are the time delay and the relative magnification factor between the NE and SW images.

Both quantities have been investigated with different methods and at different wavelengths (see Table 3). At low radio frequencies, the required angular resolution of < 1″ and the contamination from the Einstein ring make it difficult to disentangle the individual light curves from both images (van Ommen et al. 1995; Lovell et al. 1998). As mentioned above, it is possible to measure the flux-density ratio at mm wavelengths using the foreground molecular absorption and the saturation level of high-opacity lines, where absorption almost completely masks the SW image and therefore makes it possible to retrieve the flux densities of both NE and SW images even with low-angular-resolution data (Wiklind & Combes 2001; Muller et al. 2006). Time-series analyses of the γ-ray light curve benefit from long time ranges and daily monitoring¹. However, results with this method are still controversial (e.g., Barnacka et al. 2011; Abdo et al. 2015). Their interpretation might even be more complicated if individual γ-flares arise at different locations (Barnacka et al. 2015), implying slightly different time delays and magnification ratios. Overall, the results seem to converge toward a time delay in the range of 20 − 30 days and a differential magnification of ∼1.0 − 1.5, which seems consistent with a reasonable lens model (Muller et al. 2020a).

The relatively featureless light curve during our ALMA 2016 monitoring and potential effects from milli- or micro-lensing prevent us from setting new strong constraints on the time delay and differential magnification of the PKS 1830−211 system. Nonetheless, we applied a measure-of-randomness analysis (see, e.g., Chelouche et al. 2017) using the Von Neumann estimator (ℰ_VN), which is defined as:

(1)

where the datapoints ℱ(t_i) correspond to the flux-density measurements from the NE image augmented with the points corresponding to the flux-density measurements of the SW image shifted in time by the time delay, Δt, and scaled by the differential magnification, ℛ_q, that is, . By varying the two quantities Δt and ℛ_q, we can calculate the corresponding values of ℰ_VN. The assumption of minimum randomness (as in the χ² minimization) points to the most likely values of the two parameters ℛ_q and Δt. We show the results of this analysis in Fig. 3. Given the quasi-monotonic and relatively featureless light curve of PKS 1830−211 during our monitoring, there is a strong degeneracy between the two parameters. Without significant a priori information on the differential magnification, we cannot better constrain the time delay with only our 2016 flux-density measurements. In Sect. 3.6, we further present a simple parametric model of the continuum evolution during the monitoring, from which we extract best-fit values of the time delay, namely 25 ± 3 days, and of the differential magnification, 1.25 ± 0.04, which fall in the valley of lowest Von Neumann randomness in Fig. 3.

Fig. 3. Measure of randomness using the Von Neumann estimator ℰVN (defined in Eq. (1)) as a function of the time delay and differential magnification parameters. The best-fit value obtained from a parametric model of the continuum evolution (Sect. 3.6) is shown by the white mark.

3.4. Differential polarimetry

To continue our analysis of the continuum emission properties of PKS 1830−211 during our monitoring, we can also investigate the evolution of its polarization state. Although the ALMA data were not taken in full polarization mode, we can use the data obtained separately from the two orthogonal polarizers, XX and YY, in the alt-azimuth frame of the antennas to retrieve some source polarization properties using a differential polarimetry analysis (see, e.g., Martí-Vidal et al. 2015, 2016; Martí-Vidal & Muller 2019). First, we can form the difference ratio Δ_XY = (S​_XX − S​_YY)/S​_I, where S​_i are the signal retrieved from the XX and YY polarizers, and I is the Stokes I total flux. We already see in Fig. 1c that a clear polarization signal is present in the data, as Δ_XY varies within ±8% during the monitoring.

Furthermore, we can determine the flux-density ratios ℛ_XX and ℛ_YY from the XX and YY polarizers separately, and form the double polarization ratio, ℛ_pol:

(2)

as defined by Martí-Vidal et al. (2015). If the quasar is not polarized or if both lensed images have strictly the same polarization properties (polarization intensity and electric vector position angle, hereafter EVPA), then ℛ_pol = 0. However, if the two images have different polarization properties, which could be due to a faster polarization variability than the time delay, then ℛ_pol would be a nonzero quantity that would depend on the polarization intensity and EVPA angle differences between the two images as well as on the parallactic angle of the observations. With observations at different frequencies, it is furthermore possible to derive the Faraday rotation from ℛ_pol (see Martí-Vidal et al. 2015, 2016).

The variations of ℛ_pol encode information on the polarization variability of the quasar. In Fig. 1d, we see that |ℛ_pol| varies significantly between 0% and ∼5% during the first half of the monitoring, but relaxes to ∼0 after MJD = 7580. This suggests that the EVPA stabilized for a longer time interval than the time delay, after which the second image also stabilized. We note that the quasar fractional polarization eventually reaches a level of ∼5%, as shown by the nearly constant Δ_XY after MJD = 7580. This value is comparable to other polarization measurements in 2015 (Band 9, Martí-Vidal & Muller 2019) and 2019 (Band 6, Martí-Vidal et al. 2020), although the polarization variability can be significant over a shorter timescale than the time delay (Martí-Vidal et al. 2020).

3.5. Spectral index

We finally show in Fig. 1e the evolution of the spectral indices² of the two lensed images and their difference. The spectral indices were measured directly from the spectra corresponding to each lensed image after flagging channels associated with absorption features. We emphasize that any instrumental or calibration effects are removed when taking the spectral index difference Δα = α_NE − α_SW between the two lensed images.

Throughout the monitoring, the difference Δα was almost always positive, which reflects the rapid intrinsic variations of the quasar at a rate faster than the time delay, with the NE image having a steeper spectral index. Until the end of April 2016 (MJD ∼ 7500), the spectral index is ∼ − 0.7. It then gets a little steeper, except for a sudden rise of the quasar spectral index for two measurements around MJD = 7550 (15 and 22 June 2016), which corresponds to a small flux-density bump in the light curve.

It is interesting to note that the spread of Δα (i.e., between −0.05 and 0.15) in the different visits of our six-month monitoring is comparable to the change of Δα measured by Martí-Vidal et al. (2020) within a two-hour observation coincident in time with the onset of the 2019 record-breaking flare of PKS 1830−211. This suggests a considerable change of activity level between these two periods.

3.6. Simple scenario of the continuum evolution

The smooth and relatively slow flux-density evolution of the quasar suggests that the continuum spectra of the two lensed images should be similar, that is, that there are no clear features in the light curve on timescales short enough to introduce large differences between the NE and SW brightness or spectral index. Also, we do not consider the effects of micro- and milli-lensing here. If the decrease in source brightness was achromatic (i.e., if the flux density were to decay at the same rate for all frequencies), the spectral indices of the two images would be equal within the error bars. However, Fig. 1e shows that the spectral indices of the two images are systematically different at almost all epochs, the spectral index of the NE image being about ∼0.05 higher than that of SW. One natural way to explain this difference in spectral index is to have a slower brightness decrease at higher frequencies, such that the leading NE image always has a slightly flatter spectrum than the SW image.

Here we propose a simple parametric model of the blazar continuum emission able to simultaneously explain the monotonic flux-density decay and the (small) difference in spectral index between the NE and SW images. The model assumes that the source brightness consists in two components: one that remains stable over time, that is, the core, and another that cools down adiabatically, that is, a moving feature, hereafter called a plasmon, injected in a conical jet. The total flux density in the time frame of the quasar I(ν, t) is then described as

(3)

where I_c is the flux density of the core (assumed to have a flat spectrum), I_p is the flux density of the plasmon at t = t₀ and ν = ν₀, α_p is the spectral-index difference between the (optically thin) plasmon and the core, and β accounts for the brightness decrease due to the adiabatic expansion of the plasmon. For sake of simplicity, we do not consider radiative losses, but only adiabatic expansion, meaning that α_p remains constant over time (i.e., time and frequency dependencies are decoupled).

If the plasmon propagates downstream of a conical jet at a constant speed, the β index and I_p are related to the initial distance between the plasmon and the jet base, as well as to the propagation speed, because (for adiabatic losses) the energy of each synchrotron-emitting electron will decrease as (Pacholczyk 1970)

(4)

where r and v are the distance (to the jet base) and the propagation speed, respectively, r₀ is the initial distance, and θ is the jet opening angle. Integrating this equation for an initial electron population N ∝ E^−p (fixing p = 3) and computing the resulting synchrotron spectrum (with the assumption that r ≫ r₀) results in , where k depends on the choice of initial time, t₀.

Using Eq. (3), we can estimate the spectral-index difference between the NE and SW images at any epoch:

(5)

with

(6)

and

(7)

where Δt is the time delay between NE and SW images, ν₂ and ν₁ are the frequency pair used to compute the spectral index, and t is the time for the NE image. The time at the observer’s frame, t′, would be stretched by the factor (which also affects Δt) given by

(8)

where z is the redshift of the quasar and δ is the angle between the jet and the line of sight. In our parametric model (Eq. (3)), the values of v, δ, and θ are all degenerated into β, I₀, and the choice of t₀. Hence, the time evolution of the model is only affected by β and Δt (in the observer’s frame).

We show the results of this simple model simultaneously fitted to the light curve and spectral-index differences in Fig. 4. The corner plot results of a Markov chain Monte Carlo (MCMC) exploration of the plasmon injection time, I_c, I_p, β, α_p, the time delay, Δt, and the differential magnification are shown in Fig. 5. In this MCMC exploration, we have added a 10% systematic uncertainty to the flux densities and spectral indices, and have computed the error function using the flux densities of NE and SW images, the spectral index of NE, and the difference in spectral index between NE and SW. Even though the model is not able to reproduce the scatter seen in spectral index, it is remarkable that such a simple model can capture the light-curve decrease and the (overall) spectral differences between the NE and SW images in a consistent way. The values of the fitting parameters obtained from the MCMC exploration are given in Table 4.

Fig. 4. Fit of the light curve, spectral index, and spectral-index difference between the NE and SW images using the parametric model provided by Eqs. (3) and (5). The envelope of extreme model values for the differential polarization quantities is plotted against measurements using the model extension described in Sect. 3.7. The fit also provides the spectral index of the NE image (red points).

The constraint on the time delay, Δt, is highly model dependent (i.e., the choice of the time evolution in Eq. (3)) and also changes if the relative weight between flux densities and spectral indices is changed in the error function. For instance, if the weight of the spectral index is increased or decreased by a factor 2, the value of Δt decreases or increases by ∼3 or 7 days, respectively. The fitted value of Δt should therefore be taken with care. Still, the best fit value of Δt = 25 ± 3 days is consistent with lensing models (Muller et al. 2020a). The differential magnification between the NE and SW images, on the other hand, is estimated to be 1.24 ± 0.04 and (together with the estimate of Δt) falls within the valley of lowest Von Neumann randomness shown in Fig. 3. The breaking of the time delay–differential magnification degeneracy that was noticed from the light-curve analysis is mainly due to the added information from the spectral index behavior, as Δt is also present in Eq. (7).

The main assumptions made for our modeling are as follows:

• We neglect effects from micro- and milli-lensing, and consider only intrinsic variations of the blazar.

• The emission from the core, I_c, is stable in time and is considered to be optically thick (α_c = 0).

• A single event of plasmon injection.

• A straight trajectory of the plasmon (i.e., a ballistic plasmon), so that there are no changes in the Doppler boosting factor due to changes in the viewing angle.

• Radiative losses are negligible compared to adiabatic losses.

• We fix the jet orientation and opening angle (the actual values of these quantities are degenerated into the fitting parameters).

• The value of the core flux density is bound to be lower than the lowest flux density measured in the last epoch.

Our two-component model for the quasar continuum emission already provides a qualitative explanation for the evolution of the differential polarization, ℛ_pol. As the plasmon cools down as it propagates and expands in the jet, its contribution to the total flux density becomes smaller. Toward the end of our monitoring, the core (which would have a stable emission at timescales larger than the time delay) starts to dominate, and therefore produces similar polarization in the two images. Therefore, the differential polarization between the NE and SW images would vanish and ℛ_pol would tend to zero, as we indeed see in Fig. 1d.

3.7. Modeling with polarization evolution

We propose an extension to the model given in Eq. (3) to further account for the evolution of the polarization properties of PKS 1830−211 during the monitoring. The model starts with the same properties as described in the previous section, that is, a stable core and an evolving plasmon, and adds parameters for the polarization properties. We assume that the core fractional polarization, m_c, evolves slowly, that is, over a much longer timescale than the time delay, and is therefore taken as a constant in the model. On the other hand, we allow the fractional polarization of the plasmon, m_p, to vary in time. Accordingly, in addition to Eq. (3), we describe the evolution of the Stokes Q parameter (in the antenna frame) as:

(9)

where ψ is the parallactic angle, ϕ_c and ϕ_p are the core and plasmon EVPAs, respectively, and F(t) is a function that must decrease quickly with time to match the observed behavior of ℛ_pol toward the end of the monitoring without increasing to very high values early on in the life of the plasmon. All the required model observable quantities (i.e., XX/YY and ℛ_pol) can then be computed from Stokes Q. By construction, the fit of the light curve (which uses total intensity) is decoupled from the polarization part, meaning that the parameters shown in Fig. 5 will not change from the effects of Eq. (9).

Fig. 5. Results from our plasmon modeling (Eq. (3)) as obtained from MCMC (see also Table 4).

For the modeling of the data using Eq. (9), we notice that the parallactic angle coverage of our data is very limited, which prevents us from putting robust constraints on the core and plasmon EVPAs. Hence, we can only test the model of Eq. (9) and try to recover a qualitative behavior of our differential polarization measurements. For the sole purpose of a qualitative exploration, we set F(t − t₀) to be Gaussian, with a width of 100 days, m_c = 0.05, and m_p = 0.2. The value of m_c is taken to be equal to (XX − YY)/I toward the end of the monitoring, where we assume that the contribution from the plasmon polarization becomes negligible. The result of this model is shown in Fig. 4 (second and third panels). The lines shown are envelopes that account for every possible value of ϕ_c and ϕ_p, and assume that ϕ_c of the core is parallel to the parallactic angle in the last epochs of observation. From a qualitative point of view, the model recovers the basic behavior of converging toward a constant value of (XX − YY)/I at late epochs (for both NE and SW), which in turn results in a null ℛ_pol at these epochs. Unfortunately, it is not possible to constrain the model further, as (besides the limited parallactic-angle coverage of our observations) the EVPA of the plasmon may be changing at a faster rate than our time sampling.

4. Results from spectral lines

The dramatic drop in flux density of PKS 1830−211 observed during our monitoring implies that there was substantial activity in the quasar, possibly with strong modifications of its continuum morphology. Only VLBI would be able to directly pinpoint such structural changes (see, e.g., Garrett et al. 1997). Nevertheless, time variations of the absorption line profile from the foreground absorber due to a varying continuum illumination can also provide some insight into the quasar’s activity and the absorbing screen properties.

Comparisons of a collection of HCO⁺J = 2 − 1 (Muller & Guélin 2008; Muller et al. 2014) and CS J = 1 − 0 (Schulz et al. 2015) absorption spectra have already revealed drastic changes over timescales of months and years. However, those observations were irregular and with consecutive visits spread over up to several months or even years.

Our dedicated 2016 monitoring with ALMA now allows us to track and investigate potentially subtle changes over a continuous six-month period and over timescales of down to ∼10 days. We targeted the fundamental ground-state transitions of ortho-H₂O and CH (two Λ-doublets J = 3/2 − 1/2), which can be observed simultaneously in one ALMA tuning in band 7. Those lines were observed by ALMA in 2012 (Muller et al. 2014) and revisited in 2014 (Table 2). The H₂O absorption is deeply saturated toward the SW image near v ∼ 0 km s⁻¹ (Muller et al. 2014), but its large opacity is well suited to investigating the variations of the line wings. The H₂O line is also well suited to probing changes in the presumably optically thin NE absorption. On the other hand, the CH absorption reaches an apparent peak opacity of τ ∼ 1 along the SW line of sight, complementing the H₂O line in tracing variations in the line center of the SW absorption.

The collection of absorption spectra of the H₂O ground-state transition and of the two CH Λ-doublets is shown in Figs. B.1 and B.2 for each visit of our 2016 monitoring for the SW and NE lines of sight, respectively. The rms noise values for each visit –normalized to the continuum level– are given in Table 1 (last column) for the spectral window centered on the water line, which showed the lowest sensitivity. The two other spectral windows centered on the CH Λ-doublets show 10%–20% improved sensitivity. The lowest rms noise levels, of namely ∼0.5% of the flux of the SW image (at the native channel resolution), were achieved at the start of the monitoring, when PKS 1830−211 had a high flux density. The rms noise increased to ∼1% toward the end of the monitoring due to the quasar’s drop in flux density. Only for one visit (2 June 2016) was the noise level up to ∼2% per channel due to poorer observing conditions. Hence, the ALMA monitoring provides us with the opportunity to trace variations of the absorption profile between two visits down to a level of about 1% of the continuum level, at a velocity resolution of < 1 km s⁻¹.

4.1. Long-term variations (2012–2022)

Figure 6 shows all the available spectra of the H₂O (SW and NE lines of sight) and CH (SW line of sight only) lines taken with ALMA between 2012 and 2016 (observation details are listed in Tables 1 and 2). The spectra with large deviations with respect to the majority are highlighted with a special color code.

Fig. 6. Overlay of all ALMA spectra of H2O (SW and NE lines of sight) and CH (SW only) observed toward PKS 1830−211 between 2012 and 2016. Special profiles showing clear deviations are highlighted with the color code given in the upper box.

The main variations for the H₂O spectra toward the SW image do not occur in the line center near v ∼ 0 km s⁻¹, where the absorption is saturated, but in the wings of this component. The blue wing (i.e., v ∼ −20 to −60 km s⁻¹) appears to have been stable between the 2012 and 2014 spectra, but jumped to a new level between 2014 and 2016. In contrast, the red wing (i.e., v ∼ +20 to +60 km s⁻¹) shows more systematic variation, with the deepest absorption level being reached at the beginning of our survey in March 2016. This is not the case for the v ∼ +170 km s⁻¹ component, which was deepest in 2014 (Fig. 7). At the same time, the CH profile of the SW absorption does not show significant variation, except for some changes between 2012 and later observations.

Fig. 7. Zoom onto the variations of the +174 km s−1 (top) and −225 km s−1 (bottom) features in the H2O SW and NE spectra, respectively. Same as for Fig. 6; special profiles showing clear deviations are highlighted with the color code given in the upper box.

Concerning the water absorption toward the NE image, there was a drastic decrease in the v = −150 km s⁻¹ component between 2012 and 2014, whereas the subsequent spectra from 2014 to 2016 remain relatively unchanged. This major change is reminiscent of what happened between 2003 and 2006 in observations of the HCO⁺ (2−1) line by Muller & Guélin (2008), with the near disappearance of the same velocity component during this time interval. For the v = −225 km s⁻¹ feature, unnoticed prior to 2009 ATCA observations by Muller et al. (2011), the H₂O absorption varies continuously between 2012 and 2016 (Fig. 7). To these water spectra, we can add the observations of the CH⁺ (1−0) line in 2015 (reported by Muller et al. 2017). If we scale down the optical depth of CH⁺ by a factor of ∼12 (see Fig. 8), its profile matches the v = −150 km s⁻¹ feature relatively well, but is not consistent for the v = −225 km s⁻¹ feature, suggesting that one or the other significantly varied again between 2015 and 2016.

Fig. 8. Comparison of the average H2O spectrum (filled blue) observed in 2016 with the ALMA CH+ (1−0) spectrum (orange line) observed in 2015 (Muller et al. 2017) for both lines of sight. For the NE spectrum, the opacity of the CH+ line, assuming a filling factor of unity, is reduced by a factor 12.

Finally, there are a handful of observations of the rare water isotopologs HO and HO between 2014 and 2022 (Table 2 and shown in Fig. 9). Those reveal a spectacular increase in the HO absorption for the v ∼ −5 km s⁻¹ component in the July 2019 spectrum, when PKS 1830−211 was enduring a record-breaking γ- and radio flare (see Fig. 2b, and, e.g., Martí-Vidal et al. 2020). The integrated opacity of this component nearly increased by a factor four, while the absorption near v = +5 km s⁻¹ did not change significantly. In particular, this remarkable change in absorption, together with the increase in the flux density of PKS 1830−211, allowed the detection of several deuterated species (Muller et al. 2020b). Interestingly, the most recent spectrum of HO obtained in August 2022 shows a return to almost the same profile as 2014, suggesting that the 2019 change was indeed due to an unusual event in PKS 1830−211, possibly related to a compact structure reminiscent of a typical Galactic dark cloud, as previously discussed by Muller et al. (2020b).

Fig. 9. Overlay of available ALMA spectra of HO and HO 11, 0 − 10, 1 toward the SW image of PKS 1830−211.

4.2. Time-variation statistics

The extremes of the 2016 H₂O and CH absorption spectra are shown in Fig. 10. The saturated region of the water absorption toward the SW image (v ∼ 0 km s⁻¹ shows almost no variation, with the variation that we do see being due to a slight change of covering factor. It is interesting to note that the largest variations (up to ∼5% of the continuum level) occurred in the line wings, especially in the red wing (v ∼ +40 km s⁻¹), but are not symmetric with respect to the line center. On the other hand, the variations seen in the CH spectra toward the SW image, of the H₂O v ∼ +170 km s⁻¹ component toward the SW image, and of the H₂O spectrum toward the NE image (i.e., all optically thin) are limited to an rms dispersion of ≲1%.

Fig. 10. Extrema and rms dispersion of the H2O and CH absorption spectra during our 2016 monitoring (the CH profiles are not shown for the NE line of sight, because they contain only noise).

Furthermore, we show the difference between consecutive spectra divided by their respective time separation (i.e., rate of change) in Fig. 11. The strongest changes – which are of up to ∼1% per day on average – appeared at the beginning of the monitoring. Nevertheless, variations are not restricted to single time events, but apparently occurred continuously during the monitoring.

Fig. 11. Rate of the variations between consecutive observations of the H2O absorption profile along the SW line of sight of PKS 1830−211.

4.3. Principal component analysis

To extract more information from the variations of the absorption line profile, we ran a principal component analysis (PCA). This method searches for maximum variance through the data and outputs the results in new orthogonal axes, that is, principal components (hereafter PCA eigenspectra) with their associated (eigen)values in order of decreasing variance. We used the sklearn.decomposition.PCA class in the scikit-learn Python package (Pedregosa et al. 2011) to retrieve the eigenspectra (i.e., as a function of velocity) and eigenvalues (i.e., as a function of time) of the monitoring spectra. The results of the PCA decomposition for the H₂O spectra toward the SW image are shown in Figs. 12 and 13 for the eigenvectors and eigenvalues, respectively. The three most significant PCA components consist in spectral features in the blue and red wings of the absorption profile, which resemble to some extent the difference spectra shown in Fig. 11. Most interestingly, the evolution of the three first eigenvalues (all together describing more than half the variance of the data) does not appear to be random but can be well described by a sinusoidal function across the entire duration of the monitoring. Fitting such a sinusoidal function, we find periods of 272 ± 15 days, 186 ± 20 days, and 117 ± 12 days, for PC1, PC2, and PC3, respectively. These periods are comparable to the duration of our monitoring and it remains to be seen if the evolution is truly periodic with monitoring over longer timescales. Using the spectra of the CH doublets along the SW line of sight instead of H₂O, or the H₂O spectra along the NE line of sight, we also find apparent wavy behaviors for PC1, as shown in Fig. 14. In addition, we separately checked for the velocity intervals corresponding to the blue and red wings of the H₂O absorption along the SW line of sight, and, again, find wavy signals. A final test, taking a velocity range outside of any absorption features (i.e., purely noisy channels), shows no significant structure in the PCA decomposition. All these measurements (Table 5; except overlapping velocity range for measurements from H₂O spectra) are independent, but still yield a consistent signal period of ∼230 days for PC1, within error bars.

Fig. 12. Eigenspectra of the PCA decomposition of the H2O spectra toward the SW image of PKS 1830−211.

Fig. 13. Time evolution of the eigenvalues of the PCA decomposition of H2O absorption spectra toward the SW image of PKS 1830−211.

Fig. 14. Time evolution of the eigenvalues of the first component (PC1) of the PCA decomposition of various subsets of the monitoring spectra. The line of sight, species, considered velocity interval, and PC1 scores are given in each box.

In order to test the PCA performances, we ran experiments synthetic spectra with similar noise properties and time sampling as in the ALMA monitoring, in which we introduced a time varying component (with a given period) in the form of a Gaussian perturbation with a given amplitude, full width at half maximum (FWHM), and velocity. We ran those synthetic spectra into the same PCA decomposition process as for the real observed spectra and the same fitting exercise with a sinusoidal function for the time evolution of the PC eigenvalues. We find that the first PC component PC1 succeeds in correctly retrieving (i.e., within error bars of the fit) the period for perturbation signals with S/N down to a few percent. In other words, the PCA, and particularly its PC1, shows a good performance in retrieving this signal.

As it is difficult to imagine any periodic behavior in the absorbing clouds simultaneously affecting both lines of sight and both H₂O and CH spectra, it is tempting to conclude that this wavy signal originates from the background illumination. This evolution could be connected to the model of a helical jet proposed by Nair et al. (2005), for which the precession period was estimated to be of one year. This interpretation would indeed bring a natural explanation for the time variability of the absorption lines toward PKS 1830−211, but would need to be tested with new monitoring data over a longer time span. In this scenario, we would observe a slow modulation of the absorption due to the jet helicity, with occasional bursts or dips of absorption, depending on whether or not there is favorable illumination, whenever a new plasmon or a local flare occurs in the jet.

4.4. 2016 averaged spectra

We take advantage of the multiple visits of the 2016 monitoring to produce deep average spectra of H₂O and CH absorptions toward both images. To achieve the best sensitivity and quality, we removed a spectral baseline as a first-order polynomial (i.e., corresponding to the intrinsic spectral index of each lensed image) on each individual spectrum. We then averaged all spectra after weighting them individually by the square invert of the rms noise levels listed in Table 1. The resulting high-sensitivity spectra have a noise level of better than 0.1%−0.2% of the normalized continuum level, as shown in Figs. 15 and 16. These “super” spectra lead us to two discoveries along the SW line of sight: the detection of the weak absorption from the ¹³CH isotopolog and a remarkable absorption trough of 500 km s⁻¹ in width, both of which are discussed in the following subsections.

Fig. 15. Weighted-average spectra toward the SW image of PKS 1830−211, combining all visits in 2016. The hyperfine structure is shown for each CH Λ-doublet, with the strongest hfs component set at a velocity v = 0 km s−1, or as indicated otherwise. There are two boxes for each line, the first one showing the whole absorption spectrum and the second one showing a smaller range around the continuum level. Best-fit of the CH and 13CH absorption is shown in red.

Fig. 16. Weighted-average spectra toward the NE image of PKS 1830−211, combining all visits in 2016. The hyperfine structure is shown for each CH Λ-doublet, with the strongest hfs component set at a velocity v = −150 km s−1. The lower panel shows the CH hfs-deconvolved profile.

4.4.1. Detection of ¹³CH

The rare isotopolog ¹³CH was recently detected for the first time in the ISM of the Milky Way by Jacob et al. (2020). As CH is ubiquitous in the ISM, with optically thin lines, and is expected not to be affected by fractionation in its chemical formation, it appears to be an excellent species for measuring the ¹²C/¹³C elemental isotopic ratio (see a detailed discussion by Jacob et al. 2020). This isotopic ratio is an interesting probe of stellar nucleosynthesis history. Chemical-evolution models (e.g., Prantzos et al. 1996; Kobayashi et al. 2011) predict a positive ¹²C/¹³C gradient in the Milky Way, increasing with galactocentric distance and decreasing with time.

Here, we present the first detection of ¹³CH toward the SW image of PKS 1830−211 in the average spectra of the 2016 monitoring data. The detection is robust, with its typical fine structure (see Table A.1) present in both Λ-doublets, as shown in Figs. 15a and b. Nevertheless, the ¹³CH absorption signal is very weak, partly due to the optical thinness of CH, to the large ¹²C/¹³C isotopic ratio, and to the dilution of the signal among the many hyperfine components. We were able to obtain a detection in absorption toward PKS 1830−211 thanks to the brightness of the quasar and the high sensitivity of ALMA, providing a very high signal-to-noise ratio after averaging all the monitoring spectra.

In order to estimate the ¹²CH/¹³CH ratio, we performed a simultaneous fit of the two Λ-doublets of ¹²CH and ¹³CH for the SW absorption spectrum using the same intrinsic velocity profile for both CH and ¹³CH convolved to their corresponding fine structure, respectively, and with the ¹²CH/¹³CH ratio as free scaling parameter. We made several attempts with increasing numbers of Gaussian velocity components to reproduce the absorption profile of CH, and obtained satisfying results using about ten components. During the process, we found that the ratio was relatively stable around approximately 150, while the uncertainties and fit residuals decrease with an increasing number of velocity components. Our final ¹²CH/¹³CH ratio is 150 ± 10. The physical relevance of every single velocity component in the fit is not certain, but it is clear that a mix of narrow and broad components is needed to reproduce the profile, in particular with some large components for the wings of the main v = 0 km s⁻¹ absorption feature.

Compared to the ¹²C-/¹³C- isotopolog ratios from other species in the SW line of sight toward PKS 1830−211, the ¹²CH/¹³CH ratio is the highest measured so far: 97 ± 6 for ¹²CH⁺/¹³CH⁺ (Muller et al. 2017), 62 ± 3 for CH₃OH/¹³CH₃OH (Muller et al. 2021), and even lower values 20 − 50 for HCO⁺, HCN, and HNC (Muller et al. 2006, 2011). This strongly suggests that all these other species may be affected by fractionation issues. A detailed analysis of the ¹²C-/¹³C- isotopolog ratios is beyond the scope of this study and will be presented elsewhere.

A ¹²C/¹³C ratio ∼150 in the z = 0.89 absorber is higher than the values observed within the inner 16 kpc of the Milky Way (see, e.g., the measurements of ¹²C/¹³C with galactocentric distance in Fig. 4 by Jacob et al. 2020); as such, it probably reflects a true evolution effect between the z = 0.89 absorber and the Milky Way today, as also found for other isotopic ratios (e.g., from N, O, S, Si, Muller et al. 2006, 2011).

4.4.2. Weak and broad water absorption trough in the SW sightline

In addition to the long-known main absorption features at v ∼ 0 km s⁻¹ and v ∼ −150 km s⁻¹ toward the SW and NE images, respectively (Wiklind & Combes 1998; Muller et al. 2006), Muller et al. (2011) discovered narrow velocity components at velocities of v = −300, −224, and +170 km s⁻¹ in an unbiased 7 mm absorption survey with the Australian Telescope Compact Array. The −224 and +170 km s⁻¹ components were later identified in ALMA data as occurring toward the NE and SW image, respectively (Muller et al. 2014), but the v = −300 km s⁻¹ component was not confirmed. However, two narrow velocity components at v = −320 and −280 km s⁻¹ were later observed in 2015 in the CH⁺ (1−0) ALMA spectrum toward the SW image (Fig. 8 and Muller et al. 2017), although near the edge of the spectral window. Both components are now readily confirmed in the H₂O “super” spectrum toward the SW image, as can be seen in Fig. 15c.

The origin of these components with large velocity spans in the line of sight is unclear, and they could be due to a high inclination of the absorber, large velocity gradients in the disk, or an extra-planar gas component, such as high-velocity clouds in the halo as observed around the Milky Way (as also discussed previously by Muller et al. 2014). They are not explained in the absorption kinematic model by Combes et al. (2021), which reproduces the H I and OH absorption spectra seen at cm wavelengths and the molecular absorption at mm wavelengths.

In addition, the ALMA “super” spectrum of the H₂O absorption now reveals an even more surprising structure in the form of a weak (τ < 0.01) continuous absorption trough covering a wide velocity range of nearly 500 km s⁻¹ toward the SW image (Fig. 15c). This is much wider than the continuous absorption span of ∼200 km s⁻¹ along the NE line of sight. It could be that this absorption trough reveals the presence of circumgalactic gas (e.g., Bahcall & Spitzer 1969; Tumlinson et al. 2017), although with up to such high velocities, it is not clear whether or not the gas would remain tied to the galaxy, depending on the halo-escape velocity, or would then (literally) form a “water fountain”. Given the lack of constraints on the distribution of such material in the halo (probably fragmented in filamentary structures), it is not possible to get a meaningful estimate of its mass content. In any case, we again emphasize the weakness of this signal, which would be extremely difficult to trace in other galaxies unless a background continuum source as bright as PKS 1830−211 were present.

4.5. Water isotopologs and oxygen isotopic ratios

With different water isotopologs observed close in time in the 2014 and 2022 datasets (Table 2), we can derive the HO:HO:HO isotopolog ratios of the water family. Those are not expected to be affected by fractionation effects and should therefore reflect the true elemental isotopic ratios of oxygen. For a given epoch, we assumed the same opacity profile (starting with five Gaussian velocity components) for all water isotopologs, as well as for CH when available in order to make the link between the high-opacity line of HO and the optically thin rare water isotopologs. We then fitted the opacity ratios as free scaling factors. When including the main water isotopolog, we also added the continuum-covering factor as a free parameter of the fit; otherwise, we fixed it to unity for the fit of the HO and HO isotopologs alone. The results of our fitting exercises are summarized in Table 6. Due to the somewhat better signal-to-noise ratio of the data in July 2014, we used a higher number of Gaussian velocity components to obtain a better match of the profile and decrease the reduced χ-squared. However, we found that the fitted ratios were not drastically affected and remained similar within the uncertainties. The values of the ¹⁶O/¹⁸O (=65.3 ± 0.7 after combining the measurements with their statistical weights) and ¹⁸O/¹⁷O (=11.5 ± 0.5) ratios are consistent with the previous measurements made from the HCO⁺ isotopologs by Muller et al. (2006, 2011), albeit with much better accuracy. Regarding the ¹²C/¹³C ratio discussed above (Sect. 4.4.1), the oxygen isotopic ratios point to a different ISM isotopic composition from that of the solar neighborhood in the Milky Way today (e.g., Lucas & Liszt 1998), namely one that is less processed and most likely dominated by massive-star nucleosynthesis (see also, e.g., Wallström et al. 2016; Martín et al. 2019; Tang et al. 2019 and references therein).

5. Discussion

5.1. Analysis of the absorption profiles and their time variations

The absorption spectrum S​_ν resulting from an absorbing screen with opacity τ_ν(x, y) (x and y being in the plane of the sky) in front of a background continuum illumination i(x, y) and normalized to it is

(10)

where I_tot = ∬i(x, y)dxdy is the total continuum intensity. Therefore, the normalized absorption spectrum is the continuum intensity-weighted average of the exponential of the opacity distribution of the illuminated material: S = ⟨e^{−τ(x, y)}⟩_i. In the case where the continuum emission is spatially extended but unresolved by the observations, the distributions i(x, y) and τ(x, y) are unknown, and we can only measure I_tot and the resulting normalized absorption spectrum S​_ν. If we consider a representative average opacity , Eq. (10) can be expressed as:

(11)

where we introduce the source-covering factor f_c, which is defined as

(12)

where Ω_abs is the solid angle fraction of the background continuum actually covered by absorbing material (weighted by intensity); by construction: 0 ≤ f_c ≤ 1.

The most straightforward way to measure the covering factor is to use the saturated region of a very optically thick line, such as the H₂O 1₁₀ − 1₀₁ in our case (τ ≳ 10). Figure 17 shows the measurements of f_c across the different visits of our 2016 monitoring. The covering factor is ∼96% with only mild variations of the order of 1%, which contrast with the large variations of the flux density between the beginning and the end of our monitoring (Fig. 1). The values of f_c are consistent with previous measurements of saturated lines (Muller et al. 2014).

Fig. 17. Measurements of the covering factor across the different visits of our 2016 monitoring. (a) Saturation level (1 − S​ν) at the center of the H2O line for all the visits in 2016. (b) The values of all spectral channels with velocities |v|< 8 km s−1 are averaged to provide measurements of the continuum-covering factor as a function of time.

To further characterize the absorption variations, we can assume that the absorbing screen is constant in time, that is, dτ(x, y)/dt = 0. This is justified by the fact that the timescales associated with potential causes of variations in the absorber are much longer than those we are interested in, even for small-scale structures on AU- to pc-scales. For example, a structure of 10 AU with a velocity of 200 km s⁻¹ in the plane of the sky (e.g., due to rotation for a face-on spiral galaxy) would have a crossing time in front of a fixed point-like continuum of three months. For a cloud collapse, the free-fall timescale is , where R is the cloud size and M its mass. Finally, we can also exclude variations due to a change of the chemical composition of the absorbing material, because typical chemical timescales are orders of magnitude longer than years (e.g., Tielens 2005; Valdivia et al. 2017). Therefore, we need only consider time variations due to structural changes of the quasar’s morphology, that is, change of i(x, y) with time.

If the absorbing screen is uniform, where τ(x, y) = τ₀, there can be no variation of the normalized absorption spectrum, regardless of variations of i(x, y). Consequently, the mere observation of variations of the absorption spectra implies that there are opacity structures in the absorbing screen at scales smaller than the continuum size. Similarly, a simple change of i(x, y) by a multiplicative factor would not induce variations of the absorption spectrum due to the normalization by I_tot, and therefore variations of S​_ν also imply small-scale structures in the quasar’s continuum emission. Hence, the absorption variations reflect the coupling of small-scale structures between the continuum i(x, y) and absorbing screen τ(x, y).

Taking Eq. (11) as a simplified description of the absorption spectra and assuming that all continuum variations result in a change of the average opacity, then the derivative dS​_ν/dτ yields an estimate of the variation; we get

(13)

The function xe^−x peaks sharply for x = 1 and tends rapidly to zero for large x, explaining why the maximum variations of the water absorption are seen in the line wings but are strongly attenuated at the saturated line center.

5.1.1. Possible causes of time variations of the continuum distribution

In the long term (months to years), we have seen some drastic changes in the absorption profiles, such as the near disappearance of the NE absorption near −150 km s⁻¹ between 2003 and 2006 (Muller & Guélin 2008), its drastic decrease between 2012 and 2014, or the doubling of the absorption depth for the v = −5 km s⁻¹ velocity component in the HO SW profile observed in 2019 (Fig. 9 and Muller et al. 2021). This last event was occurring at the time of a record-breaking flare of PKS 1830−211 at γ-rays and radio. It is possible that one or another cloud becomes occasionally illuminated (or highlighted) by a new continuum component. This “occasional event” scenario is consistent with the fact that the last HO spectrum in 2022 shows that the absorption returned to almost the same profile as pre-2019, which we could then qualify as a “quiescent absorption profile”. Nevertheless, at this stage, we cannot exclude the possible effect of milli-lensing on timescales of months to years. In 2019, the flux-density ratio was indeed in the low range of values (i.e., either the NE image was fainter than usual with respect to the SW image, or the SW image was brighter than the NE image). We would need more statistics of drastic absorption-change events to investigate a possible connection with γ-ray and radio flares.

On shorter timescales of weeks to months, the effect of micro-lensing on the variability of absorption spectra was investigated by Lewis & Ibata (2003). The picture is highly complicated by the large parameter space defined by absorbing cloud properties (shape, size, absorption profile) and distribution with respect to the background continuum illumination (all unknown), but Lewis & Ibata (2003) validate the fact that micro-lensing can potentially induce significant modulation of the absorbing signal provided absorbing clouds have subparsec-scale structures.

An object of mass M in the lens plane can produce significant variability by (micro-)lensing for a typical source size below the Einstein radius:

(14)

where D_ij are the angular diameter distances between the observer o, the lens l, and the source s, respectively. For PKS 1830−211 and a stellar-mass object, θ_E ∼ 2 μas.

The size of the SW core was measured to be θ_SW ∼ 0.1 mas at λ ∼ 1 cm and to scale with a ∼λ² dependence at radio–cm wavelengths (Guirado et al. 1999). Extrapolating this behavior to mm wavelengths, we would obtain a size of 1 μas at λ ∼ 1 mm. However, broadening due to interstellar scattering may have a far less significant effect at mm/submm wavelengths. In fact, Muller et al. (2021) estimated a flatter dependence in λ^0.3 from the covering factors derived from a set of submm methanol lines. Extrapolating the size measurement at 1 cm with this flatter dependence, we obtain a size estimate of about 50 μas at 1 mm³, which is roughly one order of magnitude larger than θ_E. Therefore, it is possible that micro-lensing induces variability on timescales of weeks to months, which corresponds to the caustic crossing time (Lewis & Ibata 2003).

On the other hand, the scenario of a helical jet would naturally induce variations of the continuum illumination with a clear periodicity. In turn, this periodicity may be imprinted in the variations of the absorption profile, which could be what we observe in the 2016 spectra, as suggested by our PC analysis (Sect. 4.3). In Sect. 3.6, we propose a simple parametric model of a single ballistic plasmon in the jet to account for the continuum evolution during the 2016 monitoring. The ballistic assumption would not necessarily be in contradiction with the helical jet scenario. The toy model would also work even if the plasmon path is helical within the jet and its opening angle. There could be a time-changing Doppler boosting that would be added to the power-law decay of the toy model (depending on how open the helical path is and how close the jet line is to the line of sight), but the effect could be hidden below the flux-density uncertainties. Moreover, the facts that (i) we see a clear sinusoidal signal in the PCA and that (ii) the largest variations of the absorption profile happened at the beginning of the monitoring may support the single plasmon injection event. Indeed, supposing there were more plasmons at different points of the helical jet, one would see a “dirty” superposition of sinusoids with different phases. Therefore, having one simple sinusoidal signal in the absorption variation may be telling us that the illumination pattern of the absorbing clouds was simple during the 2016 monitoring.

Given the quasi-monotonic and severe decrease in the flux density of PKS 1830−211 during our six-month monitoring (by ∼80%), the short (timescale of the order of one week) and regular variations of the absorption spectra, and the apparent coherence (periodicity?) of the variations over ≳200 days seen toward both lensed images, we conclude that they are most likely due to intrinsic changes in the background continuum distribution, that is, the activity of the quasar, rather than to micro-lensing events. Eventually, the different scenarios (helical jet or lensing-induced variations) could be tested with long-term (i.e., more than one year long) monitoring of the absorption variation and its potential periodicity (e.g., with ALMA) and continuum structural changes (e.g., with VLBI), although flaring or micro- and milli-lensing events could occasionally confuse a periodic signal induced by the jet helicity.

5.1.2. Variations on short timescales

We now consider variations on short timescales – that is, between two consecutive observations at t₁ and t₂ – resulting from intrinsic morphological changes in the blazar. Accordingly, we assume that the characteristic size scale, a causal region of size c(t₂ − t₁), is small compared to the overall continuum emission.

We can then introduce the continuum distribution i₂(x, y, t₂) = i₁(x, y, t₁)+δi(x, y), where δi(x, y) = 0 outside of the small region of the continuum that was indeed affected by changes between t₁ and t₂. We do not need to assume any specific shape for this region, except that it is strictly included within a circle of radius c(t₂ − t₁) or θ_E centered toward the direction x_c, y_c in the sky. We express the new normalized absorption spectrum S​₂ from Eq. (10): . Rewriting the equation for δI/I₁ ≪ 1, and neglecting the second-order terms in δI/I₁, we get:

(15)

Furthermore, we assume that the absorbing material illuminated by the continuum region δi(x, y) has homogeneous opacity τ_c = τ(x_c, y_c). After simplifying the integral summation, we obtain:

(16)

As already mentioned, ⟨e^−τ⟩_i is the continuum intensity-weighted average of the exponential of the opacity distribution. Therefore, in order to observe variations between normalized absorption spectra, the necessary condition is that the opacity of the absorbing region behind which the continuum changed is different from the overall average opacity; in other words, that the absorbing material has substructures at scales smaller than the continuum illumination solid angle (Ω_abs).

Because the function e^−τ takes values strictly between 0 and 1, the most favorable cases to detect variations according to Eq. (16) are: when ⟨e^−τ⟩_i ∼ 1 and e^−τ_c ∼ 0, that is when a small and optically thick region is suddenly illuminated over an absorbing medium of mostly very low opacity; or when ⟨e^−τ⟩_i ∼ 0 and e^−τ_c ∼ 1, that is when the sudden change of continuum occurs behind an optically thin region, while previously the bulk of the absorption concerned optically thick material. The maximum possible change in the difference of normalized absorption spectra is δI/I₁.

5.2. Correlation H₂O – CH

Assuming that the absorbing material covers the background continuum source entirely (f_c = 1), we can directly convert the absorption profile into optical depth from Eq. (11). With this assumption, it is also possible to deconvolve the hfs pattern of the CH doublets, for example using a Clean algorithm method assuming that the relative intensities of the hfs components follow their statistical weights. Accordingly, we can investigate the optical-depth correlation between the H₂O and CH lines, as shown in Fig. 18 for both the SW and NE lines of sight.

Fig. 18. Correlation between the optical depths of the ortho-H2O 11, 0 − 10, 1 line and of the deconvolved profile of CH 3/2–1/2 (given for an equivalent hyperfine component with normalized intensity S​ul = 1) along the SW (left) and NE (right) lines of sight, respectively. In both lines of sight, we assume a source-covering factor of fc = 1. The colors encode the velocity intervals shown on the line profiles (top). The dashed lines correspond to opacity ratios τ(H2O)/τ(CH*) = 3 and 9, respectively.

In the Milky Way, both H₂O and CH are found to be excellent tracers of the H₂ along diffuse molecular sightlines. Sheffer et al. (2008) derive a quasi-linear relationship between the column densities of CH and H₂, with an average abundance ratio of [CH]/[H₂] = along diffuse molecular sightlines and dark clouds. For water, Flagey et al. (2013) also find a remarkably constant abundance ratio relative to H₂ in Galactic translucent clouds, with [H₂O]/[H₂] = (4.8 ± 0.3)×10⁻⁸. Taking these relative abundances with their uncertainties, and assuming (i) that the relative abundances of H₂O and CH for the absorbing material in the PKS 1830−211 absorber are similar to the Galactic clouds, (ii) an ortho-to-para ratio of three for H₂O, and (iii) that the molecules are at excitation equilibrium with cosmic-microwave-background photons (T_CMB = 5.14 K at z = 0.89, Muller et al. 2013), we would expect opacity ratios τ(H₂O)/τ(CH^*) of between 3 and 9 (at the 1σ confidence level), CH^* being the normalized CH spectrum (equivalent to a hyperfine structure component of relative strength S​_ul = 1) obtained after deconvolution of the hyperfine structure of the two CH Λ-doublets.

Discarding the center of the SW line profile (velocities between ∼ − 10 and +10 km s⁻¹), for which the water line is heavily saturated and has f_c ∼ 95% (which is therefore not consistent with our working assumption f_c = 1), we find that the H₂O/CH opacity ratios along the SW line of sight are roughly around 9 (see also fit results in Table 6), that is, in the upper part of Galactic ratios. Along the NE line of sight, the τ(H₂O)/τ(CH^*) opacity ratios are somewhat lower than in the SW line of sight, although the spread may partly be due to lower signal-to-noise ratios. These ratios are nevertheless still within the range of values observed in Milky Way translucent clouds. Based on the overall similarity between the τ(H₂O)/τ(CH^*) ratios in the PKS 1830−211 absorber and in Galactic clouds, and except for the heavily saturated part of the SW water line, we do not find evidence for large opacity effects in other velocity components. Should the source-covering factor be lower than unity for some velocity components, the true optical depth of the water line would rise faster than that of CH, and the points in Fig. 18 would deviate from the straight correlation in the same way as for the saturated region near v = 0 km s⁻¹. Therefore, regarding H₂O and CH, it seems that the chemical conditions in the different velocity components along the SW line of sight are consistent with a relatively constant H₂O/CH abundance ratio.

Moreover, the differences in H₂O and CH^* opacity ratios between the SW and NE lines of sight follow a similar trend to that observed in the Milky Way between clouds in the central molecular zone (CMZ) and the Galactic disk, for which Sonnentrucker et al. (2013) found a water abundance three to five times higher in the CMZ. The cosmic-ray (CR) ionization rate deduced from a simple analysis of the OH⁺ and H₂O⁺ ions also follows a similar trend between the CMZ and disk of the Milky Way (Indriolo et al. 2015) on one hand, and the SW and NE lines of sight toward PKS 1830−211 (Muller et al. 2016) on the other. As these ions are precursors of H₂O, the somewhat larger H₂O abundance could be related to a higher CR ionization rate.

6. Summary and conclusions

We monitored the lensed quasar PKS 1830−211 with ALMA over a time span of six months in order to investigate the variability of the submm continuum emission of the blazar and of its foreground molecular absorption. We summarize our results as follows. First, for the continuum activity:

• The submm light curve shows a dramatic and nearly monotonic decrease of 80% in the flux density during our monitoring.

• In contrast, the flux-density ratios between the two lensed images of the quasar show only mild variations of less than 10%. There are signs of variations of the polarization properties, which manifest as fluctuations in the differential polarization between the two images. Toward the end of the monitoring, the quasar seems to have reached a stable stage in its polarization variability.

• The nearly featureless light curve prevents us from obtaining a strong constraint on the time delay between the two images. Looking backwards at the previous measurements of the flux-density ratios, it appears that it is also difficult to precisely determine the differential magnification ratio of the system, as the long-term evolution of the flux-density ratios suggests milli-lensing events on timescales of several years.

• We propose a relatively simple parametric model of a single ballistic synchrotron-cooling plasmon to account for the evolution of the continuum emission during the monitoring.

• From this parametric model, we also obtain a new measurement of the time delay between the NE and SW images of Δt = 25 ± 3 days, and of their differential magnification, namely 1.24 ± 0.04, both of which are consistent with values provided by modeling of the system and its lensing geometry.

For the absorption lines arising in the z = 0.89 lens–absorber:

• The absorption line profiles of H₂O (saturated with optical depth > 1) and CH (peak opacity τ ∼ 1) reveal clear time variations (on average up to ∼1% of the continuum level per day for the water line) between spectra taken 1–3 weeks apart.

• A principal component analysis of the absorption spectra toward both images of PKS 1830−211 reveals a seemingly wavy behavior of the variations during the monitoring, with an apparent period of ∼230 days, although this timescale is longer than the time span of our survey and would therefore need to be confirmed with a longer monitoring. It is tempting to relate this period to the jet precession period of approximately one year in the scenario of a helical jet proposed by Nair et al. (2005). This would provide a natural explanation for the time variability of the absorption line.

• Combining all the spectra along the SW image, we obtain a detection of the rare isotopolog ¹³CH, with a ¹²CH/¹³CH ratio of ∼150. From complementary data, we also measure oxygen isotopic ratios and derive ¹⁶O/¹⁸O = 65.3 ± 0.7 and ¹⁸O/¹⁷O = 11.5 ± 0.5.

• We discover a remarkable, broad but shallow absorption trough of 500 km s⁻¹ in width observed in the water line toward the SW image, which could be the signature of extra-planar molecular gas.

• We explore the chemical correlation between CH and H₂O using their optical depth ratio and find a H₂O/CH ratio comparable to that observed in Milky Way clouds.

The lens–absorber system toward PKS 1830−211 acts as a powerful cosmological microscope, allowing us to investigate the coupling between the quasar activity and the variations of the absorption pattern. In particular, whether the absorption variability is periodic – and is therefore most likely connected with the precession of the quasar’s helical jet – or irregular – and is connected to occasional flaring or micro-/milli-lensing events – remains to be confirmed.

Acknowledgments

This paper makes use of the following ALMA data: ADS/JAO.ALMA#2011.0.00405.S (observations of H₂O and CH in April 2012), ADS/JAO.ALMA#2012.1.00056.S (observations of H₂O and CH in May 2014), ADS/JAO.ALMA#2013.1.00020.S (observations of H₂O, CH, HO, and HO in July 2014), ADS/JAO.ALMA#2015.1.00075.S (2016 monitoring of H₂O and CH), ADS/JAO.ALMA#2018.1.00051.S (observations of HO in 2019), ADS/JAO.ALMA#2021.A.00028.S (observations of HO and HO in 2022). ALMA is a partnership of ESO (representing its member states), NSF (USA) and NINS (Japan), together with NRC (Canada), NSC and ASIAA (Taiwan), and KASI (Republic of Korea), in cooperation with the Republic of Chile. The Joint ALMA Observatory is operated by ESO, AUI/NRAO and NAOJ. This work has been partially supported by Generalitat Valenciana (GenT Project CIDEGENT/2018/021 and grant ASFAE/2022/018) and the MICINN (Research Project PID2019-108995GB-C22 and the Astrophysics and High Energy Physics programme, with funding from European Union NextGenerationEU, PRTR-C17I1). This research has made use of NASA’s Astrophysics Data System.

References

Appendix A: Spectroscopic data

Spectroscopic data for the transitions of H₂O and CH isotopologs observed in this study are given in Tab. A.1.

Appendix B: Additional material

We provide an overview of the ALMA continuum-normalized absorption spectra of H₂O and CH toward the SW and NE image of PKS 1830−211 during our monitoring in Figs. B.1 and B.2, respectively.

Fig. B.1. Overview of the continuum-normalized absorption spectra of the CH Jp = 3/2+ − 1/2− (left), CH Jp = 3/2− − 1/2+ (middle), and ortho-H2O JKa, Kc = 11, 0 − 10, 1 (right) lines toward the SW image of PKS 1830−211, for each visit of the 2016 monitoring. The visit numbers #i correspond to the entries in Table 1.

Fig. B.2. Same as Fig. B.1, but toward the NE image of PKS 1830−211.

All Tables

All Figures

Fig. 1. Time evolution of the continuum properties of PKS 1830−211 during our monitoring in 2016 (see Table 1): (a) Flux density of the NE and SW images (error bars of 20% are indicated for the NE image data points only; the color code indicates the source used for flux density calibration). (b) NE/SW flux density ratio. (c) Difference in flux density between the XX and YY polarizers normalized by the Stokes I flux density. (d) Polarization ratio ℛpol (values are encoded as a function of the PA of the observations). (e) Difference in spectral index. The expected time delay of ∼27 days between the two lensed images of the quasar is indicated in panel a.

In the text

Fig. 2. History of ALMA measurements toward PKS 1830−211. (a) Flux-density ratios of PKS 1830−211 NE and SW images (with the color code corresponding to ALMA bands, as indicated at the top of the figure). (b) Flux densities of the NE image of PKS 1830−211, scaled to an equivalent frequency of 300 GHz such as S​300 GHz = S​ν × (ν/300 GHz)α with a spectral index α = −0.7 (error bars indicate a nominal flux uncertainty of 20%) overlaid on top of the Fermi-LAT γ-ray daily-monitoring light curve (in green).

In the text

	Fig. 3. Measure of randomness using the Von Neumann estimator ℰVN (defined in Eq. (1)) as a function of the time delay and differential magnification parameters. The best-fit value obtained from a parametric model of the continuum evolution (Sect. 3.6) is shown by the white mark.
In the text

	Fig. 4. Fit of the light curve, spectral index, and spectral-index difference between the NE and SW images using the parametric model provided by Eqs. (3) and (5). The envelope of extreme model values for the differential polarization quantities is plotted against measurements using the model extension described in Sect. 3.7. The fit also provides the spectral index of the NE image (red points).
In the text

	Fig. 5. Results from our plasmon modeling (Eq. (3)) as obtained from MCMC (see also Table 4).
In the text

	Fig. 6. Overlay of all ALMA spectra of H2O (SW and NE lines of sight) and CH (SW only) observed toward PKS 1830−211 between 2012 and 2016. Special profiles showing clear deviations are highlighted with the color code given in the upper box.
In the text

	Fig. 7. Zoom onto the variations of the +174 km s−1 (top) and −225 km s−1 (bottom) features in the H2O SW and NE spectra, respectively. Same as for Fig. 6; special profiles showing clear deviations are highlighted with the color code given in the upper box.
In the text

	Fig. 8. Comparison of the average H2O spectrum (filled blue) observed in 2016 with the ALMA CH+ (1−0) spectrum (orange line) observed in 2015 (Muller et al. 2017) for both lines of sight. For the NE spectrum, the opacity of the CH+ line, assuming a filling factor of unity, is reduced by a factor 12.
In the text

	Fig. 9. Overlay of available ALMA spectra of HO and HO 11, 0 − 10, 1 toward the SW image of PKS 1830−211.
In the text

	Fig. 10. Extrema and rms dispersion of the H2O and CH absorption spectra during our 2016 monitoring (the CH profiles are not shown for the NE line of sight, because they contain only noise).
In the text

	Fig. 11. Rate of the variations between consecutive observations of the H2O absorption profile along the SW line of sight of PKS 1830−211.
In the text

	Fig. 12. Eigenspectra of the PCA decomposition of the H2O spectra toward the SW image of PKS 1830−211.
In the text

	Fig. 13. Time evolution of the eigenvalues of the PCA decomposition of H2O absorption spectra toward the SW image of PKS 1830−211.
In the text

	Fig. 14. Time evolution of the eigenvalues of the first component (PC1) of the PCA decomposition of various subsets of the monitoring spectra. The line of sight, species, considered velocity interval, and PC1 scores are given in each box.
In the text

Fig. 15. Weighted-average spectra toward the SW image of PKS 1830−211, combining all visits in 2016. The hyperfine structure is shown for each CH Λ-doublet, with the strongest hfs component set at a velocity v = 0 km s−1, or as indicated otherwise. There are two boxes for each line, the first one showing the whole absorption spectrum and the second one showing a smaller range around the continuum level. Best-fit of the CH and 13CH absorption is shown in red.

In the text

	Fig. 16. Weighted-average spectra toward the NE image of PKS 1830−211, combining all visits in 2016. The hyperfine structure is shown for each CH Λ-doublet, with the strongest hfs component set at a velocity v = −150 km s−1. The lower panel shows the CH hfs-deconvolved profile.
In the text

	Fig. 17. Measurements of the covering factor across the different visits of our 2016 monitoring. (a) Saturation level (1 − S​ν) at the center of the H2O line for all the visits in 2016. (b) The values of all spectral channels with velocities |v|< 8 km s−1 are averaged to provide measurements of the continuum-covering factor as a function of time.
In the text

Fig. 18. Correlation between the optical depths of the ortho-H2O 11, 0 − 10, 1 line and of the deconvolved profile of CH 3/2–1/2 (given for an equivalent hyperfine component with normalized intensity S​ul = 1) along the SW (left) and NE (right) lines of sight, respectively. In both lines of sight, we assume a source-covering factor of fc = 1. The colors encode the velocity intervals shown on the line profiles (top). The dashed lines correspond to opacity ratios τ(H2O)/τ(CH*) = 3 and 9, respectively.

In the text

	Fig. B.1. Overview of the continuum-normalized absorption spectra of the CH Jp = 3/2+ − 1/2− (left), CH Jp = 3/2− − 1/2+ (middle), and ortho-H2O JKa, Kc = 11, 0 − 10, 1 (right) lines toward the SW image of PKS 1830−211, for each visit of the 2016 monitoring. The visit numbers #i correspond to the entries in Table 1.
In the text

	Fig. B.2. Same as Fig. B.1, but toward the NE image of PKS 1830−211.
In the text