Gravitational Lens System PS J0147+4630 (Andromeda's Parachute): Main Lensing Galaxy and Optical Variability of the Quasar Images

Because follow-up observations of quadruple gravitational lens systems are of extraordinary importance for astrophysics and cosmology, we present single-epoch optical spectra and $r$--band light curves of PS J0147+4630. This recently discovered system mainly consists of four images ABCD of a background quasar around a foreground galaxy G that acts as a gravitational lens. First, we use long--slit spectroscopic data in the Gemini Observatory Archive and a multi-component fitting to accurately resolve the spectra of A, D, and G. The spectral profile of G resembles that of an early-type galaxy at a redshift of 0.678 $\pm$ 0.001, which is about 20\% higher than the previous estimate. Additionally, the stellar velocity dispersion is measured to $\sim$5\% precision. Second, our early $r$-band monitoring with the Liverpool Telescope leads to accurate light curves of the four quasar images. Adopting time delays predicted by the lens model, the new lens redshift, and a standard cosmology, we report the detection of microlensing variations in C and D as large as $\sim$0.1 mag on timescales of a few hundred days. We also estimate an actual delay between A and B of a few days (B is leading), which demonstrates the big potential of optical monitoring campaigns of PS J0147+4630.


Introduction
Optical light curves of the four images of a quadruple gravitationally lensed quasar (quad) have a great potential to reveal physical properties of the lensed quasar, the main lensing galaxy, and the universe as a whole. Thus, monitoring campaigns of QSO 2237+0305 (Einstein Cross; Woźniak et al. 2000;Alcalde et al. 2002;Eigenbrod et al. 2008b) have proved to be invaluable tools for unveiling the structure of the quasar accretion disk and the composition of the lensing galaxy (e.g., Shalyapin et al. 2002;Kochanek 2004;Eigenbrod et al. 2008a). Most recently, light curves of other quads also led to important astrophysical results, e.g., an estimate of the accretion disk size in HE 0435−1223 and WFI 2033−4723 using microlensinginduced variations (Fian et al. 2018;Morgan et al. 2018), and a robust measurement of the Hubble constant H 0 from the time delays of HE 0435−1223, RXJ1131−1231, and B1608+656 (assuming a flat ΛCDM cosmology with free energy density Ω Λ ; Bonvin et al. 2017). For a given quad, in addition to accurate time delays, the redshift and stellar velocity dispersion of the main lensing galaxy (and the properties of the lens environment) are required to obtain strong constraints on cosmological parameters (e.g., Suyu et al. 2017). Berghea et al. (2017) reported the serendipitous discovery of the bright quad PS J0147+4630 using multiband frames from the Panoramic Survey Telescope and Rapid Response System (Pan-STARRS; Chambers et al. 2016). The four images of this first Pan-STARRS lensed quasar are not arranged like a cross, but rather resembling the shape of an open parachute (Andromeda's Parachute). An arc-like structure contains the three brightest images (A, B, and C; r∼16−16.5 mag), while the fourth component D (r∼18 mag) is located about 3″ from such arc. The main lensing galaxy (G;r 20 mag) and D are only ∼1″ apart. Additionally, spectroscopic observations indicated the broad absorption line nature of the quasar, which has a redshift of 2.341±0.001 (Lee 2017) and is absorbed by several intervening metal systems (Rubin et al. 2018). From long-slit optical spectroscopy with the Gemini Multi-Object Spectrograph (GMOS) on the 8.1 m Gemini North Telescope, Lee (2018) also found a lens redshift z G =0.5716±0.0004. These redshifts and Hubble Space Telescope (HST) imaging of the lens system allowed Shajib et al. (2019) to predict time delays between the quasar images for a flat ΛCDM cosmology with H 0 =70 km s −1 Mpc −1 and Ω Λ =0.7 (Ω m =0.3).
The Gravitational LENses and DArk MAtter (GLENDAMA) project is conducting optical observations of a sample of ten lensed quasars with bright images (r < 20 mag) at < < z 1 3 (Gil- Merino et al. 2018). This representative sample includes the quad PS J0147+4630, which is being monitored with the 2.0 m Liverpool Telescope (LT) since 2017 August. In Section 2, we reanalyze the GNT-GMOS data of PS J0147+4630 to accurately extract spectra for the three close sources within the slit (A, D, and G). We then focus on the new spectrum of G, measuring its redshift and the stellar velocity dispersion. Implications of this data reanalysis are also discussed in Section 2. In Section 3, we present LT r-band light curves of the four images of PS J0147 +4630 spanning two complete observing seasons from 2017 to 2019, and show the potential of optical monitoring programs of this quad. We summarize the paper in Section 4. seeing conditions. The 0 5-width slit with a spatial pixel scale of 0 1614 was oriented along the line joining A and D (and crossing G). We downloaded these publicly available observations 4 to reanalyze them. Our aim was to accurately resolve the individual spectra of the three sources in the crowded region through a well-tested state-of-the-art technique (see below). Before doing the extraction of spectra for A, D, and G, usual data reductions were performed using the Gemini IRAF package. 5 Regarding the wavelength range and dispersion, they were 4403-7605 Å and 1.02 Å pix −1 , respectively. We also inferred a resolving power of ∼2500 from the 2.24 Å width of the 5577 Å [O I] line in the night airglow.
We extracted the instrumental spectra of A, D, and G by fitting three 1D Moffat profiles in the spatial direction for each wavelength bin (e.g., Sluse et al. 2007;Shalyapin & Goicoechea 2017). For an extended source, a 1D point-spread function (PSF) model does not describe its total light. Therefore, because G was treated as a point-like object, we actually did not derive the overall spectrum of the galaxy. However, the flux ratio A/G is ∼250, and the very faint galaxy does not play a relevant role when modeling the light distribution along the slit. This distribution is dominated by emission from A and D, so the critical issue is our ability to accurately account for the flux of both quasar images, minimizing the contamination of G. We used HST astrometry (Shajib et al. 2019) to set the positions of D and G with respect to A. In regard to the 1D PSF model, we checked that a skewed Moffat profile (Schönebeck et al. 2014) works better than the Gaussian or the symmetric Moffat one (see Appendix A; Moffat 1969). Thus, we adopted a skewed Moffat function with a 6 =a 7 =0, i.e., only the asymmetry parameter a 5 was considered (see Equation (10) of Schönebeck et al. 2014). We also set the Moffat power-law index to an optimal value of 2. In a first fit to the multiwavelength 1D flux distribution, the position (centroid) of A, the width and asymmetry of the Moffat function, and the amplitudes of the three components were allowed to vary at each wavelength. We then fitted A positions and Moffat structure parameters to smooth polynomial functions of the observed wavelength, fixing positionstructure parameters to their polynomial values and leaving only the three amplitudes as free parameters in a second iteration.
The four instrumental spectra of each source (one per exposure) were calibrated in flux and combined into a single spectral energy distribution. To carry out the flux calibrations, we used GNT-GMOS-B600 0 5-width slit observations of the standard star EG 131 (Bessell 1999) on 2017 September 4. Final spectra of A, D, and G are included in Table 1 and plotted in Figure 1. A comprehensive analysis of the quasar spectra will be presented in a forthcoming paper. Here, we concentrate on the optical spectrum of the main lensing galaxy. The spectral slopes of G, as well as its Ca II H&K doublet at about 6600 and 6660 Å and its G-band around 7225 Å, are in very good agreement with an early-type galaxy template 6 at z = 0.678 (see Figure 1). After doing an initial estimate for z G from the observed absorption features, we accurately measured the lens redshift using the Penalized PiXel-Fitting (pPXF) package 7 (Cappellari & Emsellem 2004;Cappellari 2017). This software compares the spectrum of G and stellar spectra, providing a correction to the initial value of z G and measuring the stellar velocity dispersion. From the MILES Library of Stellar Spectra 8 (Sánchez-Blázquez et al. 2006;Falcón-Barroso et al. 2011), we obtained z G = 0.678 ± 0.001 and σ G = 313 ± 14 km s −1 (see Figure 2).
We note that our G spectrum and the redshifted template in Figure 1 are very similar, while the G spectrum in Figure 2 of Lee (2018) does not resemble those of typical galaxies and seems to be heavily contaminated by quasar light. Thus, the previous identification of absorption lines of G is unreliable, and we adopt the new value of z G as the lens redshift. For a flat ΛCDM cosmology with H 0 =70 km s −1 Mpc −1 and Ω Λ = 0.7, an 18.5% increase in the lens redshift, from 0.572 to 0.678, results in predicted time delays that are ∼27% longer than those in Table  C1 of Shajib et al. (2019). These new delays are Dt AB = −2.7, Table 1 GNT Notes. The spectrum of G is not properly calibrated, because fluxes are underestimated in a factor ∼5. Table 1 is published in its entirety in the machine-readable format. A portion is shown here for guidance regarding its form and content. a Observed wavelength in Å. b Flux in 10 −17 erg cm −2 s −1 Å −1 . c Flux error in 10 −17 erg cm −2 s −1 Å −1 .
(This table is available in its entirety in machine-readable form.) Figure 1. GNT-GMOS-B600 spectra of PS J0147+4630ADG in 2017 September. Flux is shown in a logarithmic scale to improve visibility, and gray highlighted bars are associated with gaps between the detectors. We note that the quasar spectra are reasonably well calibrated in flux, while the galaxy spectrum is not. The redshifted (z = 0.678) template of an early-type galaxy is also displayed for comparison purposes (see main text).
Δt AC =−9, and Δt AD =245 days (Δt XY = t Y − t X ). We also remark that the value of z G plays a critical role in determining H 0 from measured time delays of the system. If we assume that

LT-IO:O Light Curves
We monitored PS J0147+4630 with the LT from 2017 August to 2018 January and from 2018 July to 2019 February, i.e., during the first two visibility periods after its discovery. All observations were made using the Sloan r-band filter on the IO: O optical camera (pixel scale of 0 30), and a single 120 s exposure was taken each observing night. Although we obtained 84 frames, six of them have poor quality and are not usable for doing photometry. Thus, in addition to basic instrumental reductions from the LT-IO:O pipeline, we used IRAF software 9 (Tody 1986(Tody , 1993 to remove cosmic-rays and bad pixels from the remaining 78 frames. The central region of one of these usable frames is shown in Figure 3. This includes the four quasar images, as well as an isolated and unsaturated star at R.A. (J2000)=26°.773246 and decl. (J2000)=+46°. 506670 (r=16.606 mag) that allows us to build up an empirical 2D PSF. A bright field star at R.A. (J2000)=26°. 746290 and decl. (J2000)=+46°. 504028 is also used as a control object having constant brightness r=15.421 mag (see the middle right side of Figure 3).
For each of the 78 frames, the brightness of A, B, C, and D, and the control star were extracted through PSF fitting, using the IMFITFITS software (McLeod et al. 1998) and the PSF star located in the center of the field of view in Figure 3. In the strong lensing region, ABCD and G were modeled as four point-like sources and a de Vaucouleurs profile convolved with the empirical 2D PSF, respectively. Our realistic overall model also incorporated several HST constraints: positions of B, C, D, and G with respect to A, and structure parameters of G (Shajib et al. 2019). We applied the IMFITFITS code to the best frames to obtain the constant galaxy flux, and then to all frames (whatever their quality), fixing the galaxy flux and allowing the remaining parameters to vary, i.e., the absolute position of A and the four quasar fluxes (see Appendix B). In Table 2, we present values of the FWHM of the seeing disk and the signalto-noise ratio (S/N) for the PSF star, along with the r-band magnitudes of the quasar images and the control star.
To produce high-quality early light curves of PS J0147 +4630, we selected the 70 observing epochs in Table 2 with FWHM1 6, and removed two additional epochs in which the B magnitude deviates significantly from adjacent values (2018 August 16 and 28). In Figure 4, we display our final light curves of the quasar images (filled symbols). From the data at these 68 epochs, we also estimated typical magnitude errors. For each image, relying on theoretical grounds, we calculated the standard deviation between magnitudes having time separations 4.5 days (i.e., using 31 pairs of consecutive magnitudes), and then divided it by the square root of 2. The resulting errors are 0.0053 (A), 0.0058 (B), 0.0093 (C), and 0.0154 (D) mag. These typical errors were multiplied by the relative S/N at each epoch, ( ) á ñ S N S N epoch , to obtain individual photometric uncertainties (á ñ S N is the average S/N; Howell 2006). We achieve ∼0.5%-1.5% photometry over a period of about 1.6 yr, which includes an unavoidable visibility gap of more than 5 months. The effective sampling rate (excluding the long gap) is ∼5 data per month.  Notes. Table 2 is published in its entirety in the machine-readable format. A portion is shown here for guidance regarding its form and content. a yymmdd. b MJD-50000. c FWHM of the seeing disk in ″. d S/N for the PSF star within a circle of radius FWHM. e r-SDSS magnitude.
(This table is available in its entirety in machine-readable form.)

What Can We Learn from Early Light Curves of the Lensed Quasar?
The expected time delays in Section 2 and the early light curves of ABCD in Section 3.1 are useful tools to analyze the origin and properties of the quasar variability in the r band.  Figure 5. Second, we computed difference light curves in a standard way, i.e., for Y=B, C, and D. For a given image Y, the dates in the magnitude-and time-shifted light curve of A, m AY (t), were taken as reference epochs. Values of m AY (t) were subtracted from averaged magnitudes of Y in bins with semiwidth α centered on the reference dates. We probed several reasonable values of α, obtaining difference curves consistent with each other, and then setting α=5 days (see the right panels of Figure 5). The AB comparison (top panels) indicates that both images exhibit almost parallel variations, so the difference light curve has a noisy behavior around the zero level. Despite how extrinsic (microlensing) variability can be present, it should consist of short timescale fluctuations with amplitudes 0.02 mag. However, the situation is quite different from the other two AC (middle panels) and AD (bottom panels) comparisons. In both cases, we detect significant microlensing variations. The difference light curves for C and D include ∼0.1 mag fluctuations over timescales from 100 to 400 days. Some other quads show similar levels of microlensing activity (e.g., Fian et al. 2018;Gil-Merino et al. 2018). We remark that microlensing signals in the D image rely on time delays Δt AD predicted by Shajib et al.ʼs (2019) mass model and plausible values of H 0 . Despite the true value of Dt AD maybe being slightly out of the delay range used in this paper, this would not produce noticeable changes with respect to the similar signals in the bottom right panel of Figure 5.
The LT-IO:O light curves of A and B over the two first monitoring seasons have similar shapes, providing evidence for intrinsic variability and enabling us to estimate Dt AB . As we only try to show the potential of monitoring campaigns of PS J0147+4630, instead of an exhaustive analysis from several cross-correlation techniques, we exclusively used the D 4,2 2 dispersion estimator (Pelt et al. 1996) to match both light curves. This method is simple and very popular, and it works in   Table 2) are marked with filled and open symbols, while filled symbols trace the final light curves after removing 10 epochs (error bars of A, B, and C have sizes similar to those of symbols we use; see main text). the presence of short timescale microlensing (e.g., Goicoechea & Shalyapin 2016). A comprehensive study of all time delays will be done when much more extended light curves of the quad are available. First, in order to obtain Δt AB and a constant magnitude offset Δm AB =m B (t)−m A (t + Δt AB ), we focused on a biparametric D 4,2 2 . Second, we checked that reasonable values of the decorrelation length (d  10 days) produce smooth dispersion spectra having minima at D < t AB 0, and then chose δ=10 days. Third, we generated 3000 simulated light curves of each quasar image at epochs equal to those of observation. Each observed magnitude was modified by adding normally distributed random numbers with zero mean and standard deviation equal to the measured uncertainty. Fourth, the D 4,2 2 estimator was applied to all (A, B) pairs of simulated curves. The distribution of magnitude offsets led to Δm AB =0.248±0.001 mag (1σ confidence interval), and the time delay histogram is depicted in Figure 6. This yields an 1σ interval Δt AB =−2.6 -+ 3.2 1.1 days, indicating that B is leading. Additionally, the median delay is in very good agreement with the expected value for H 0 =70 km s −1 Mpc −1 (vertical solid line in Figure 6; see the end of Section 2). Figure 6 also displays a vertical "zipper" indicating the narrow range of expected delays for H 0 =65−75 km s −1 Mpc −1 .

Summary
Within the framework of the GLENDAMA project, we are analyzing archive data and conducting follow-up observations of a sample of ten gravitationally lensed quasars (Gil-Merino et al. 2018), and this paper focuses on the recently discovered Pan-STARRS quad PS J0147+4630 (Berghea et al. 2017). After retrieving long-slit spectroscopic data in the Gemini Observatory Archive (Lee 2018), we use a robust multicomponent fitting (e.g., Sluse et al. 2007;Shalyapin & Goicoechea 2017) to extract individual spectra of the main lensing galaxy and two quasar images in the observed crowded region. Despite both quasar spectra containing valuable imprints of intervening objects, we only study in detail the optical spectrum of G. This agrees very well with the spectral profile of an early-type galaxy at a redshift of 0.678, so the previous redshift determination through a heavily contaminated spectral energy distribution is very likely to be biased. An accurate estimate of z G is crucial to predict time delays and carry out cosmological studies (e.g., Jackson 2015). We also measure a stellar velocity dispersion σ G =313±14 km s −1 , which is relevant, along with other results from ongoing observational efforts (e.g., macrolens flux ratios could be revealed by radio imaging with the Very Large Array; Berghea et al. 2017), to better constrain the lensing mass distribution.
We have also conducted an early r-band monitoring of the four quasar images ABCD with the LT, and the corresponding light curves are used to probe the potential of optical monitoring campaigns. We take the new redshift z G =0.678 to properly modify the time delays predicted by Shajib et al. (2019), which permits us to construct reliable combined and difference light curves. These curves indicate that C and D images are affected by microlensing effects presumably produced by stars within G, and the microlensing-induced variations are promising tools to constrain the accretion disk size in PS J0147+4630 (e.g., Fian et al. 2018;Morgan et al. 2018). Lee (2018) also pointed out the presence of microlensing when comparing his single-epoch Gemini spectra of A and D. In addition, we find that A and B vary in an almost parallel way, suggesting that we basically detect intrinsic activity of the source quasar in both images. Using a single magnitude offset and a time delay to cross-correlate the light curves of A and B, we obtain Δt AB =−2.6 -+ 3.2 1.1 days. Hence, as expected from lens modeling, A arrives a few days later than B. Although the Δt AB estimation might be improved by considering two magnitude offsets or other sophisticated approaches, more extended brightness records are required to disentangle intrinsic from extrinsic (microlensing) signals in ABCD, and thus, to accurately measure the three independent delays Δt AB , Δt AC , and Δt AD . These delays are essential pieces for a lens-based cosmology (e.g., Bonvin et al. 2017).
We thank the anonymous referee for helpful comments and suggestions to improve the presentation of this paper. The Liverpool Telescope 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 Astrofisica de Canarias with financial support from the UK Science and Technology Facilities Council. This article is also based on observations obtained at the Gemini Observatory (acquired through the Gemini Observatory Archive and processed using the Gemini IRAF package), which is operated by the Association of Universities for Research in Astronomy, Inc., under a cooperative agreement with the NSF on behalf of the Gemini partnership: the National Science Foundation Software: IRAF (Tody 1986(Tody , 1993, pPXF (Cappellari & Emsellem 2004;Cappellari 2017), IMFITFITS (McLeod et al. 1998).

Appendix A Extracting Individual Spectra of A, D, and G
To illustrate how the multi-component fitting works on instrumental flux distributions along the slit, we considered spatial flux distributions at wavelengths around 6800 Å for the first 1200 s spectroscopic exposure (see Section 2). By averaging over an 100 Å wavelength interval, we then obtained the spatial profile that appears in Figure 7 (black filled circles). This 1D flux distribution was fitted to three skewed Moffat profiles with the same structure parameters (contributions of the three close sources A, D, and G; solid lines with different colors), resulting in a very good global solution (see A+D+G in Figure 7). We also show a fit to three Gaussian functions (dashed lines) for comparison purposes. These symmetric profiles do not reproduce the data in an accurate way, leading to significant fit residuals.

Appendix B Extracting Individual r-band Fluxes of A, B, C, and D
As an example of the photometric fitting procedure, we first selected the r-band frame on 2018 August 13 (see Figure 3). The subframe of interest (containing the lens system) is shown in the left panel of Figure 8. Second, the 2D flux distribution in this subframe was modeled as four empirical PSFs (quasar images ABCD) plus a de Vaucouleurs profile convolved with the empirical PSF (main lensing galaxy G). The best-fit model from the IMFITFITS software package (McLeod et al. 1998) has a reduced χ 2 of 1.28, which is a typical value when fitting system subframes at other epochs. The middle panel of Figure 8 displays the residual signal after subtracting the four quasar images. Regarding this panel, it is worth noting that the fluxes of ABCD are one or two orders of magnitude larger than the flux of G, and the galaxy is an extended source. Additionally, residuals correspond to a relatively short exposure with a 2 m telescope. Hence, the light distribution in the middle panel is not dominated by emission from G. In the right panel of Figure 8, we also show the residual light after subtracting the full best-fit model (all sources). Pixels in this last subframe only contain the expected noisy signal.
It is worthy to mention that we analyze a crowded system using a well-tested photometric method. This worked quite well in other gravitational lens systems, since light curves from the IMFITFITS method basically agreed with those from alternative techniques, even using different telescopes Ullán et al. 2006;Hainline et al. 2013;Gil-Merino et al. 2018). We also note that cross-talk between two close sources depends on seeing, and sometimes it produces evident deviations (with respect to adjacent magnitudes) in the light curves of both sources (e.g., see discussion on systematic effects and the corresponding outliers in Section 4 of Goicoechea & Shalyapin 2016). However, after removing data from frames with relatively poor seeing and two additional epochs producing evident outliers in the brightness record of the B image (see Section 3.1), our final light curves of PS J0147+4630 are not affected by significant systematic effects.