{'fraction_non_alphanumeric': 0.03912263210368893, 'fraction_numerical': 0.027437686939182453, 'mean_word_length': 4.342138900724329, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Discovery of the First Quadruple Gravitationally Lensed Quasar Candidate with Pan-STARRS

We report the serendipitous discovery of the first gravitationally lensed quasar candidate from Pan-STARRS. The grizy images reveal four point-like images with magnitudes between 14.9 and 18.1 mag. The colors of the point sources are similar, and they are more consistent with quasars than with stars or galaxies. The lensing galaxy is detected in the izy bands, with an inferred photometric redshift of ∼0.6, lower than that of the point sources. We successfully model the system with a singular isothermal ellipsoid with shear, using the relative positions of the five objects as constraints. While the brightness ranking of the point sources is consistent with that of the model, we find discrepancies between the model-predicted and observed fluxes, likely due to microlensing by stars and millilensing due to the dark matter substructure. In order to fully confirm the gravitational lens nature of this system and add it to the small but growing number of the powerful probes of cosmology and astrophysics represented by quadruply lensed quasars, we require further spectroscopy and high-resolution imaging.


INTRODUCTION
Since the serendipitous discovery of the first gravitationally lensed quasar (Walsh, Carswell, & Weymann 1979), these systems have become powerful probes of astrophysics and cosmology, as illustrated by the wealth of science they have provided over the years. For example, Peng et al. (2006) and Ding et al. (2017) have studied the coevolution of supermassive black holes and the quasar hosts harboring them up to high redshift (Claeskens & Surdej 2002;Treu 2010). Oguri et al. (2012), Bonvin et al. (2017) and others have constrained the cosmological constant and the Hubble constant from samples of lensed quasars (see the recent review by Treu & Marshall 2016). Flux ratio anomalies have revealed luminous satellites or set constraints on the dark matter substructure in the lensing galaxies (e.g., Chiba et al. 2005;McKean et al. 2007;Fadely & Keeton 2012).
These applications all require increasingly large samples of lensed quasars, and in particular, quadruply lensed quasars (quads), due to the increased number of modeling constraints they provide. To date, there are only about three dozen known quads over the whole sky 1 . Due to the rare nature of lensed quasars, large-scale surveys of sufficient depth and resolution are required to significantly increase the present sample (Oguri & Marshall 2010). As such surveys have become available, pioneering searches for lensed quasars, such as the Cosmic Lens All Sky Survey (CLASS; Myers et al. 2003) have been succeeded by the SDSS Quasar Lens Search (SQLS; Oguri et al. 2006), and more recently by the STRong lensing Insight in the Dark Energy Survey (STRIDES; Agnello et al. 2015) and the extension to SQLS using the Baryon Oscillation Spectroscopic Survey (BQLS; More et al. 2016), amongst others.
In this paper, we report the serendipitous discovery of a quad lensed quasar candidate from the Panoramic Survey Telescope and Rapid Response System (Pan-STARRS1, hereafter PS1) released images. The lensed quasar was discovered as part of a variability study of active galactic nuclei selected based on Mid-Infrared characteristics (Secrest et al. 2015), and data from the USNO Robotic Astrometric Telescope (URAT; Zacharias et al. 2015). To our knowledge, if proven with spectroscopic observations, this would be the first published gravitational lens discovered in the PS1 data. We note that this system is very similar to the well-studied quad lens RX J1131-1231 (Sluse et al. 2003), in terms of image separation and overall configuration. We model the system to provide evidence in support of its lensing nature, and we conclude with the necessity of obtaining spectroscopic confirmation. Due to its position, from the ground this is not possible until the end of the year.
PS1 is a wide-field imaging system, with a 1.8 m telescope and 7.7 deg 2 field of view (FOV), located on the summit of Haleakala in the Hawaiian island of Maui. The first PS1 data was released in 2016 December, including both images and catalogs (see Chambers et al. 2016). The 1.4 Gpixel camera consists of 60 CCDs with pixel a size of 0.256 arcsec (Onaka & al. 2008;Tonry & Onaka 2008). It uses five filters (g P 1 , r P 1 , i P 1 , z P 1 , y P 1 , hereafter grizy), similar to the ones used by the Sloan Digital Sky Survey (SDSS; York et al. 2000). The largest survey PS1 performs is the 3π survey, covering the entire sky north of −30 deg decl. Given the large sky coverage, resolution, and depth, PS1 is expected to contain nearly 2000 gravitationally lensed quasars, of which about 300 will be quad lensed quasars (Oguri & Marshall 2010).
The structure of the paper is as follows. In Section 2, we measure the relative astrometry, photometry, and morphology of this system. In Section 3, we infer photometric redshifts, and in Section 4, we conduct a photometric variability study to support the photometric redshifts. In Section 5, we fit a lensing mass model to the system. We present our conclusions and future works in Section 6. We assume a cosmological model with Ω M = 0.274, Ω L = 0.726, and h = 0.71. All magnitudes are in the AB system except for the WISE ones, which are in the Vega system.

ASTROMETRY, PHOTOMETRY, AND MORPHOLOGY
We used the stacked images from PS1 in all five filters in order to perform a morphological modeling. We present a close-up color image of the system in Fig. 1, which clearly shows four objects within 3.8 ′′ of each other. The three brightest of these (A, B and C) are arranged in an arc-like configuration. The PS1 catalog identifies only the sources A and D, with B and C being blended with A.
The simple modeling with a point spread function (PSF) constructed from stars in the PS1 FOV left large residuals at the locations of the four objects, which were significant enough to affect the derived parameters. We suspect that this is due to spatial variation in the PSF across the PS1 FOV, and have therefore adopted a technique of fitting the system with an analytical PSF, which, it is assumed will not change on the small scale represented by the system. We did this using Hostlens (Rusu et al. 2016), a variant of glafic (Oguri 2010) that uses χ 2 minimization in order to fit point sources as well as Sérsic (Sérsic 1963) profiles convolved with an analytical PSF. The analytical PSF comprises two concentric Moffat (1969) profiles, each one with its own FWHM, ellipticity, orientation, and shape parameter. The two profiles are also characterized by the relative flux of one to the other. We successfully modeled the four objects as point sources convolved with the PSF, therefore showing that they are point like. Our strategy was to run Hostlens starting from 100 different positions in the parameter space, select the best of the resulting models, and further run 10 Markov Chain Monte Carlo (MCMC) chains using the Metropolis-Hastings algorithm around it. The chains consist of one million steps, with an acceptance rate of ∼ 0.3, and we removed the first one-fifth of these (the "burn-in" steps). We also checked that the chains have converged, using the method of Gelman (1995). The MCMC-derived uncertainties between the various analytical parameters are shown in Fig. 2, and generated using corner.py (Foreman-Mackey 2016). In order to quantify the additional uncertainties inherent in our method due to the non-analytical nature of the "true" PSF, we ran 100 simulations where we added noise of similar properties to the real images on top of the best-fit model. Here, we created the best-fit model by using a PSF constructed from nearby stars, and we ran Hostlens on each of the 100 simulations, using an analytical PSF. As a final estimate of our uncertainties, we used the maximum between the MCMC-derived uncertainties and those from the simulations. We were able to obtain much improved residuals in all bands, and we show these in Fig. 3. Fig. 3 shows conspicuous residuals in the i, z and y bands (middle row). These are consistent with the discovery of a lensing galaxy, in support of the lensing nature of this system. Due to the relatively low signal-to-noise ratio and the proximity of the point sources, we were unable to fit the morphology of this galaxy with free parameters. The red colors, suggested by the non-detection in the bluer bands, suggest that it is a red, early-type galaxy. As a result, we modeled it with a Sérsic index of 4, typical of early-type galaxies.
We checked that this produces a better fit than a Sérsic index of 1. For the Sérsic profile, we used an axis ratio of 1 and we fixed the effective radius at 0.5 ′′ . As the galaxy is most conspicuous in the i and y bands, we consider the most reliable estimate of the relative positions of all objects to be derived from the weighted average of the positions measured in these two bands, where we weight by the inverse of the measured uncertainties. The resulting astrometry and photometry in each band are shown in Table 1. We have corrected the magnitudes for galactic extinction using the maps of Schlafly & Finkbeiner (2011).
In addition to the PS1 data, we have also looked for archival data of this system. We use the Gaia position for component D as the absolute astrometric reference position for this system: 26.792307, 46.511273. The errors for these coordinates given in the Gaia catalog are 15.6 and 8.9 mas, respectively.
The system is also bright in the infrared, with the unresolved magnitudes being 11.524±0.022, 10.434±0.020, 6.769±0.015, and 4.518±0.023 in the WISE bands − W1, W2, W3, and W4, respectively. Finally, the system is very likely also a radio source. The NRAO VLA Sky Survey (NVSS) images show a source, 2MASX J01471020+4630433, with flux 12.6 mJy at 1.4 GHz. The listed position for this source is only 3 ′′ away from the center of the system in the optical bands, which is comparable with NVSS astrometric errors (Condon et al. 1998).
We are currently in the process of further observing the system in the radio with the Very Large Array (VLA). Radio data are useful for the study of lensed quasars, as it provides high-resolution imaging of the source, as well as fluxes not affected by microlensing and extinction.

PHOTOMETRIC REDSHIFTS
We fit the photometry measured in Section 2 for the point-like sources and the faint red galaxy using spectral energy distribution (SED) templates. The colors of the point-like images are very similar (Table 1), consistent with being multiple images of the same source.
Due to their point-like morphology, we suspect these to be quasars, and have therefore fit them with quasar templates in order to check the quality of the fit and to infer photometric redshifts. Photometric redshifts of quasars are notoriously difficult to constrain, due to the relatively featureless continua (Richards et al. 2009) and the variable equivalent widths of the broad emission lines. We measured redshifts using several methods: • The method of Wu & Jia (2010) uses derived quasar colors as a function of redshifts.
This is based on SDSS colors. PS1 bands are similar but not identical to the SDSS bands, and we used the corrections of Finkbeiner et al. (2016) to obtain magnitudes in the SDSS system. Following Wu & Jia (2010), we minimized χ 2 to obtain redshifts for each of the point sources. We only used three of the colors as the Y band from UKIRT is very different from the PS1 y band. We plot χ 2 ν in Figure 4 for each quasar image, and also by fitting all sources together. We notice that the curves are quite consistent with each other, as expected from the similar colors of the sources. In the latter case, we obtain the best fit at z = 0.820 +0.018 −0.014 . We notice a second minimum at z ∼ 2.6. The χ 2 ν for the minima is large, which could be due to the large spread seen in the colors of Wu & Jia (2010). We note that no numerical uncertainties are listed in that paper for the derived colors. We have also explored incorporating a prior based on the quasar luminosity function, given the unmagnified observed magnitudes of the source (see Section 5). While the prior prefers low-redshift values, its effect on Fig. 4 is negligible.
• A second method uses the photometric redshift code LePhare (Arnouts et al. 1999;Ilbert et al. 2006) to fit a quasar template to the observed colors. We obtain the best fit at a redshift of ∼ 2.6 − 2.8 for all four point sources. These results agree with the second peak of the method above. The χ 2 is between 13 and 125 for 4 degrees of freedom (5 filters). We find poorer fits (larger χ 2 ) when using stellar or galactic SED templates, providing further support that the sources are quasars. In particular, the best-fit galaxy templates correspond to starburst galaxies at low redshift 0.1. At such low redshifts the sources would likely be resolved, not point like.
• A combination of SDSS and WISE photometry provides good discriminators between high-and low-redshift quasars. DiPompeo et al. (2015) used the ratio between WISE W1 and W2 bands to the i-band flux, and found that high values are expected for quasars below a redshift of 1. We added the fluxes for A, B, C, and D in the i band (to account for the fact that the system is not resolved in WISE; G is of negligible flux), and we obtained a total magnitude of 14.4. We obtain ratios W1/i = 14.4 and W2/i = 39.0, which both indicate a redshift > 2, and are not compatible with the low value of ∼ 0.8, which would require ratios of W1/i = 50 and W2/i = 120. We also used the diagnostics from Wu et al. (2012). As seen in Fig. 13 of that paper, the colors z -W1 and W2 -W3 discriminate well between low and high redshift. We obtain z -W1 = 2.8 and W2 -W3 = 3.7, which both suggest a redshift larger than 2.
The lens galaxy G is very faint but we managed to fit a galaxy SED with an early-type template, resulting in a best-fit z = 0.57 +0.20 −0.13 . For this, we used the magnitudes measured in the izy bands and assumed detection limits of 21 mag in the gr bands, ∼ 2 mag fainter than the rest. We independently checked that the observed magnitudes are reasonable for a redshift of z = 0.57: for this redshift, the observed magnitudes imply an absolute magnitude of about −23.0 ∼ −23.5 in the R-band. From the Faber-Jackson relation (Faber & Jackson 1976), using the velocity dispersion we estimate from the lens modeling in Section 5, we obtain R ∼ −23.6, therefore not inconsistent.

URAT VARIABILITY
Here we undertake a quasar variability study as an independent check of the redshift estimate we derived in the previous section. URAT1 was a U. S. Naval Observatory project designed to provide accurate astrometry in the northern hemisphere down to about 18 mag.
A single filter centered at ∼ 720 nm was used, which is almost as red as the PS1 i filter.
Observations at the Flagstaff Station over more than 3 years provided more than 50 epochs for a large number of objects.
To characterize the variability of this system, we used a typical method used for quasars, the structure function V, as presented in Vanden Berk et al. (2004). It uses the magnitude differences between different epochs of observation to estimate the variability at different time intervals (time lag). It is defined as where ∆m is the measured magnitude difference for each pair of epochs in the light curve, σ is the statistical measurement uncertainty of ∆m and the brackets denote average quantities.
Vanden Berk et al. (2004) used the structure function to characterize the variability of a large sample of SDSS quasars. We note that their study is statistical in nature and therefore only provides general properties which we use to try to discriminate between the two photometric redshifts we obtained in the previous section.
In order to improve the results, we performed relative photometry of URAT epoch data using nearby stars in the 12 to 14 magnitude range. We removed stars that show variability and also removed some poor quality observations. We used more than 50 stars and 36 epochs spanning over two years. We estimate errors using the nearby stars with similar magnitudes. We obtained V = 0.074 and an average time lag of 297 days. According to Fig. 16 of Vanden Berk et al. (2004), this V value matches quasars at about z = 2 and therefore our photometric value of z = 2.6 is a better match than the low value of z = 0.8, which would require a V value of 0.05. We also used the absolute magnitude dependence in Fig. 11 of Vanden Berk et al. (2004). We estimate an absolute magnitude for our source quasar using the demagnified average of the observed magnitudes of A, B, and C (see Section 5). Assuming a composite quasar spectrum, we estimate that the i-band rest-frame absolute magnitude ranges between −24.8 and −25.2 at redshift 0.8, and between −27.6 and −27.9 at redshift 2.6. Again, we see that the V value better matches a high luminosity quasar, and therefore the higher redshift is preferred. We note that in this study we have ignored the variability caused by microlensing (see the next section).

LENS MODELING
The results of the previous sections show that the source is likely a quadruply imaged quasar, lensed by a foreground early-type galaxy. We subsequently modeled the observed configuration (the relative positions of the objects A,B,C,D, and G) using the lensmodel software from the gravlens package (Keeton 2001). We chose not to use the observed relative fluxes as constraints, as it is well known that at these wavelengths they usually show large discrepancies due to microlensing, millilensing, intrinsic variability, or extinction (e.g., Sluse et al. 2008). For the mass profile of the lensing galaxy, we use a Singular Isothermal Ellipsoid (SIE) with external shear. There are eight constraints from the four image positions and two constraints from the galaxy position, for a total of 10. Also, there are nine free parameters: two for the galaxy position, one for the Einstein radius, two for ellipticity, two for shear, and two for the source position. The model therefore has one degree of freedom.
For the best-fit mass model parameters, we find an Einstein radius b = 1.932 +0.008 −0.011 , ellipticity e c = 0.163 +0.063 −0.070 , e s = −0.035 +0.023 −0.024 , and shear γ c = 0.122 +0.015 −0.016 , γ s = 0.065 ± 0.005. Here, we expressed the ellipticity and shear in quasi-Cartesian coordinates, rather than polar coordinates, as the distribution of uncertainties is more gaussian. The two coordinate systems are related by (e c , γ c ) = (e, γ) cos 2(θ e , θ γ ) and (e s , γ s ) = (e, γ) sin 2(θ e , θ γ ). The derived image magnifications are approximately 19, 11, 11, and 0.7 for A, B, C and D, respectively. The uncertainties were determined using MCMC as plotted in Fig. 5, where the convergence of the chains was again checked using the method of Gelman (1995). The Ds , a velocity dispersion of σ ∼ 500 km/s. Here, c is the speed of light, and D s , D ls are the angular diameter distances to the source and between the lens and source, respectively. This value is too large for a single galaxy, but if we employ the second-best source redshift estimate, z ∼ 2.6, which is preferred by LePhare, the combination of WISE and PS1 colors (Section 3) and the variability results in Section 4, we obtain σ ∼ 315 km/s, a more common value for massive early-type galaxies. We therefore assume this redshift to infer the predicted time delays, and we obtain ∆t BC = 0.1 days, ∆t AC = 1.7 days and ∆t DC = 226.2 days, with image C (a minimum of the time delay surface) leading. Image B is another minimum, and images A and D are saddle points.
In Figure 6 we show the critical curve and the caustics, as well as the positions of the source and the observed images. The source is crossing a cusp, which results in the production of three bright images in close proximity and one isolated faint image. This is a classic configuration for quad lenses. The predicted magnitude differences are shown in Table 1, and a comparison with the observed values are shown in Figure 7.
We find that the predicted brightness ranking of the images is consistent with the observations, with image A being brightest and image D being faintest. However, there are still large discrepancies between the observed and predicted fluxes. For example, the fluxes of images B and C are different, whereas they are expected to be the same for any smooth mass profile of the lensing galaxy. Such "flux anomalies" have been widely encountered for cusp lensed quasars and can be explained by small-scale structure in the lensing galaxy (e.g., Keeton, Gaudi, & Petters 2003). We also remark that the colors of the four images are very similar, within 0.2 mag across all filters. Small chromatic variations between multiple quasar images are known to be caused by intrinsic variability, microlensing, or differential extinction (Yonehara et al. 2008). In any case, differential extinction would not explain the flux discrepancy in image D, as this is ∼ 1 mag brighter than expected.
Since image D is closest to the lensing galaxy, the stellar density is comparatively high at its location, meaning that microlensing due to stars in the lensing galaxy is a plausible explanation for the discrepancy. Second, we remark that the time delays between image D and the other images is large, which means that while intrinsic flux variations between A, B and C will be washed out, they could stand out in image D, and therefore be responsible for the discrepancy. Lastly, small-scale substructure (millilensing) could cause contribute to the discrepancy.

CONCLUSIONS AND FUTURE WORK
We have presented the first quadruply imaged gravitationally lensed quasar candidate in the Pan-STARRS1 Survey, discovered via visual inspection of the multiband images.
We find that the evidence supporting the gravitational lens nature of this system is overwhelming and consists of the following: • The presence of four point-like images within 3.8 ′′ of each other, with similar colors, more consistent with quasar templates, but less so with stars and galaxies.
• The detection of a red galaxy between the point source, such that the relative positions of all five objects are fully consistent with a well-known "cusp" configuration. We successfully reproduce this configuration using a SIE + shear model for the lensing galaxy. Furthermore, the brightness ranking of the point-like images is consistent with the one predicted by the model.
• The photometric redshift analysis shows that the galaxy is of elliptical template, typical for lensing galaxies, and of redshift smaller than the point sources. For the inferred redshifts (assuming a source redshift of ∼ 2.6), the velocity dispersion of the lensing galaxy is consistent with those of massive early-type galaxies.
Nevertheless, we stop short of claiming that this is a confirmed lens until we can obtain spectroscopy of this system. Spectroscopy of the four point-like images would unequivocally demonstrate whether they are multiple images of the same background quasar, as well as allow one to infer its redshift and physical properties. We also note that the system is accessible to ground-based high-resolution observations at most adaptive optics-capable facilities in the northern hemisphere, due to the its proximity to a R ∼ 12 mag star ∼ 18 ′′ away.
Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement.

3.55±0.06
Note. -The first four columns are the four quasar images and the fifth column is the lens. We also label the images according to the arrival time. The first five rows are the magnitude measurements derived using Hostlens, corrected for extinction. The next two rows are the relative positions of each object relative to A. We assume a pixel scale of 0.256 ′′ . ∆α is positive toward the west, and ∆δ toward the north. The last row shows the model prediction for the relative magnitudes relative to image A, from the best-fit mass model.  (1,2,3,4), the same for the second Moffat profile (5,6,7,8), the relative flux as flux1/(flux1 + flux2) (9), sky level in counts (10), positions in pixels on the x and y axes of point source A (11,12) and its flux in counts (13), and the same for objects B,C and D (14-25).
-25 - Fig. 3.-Original images (10 ′′ × 10 ′′ , first row) and residuals after subtracting the best-fit model with Hostlens. From left to right: g, r, i, z, and y bands. In the second row, the lensing galaxy is not included in the modeling, whereas in the third row it is. The images are in linear scale and cover the full dynamic range. The remaining structure in the residuals amounts to less than 6% of the peak fluxes and can be attributed to the large Poisson noise at the center of the point-like images, as well as to the faint arcs expected from the lensed quasar host galaxy, which are not accounted for in the modeling.