Andromeda's Parachute: Time Delays and Hubble Constant

Optical frames from the Panoramic Survey Telescope and Rapid Response System (Pan-STARRS; Chambers et al. 2019) led to the serendipitous discovery of the strong gravitational lens system with a quadruply imaged quasar (quad) PS J0147+4630 (Berghea et al. 2017). Due to its position in the sky and the spatial arrangement of the four quasar images, this quad is also called Andromeda's Parachute (e.g., Rubin et al. 2018). The three brightest images (A, B, and C) form an arc that is about 3'' from the faintest image D, and the main lens galaxy G is located between the bright arc and D. This configuration is clearly seen in the left panel of Figure 1, which is based on Hubble Space Telescope (HST) data.

Zoom In Zoom Out Reset image size

Early optical spectra of the system confirmed the gravitational lensing phenomenon and revealed the broad absorption-line nature of the quasar, which has a redshift z_s ∼ 2.36 (Lee 2017; Rubin et al. 2018). Lee (2018) also performed the first attempt to determine the redshift of G from spectroscopic observations with the 8.1 m Gemini North Telescope (GNT). An accurate reanalysis of these GNT data showed that the first estimate of the lens redshift was biased, by enabling better identification of G as an early-type galaxy at z_l = 0.678 ± 0.001 with stellar velocity dispersion σ_l = 313 ± 14 km s⁻¹ (Goicoechea & Shalyapin 2019), in good agreement with the recent measurements of z_l and σ_l by Mozumdar et al. (2023).

As far as we know, the quasar PS J0147+4630 is the brightest source in the sky at redshifts z > 1.4 (apart from transient events such as gamma-ray bursts), and its four optical images can be easily resolved with a ground-based telescope in normal seeing conditions. Thus, it is a compelling target for various physical studies based on high-resolution spectroscopy (e.g., Rubin et al. 2018) and detailed photometric monitoring (e.g., Lee 2018). Early two-season monitoring campaigns with the 2.0 m Liverpool Telescope (LT; Goicoechea & Shalyapin 2019) and the 2.5 m Nordic Optical Telescope (NOT; Dyrland 2019) provided accurate optical light curves of all quasar images, as well as preliminary time delays and evidence of microlensing-induced variations. A deeper look at the optical variability of Andromeda's Parachute is of great importance, since robust time delays and well-observed microlensing variations can be used to determine cosmological parameters (e.g., Treu & Marshall 2016) and the structure of the quasar accretion disk (e.g., Schmidt & Wambsganss 2010).

This paper is organized as follows. In Section 2, we present combined LT and NOT light curves of the four images of PS J0147+4630 spanning six observing seasons from 2017 to 2022. In Section 3, using these optical light curves, we carefully analyze the time delays between images and the quasar microlensing variability. In Section 4, we discuss the Hubble constant (H₀) value from the measured time delays and lens mass models. Our main conclusions are included in Section 5.

We monitored PS J0147+4630 with the LT from 2017 August to 2022 October using the IO:O optical camera with a pixel scale of ∼030. Each observing night, a single 120 s exposure was taken in the Sloan r-band filter, and over the full monitoring period, 212 r-band frames were obtained. The LT data reduction pipeline carried out three basic tasks: bias subtraction, overscan trimming, and flat-fielding. Additionally, the IRAF software ⁶ (Tody 1986, 1993) allowed us to remove cosmic rays and bad pixels from all frames. We extracted the brightness of the four quasar images ABCD through point-spread function (PSF) fitting, using the IMFITFITS software (McLeod et al. 1998) and following the scheme described by Goicoechea & Shalyapin (2019). Table 1 includes the position and magnitudes of the PSF star, as well as of other relevant field stars. These data are taken from the Data Release 1 of Pan-STARRS ⁷ (Flewelling et al. 2020). Our photometric model consisted of four point-like sources (ABCD) and a de Vaucouleurs profile convolved with the empirical PSF (lensing galaxy G). Positions of components with respect to A and structure parameters of G were constrained from HST data (Shajib et al. 2019).

We also selected six nonvariable blue stars in the field of PS J0147+4630 and performed PSF photometry on five of them (see the calibration stars Cal1–Cal5 in Table 1; Cal6 is a saturated star in LT frames). For each of the five calibration stars, we calculated its average magnitude within the monitoring period and magnitude deviations in individual frames (by subtracting the average). In each individual frame, the five stellar magnitude deviations were averaged together to calculate a single magnitude offset, which was then subtracted from the magnitudes of quasar images. After this photometric calibration, we removed 22 observing epochs in which quasar magnitudes deviate appreciably from adjacent values. Thus, the final LT r-band light curves are based on 190 frames (epochs), and the typical uncertainties in the light curves of the quasar images and control star (see Table 1) were estimated from magnitude differences between adjacent epochs separated by no more than 4.5 days (Goicoechea & Shalyapin 2019). We derived typical errors of 0.0062 (A), 0.0077 (B), 0.0097 (C), 0.0197 (D), and 0.0058 (control star) mag. For the control star, we have also verified that its typical error practically coincides with the standard deviation of all measures (0.0055 mag). To obtain photometric uncertainties at each observing epoch, the typical errors were scaled by the relative signal-to-noise ratio of the PSF star (Howell 2006).

The optical monitoring of PS J0147+4630 with the NOT spanned from 2017 August to 2019 December. We used the ALFOSC camera with a pixel scale of ∼021 and the R-Bessel filter. This passband is slightly redder than the Sloan r band. Each observing night, we mainly took three exposures of 30 s each under good seeing conditions. The FWHM of the seeing disk was about 10 (we also estimated FWHM seeing = 135 ± 015 from LT frames), and we collected 298 individual frames over the entire monitoring campaign. After a standard data reduction, IMFITFITS PSF photometry yielded magnitudes for the quasar images (see above for details on the photometric model). To avoid biases in the combined LT-NOT light curves, the same photometric method was applied to LT and NOT frames. This method differs from that of Dyrland (2019), who used the DAOPHOT package in IRAF (Stetson 1987; Massey & Davis 1992) to extract magnitudes from NOT frames.

The six calibration stars in Table 1 were used to adequately correct quasar magnitudes (see above), and we were forced to remove 17 individual frames leading to magnitude outliers. We then combined R-band magnitudes measured on the same night to obtain photometric data of the lensed quasar and control star at 77 epochs. Again, typical errors were derived from magnitudes at adjacent epochs that are separated <4.5 days. This procedure led to uncertainties of 0.0122 (A), 0.0122 (B), 0.0144 (C), 0.0197 (D), and 0.0170 (control star) mag. Errors at each observing epoch were calculated in the same way as for the LT light curves.

As a last step, we combined the r-band LT and R-band NOT light curves. If we focus on the quasar images and consider rR pairs separated by no more than 2.5 days, the values of the average color 〈r − R〉 are 0.0565 (A), 0.0616 (B), 0.0546 (C), and 0.0652 (D). Brightness records of the ABC images are more accurate than those of D, and thus we reasonably take the average colors of ABC to estimate a mean r − R offset of 0.0576 mag. After correcting the R-band curves of the quasar for this offset, we obtain the new records in Table 2. Table 2 contains r-band magnitudes of the quasar images and the control star at 267 observing epochs (MJD−50,000). In Figure 1, we also display our new 5.2 yr light curves.

Previous efforts focused on early monitorings with a single telescope, trying to estimate delays between the image A and the other quasar images, Δt_AX = t_A − t_X (X = B, C, D), and find microlensing signals (Dyrland 2019; Goicoechea & Shalyapin 2019). ⁸ Here, we use the new light curves in Section 2 along with state-of-the-art curve-shifting algorithms to try to robustly measure time delays between images. At the end of this section, we also discuss the extrinsic (microlensing) variability of the quasar.

As is clear from Figure 1, there are short time delays between images ABC, while it is hard to get an idea about the Δt_AD value by eye. Fortunately, there are several cross-correlation techniques to measure time delays between light curves containing microlensing variations (e.g., Liao et al. 2015, and references therein), and thus we considered PyCS3 curve-shifting algorithms ⁹ (Tewes et al. 2013; Millon et al. 2020a, 2020b) to obtain reliable time delays of PS J0147+4630. PyCS3 is a well-tested software toolbox to estimate time delays between images of gravitationally lensed quasars, and we focused on the χ² technique, assuming that the intrinsic signal and the extrinsic ones can be modeled as a free-knot spline (FKS). This technique shifts the four light curves simultaneously (ABCD comparison) to better constrain the intrinsic variability, and relies on an iterative nonlinear procedure to fit the four time shifts and splines that minimize the χ² between the data and model (Tewes et al. 2013). Results depend on the initial guesses for the time shifts, so it is necessary to estimate the intrinsic variance of the method using a few hundred initial shifts randomly distributed within reasonable time intervals. In addition, an FKS is characterized by a knot step, which represents the initial spacing between knots. The model consists of an intrinsic spline with a knot step η and four independent extrinsic splines with η_ml that account for the microlensing variations in each quasar image (Millon et al. 2020a).

To address the intrinsic variability, we considered three η values of 30, 50, and 70 days. Intrinsic knot steps shorter than 30 days fit the observational noise, whereas η values longer than 70 days do not fit the most rapid variations of the source quasar. Intrinsic variations are usually faster than extrinsic ones, and additionally, the software works fine when the microlensing knot step is significantly longer than η. Therefore, the microlensing signals were modeled as FKSs with η_ml = 350−400 days (i.e., values intermediate between those shown in Table 2 of Millon et al. 2020a). We also generated 500 synthetic (mock) light curves of each quasar image, optimized every mock ABCD data set, and checked the similarity between residuals from the fits to the observed curves and residuals from the fits to mock curves. The comparison of residuals was made by means of two statistics: standard deviation and normalized number of runs Z_r (see details in Tewes et al. 2013). For η = 50 days and η_ml = 400 days, histograms of residuals derived from mock curves (gray) and from the LT-NOT light curves of PS J0147+4630 are included in the top panels of Figure 2. It is apparent that the standard deviations through the synthetic and the observed curves match very well. Additionally, the bottom panels of Figure 2 show distributions of Z_r from synthetic light curves (gray) for η = 50 days and η_ml = 400 days. These bottom panels also display the Z_r values from the observations (vertical lines), which are typically located at ∼0.4σ of the mean values of the synthetic distributions.

Zoom In Zoom Out Reset image size

Four pairs of (η, η_ml) values (see above) led to the set of time delays in Figure 3. We have verified that other feasible choices for η_ml (e.g., η_ml = 200 days) do not substantially modify the results in this figure. The black horizontal bars correspond to 1σ confidence intervals after a marginalization over results for all pairs of knot steps, and those in the left panels and bottom panels of Figure 3 are included in Table 3. We finally adopted the time delays in Table 3, which are symmetric about central values and useful for subsequent studies.

Zoom In Zoom Out Reset image size

It seems to be difficult to accurately determine delays between the brightest images ABC because they are really short. To robustly measure Δt_AC in a near future, we will most likely need to follow a nonstandard strategy focused on several time segments associated with strong intrinsic variations and weak extrinsic signals. Fortunately, we find an accurate and reliable value of Δt_AD (uncertainty of about 4%), confirming the early result from two monitoring seasons with the NOT and a technique different to that we used in this paper (Dyrland 2019). It is also worth mentioning that the dispersion method ignoring microlensing variations (the simplest approximation with fewer free parameters; Pelt et al. 1996) produces an optimal AD delay separated by only 10 days from that obtained with PyCS3. We also note that Δt_BD and Δt_CD have errors of 3.5%–3.7%, and thus we present accurate values of the three independent time delays relative to image D.

An image comparison spanning 13 yr is also depicted in Figure 4. We have downloaded five r-band warp frames of PS J0147+4630 that are included in the Data Release 2 of Pan-STARRS. These Pan-STARRS frames were obtained on three nights in the 2010−2013 period, i.e., a few years before the discovery of the lens system. Two frames are available on two of the three nights, so rough photometric uncertainties through average intranight variations are 0.012 (A), 0.008 (B), 0.019 (C), and 0.033 (D) mag. To discuss the differential microlensing variability of the images BCD with respect to A, Figure 4 shows the original curve of A along with shifted curves of BCD. We used the central values of the delays relative to image A and constant magnitude offsets to shift curves. The offsets Δm_AB, Δm_AC, and Δm_AD are average differences between magnitudes of A and those of B, C, and D, respectively. The global shapes of the four brightness records indicate the presence of long-term microlensing effects and suggest that PS J0147+4630 is a suitable target for a deeper analysis of its microlensing signals. Over the last six years, it is noteworthy that there is good overlap between the original curve of A and the shifted curve of D. In addition, the differential microlensing variation of C is particularly prominent, showing a microlensing episode with a total amplitude greater than 0.1 mag.

Zoom In Zoom Out Reset image size

Berghea et al. (2017) presented preliminary modeling of the lens mass of PS J0147+4630 from Pan-STARRS data, whereas Shajib et al. (2019) have recently modeled the lens system using HST imaging. To model the HST images, Shajib et al. (2019) have considered the distributions of light of the lens and source, and the lens mass distribution. Their solution for the lensing mass relies on a lens scenario consisting of a singular power-law ellipsoid (SPLE; describing the gravitational field of the main lens galaxy G) and an external shear (ES; accounting for the gravitational action of other galaxies). The dimensionless surface mass density (convergence) profile of the SPLE was characterized by a power-law index β = 2.00 ± 0.05, where β = 2 for an isothermal ellipsoid. ¹⁰

We first considered Shajib et al.'s (2019) solution, a flat ΛCDM (standard) cosmology with matter and dark energy densities of Ω_M = 0.3 and Ω_Λ = 0.7, respectively, ¹¹ updated redshifts z_l = 0.678 (Goicoechea & Shalyapin 2019) and z_s =2.357 (based on emission lines that are observed at near-IR wavelengths), and the time delay in the third column of Table 3 to calculate ${H}_{0}^{\mathrm{model}}$ and put it into perspective. We obtained ${H}_{0}^{\mathrm{model}}$ = 100 ± 10 km s⁻¹ Mpc⁻¹, which significantly exceeds a concordance value of ∼70 km s⁻¹ Mpc⁻¹ (e.g., Jackson 2015). If additional mass along the line of sight is modeled as an external convergence κ_ext, then ${H}_{0}^{\mathrm{true}}={H}_{0}^{\mathrm{model}}(1-{\kappa }_{\mathrm{ext}})$ (e.g., Rusu et al. 2017). The factor 1 − κ_ext should be ∼0.7 (κ_ext ∼ 0.3) to decrease H₀ until accepted values. Therefore, the external convergence required to solve the H₀ crisis is one order of magnitude higher than typical values of κ_ext (e.g., Rusu et al. 2017; Birrer et al. 2020).

The Hubble constant can be also inferred from another lens mass solution based on approaches similar to those of Shajib et al. (2019). Adopting a standard cosmology and updated redshifts (see above), the solution of Schmidt et al. (2023) and the three time delays relative to image D (last three columns in Table 3) led to ${H}_{0}^{\mathrm{model}}$ values in the range 116–131 km s⁻¹ Mpc⁻¹. Thus, Schmidt et al.'s (2023) solution with power-law index β = 2.08 ± 0.02 produces even higher ${H}_{0}^{\mathrm{model}}$ values than those from Shajib et al.'s (2019) solution. Although the H₀ crisis may be related to an inappropriate (SPLE + ES) lens scenario or a very high external convergence, we have sought for a new mass reconstruction using astrometric and time-delay constraints, an SPLE + ES scenario, updated redshifts, a standard cosmology, and ${H}_{0}^{\mathrm{model}}$ = 70 km s⁻¹ Mpc⁻¹. In presence of a typical (weak) external convergence, the ${H}_{0}^{\mathrm{true}}$ value would be consistent with accepted ones.

Our standard astrometric constraints consisted of the HST positions of ABCD (with respect to G at the origin of coordinates; Shajib et al. 2019). SPLE + ES mass models of quads usually indicate the existence of an offset between the center of the SPLE and the light centroid of the galaxy (e.g., Sluse et al. 2012; Shajib et al. 2019). Hence, instead of formal astrometric errors for G, we adopted σ_x = σ_y = 004. This uncertainty level equals the rms of mass/light positional offsets for most quads in the sample of Shajib et al. (2019). In addition to astrometric data, the set of constraints incorporated the LT-NOT time delays relative to image D (see Table 3). The number of observational constraints and the number of model parameters were 13 and 10, respectively. For 3 degrees of freedom, the GRAVLENS/LENSMODEL software ¹² (Keeton 2001, 2004) led to the 1σ intervals in Table 4 (χ² = 3.56 for the best fit).

While our solution for the mass of the early-type galaxy G is characterized by a convergence a little shallower than isothermal (β < 2; see Table 4), Shajib et al.'s (2019) and Schmidt et al.'s (2023) solutions for the surface mass density are more centrally concentrated (β ≥ 2), suggesting this is a key reason to infer such high values of ${H}_{0}^{\mathrm{model}}$ from previous models (e.g., Refsdal & Surdej 1994; Kochanek & Schechter 2004; Jackson 2015). The only issue with all SPLE + ES mass models is the existence of a significant mass/light misalignment, i.e., the light and mass distributions of the lens galaxy do not match. This misalignment could be genuine or due to an oversimplification of the lens scenario (e.g., Sluse et al. 2012; Shu et al. 2016; Gomer & Williams 2021). Most early-type galaxies reside in overdense regions, so external tidal fields in their vicinity are expected to have relatively high amplitudes. External shear strengths for quads exceeding 0.1 are consistent with N-body simulations and semianalytic models of galaxy formation (Holder & Schechter 2003). Using a model consisting of a singular isothermal elliptical potential and external shear, Luhtaru et al. (2021) have also shown that PS J0147+4630 is a shear-dominated system.

In this paper, we performed a comprehensive analysis of the optical variability of the quadruply imaged quasar PS J0147+4630. Well-sampled light curves from its discovery in 2017 to 2022 were used to robustly measure the three time delays relative to image D. However, these light curves did not allow us to accurately (in terms of fractional uncertainty) determine the very short time delays between the three bright images ABC forming a compact arc. Additionally, the microlensing-induced variation of the C image (with respect to A) was particularly large in the period 2017−2022. Combining our new brightness records with quasar fluxes from Pan-STARRS imaging in 2010−2013, the extended light curves also revealed significant long-term microlensing effects. A microlensing analysis of current data and future light curves from a planned optical multiband monitoring is expected to lead to important constraints on the spatial structure of the quasar accretion disk (Eigenbrod et al. 2008; Poindexter et al. 2008; Cornachione et al. 2020; Goicoechea et al. 2020).

From HST imaging of the quad, Shajib et al. (2019) and Schmidt et al. (2023) have carried out reconstruction of the lensing mass from an SPLE + ES scenario. However, using updated redshifts of the source and lens (and assuming a standard cosmology), their mass reconstructions along with measured delays relative to image D led to an unacceptably large value of the Hubble constant. Although the integrated mass from objects along the line of sight to PS J0147+4630 is still unknown, an unexpected (unusually high) external convergence is required to fix this H₀ issue. To try to overcome the H₀ crisis, we have sought and found a new mass model that is consistent with astrometric and time-delay constraints, a typical external convergence, and currently accepted values for H₀ around 70 km s⁻¹ Mpc⁻¹ (e.g., see Figure 2 of Di Valentino et al. 2021). Time delays are very sensitive to the slope of the mass profile of the main lens galaxy G (e.g., Kochanek & Schechter 2004), and the new model incorporates a surface mass density less centrally concentrated than previous ones.

Alternatively, the SPLE + ES lens scenario might be an oversimplification of the actual one, since all SPLE + ES models indicate that there is a mass/light misalignment. While this misalignment may be true, it could also be due to the presence of nonmodeled components such as substructures and/or companions of G (e.g., Sluse et al. 2012; Gomer & Williams 2021). Further refinement of the lens scenario along with an extension and improvement of the set of observational constraints (future deep photometry and spectroscopy is a pending task of special relevance) will contribute to an accurate determination of H₀ and other cosmological parameters (e.g., Bonvin et al. 2017; Birrer et al. 2020). The forthcoming Legacy Survey of Space and Time at the Vera C. Rubin Observatory should provide the strong lens community with a strong increase in the number of known lensed quasars with measured time delays. To be able to utilize such a large increase in the statistical sample to provide correspondingly precise and accurate measurements of H₀, it is crucial to reliably identify the systems with more complex lens scenarios that could otherwise bias the H₀ measurement. PS J0147+4630 provides an interesting case study in this respect.

We thank Martin Millon for making publicly available a Jupiter notebook that has greatly facilitated the use of the PyCS3 software. We also thank anonymous comments and suggestions to a preliminary version of this manuscript, which have helped us to build the current version. This paper is based on observations made with the Liverpool Telescope (LT) and the Nordic Optical Telescope (NOT). The LT is operated on the island of La Palma by Liverpool John Moores University in the Spanish Observatorio del Roque de los Muchachos of the Instituto de Astrofísica de Canarias with financial support from the UK Science and Technology Facilities Council. The NOT is operated by the Nordic Optical Telescope Scientific Association at the Observatorio del Roque de los Muchachos, La Palma, Spain, of the Instituto de Astrofísica de Canarias. The data presented here were in part obtained with ALFOSC, which is provided by the Instituto de Astrofísica de Andalucia (IAA) under a joint agreement with the University of Copenhagen and NOTSA. We thank the staff of both telescopes for a kind interaction. We have also used imaging data taken from the Pan-STARRS archive and the Barbara A. Mikulski archive for the NASA/ESA Hubble Space Telescope, and we are grateful to all institutions developing and funding such public databases. V.N.S. would like to thank the Universidad de Cantabria (UC) and the Spanish AEI for financial support for a long stay at the UC in the period 2022–2023. H.D. acknowledges support from the Research Council of Norway. This research has been supported by the grant PID2020-118990GB-I00 funded by MCIN/AEI/10.13039/501100011033.

Facilities: Liverpool:2m (IO:O) - Liverpool 2 meter telesope at Observatorio del Roque de los Muchachos, NOT (ALFOSC) - Nordic Optical Telescope, PS1 - Panoramic Survey Telescope and Rapid Response System Telescope #1 (Pan-STARRS), HST (WFC3) - Hubble Space Telescope satellite.

Software: IRAF (Tody 1986, 1993), IMFITFITS (McLeod et al. 1998), Python (https://www.python.org/), PyCS3 (https://gitlab.com/cosmograil/PyCS3), GRAVLENS/LENSMODEL (http://www.physics.rutgers.edu/keeton/gravlens/).