# Coherent and incoherent $\mathrm{J/}\psi$ photonuclear production in an
energy-dependent hot-spot model

J. Cepila  Faculty of Nuclear Sciences and Physical Engineering, Czech
Technical University in Prague, Czech Republic  J. G. Contreras  Faculty of
Nuclear Sciences and Physical Engineering, Czech Technical University in
Prague, Czech Republic  M. Krelina  Departamento de Física, Universidad
Técnica Federico Santa María; Centro Científico-Tecnológico de Valparaíso-
CCTVal, Casilla 110-V, Valparaíso, Chile  Faculty of Nuclear Sciences and
Physical Engineering, Czech Technical University in Prague, Czech Republic

###### Abstract

In a previous publication, we have presented a model for the photoproduction
of $\mathrm{J/}\psi$ vector mesons off protons, where the proton structure in
the impact-parameter plane is described by an energy-dependent hot-spot
profile. Here we extend this model to study the photonuclear production of
$\mathrm{J/}\psi$ vector mesons in coherent and incoherent interactions of
heavy nuclei. We study two methods to extend the model to the nuclear case:
using the standard Glauber-Gribov formalism and using geometric scaling to
obtain the nuclear saturation scale. We find that the incoherent cross section
changes sizably with the inclusion of subnucleonic hot spots, and that this
change is energy dependent. We propose to search for this behavior by
measuring the ratio of the incoherent to coherent cross section at different
energies. We compare the results of our model to results from RHIC and from
the Run 1 at the LHC finding a satisfactory agreement. We also present
predictions for the LHC at the new energies reached in Run 2. The predictions
include $\mathrm{J/}\psi$ production in ultra-peripheral collisions, as well
as the recently observed photonuclear production in peripheral collisions.

###### pacs:

24.85.+p,25.20.Lj,14.40.Pq

##  I Introduction

Diffractive photoproduction of $\mathrm{J/}\psi$ vector mesons in high-energy
interactions is a sensitive probe of the gluon distribution of hadrons in the
small-$x$ region, where saturation effects are expected to be important
Albacete and Marquet (2014). As such, it has been extensively studied
experimentally at the LHC Contreras and Tapia Takaki (2015); Andronic  _et
al._ (2016) and it is an essential component of the research program envisaged
for future facilities Accardi  _et al._ (2016); Abelleira Fernandez  _et al._
(2012).

The idea of describing the structure of protons with subnucleonic degrees of
freedoms representing regions of high-gluon density, so-called hot spots, has
recently yielded very interesting results. In particular, it has been found
that, using a Good-Walker formalism Good and Walker (1960), the dissociative
photoproduction of $\mathrm{J/}\psi$ off proton targets was sensitive to
geometric fluctuations of the positions of hot spots in the impact-parameter
plane of the interaction Mäntysaari and Schenke (2016a, b). In this model, the
proton is made up of three hot spots, and the comparisons to experimental data
is performed at a fixed center-of-mass energy ($W_{\gamma\mathrm{p}}$) of the
photon-proton system. Recently, this model has also been extended to the case
of photonuclear production at the LHC Mäntysaari and Schenke (2017).

In a further development, a model in which the number of hot spots increased
with decreasing $x$ was presented in Cepila  _et al._ (2017). Such a change in
the transverse profile of the target is qualitatively similar to the formal
results obtained by evolving the JIMWLK equation and depicted in Fig. 1 of
Schlichting and Schenke (2014). The model presented in Cepila  _et al._ (2017)
was able to correctly describe experimental data on the $W_{\gamma\mathrm{p}}$
dependence of both, the exclusive and the dissociative production of
$\mathrm{J/}\psi$. In this model, the energy dependence of the dissociative
process shows a sharp decrease, quantitatively different than expectations
from extrapolations of HERA data, providing a new signature of saturation
effects, which could be measurable at current LHC energies Cepila  _et al._
(2017).

Here, we extend our model from Cepila  _et al._ (2017) to the case of
photonuclear interactions. We follow two different approaches to go from
proton to nuclear targets. One uses the standard Glauber-Gribov formalism as
proposed in Armesto (2002), while the other relies on geometric scaling to
compute the nuclear saturation scale using as an input the saturation scale in
the proton Armesto  _et al._ (2005).

In agreement with the results of Mäntysaari and Schenke (2017), we show that
the incoherent photonuclear production of $\mathrm{J/}\psi$ is sensitive to
the hot spot structure of nucleons. In addition, we predict that this effect
is energy-dependent and that this dependence could be observed by measuring
the ratio of the incoherent to the coherent photonuclear cross section as a
function of the energy of the interaction.

We compare the results of our model with data from RHIC and from the Run 1 at
the LHC and make predictions for measurements to be performed at the LHC using
Run 2 data, including the prediction of $\mathrm{J/}\psi$ photonuclear
production in both, peripheral and ultra-peripheral interactions.

These studies of fluctuations in the transverse structure of hadronic targets
have been shown to be relevant in other contexts besides the diffractive
photonuclear production of vector mesons. For example, similar ideas have been
applied to explain the hollowness effect in proton-proton interactions at high
energies Albacete and Soto-Ontoso (2017). Another area where such models have
been applied is the initial state of collisions involving relativistic nuclei
where effects of initial spatial asymmetries may impact on key measurements
and their interpretations, e.g. Welsh  _et al._ (2016); M ntysaari  _et al._
(2017).

##  II Description of the formalism

Figure 1:  Diffractive photonuclear production of $\mathrm{J/}\psi$ in A-A
collisions. In diagram (a) the photon interacts with the full nucleus A, while
in (b) it scatters off one of the nucleons in A. Accordingly these processes
are called coherent and incoherent $\mathrm{J/}\psi$ production, respectively.

The processes which we study here are shown in Fig 1. Two heavy nuclei A
approach at high energies. One of the nuclei coherently emits a quasi-real
photon, which collides with the other nucleus. We are interested in the
following two processes: in the coherent case the photon interacts with the
full nucleus, while in the incoherent case it scatters off a single nucleon in
the nucleus. In both cases a $\mathrm{J/}\psi$ vector meson is produced.

The center-of-mass energy of the photon-nucleus system is denoted by
$W_{\gamma\mathrm{A}}$. The square of the four-momentum transfer in the
nucleus vertex is denoted by $t$ and it is related to the transverse momentum
($p_{T}$) of the produced $\mathrm{J/}\psi$, in the laboratory frame, by
$-t=p^{2}_{T}$. The $t$ distribution is related to the transverse distribution
of matter in the target. For the case of a large nucleus the $\mathrm{J/}\psi$
produced in a coherent process has a transverse momentum of just a few tens of
MeV/$c$, while in the incoherent process it reaches a few hundred of MeV/$c$,
reflecting the difference in size between the nucleon and the nucleus.

###  II.1 The nucleus-nucleus cross section

The cross section $d\sigma_{\rm AA}/dy$ for the photoproduction of a
$\mathrm{J/}\psi$ at rapidity $y$ in collisions of nuclei can be factorized as
the product of the photon flux $n_{\gamma}(y,\\{\vec{b}\\})$ and the
photonuclear cross section $\sigma_{\gamma\mathrm{A}}(y)$ as:

$\frac{d\sigma_{\rm
AA}}{dy}=n_{\gamma}(y;\\{\vec{b}\\})\sigma_{\gamma\mathrm{A}}(y)+n_{\gamma}(-y;\\{\vec{b}\\})\sigma_{\gamma\mathrm{A}}(-y),$
(1)

where $\\{\vec{b}\\}$ delimits the range in impact-parameter $\vec{b}$ taken
into account in the interaction. The two terms in this equation reflect the
fact that any of the two nuclei can act as the source of the photon.

In the following, we will consider the rapidity $y$ in Eq. (1) as given in the
laboratory frame. It is defined with respect to the direction of the target.
The rapidity of the $\mathrm{J/}\psi$ is related to the center-of-mass energy
of the photon-nucleus system through

$W_{\gamma\mathrm{A}}^{2}=\sqrt{s_{\rm NN}}M_{\mathrm{J/}\psi}e^{-y},$ (2)

where $M_{\mathrm{J/}\psi}$ is the mass of the $\mathrm{J/}\psi$ vector meson
and $\sqrt{s_{\rm NN}}$ is the center-of-mass energy per nucleon pair in the
A-A system.

The formalism for the computation of the flux has been described in detail in
Contreras (2017). The case for ultra-peripheral collisions (UPC) was
originally proposed in Klein and Nystrand (1999). It is well known and it has
been used extensively. For the case of peripheral collisions and for the
centrality class we are interested in (70–90% centrality class) the formalism
presented in Contreras (2017) produces numerically similar results as other
proposals to compute the flux in peripheral collisions as discussed in Klusek-
Gawenda and Szczurek (2016); Gay Ducati and Martins (2017).

###  II.2 The photon-nucleus cross section

The photon-nucleus cross section appearing in Eq. (1) is given by the integral
over $t$ of the following cross sections:

$\left.\frac{d\sigma_{\gamma\mathrm{A}}}{dt}\right|^{\rm
coh}_{T,L}=\frac{(R^{T,L}_{g})^{2}}{16\pi}\left|\left<A^{j}(x,Q^{2},\vec{\Delta})_{T,L}\right>_{j}\right|^{2}$
(3)

for the coherent process, and

$\displaystyle\left.\frac{d\sigma_{\gamma\mathrm{A}}}{dt}\right|^{\rm
inc}_{T,L}$ $\displaystyle=$
$\displaystyle\frac{(R^{T,L}_{g})^{2}}{16\pi}\left(\left<\left|A^{j}(x,Q^{2},\vec{\Delta})_{T,L}\right|^{2}\right>_{j}\right.$
(4)
$\displaystyle-\left.\left|\left<A^{j}(x,Q^{2},\vec{\Delta})_{T,L}\right>_{j}\right|^{2}\right)$

for incoherent production, where the average is over different geometrical
configurations $j$ of the nucleons, respectively hot spots, inside the
nucleus, as detailed below.

The amplitude is given by

$\displaystyle A^{j}(x,Q^{2},\vec{\Delta})_{T,L}$ $\displaystyle=$
$\displaystyle i\int d\vec{r}\int^{1}_{0}\frac{dz}{4\pi}(\Psi^{*}\Psi_{\rm
V})_{T,L}$ (5) $\displaystyle\int
d\vec{b}\;e^{-i(\vec{b}-(1-z)\vec{r})\cdot\vec{\Delta}}\left(\frac{d\sigma_{\rm
dA}}{d\vec{b}}\right)_{j},$

where $\Psi_{\rm V}$ is the wave function of the vector meson, $\Psi$ is the
wave function of a virtual photon fluctuating into a quark-antiquark color
dipole, $\vec{r}$ is the transverse distance between the quark and the
antiquark, and $z$ is the fraction of the longitudinal momentum of the dipole
carried by the quark.

The amplitude depends on the virtuality of the quasi-real photon $Q^{2}$,
which we take to be $Q^{2}=0.05$ GeV2. It also depends on
$\vec{\Delta}^{2}=-t$ and on
$x=(M_{\mathrm{J/}\psi}/W_{\gamma\mathrm{A}})^{2}$. The indices $T$ and $L$
represent the contributions of transversely and longitudinally polarized
photons, respectively. The total cross section is the sum of the $T$ and $L$
contributions. The QCD part of the model is contained in the dipole-target
cross section $d\sigma_{\rm dA}/d\vec{b}$, which is discussed below.

The so-called skewedness correction Shuvaev  _et al._ (1999) is given by

$R^{T,L}_{g}(\lambda^{T,L}_{g})=\frac{2^{2\lambda^{T,L}_{g}+3}}{\sqrt{\pi}}\frac{\Gamma(\lambda^{T,L}_{g}+5/2)}{\Gamma(\lambda^{T,L}_{g}+4)},$
(6)

with $\lambda^{T,L}_{g}$ defined as

$\lambda^{T,L}_{g}\equiv\frac{\partial\ln(A_{T,L})}{\partial\ln(1/x)}.$ (7)

As we did in Cepila  _et al._ (2017), for $\Psi_{T,L}$ we use the definitions
and parameter values of Golec-Biernat and Wusthoff (1998) and for the wave
function of the vector meson, $\Psi_{\rm V}$, we use the boosted-Gaussian
model Nemchik  _et al._ (1994, 1997), with the numerical values of the
parameters as in Kowalski  _et al._ (2006).

###  II.3 The dipole-nucleus cross section

We have followed two approaches to use the dipole-nucleon cross section
determined in Cepila  _et al._ (2017) to define the dipole-nucleus cross
section without introducing any new additional free parameters.

In one approach, denoted GG in the figures, we chose the Glauber-Gribov
methodology proposed in Armesto (2002), which relates the dipole-proton (dp)
to the dipole-nucleus cross section via the nuclear profile $T_{\rm
A}(\vec{b})$:

$\left(\frac{d\sigma_{\rm
dA}}{d\vec{b}}\right)_{j}=2\left[1-\exp\left(-\frac{1}{2}\sigma_{\rm
dp}(x,r)T^{j}_{\rm A}(\vec{b})\right)\right],$ (8)

where we use the Golec-Biernat and Wusthoff Golec-Biernat and Wusthoff (1998)
model for the dipole–proton cross section:

$\sigma_{\rm
dp}(x,r)=\sigma_{0}\left[1-\exp\left(-r^{2}Q^{2}_{s}(x)/4\right)\right],$ (9)

with $\sigma_{0}=\pi R^{2}_{\rm p}$, $R_{\rm p}$ the proton radius, and the
saturation scale given by

$Q^{2}_{s}(x)=Q^{2}_{0}(x_{0}/x)^{\lambda},$ (10)

where the values of the parameters are as in Cepila  _et al._ (2017).

In the second approach we use a model that factorizes the $\vec{b}$ and $x$
dependences as proposed in Cepila  _et al._ (2017):

$\left(\frac{d\sigma_{\rm dA}}{d\vec{b}}\right)_{j}=\sigma^{\rm
A}_{0}\left[1-\exp\left(-r^{2}Q^{2}_{A,s}(x)/4\right)\right]T^{j}_{\rm
A}(\vec{b}).$ (11)

Here, as in Cepila  _et al._ (2017), $\sigma^{\rm A}_{0}$ is related to the
area of the target by $\sigma^{\rm A}_{0}=\pi R^{2}_{\rm A}$, where $R_{\rm
A}$ is the ratio of the nucleus as defined in the corresponding Woods-Saxon
distribution. The saturation scale of the nucleon is related, through
geometric scaling, to the saturation scale of the nucleus as proposed in
Armesto  _et al._ (2005):

$Q^{2}_{s,{\rm A}}(x)=Q^{2}_{s}(x)\left(\frac{{\rm A}\pi R^{2}_{\rm p}}{\pi
R^{2}_{\rm A}}\right)^{\frac{1}{\delta}},$ (12)

with $\delta=0.8$ as found in Armesto  _et al._ (2005) and $\sigma_{0}$ fixed
from the analysis in Cepila  _et al._ (2017). This model is denoted GS in the
figures.

We explore two models for the nuclear profile, one considers the nucleus to be
made up of nucleons and the second includes subnucleonic degrees of freedom in
the form of hot spots. (In the figures these approaches are denoted by n and
hs, respectively.) The key distinguishing feature of our model, with respect
to that presented in Mäntysaari and Schenke (2017), is that the number of hot
spots in a given nucleon increases with decreasing $x$ as proposed in Cepila
_et al._ (2017).

In detail, we have for the first case that each nucleon has a Gaussian profile
of width $B_{\rm p}$, centered at random positions $\\{\vec{b}_{i}\\}^{j}$
sampled from a Woods-Saxon nuclear profile for each configuration $j$:

$T^{j}_{\rm A}(\vec{b})=\frac{1}{2\pi B_{\rm p}}\sum^{\rm
A}_{i=1}\exp\left(-\frac{(\vec{b}-\vec{b}^{j}_{i})^{2}}{2B_{\rm p}}\right).$
(13)

For the hot-spot model we have that the nucleons are themselves made up of hot
spots, distributed according to a Gaussian of width $B_{\rm hs}$.

$T^{j}_{\rm A}(\vec{b})=\frac{1}{2\pi B_{\rm hs}}\sum^{\rm
A}_{i=1}\frac{1}{N_{\rm hs}}\sum^{N_{\rm
hs}}_{k=1}\exp\left(-\frac{(\vec{b}-\vec{b}^{j}_{i}-\vec{b}^{j}_{k})^{2}}{2B_{\rm
hs}}\right),$ (14)

where in this case $N_{\rm hs}$ is a random number drawn from a zero-truncated
Poisson distribution, where the Poisson distribution has a mean value

$\langle N_{\rm hs}(x)\rangle=p_{0}x^{p_{1}}(1+p_{2}\sqrt{x}).$ (15)

All the parameters of the model have been fixed from a comparison to data on
$\mathrm{J/}\psi$ photoproduction off protons Cepila  _et al._ (2017). The
values of the parameters and the associated discussion can be found in Cepila
_et al._ (2017) and will not be repeated here.

In the GS model, the integrals over impact parameter $\vec{b}$ are done
analytically and factorize from the rest of the integrals. This allowed us to
use 10000 different profile configurations, which increased the numerical
precision of the computation to the percent level.

On the other hand, the GG model has to be integrated numerically also over
impact parameter. We have used globally adaptive subdivision with importance
sampling as implemented in the Suave method Hahn (2005). The program has been
configured to reach either 2% precision or to stop after a given number of
maximum tries. This last criteria was needed because the convergence was slow.
It has been checked that in no case the numerical precision of the method
exceeded 4%. The previous discussion corresponds to the numerical integration
of one configuration. But the cross sections are proportional to the average
(the variance) over configurations for the coherent (incoherent) process. Due
to the amount of computer resources needed in the GG model, we have only used
200 configurations for this case. This translates into a $\approx 2\%$
($\approx 13\%$) uncertainty for the coherent (incoherent) results shown in
Figs. 2 and 3. Similar observations for the fluctuations in the variance have
already been noted in Toll and Ullrich (2013). There is another source of
numerical uncertainty related to the fact that, to speed up the computation,
we used a grid in impact parameter space to store the profiles and then
interpolations have been used to obtain the specific values of the profile at
the particular value of $\vec{b}$ needed at a given point. The effect of this
approximation can be seen in Fig. 2 as a difference between the computation
for the hot-spot and the nucleon profiles for the coherent cross section.

##  III Predictions of the model and comparison to data

###  III.1 Results for the photonuclear cross sections

Figure 2:  Predictions for the coherent photonuclear production of $\mathrm{J/}\psi$ off lead targets as a function of $x$ ($W_{\gamma\mathrm{Pb}}$) compared with data extracted from ALICE measurements Abelev  _et al._ (2013a); Abbas  _et al._ (2013); Adam  _et al._ (2016), as explained in Contreras (2017).  Figure 3:  Predictions for the incoherent photonuclear production of $\mathrm{J/}\psi$ off lead targets as a function of $x$ ($W_{\gamma\mathrm{Pb}}$) compared with data extracted from ALICE measurements Abbas  _et al._ (2013), as explained in Contreras (2017). Figure 4:  Ratio of the incoherent to the coherent cross section for the photonuclear production of $\mathrm{J/}\psi$ off lead targets as a function of $x$ ($W_{\gamma\mathrm{Pb}}$) compared with ALICE Abbas  _et al._ (2013).  Figure 5:  Predictions for the $t$ dependence at mid-rapidity for the LHC Run 1 of the coherent and incoherent photonuclear cross section using the GG approach with (hs) or without (n) subnucleonic degrees of freedom.  Table 1:  Predictions of our model for the photonuclear production of $\mathrm{J/}\psi$ and comparison with data from Au-Au collisions measured by PHENIX at $x=0.015$ Afanasiev  _et al._ (2009).  Source |  $\sigma_{\gamma\rm Au}$ ($\mu$b)
---|---
PHENIX coherent Afanasiev  _et al._ (2009) | 5.7 $\pm$ 2.3 (stat) $\pm$ 1.2 (syst)
PHENIX incoherent Afanasiev  _et al._ (2009) | 3.6 $\pm$ 1.4 (stat) $\pm$ 0.7 (syst)
GS-hs coherent | 6.9
GS-hs incoherent | 4.5
GS-n coherent | 6.9
GS-n incoherent | 1.8
GG-hs coherent | 10.4
GG-hs incoherent | 2.4
GG-n coherent | 11.2
GG-n incoherent | 1.2

The predictions of our model for coherent photonuclear production of
$\mathrm{J/}\psi$ off lead are shown in Fig. 2. We present results for the
case where the dipole-nucleus cross section was computed with the Glauber-
Gribov approach using either a nuclear profile made up of nucleons or of hot
spots (GG-n and GG-hs, respectively), as well as for the equivalent
calculation performed within the geometric-scaling methodology (denoted by
GS-n and GS-hs, respectively). These results are compared to data extracted
from ALICE measurements in peripheral Adam  _et al._ (2016) and ultra-
peripheral collisions Abelev  _et al._ (2013a); Abbas  _et al._ (2013) by
taking into account the photon flux as explained in Contreras (2017).

The predictions from the GS-hs and GS-n cases are very similar for coherent
production and cannot be distinguished in the figure, while in the GG approach
there is a small difference in the prediction when considering a nuclear
structure made of hot spots or made of nucleons. The difference can be traced
back to the limited number of con

gurations explored and to the granularity of the grid used to store the
different profiles shown in Eqs. (13) and (14) as explained above. The
predictions of the GS approach give an excellent description of data, while
the GG describe the lowest $x$ data, but overestimate the measurements at
larger values of $x$.

The predictions of our model for incoherent photonuclear production of
$\mathrm{J/}\psi$ off lead are shown in Fig. 3. In this case, there is a clear
difference in the predictions when using a profile with subnucleonic degrees
of freedom, with respect to the standard nuclear profile made of nucleons. The
comparison with the ALICE measurement Abbas  _et al._ (2013) suggests that
data seems to prefer the hs models, specially for the GG computation. The same
impression is obtained by looking at the ratio of the incoherent to coherent
cross sections, as shown in Fig. 4. In addition, it can be seen in this
figure, that the subnucleonic degrees of freedom introduce a dependence on $x$
for this ratio.

In order to exploit this $x$ dependence, we compare in Tab. 1 the results of
our model for Au-Au collisions with data measured by PHENIX Afanasiev  _et
al._ (2009) at RHIC. This measurement was performed at mid-rapidity for a
center-of-mass energy of the Au-Au system $\sqrt{s_{\rm NN}}=0.2$ TeV, which
corresponds to $x=0.015$. Even though the measurement has large uncertainties,
PHENIX result, in conjunction with ALICE data, supports the prediction that
the ratio of the incoherent to the coherent cross section for photonuclear
production of $\mathrm{J/}\psi$ depends on $x$. PHENIX data, as already seen
in the case of ALICE, also seems to prefer the hot-spot models.

Figure 5 shows the dependence of the coherent and incoherent cross sections on
$t$ at mid-rapidity for the LHC Run 1 energies for the GG approach. A similar
picture is obtained in the GS methodology. For the coherent case, the
diffractive structure is clearly seen, while the incoherent cross section
shows a large increase when including subnucleonic degrees of freedom in the
nuclear profile, with respect to the case of a nucleon-only composition.
Similar results have been shown in Mäntysaari and Schenke (2017).

###  III.2 Results for the Pb-Pb cross section

Figure 6:  Cross section for coherent photonuclear $\mathrm{J/}\psi$
production in Pb-Pb collisions at $\sqrt{s_{\rm NN}}=2.76$ TeV, corresponding
to the LHC Run 1 energies, as a function of rapidity. Predictions of the model
are compared to data measured in UPC by ALICE Abelev  _et al._ (2013a); Abbas
_et al._ (2013) and CMS Khachatryan  _et al._ (2017), as well as data measured
by ALICE Adam  _et al._ (2016) in peripheral collisions in the 70-90%
centrality class. Figure 7:  Cross section for incoherent photonuclear
$\mathrm{J/}\psi$ production in Pb-Pb collisions at $\sqrt{s_{\rm NN}}=2.76$
TeV, corresponding to the LHC Run 1 energies, as a function of rapidity.
Predictions of the model are compared to data measured in UPC by ALICE Abbas
_et al._ (2013). Figure 8:  Cross section for coherent photonuclear
$\mathrm{J/}\psi$ production in Pb-Pb collisions at $\sqrt{s_{\rm NN}}=5.02$
TeV, corresponding to the LHC Run 2 energies, as a function of rapidity.
Figure 9:  Cross section for coherent photonuclear $\mathrm{J/}\psi$
production in Pb-Pb collisions at $\sqrt{s_{\rm NN}}=5.02$ TeV, corresponding
to the LHC Run 2 energies, as a function of rapidity.

The predictions of our model for nucleus-nucleus cross sections, see Eq. (1),
are shown in Figs. 6 to 9. The predictions for energies corresponding to the
LHC Run 1 are compared to the available data. As of now there is no published
data from the LHC Run 2 period for these observables.

In peripheral collisions, we have used the results from ALICE Abelev  _et al._
(2013b) to convert the centrality class to a range in impact parameter.
According to which, the 70-90% centrality class corresponds to a range
(13.05,14.96) fm in impact parameter at Run 1 energies. The same formalism
applied to LHC Run 2 energies yields a range (13.1,15.0) fm.

Figure 6, shows that the GS-hs version of our model matches quite well the
data for coherent production available from ALICE Abelev  _et al._ (2013a);
Abbas  _et al._ (2013) and CMS Khachatryan  _et al._ (2017) in UPC, as well as
the data measured in peripheral interactions by ALICE Abbas  _et al._ (2013).
The GG-hs version of the model slightly overestimates the data. A similar
picture is obtained when comparing the prediction and the measurement for
incoherent production as shown in Fig. 7.

Finally, Figs. 8 and 9 show the prediction of our model for Pb-Pb collisions
at an energy corresponding to the LHC Run 2. All LHC collaborations already
took Run 2 data in 2015 and a second data taking period is expected towards
the end of 2018. These data, will produce measurements which are expected to
have small experimental uncertainties.

##  IV Discussion

The key new element of our computations, with respect to those presented in
Mäntysaari and Schenke (2017) is the energy dependence of the hot-spot
structure of the target in the impact-parameter plane. The parameters defining
this dependence have been fixed in Cepila  _et al._ (2017). The most striking
product of this difference is that the ratio of the incoherent to coherent
cross section develops a dependence on energy, decreasing as the energy ($x$)
is increase (decreased) as shown in Fig. 4.

Another difference of our model, with respect to that of Mäntysaari and
Schenke (2017) is the implementation of the transition from proton to nuclear
targets. The approach in Mäntysaari and Schenke (2017) is similar in spirit to
the Glauber-Gribov model used here, but they use a different prescription to
compute the dipole-proton interaction. Numerically the predictions from
Mäntysaari and Schenke (2017) differ from ours. In particular, for the
comparison with UPC data at mid rapidity, they overestimate the measurement
for both, the coherent and the incoherent production, while in our case the
GG-hs model overestimates the coherent process, but underestimates the
incoherent production.

At larger rapidites, the model from Mäntysaari and Schenke (2017) agrees with
the data from CMS in UPC slightly better than our predictions. They do not
offer predictions at the largest rapidities measured at the LHC, because this
region covers large values of $x$.

In phenomenological studies, it is customary to define the upper boundary of
the small $x$ region as $x\approx 0.01$. There is no strong quantitative
argument to chose exactly $x=0.01$ as the upper limit for these studies. Here,
we present results including values of $x$ up to 0.034, which corresponds to a
mild 25% change in $\ln(1/x)$ with respect to $x=0.01$. This value was chosen
to match the maximum rapidity coverage of the ALICE detector for the LHC
energies available in Run 2 and to be able to compare to measurements from
RHIC.

Another observation made in Mäntysaari and Schenke (2017) concerns the
dependence of the results on the chosen model for the wave function of the
$\mathrm{J/}\psi$, where the main effect between different models is an
overall change in normalization. In our model the normalization is given by
the $\sigma_{0}$ parameter which has been fixed by construction to the
effective area of the target. Once the normalization was fixed, we checked in
Cepila  _et al._ (2017) that the same parameters described correctly the
behavior of inclusive deeply inelastic scattering - which do not depend on the
wave function of the vector meson - at the appropriate scale. This means that
with the current set of parameters in our model, a change of model for the
wave function of the $\mathrm{J/}\psi$ would spoil the agreement with data, so
in that sense all the components of our model are fixed and we cannot change
them. Note that it is expected that the ratio of cross sections is less
susceptible to these effects, so that it is a more solid computation from the
phenomenology side, as well as having smaller experimental uncertainties.

The agreement between the predictions of our model and the experimental data
is noteworthy, because it has been reached without adding any new free
parameter to the model that successfully described the equivalent
photoproduction processes off proton targets. In particular, the comparison of
data with the predictions of the hot-spot models supports the existence of
subnucleonic degrees of freedom in the nucleus and that this hot-spot
structure evolves with $x$, respectively with $W_{\gamma\mathrm{A}}$.

Future measurements at RHIC and at the LHC are expected to have substantially
lower experimental uncertainties providing valuable constraints to improve our
picture of the nuclear structure. In particular, already existing data from
UPC at mid rapidity from RHIC and at forward rapidities from the LHC-Run 2
energies, are dominated by contributions at large $x$. The measurement of the
incoherent photonuclear production of $\mathrm{J/}\psi$ in these kinematic
domains, will confirm or disprove the $x$ ($W_{\gamma\mathrm{A}}$) dependence
of the incoherent to coherent cross section ratio.

##  V Summary, conclusions and outlook

We have presented a model of the coherent and incoherent photoproduction of
$\mathrm{J/}\psi$ in hadronic targets, which includes an energy evolution of
the QCD structure of the target in the impact-parameter plane. The parameters
of the model had been fixed before by data from HERA and the LHC on proton
targets. The model has been extended here to nuclear targets and compared to
existing data from RHIC and LHC. The agreement between data and the model
predictions is noteworthy, because it has been achieved without the addition
of any new free parameters.

The main new ingredient of our approach, namely the inclusion of an energy
evolution of the number of hot spots in the hadronic target, translates in
this case in an energy dependence of the ratio of the incoherent to coherent
cross sections. This prediction can be tested with already existing data from
RHIC and the LHC Run 2. Future data from the LHC will provide new and stronger
constraints to our model and help to understand better the subnucleonic
structure of hadronic targets and its energy evolution in the small $x$
regime.

## Acknowledgements

We would like to thank Daniel Tapia-Takaki for useful discussions. This work
has been partially supported by the following grants: Conicyt PIA/ACT 1406
(Chile), Conicyt PIA/Basal FB0821 (Chile) and LTC17038 of the INTER-EXCELLENCE
program at the Ministry of Education, Youth and Sports (Czech Republic).

## References

  * Albacete and Marquet (2014) J. L. Albacete and C. Marquet, Prog. Part. Nucl. Phys. 76, 1 (2014), arXiv:1401.4866 [hep-ph] .
  * Contreras and Tapia Takaki (2015) J. G. Contreras and J. D. Tapia Takaki, Int. J. Mod. Phys. A30, 1542012 (2015).
  * Andronic  _et al._ (2016) A. Andronic _et al._ , Eur. Phys. J. C76, 107 (2016), arXiv:1506.03981 [nucl-ex] .
  * Accardi  _et al._ (2016) A. Accardi _et al._ , Eur. Phys. J. A52, 268 (2016), arXiv:1212.1701 [nucl-ex] .
  * Abelleira Fernandez  _et al._ (2012) J. Abelleira Fernandez _et al._ (LHeC Study Group), J.Phys. G39, 075001 (2012), arXiv:1206.2913 [physics.acc-ph] .
  * Good and Walker (1960) M. L. Good and W. D. Walker, Phys. Rev. 120, 1857 (1960).
  * Mäntysaari and Schenke (2016a) H. Mäntysaari and B. Schenke, Phys. Rev. Lett. 117, 052301 (2016a), arXiv:1603.04349 [hep-ph] .
  * Mäntysaari and Schenke (2016b) H. Mäntysaari and B. Schenke, Phys. Rev. D94, 034042 (2016b), arXiv:1607.01711 [hep-ph] .
  * Mäntysaari and Schenke (2017) H. Mäntysaari and B. Schenke, Phys. Lett. B772, 832 (2017), arXiv:1703.09256 [hep-ph] .
  * Cepila  _et al._ (2017) J. Cepila, J. G. Contreras, and J. D. Tapia Takaki, Phys. Lett. B766, 186 (2017), arXiv:1608.07559 [hep-ph] .
  * Schlichting and Schenke (2014) S. Schlichting and B. Schenke, Phys. Lett. B739, 313 (2014), arXiv:1407.8458 [hep-ph] .
  * Armesto (2002) N. Armesto, Eur. Phys. J. C26, 35 (2002), arXiv:hep-ph/0206017 [hep-ph] .
  * Armesto  _et al._ (2005) N. Armesto, C. A. Salgado, and U. A. Wiedemann, Phys. Rev. Lett. 94, 022002 (2005), arXiv:hep-ph/0407018 [hep-ph] .
  * Albacete and Soto-Ontoso (2017) J. L. Albacete and A. Soto-Ontoso, Phys. Lett. B770, 149 (2017), arXiv:1605.09176 [hep-ph] .
  * Welsh  _et al._ (2016) K. Welsh, J. Singer, and U. W. Heinz, Phys. Rev. C94, 024919 (2016), arXiv:1605.09418 [nucl-th] .
  * M ntysaari  _et al._ (2017) H. M ntysaari, B. Schenke, C. Shen, and P. Tribedy, Phys. Lett. B772, 681 (2017), arXiv:1705.03177 [nucl-th] .
  * Contreras (2017) J. G. Contreras, Phys. Rev. C96, 015203 (2017), arXiv:1610.03350 [nucl-ex] .
  * Klein and Nystrand (1999) S. R. Klein and J. Nystrand, Phys. Rev. C60, 014903 (1999), arXiv:hep-ph/9902259 [hep-ph] .
  * Klusek-Gawenda and Szczurek (2016) M. Klusek-Gawenda and A. Szczurek, Phys. Rev. C93, 044912 (2016), arXiv:1509.03173 [nucl-th] .
  * Gay Ducati and Martins (2017) M. Gay Ducati and S. Martins, Phys. Rev. D96, 056014 (2017), arXiv:1705.05053 [hep-ph] .
  * Shuvaev  _et al._ (1999) A. G. Shuvaev, K. J. Golec-Biernat, A. D. Martin, and M. G. Ryskin, Phys. Rev. D60, 014015 (1999), arXiv:hep-ph/9902410 [hep-ph] .
  * Golec-Biernat and Wusthoff (1998) K. J. Golec-Biernat and M. Wusthoff, Phys. Rev. D59, 014017 (1998), arXiv:hep-ph/9807513 [hep-ph] .
  * Nemchik  _et al._ (1994) J. Nemchik, N. N. Nikolaev, and B. G. Zakharov, Phys. Lett. B341, 228 (1994), arXiv:hep-ph/9405355 [hep-ph] .
  * Nemchik  _et al._ (1997) J. Nemchik, N. N. Nikolaev, E. Predazzi, and B. G. Zakharov, Z. Phys. C75, 71 (1997), arXiv:hep-ph/9605231 [hep-ph] .
  * Kowalski  _et al._ (2006) H. Kowalski, L. Motyka, and G. Watt, Phys. Rev. D74, 074016 (2006), arXiv:hep-ph/0606272 [hep-ph] .
  * Hahn (2005) T. Hahn, Comput. Phys. Commun. 168, 78 (2005), arXiv:hep-ph/0404043 [hep-ph] .
  * Toll and Ullrich (2013) T. Toll and T. Ullrich, Phys. Rev. C87, 024913 (2013), arXiv:1211.3048 [hep-ph] .
  * Abelev  _et al._ (2013a) B. Abelev _et al._ (ALICE), Phys. Lett. B718, 1273 (2013a), arXiv:1209.3715 [nucl-ex] .
  * Abbas  _et al._ (2013) E. Abbas _et al._ (ALICE), Eur. Phys. J. C73, 2617 (2013), arXiv:1305.1467 [nucl-ex] .
  * Adam  _et al._ (2016) J. Adam _et al._ (ALICE), Phys. Rev. Lett. 116, 222301 (2016), arXiv:1509.08802 [nucl-ex] .
  * Afanasiev  _et al._ (2009) S. Afanasiev _et al._ (PHENIX), Phys. Lett. B679, 321 (2009), arXiv:0903.2041 [nucl-ex] .
  * Khachatryan  _et al._ (2017) V. Khachatryan _et al._ (CMS), Phys. Lett. B772, 489 (2017), arXiv:1605.06966 [nucl-ex] .
  * Abelev  _et al._ (2013b) B. Abelev _et al._ (ALICE), Phys. Rev. C88, 044909 (2013b), arXiv:1301.4361 [nucl-ex] .


# Two-component mass models of the lensing galaxy in the quadruply imaged
supernova iPTF16geu

Liliya L.R. Williams  School of Physics and Astronomy, University of
Minnesota, 116 Church Street, Minneapolis, MN 55455, USA<EMAIL_ADDRESS>(LLRW)
David Zegeye  Department of Astronomy and Astrophysics, University of Chicago,
5640 South Ellis Avenue, Chicago, IL 60637, USA

###### Abstract

The first resolved, multiply imaged supernova Type Ia, iPTF16geu, was observed
4 years ago, five decades after such systems were first envisioned. Because of
the unique properties of the source, these systems hold a lot of promise for
the study of galaxy structure and cosmological parameters. However, this very
first example presented modelers with a few puzzles. It was expected that to
explain image fluxes a contribution from microlensing by stars would be
required, but to accommodate the magnitude of microlensing, the density slope
of the elliptical power law lens model had to be quite shallow,
$\rho_{2D}\propto r^{-0.7}$. Furthermore, the center of mass had to be
displaced from that of observed light by $\sim 0.1$kpc, and the position angle
of light distribution was misaligned with that of mass by $\sim 40^{o}$. In
this paper we present mass models that resolve the first two problems, and
suggest a resolution of the third. Motivated by observations of local
ellipticals, and some recent analysis of galaxy-scale lenses, our mass models
consist of two offset (baryonic) mass components. The resulting mass
distributions have a single centroid, but are lopsided, and have isodensity
contours that are not purely elliptical and not self-similar with radius. For
many of our models the microlensing requirements are modest, and the ring
formed by the extended supernova host galaxy resembles the observed one.

###### Subject headings:

gravitational lensing: strong – dark matter – galaxies:

††slugcomment: Submitted to The Open Journal of Astrophysics

##  1\. Introduction

The pioneering work of Refsdal (1964) described a new method to measure the
Hubble constant, using a phenomenon—gravitational lensing—that would not be
observed for another 15 years (Walsh et al., 1979). Though the lensed source
was envisioned to be a supernova, multiply imaged supernova were not detected
until only a few years ago (Quimby et al., 2014; Kelly et al., 2015; Goobar et
al., 2017).

There is considerable interest in multiply imaged supernova, and especially of
Type Ia, because of the wealth of information that can be extracted from them
(e.g., Oguri, 2019; Wagner & Meyer, 2019), including mass distribution in
galaxies (Oguri & Kawano, 2003; Jönsson et al., 2010), global cosmology (Qi et
al., 2019), and constraints on the nature of dark matter (Zumalacárregui &
Seljak, 2018; Pandolfi et al., 2014). Systems with lensed supernova present
some challenges as well, for example, finding these in surveys (Wojtak et al.,
2019; Porciani & Madau, 2000), and accounting for microlensing of the multiple
images (Dobler & Keeton, 2006; Yahalomi et al., 2017; Foxley-Marrable et al.,
2018; Pierel & Rodney, 2019; Goldstein et al., 2018).

Supernova iPTF16geu was the first Type Ia supernova with resolved multiple
images (Goobar et al., 2017). The source and lens are at redshifts $z=0.409$
and $z=0.2163$, respectively. Using supernova lightcurve fitting, Dhawan et
al. (2020) estimate that the total magnification is $\mu=67.8\pm 2.7$ at 68%
confidence level. Microlensing by stars in the lensing galaxy was recognized
to be important in this system since it was first detected (Foxley-Marrable et
al., 2018; More et al., 2017). However, it was soon realized that the
combination of proposed macro-lensing (i.e., due to the main lensing galaxy)
and microlensing models are not able to fully explain the system.
Specifically, assuming that the galaxy lens is described by a power law
projected density profile with an “isothermal" slope, $\Sigma\propto r^{-1}$,
would require a considerable, even unrealistic amount of microlensing
(de)magnification of the images (Yahalomi et al., 2017). Varying the fraction
of dark matter compared to total mass, or including external shear does not
help (Yahalomi et al., 2017).

This conclusion suggests that an isothermal power law lens is unlikely.
Mörtsell et al. (2019) considered a range of density profile slopes,
$\Sigma\propto r^{\alpha-2}$, and conclude that a shallower galaxy lens will
result in predicted macro-magnifications that will not put too much strain on
microlensing contribution. The best fitting slope from microlensing
considerations is $\alpha\approx 1.3$, significantly shallower than isothermal
($\alpha=1$). Such shallow slopes are far from SLACS findings (Koopmans et
al., 2006), of $\alpha=0.99\pm 0.12$, but in the range of profile slopes found
in Read et al. (2007) for lensing galaxies.

While the central regions of some galaxies may have shallow density profiles,
the profiles must steepen at larger radii. In addition to this, a further
feature of their model suggests that it might be a simplification of the
actual lens: the center of observed light is displaced from the center of
recovered mass by $\sim 0.03^{\prime\prime}$, or about $0.1$kpc, much larger
than the astrometric uncertainties of $0.002^{\prime\prime}$. Finally, the
position angle of the fitted elliptical mass distribution is misaligned with
the observed light orientation by $\sim 40^{\circ}$.

Taken together, these properties of the lens model lead us to conclude that
the actual mass distribution is more complicated than an elliptical density
power law. What form could these complications take? Addition of dark matter
substructure, as predicted by $\Lambda$CDM cosmological model (Klypin et al.,
1999; Moore et al., 1999) is unlikely to reconcile the simple mass model with
observations because substructure cannot change profile slope, or affect image
positions significantly enough to reconcile light and mass centroids.
Microlensing by stars can significantly affect image fluxes, but not image
positions. The additional complications must be due to the macro-scale mass
model of the galaxy lens.

Wagner (2018) has shown that images provide only local information about the
lensing potential; the rest is usually filled in by model assumptions. This
motivates us to explore other possibilities for what projected lens looks like
in the case of iPTF16geu. Specifically, we extend the published lens models
beyond a single power law; our models consist of two superimposed, cored power
law mass components, each with its own set of parameters (§ 2). Our recovered
mass distributions, generated using the method described in § 3, account for
image positions, and produce coincident light and mass centroids within
astrometric error (§ 4). We do not use time delays, because of the large
uncertainties, or the extended image of the supernova host galaxy, which was
found by Mörtsell et al. (2019) not to be very useful. We do compare extended
rings and time delays generated by our models with observations, in § 4.3 and
§ 4.4.

##  2\. Galaxy model

To motivate our mass model we assume that the lens galaxy at $z\approx 0.2$
must be analogous to nearby ellipticals of similar dynamical mass and
effective radius of the light distribution, $R_{e}$. Its velocity dispersion
is 129 km/s (Johansson et al., 2020)111An earlier estimate by Mörtsell et al.
(2019) was 164 km/s., and Dhawan et al. (2020) measure $R_{e}\approx
0.55^{\prime\prime}$ (depending on the filter), or 1.83 kpc, at the redshift
of the lens.

Virgo cluster ellipticals are well studied, and some are analogues to the
iPTF16geu lensing galaxy. For example, NGC 4494, has $\sigma\sim 150$ km/s
(Romanowsky et al., 2003), and $R_{e}\sim 3.8$ kpc (Foster et al., 2011), and
NGC 4434 has $\sigma\sim 100$ km/s (Dressler, 1984), and $R_{e}\sim 1.5$ kpc
(Kormendy et al., 2009; Krajnović et al., 2018). Central regions of galaxies
inside the light’s effective radius are dominated by stars. This is especially
true of lower mass ellipticals (Ferreras et al., 2005). In iPTF16geu, quad
images are formed at $\sim 1$kpc, or $\sim 0.55R_{e}$, where the bulk of the
mass is in the form of stars. Therefore, to motivate our mass model we need to
look at the stellar mass distribution in low mass ellipticals.

Dhar & Williams (2012) used data from Kormendy et al. (2009) to fit sky
projected radial light distribution of Virgo ellipticals over $4-5$ decades in
radius. A single component, either a power law or a Sersic profile does not
provide a good fit, while fits using 2 or 3 projected Einasto components have
very small residuals, with rms of $0.025-0.065$ magnitudes for NGC 4494 and
4434. When fitted with 2 components, the ratio of the two scale-lengths is
about $10-25$, for both of these galaxies. If these Virgo galaxies were the
lensing galaxies in iPTF16geu, images would form somewhat interior to the
scale-length of the larger component, at the radius where the total light
distribution undergoes the most rapid change in (projected) slope.

Our mass models represent the lensing galaxy as a superposition of two mass
components, but instead of using projected Einasto profiles, which do not have
analytical expressions for lensing potential and deflection angles, we use
alphapot analytical lensing potential (Keeton, 2011), whose expressions for
deflection angles and projected density, or convergence are presented, for
example, in Ghosh et al. (2020). Each alphapot is characterized by a center,
normalization $b$, potential profile slope $\alpha$, core radius $s$, and
ellipticity parameters $q$, and $K$:

$\Psi=b(s^{2}+x^{2}+y^{2}/q^{2}+K^{2}xy)^{\alpha/2}$ (1)

If $K=0$, the elongation of the iso-potential contours is along the $x$ or $y$
axis, and $q$ gives the axis ratio. If $K\neq 0$, the ellipticity position
angle is no longer aligned with the principal axes. Because of adjustable core
size and power law density slope, this profile captures the most important
feature of Einasto profile over the radial range relevant to lensing, namely
the ability of have constant or changing density slope.

We allow the centers of the two mass components to be offset from each other,
which amounts to having a dipole moment in the mass distribution around the
center of light. Gomer & Williams (2018) showed that to reproduce the
statistical properties of image distribution around the lens center of the
sample of $\sim 40$ galaxy-scale quad lenses, one needs to introduce small,
dipolar-like azimuthal mass perturbations around the lens center. The authors
implemented these by offsetting the centers of the two mass components, while
the total mass distribution still had a single centroid. Further evidence for
such complexity in mass distribution comes from the modeling presented in
Nightingale et al. (2019): all three of their massive elliptical lenses
consist of two baryonic components, and show lopsidedness in their mass
distribution.

The more compact of the two mass components in our models is designated as the
main one, while the secondary is more extended. The two need not have equal
mass, or brightness. In fact, the central brightness of the secondary
component tends to be fainter by about a magnitude (Dhar & Williams, 2012). In
principle, the second component in our models could also represent dark
matter, with considerable or minimal stellar contribution. For lensing, all
that matters is their mass content.

Figure 1.— Image magnifications, in magnitudes. Left panel: Minima point
images: 1st and 2nd arriving (labelled $\\#2$ and $\\#4$ in previous works).
Right panel: Saddle point images: 3rd and 4th arriving (labelled $\\#1$ and
$\\#3$ in previous works). In both panels, the observed magnifications of
iPTF16geu images (Dhawan et al., 2020) are denoted by a magenta square. Empty
black symbols represent published models: circles are More et al. (2017), and
triangles are Mörtsell et al. (2019) (labeled by the density profile slope).
Orange points are our 458 models, all satisfying $\chi^{2}<1$ for image
positions. Green points are a subset of these, chosen based on criteria
described at the end of § 3.

##  3\. Model fitting

The positions of the 4 images of the quad give us 8 data constraints. We do
not use the ring formed by the extended light from the supernova host galaxy,
as Mörtsell et al. (2019) found it to be not very useful; the ring did not
improve the fit or change if from the one found using point images alone. We
do not use time delays either, because they are expected to be short in this
cross-type quad, and because of their rather large uncertainties. Later in
this Section we present $\chi^{2}$ values calculated without the time delay
contribution, and in § 4.4 we include time delay information, and show that in
this system time delays do not help discrimanate between models. We present
images of extended host galaxy for some of our final models in § 4.3 and time
delays in § 4.4.

Our two mass components, $A$ and $B$, have 5 free parameters each, $b$, $s$,
$q$, $K$, and $\alpha$, defined in eq. (1). The offset between the two,
$\Delta x,\Delta y$, adds 2 more free parameters, and 2 additional ones are
due to the unknown location of the source. From these 14, we eliminate 2 that
correspond to the mass normalization of the lens and its orientation
(rotation) on the plane of the sky. In practice that means we set $K_{A}=0$,
and only the ratio $b_{B}/b_{A}$ matters, not their individual values. In
other words, we find fits to the quad characterized by 3 image distance ratios
with respect to the lens center, and 3 relative polar image angles. Having
found a fit, the rotation of the whole system with respect to the center of
light, and the mass normalization are adjusted to correspond to the images of
iPTF16geu. Specifically, the model quad system is shrunk or expanded, and
rotated, such that the model’s first arriving image is exactly at the distance
and position angle of the observed first arriving image.

$\Delta x,y$ | $b_{B}/b_{A}$ | $\alpha_{A}$ | $q_{A}^{2}$ | $K_{A}^{2}$ | $\alpha_{B}$ | $q_{B}^{2}$ | $K_{B}^{2}$
---|---|---|---|---|---|---|---
$\pm$0.075” | 1-100 | 0.6–1.4 | 0.7–1.0 | 0 | 0.6–1.4 | 0.7–1.3 | -0.4–0.4
Table 1Ranges of parameter values used when randomly searching the 12
dimensional parameter space of the 2 component mass model.

This means that our merit function is

$\chi_{p}^{2}=\frac{1}{6}\,\Big{(}\sum_{i=2}^{4}\,\Big{[}\frac{(d_{i,m}-d_{i,o})^{2}}{(\sigma_{d,i})^{2}}\Big{]}+\sum_{i=2}^{4}\,\Big{[}\frac{(\theta_{i,m}-\theta_{i,o})^{2}}{(\sigma_{\theta,i})^{2}}\Big{]}\Big{)},$
(2)

where $d_{i}$ is the ratio of the distance of image $i$ from the lens center
to the distance of the first arriving image, and similarly for image angle
$\theta_{i}$. Subscripts $m$ and $o$ refer to the model and observed image
properties, respectively. Uncertainties $\sigma$’s were obtained from the
astrometric uncertainties in Mörtsell et al. (2019). Subscript $p$ on
$\chi^{2}$ means that it is based on positional data only. In § 4.4 we
consider time delays as well, which modifies the merit function to

$\chi_{p,t}^{2}=\frac{1}{9}\,\Big{(}6\chi_{p}^{2}+\sum_{j=2}^{4}\,\Big{[}\frac{(t_{j,m}-t_{j,o})^{2}}{(\sigma_{t,j})^{2}}\Big{]}\Big{)},$
(3)

where the time delays given in Mörtsell et al. (2019) are with respect to the
third arriving image.

We search over a ranges of parameters, given in Table 1. Ellipticity
parameters were not allowed to be too far from circular, in keeping with the
appearance of the light of the lensing galaxy. Additionally, the core radii of
the two components, $s_{A}$ and $s_{B}$, were constrained to be smaller than
the image radius. The ratio of normalization parameters, $b_{B}/b_{A}$ was
picked from a wide range, while the offset of the second component was
confined to about a quarter of the Einstein radius. While these starting
parameter ranges are somewhat arbitrary, we note that the second stage of our
two-step modeling procedure (see below) allows parameters to wonder well
outside of the ranges set in this first step.

We use the following method to search the 12 dimensional parameter space for
solutions. Since the model parameters outnumber lensing constraints, we do not
expect a unique best fit solution. Instead, many solutions will satisfy the
constraints of iPTF16geu. Our modeling consists of two stages. First, we
randomly search a wide range of parameters, looking for approximate fits to
the lens system. It is not hard to find many fits with
$\chi_{p}^{2}\hbox{\,\hbox to0.0pt{\raise 1.82986pt\hbox{$<$}\hss}\lower
2.79857pt\hbox{$\sim$}\,}9$ using a random search. We generated 300 of these.
As expected, the distribution of $\chi_{p}^{2}$ values in this sample is
heavily biased towards high values. Then, we use randomly selected pairs of
these bad solutions to define pairs of points in the 12 dimensional parameter
space. The two solutions in a pair serve as “opposite” vertices of a simplex.
In 2D such a simplex will be a square, in 3D, a cube, etc. Each simplex is
input into downhill simplex method to find better solutions. Of the 2000
downhill simplex runs, 458, or about $23\%$ had $\chi_{p}^{2}<1$. We checked
that these solutions are in fact separate local minima in the 12 dimensional
parameter space, because a straight line between any two contains worse fits.
Our method, which combines a random search followed by downhill simplex
search, is an easily implemented and effective way of finding solutions in
multidimensional parameter space that has many local minima. The two steps of
our search method resemble those of the particle swarm optimization, though
the latter looks for a single solution, while we look for many. Note that
there is no need to find the global minimum; all local minima have
$\chi_{p}^{2}<1$, and many have $\chi_{p}^{2}<0.1$.

Because our mass models contain more parameters than data constraints, we need
to check that these extra parameters are statistically justified. To do that,
we compare our models to commonly used, single mass component models. To allow
for better fits we include external shear, which results in $k=5$ total free
parameters. In our two-component models $k=12$ (see above in this section).
When the center of the one-component lens is fixed at the center of the
observed light, as it is in the two-component models, a typical
$\chi_{p}^{2}\sim 27$. Both models have $n=6$ data points (three relative
polar image angles and distance ratios, with respect to the lens center).
Using Bayesian Information Criterion, BIC $=\chi^{2}+k\ln n$, we get $\approx
36$ for one-component models, and $\approx 22.5$ for our two-component models.
BIC differences of $>10$ are considered very strong evidence for models with
smaller BIC values. Thus, two-component models are well justified in this
system.

Each solution is a lens model of iPTF16geu, and has associated image
magnifications. The two panels of Figure 1 show magnifications of the two
images formed at the minima of the arrival time surface (left panel), and two
images formed at the saddle points (right panel). We designate the images in
two ways: the published papers label them as #1, …#4, and we also add
designation based on the arrival time, as 1st, …, 4th (Saha & Williams, 2003).

We show these 458 solutions as small orange filled points in these plots. (We
checked that the distribution of solutions does not depend on their
$\chi_{p}^{2}$ values, so a somewhat different value for the $\chi_{p}^{2}$
cut would have generated a statistically similar distribution of solutions.)
Most of them lie near a straight line, running approximately diagonally from
low to high magnifications in each of the panels. The models from published
literature are denoted by empty black triangles (Mörtsell et al., 2019, Table
9) and circles (More et al., 2017, Table 1), and also lie approximately along
the same diagonal lines. Solutions along these two lines share similar
geometry of their mass distributions, specifically, they are connected by
approximate mass sheet degeneracy (MSD), which is more aptly called the
steepness degeneracy (Gorenstein et al., 1988; Saha, 2000); flatter density
profiles correspond to higher image magnifications. Note that MSD here is not
exact: rescaling of $\kappa$ of one solution by some $\lambda$ and adding a
flat mass sheet of density $1-\lambda$ (in units of critical lensing density)
will generate a solution which is only approximately similar to the ones
plotted here along these diagonals. The actual relation between these
solutions (both orange and empty black symbols) are probably better described
by non-trivial source plane transformations (SPT; Schneider & Sluse, 2014).

Not all of our 458 solutions are equally good representations of galaxies.
Because models consist of two offset mass components, many total mass models
will have two separated centroids. These are probably un-astrophysical, and we
do not consider them further. In the remaining $\sim 30\%$ of the cases the
offset centers merely result in lopsidedness, or dipole-like feature in the
mass distribution. The total mass distribution in these has a single centroid;
in other words, the center of light and mass coincide. Aiming to isolate
models that are closer to single component ones, we further select solutions
where the main (compact) component contains at least half of the total mass
within the image ring. We end up with a subset of 57 solutions, plotted as
green dots in Figure 1. By selecting these 57, we do not preferentially select
better or worse fits: we checked that $\chi_{p}^{2}$ does not correlate with
the fraction of mass in the main component, or how far the second centroid is
from the light center.

Though our imposed cuts do not correlate with the goodness of fit, they do
correlate with other lens properties. Most of these 57 solutions deviate from
the diagonal lines, especially in the right panel, and many of these lie
closer to iPTF16geu, than further from it, meaning that these solutions
require less microlensing (de)magnification.

We cannot examine all 57 models in detail. Instead, we pick a small handful of
examples, based on their magnification properties, as these are important for
microlensing. Models M1 and M3 are on opposite sides of the main diagonal in
Figure 1, and are closer to the iPTF16geu than other green points within their
respective islands of solutions. Additionally, we selected M3 because it is
close to iPTF16geu in both panels of Figure 1. Thus M3 would require minimal
contribution from microlensing. M2 was selected because it is on the main
diagonal. We note that we could have picked a different set of models, or
models based on physical criteria other than magnification. But the wide range
of models allowed by the data requires us to make choices.

##  4\. Results: examining models M1-M3

###  4.1. Mass distributions

In this section we examine the three solutions selected in the previous
section, and depicted as large green filled symbols in Figure 1. They span a
range of magnifications, and have coincident mass and light centroids. We
present their mass distributions in Figure 2. The magenta is the $\kappa=1$
contour of the main (compact) mass component, while the green curves represent
the isodensity contours of the total mass distribution, logarithmically spaced
with an interval of 0.05. The thick green contour denotes $\kappa=1$. The
green isodensity contours are not exactly elliptical because they represent a
sum of two offset mass components. The fact that the offset results in a
dipole-like mass distribution around the lens center is evident in these
plots.

Figure 2.— Three mass models M1-M3, shown in Figure 1. In all 3, the center of
light is coincident with the center of mass, even though the mass contours are
not purely elliptical. The thick green and magenta contours represent
$\kappa=1$ levels of the total mass and the main (compact) mass component,
respectively. The green thin contours of the total mass are spaced by 0.05 in
the log. Images are labeled as in Mörtsell et al. (2019), and the arrival
sequence is given in parentheses. Thin blue line is the orientation of the
elliptical mass models in Mörtsell et al. (2019), while the two thicker blue
lines give the range of the position angle of the observed light. (North is up
and East to the left.)

Figure 3.— Radial density profiles of models M1-M3, shown in Figure 1. Solid
and short-dashed blue lines are profiles taken along the positive and negative
directions of the major axis, while the long-dashed line is along the positive
direction of the axis perpendicular to the major axis. The locations of the 4
images are marked with short vertical magenta lines. Two reference slopes, at
$\alpha=1.0$ and $\alpha=1.4$ are also indicated.

Figure 4.— Images of extended host galaxy of iPTF16geu, lensed by models
M1-M3. Images of supernova are shown as black dots, and their IDs are the same
as in Figure 2. The center of the host galaxy and iPTF16geu are displaced by
varying amounts, $<0.014^{\prime\prime}$.

Published works agree on the position angle of the recovered mass
distribution, depicted with the thin blue line going through the lens center,
pointing at image #3 to within a few degrees. This elongation direction can be
understood even without modeling, from the consideration of image
distribution: image #3 is the 4th arriving, and is usually expected to be
closest to the lens center. Since this is not the case here, there has to be
mass extention in that direction, to “push” the last arriving saddle further
from center. Not surprisingly, our models also show elongation along that
axis, though because the total mass distribution is not elliptical, the
position angle depends somewhat on the radius.

The axis of the observed light distribution is within the cone outlined by two
blue thicker lines, going through the center of the lens in Figure 2. It is
not aligned with the elongation of our mass distributions, echoing the
conclusion of Mörtsell et al. (2019). We interpret this misalignment as
further indication that the mass, and light, distribution in the lensing
galaxy cannot be represented by simple mass models. Given that the reason for
the orientation of the mass ellipticity (see the preceding paragraph) is
largely model-independent, reconciling it with the position angle of the light
could require radial variations in the stellar mass-to-light ratio, or local
elongations in the light distribution, such as the ones in the green contours
in the southern direction of model M3 in Figure 2.

The density profiles along three radial spokes emanating from the lens center
are shown in Figure 3. The solid and short-dashed blue lines are profiles
taken along the positive and negative directions of the major axis, while the
long-dashed line is along the positive direction of the axis perpendicular to
the major axis. The major axis was determined as the axis of the largest first
moment of the total mass distribution. Two reference slopes, of $\alpha=1.0$
and $1.4$ are shown for comparison. Image locations are depicted by vertical
magenta line segments. At the location of the images the profile slopes range
between $\alpha\sim 1.0$ and $\sim 1.4$, with the average being somewhat
steeper than $\alpha=1.3$ slope obtained in Mörtsell et al. (2019).

The fact that our mass maps show a variety of isodensity contour shapes means
that they are not related by MSD, but instead, by more complex degeneracies
that involve transformations of shape (Saha & Williams, 2006; Read et al.,
2007), like source plane transformations (Schneider & Sluse, 2014).

###  4.2. Microlensing

Microlensing (de)magnifications implied by our mass models can be deduced from
Figure 1, as the displacement of the magenta point, labeled iPTF16geu, from
the three large green dots, labeled M1-M3. One of the features of our models
that make them different from the published ones is that we predict solutions
whose macro magnifications do not lie on the same diagonal relation on that
figure, as those of ellipsoidal power law density profiles. Most
interestingly, there is an island of solutions in the right panel of Figure 1
that is closer to iPTF16geu than ellipsoidal power laws; our model M3 is an
example. Of the 3 models we consider, changes in flux due to microlensing are
smallest for model M3: $\hbox{\,\hbox to0.0pt{\raise
1.82986pt\hbox{$<$}\hss}\lower 2.79857pt\hbox{$\sim$}\,}0.4$ magnitudes for 3
images, and $\sim 0.75$ for the brightest, 1st arriving image. Model M1 would
require largest microlensing contribution, with the 2nd arriving image
requiring a demagnification of $\sim 1.1$ magnitudes.

Though our models are different from the elliptical power law models, we agree
with the latter that image fluxes cannot be explained without microlensing.
Microlensing is almost definitely present in this system, and most likely
affects the fluxes of at least three, or all four images. In that sense we
agree with the conclusion of Yahalomi et al. (2017) that microlensing makes it
hard to accurately determine the standard candle brightness of the supernova,
and therefore use it to help break the mass sheet degeneracy in iPTF16geu.

###  4.3. Ring image of the extended galaxy host

Extended images of host galaxies, i.e., rings, have been used extensively in
the literature for nearly three decades (e.g., Kochanek et al., 1989, 2001);
most recently to help constrain $H_{0}$ (Suyu et al., 2017), and to uncover
mass substructure in the lens galaxy or along the line of sight (e.g., Vegetti
et al., 2012). However, just like point images, they are subject to
degeneracies (Saha & Williams, 2001; Walls & Williams, 2018), and so may be of
limited value in some circumstances. The example in Figure 1 of Denzel et al.
(2020), which shows two completely different lens mass distributions producing
indistinguishable rings, is a sobering reminder that we do not yet have a full
understanding of how lensing degeneracies operate.

Rings generated from our M1-M3 mass models are shown in Figure 4. The host
galaxy source in all 3 cases was assumed to be circular, with a Gaussian light
profile. The center of the host galaxy is offset from the supernova by various
amounts, $<0.014^{\prime\prime}$, in these 3 models, comparable to the
estimate presented in Mörtsell et al. (2019), $\sim 0.02^{\prime\prime}$. The
offset was chosen to better resemble the observed ring, but we did not fine-
tune the offset (or other parameters of the source) to produce as a close a
match to the ring as possible. Given the appearance of the rings presented in
Dhawan et al. (2020), the model that appears closest to observations is M3,
because its brightest ring segment is between images $\\#1$ and $\\#4$. In M1
and M2, on the other hand, the segment between $\\#3$ and $\\#4$ is equally
bright, in contrast to the observations. This suggests that of the 3 models,
M3 is the best representation of the lensing galaxy.

###  4.4. Time Delays

We did not consider time delay information when generating our mass models in
§ 3. Had we included observed time delays the model selection would not have
been affected. To check that, we calculate $\chi^{2}$ for all our 458 models,
but now using both versions, eq. 2 and eq. 3. The results are shown in Figure
5. Because of their large uncertainties, time delays add an approximately
constant contribution to each model, and therefore have little discriminating
power compared to image positions.

Figure 6 shows the distributions of time delays predicted by our lens mass
models, orange points in Figure 1. The three histograms are labeled by the
corresponding time delay, using the same image naming convention as in
Mörtsell et al. (2019). The three sets of three cross symbols show the values
from models M1, M2 and M3. The horizontal bars at the top of the plot show
Mörtsell et al. (2019) values, for lens potential profile slopes $\alpha=1$
and $\alpha=1.4$, as labeled. The thin lines show the range of values
resulting from the continuous range of $\alpha$’s. While the time delay ranges
predicted by our models and these models agree, both are somewhat at odds with
observations for $\tau_{\\#1\\#4}$ and $\tau_{\\#1\\#3}$, but only at
$\hbox{\,\hbox to0.0pt{\raise 1.82986pt\hbox{$<$}\hss}\lower
2.79857pt\hbox{$\sim$}\,}1.5\sigma$, given the uncertainties of $\sim 1$ day.

Figure 5.— Distribution of $\chi^{2}$ values for our 458 models. The symbols
are the same as in Figure 1. The horizontal and vertical axes contain
$\chi^{2}$ evaluated based on image positions, eq. 2, and image positions
combined with time delays, eq. 3, respectively. Introducing time delay
information does not help constrain models.  Figure 6.— The three histograms
show the distributions of time delays (as labeled) of our models represented
by orange points in Figure 1. The three sets of cross symbols at the bottom
show the predictions of the 3 models M1, M2 and M3. The top of the plot shows
model predictions of Mörtsell et al. (2019). Thick line segments represent
$1\sigma$ ranges (their Table 11) for models with potential slopes $\alpha$,
as labeled. The thin horizontal lines connect time delays for the whole range
of slopes $\alpha=1\rightarrow 1.4$. The $1\sigma$ uncertainties in the
observed time delays are represented by dashed horizontal lines. Note that the
central values for $\tau_{\\#1\\#4}$ and $\tau_{\\#1\\#3}$ are outside the
limits of the plot.

##  5\. Conclusions

Multiply imaged supernova, and especially Type Ia’s, hold a lot of promise for
the study of galaxies and cosmology. They are superior to quasars as sources
because their time delays are relatively easy to measure, their absolute
fluxes are known, and their light eventually fades, allowing observations of
the lens galaxy without the glare of the bright images. All these reasons make
the very first multiply imaged supernova, iPTF16geu, very exciting.

As is often the case, the first of its kind object had some surprizes for
modelers. If the mass model is assumed to be an elliptical power law, the mass
and light centers are displaced by $\sim 0.1$kpc, which is unlikely to be the
case in an equilibrium galaxy. Furthermore, to accommodate microlensing
constraints, the power law density slope has to be rather shallow, and so
cannot hold at larger radii.

The main motivation for the present paper was to generate lens models whose
mass distribution resembles that of local ellipticals, and where the center of
observed light coincides with the center of recovered mass. Similar to local
ellipticals whose light is represented by more than one component, our mass
models consist of two mass components. The two components are offset,
reflecting conclusions of some recent lensing reconstructions. Our final
selected models have a single mass peak, and hence coincident mass and light
centers. The consequence of using two offset components was to produce dipole-
like lopsidedness in the total mass distribution. Our modeling assumptions,
namely two offset elliptical potentials, result in a range of mass models—not
a single model—all satisfying positional and time delay image information. For
illustration purposes we concentrated on three of these models, based on their
magnification properties, but could have chosen a different set.

Though our mass distributions are not elliptical power laws, we do agree with
the published models on some key features of the lens: (1) To explain the
fluxes of the iPTF16geu, microlensing by stars in the lens is required for at
least 3, or all 4 images of the quad; (2) The total mass density profile at
the location of the images is significantly shallower than isothermal; (3) The
position angle of galaxy mass elongation is misaligned with the light, which
is probably an indication that the mass, and light, distribution in the lens
is not simple.

We note that our two-component models may not be the only way to account for
the observations of iPTF16geu. It is possible that the main lensing galaxy has
a simple one-component structure, but there are additional significant
secondary galaxies, either at the same redshift or along the line of sight.
However, there are no observed nearby galaxies, making this scenario less
likely. The line-of-sight structure (LoS) is also unlikely because the source
is very nearby, $z=0.409$, so the line of sight dimension is short. Despali et
al. (2018) estimate that the typical number of LoS subhalos along a given
direction is much less than 1, and their mass is a small fraction of the mass
enclosed by the images.

## Acknowledgments

The authors would like to thank Edvard Mörtsell for kindly providing addition
information about their modeling, and the referee for helping us improve the
paper. This work was supported by the NSF through grant #1757388.

## References

  * Denzel et al. (2020) Denzel P., Mukherjee S., Coles J. P., Saha P., 2020, MNRAS, 492, 3, 3885
  * Despali et al. (2018) Despali G., Vegetti S., White S. D. M., Giocoli C., van den Bosch F. C., 2018, MNRAS, 475, 4, 5424
  * Dhar & Williams (2012) Dhar B. K., Williams L. L. R., 2012, MNRAS, 427, 1, 204
  * Dhawan et al. (2020) Dhawan S., Johansson J., Goobar A., et al., 2020, MNRAS, 491, 2, 2639
  * Dobler & Keeton (2006) Dobler G., Keeton C. R., 2006, ApJ, 653, 2, 1391
  * Dressler (1984) Dressler A., 1984, ApJ, 281, 512
  * Ferreras et al. (2005) Ferreras I., Saha P., Williams L. L. R., 2005, ApJ, 623, 1, L5
  * Foster et al. (2011) Foster C., Spitler L. R., Romanowsky A. J., et al., 2011, MNRAS, 415, 4, 3393
  * Foxley-Marrable et al. (2018) Foxley-Marrable M., Collett T. E., Vernardos G., Goldstein D. A., Bacon D., 2018, MNRAS, 478, 4, 5081
  * Ghosh et al. (2020) Ghosh A., Williams L. L. R., Liesenborgs J., 2020, MNRAS, 494, 3, 3998
  * Goldstein et al. (2018) Goldstein D. A., Nugent P. E., Kasen D. N., Collett T. E., 2018, ApJ, 855, 1, 22
  * Gomer & Williams (2018) Gomer M. R., Williams L. L. R., 2018, MNRAS, 475, 2, 1987
  * Goobar et al. (2017) Goobar A., Amanullah R., Kulkarni S. R., et al., 2017, Science, 356, 6335, 291
  * Gorenstein et al. (1988) Gorenstein M. V., Falco E. E., Shapiro I. I., 1988, ApJ, 327, 693
  * Johansson et al. (2020) Johansson J., Goobar A., Price S. H., et al., 2020, arXiv e-prints, arXiv:2004.10164
  * Jönsson et al. (2010) Jönsson J., Sullivan M., Hook I., et al., 2010, MNRAS, 405, 1, 535
  * Keeton (2011) Keeton C. R., 2011, GRAVLENS: Computational Methods for Gravitational Lensing
  * Kelly et al. (2015) Kelly P. L., Rodney S. A., Treu T., et al., 2015, Science, 347, 6226, 1123
  * Klypin et al. (1999) Klypin A., Kravtsov A. V., Valenzuela O., Prada F., 1999, ApJ, 522, 1, 82
  * Kochanek et al. (1989) Kochanek C. S., Blandford R. D., Lawrence C. R., Narayan R., 1989, MNRAS, 238, 43
  * Kochanek et al. (2001) Kochanek C. S., Keeton C. R., McLeod B. A., 2001, ApJ, 547, 1, 50
  * Koopmans et al. (2006) Koopmans L. V. E., Treu T., Bolton A. S., Burles S., Moustakas L. A., 2006, ApJ, 649, 2, 599
  * Kormendy et al. (2009) Kormendy J., Fisher D. B., Cornell M. E., Bender R., 2009, ApJS, 182, 1, 216
  * Krajnović et al. (2018) Krajnović D., Cappellari M., McDermid R. M., et al., 2018, MNRAS, 477, 3, 3030
  * Moore et al. (1999) Moore B., Ghigna S., Governato F., et al., 1999, ApJ, 524, 1, L19
  * More et al. (2017) More A., Suyu S. H., Oguri M., More S., Lee C.-H., 2017, ApJ, 835, 2, L25
  * Mörtsell et al. (2019) Mörtsell E., Johansson J., Dhawan S., Goobar A., Amanullah R., Goldstein D. A., 2019, arXiv e-prints, arXiv:1907.06609
  * Nightingale et al. (2019) Nightingale J. W., Massey R. J., Harvey D. R., et al., 2019, MNRAS, 489, 2, 2049
  * Oguri (2019) Oguri M., 2019, Reports on Progress in Physics, 82, 12, 126901
  * Oguri & Kawano (2003) Oguri M., Kawano Y., 2003, MNRAS, 338, 4, L25
  * Pandolfi et al. (2014) Pandolfi S., Evoli C., Ferrara A., Villaescusa-Navarro F., 2014, MNRAS, 442, 1, 13
  * Pierel & Rodney (2019) Pierel J. D. R., Rodney S., 2019, ApJ, 876, 2, 107
  * Porciani & Madau (2000) Porciani C., Madau P., 2000, ApJ, 532, 2, 679
  * Qi et al. (2019) Qi J., Cao S., Biesiada M., Zheng X., Ding X., Zhu Z.-H., 2019, Phys. Rev. D, 100, 2, 023530
  * Quimby et al. (2014) Quimby R. M., Oguri M., More A., et al., 2014, Science, 344, 6182, 396
  * Read et al. (2007) Read J. I., Saha P., Macciò A. V., 2007, ApJ, 667, 2, 645
  * Refsdal (1964) Refsdal S., 1964, MNRAS, 128, 307
  * Romanowsky et al. (2003) Romanowsky A. J., Douglas N. G., Arnaboldi M., et al., 2003, Science, 301, 5640, 1696
  * Saha (2000) Saha P., 2000, AJ, 120, 4, 1654
  * Saha & Williams (2001) Saha P., Williams L. L. R., 2001, AJ, 122, 2, 585
  * Saha & Williams (2003) Saha P., Williams L. L. R., 2003, AJ, 125, 6, 2769
  * Saha & Williams (2006) Saha P., Williams L. L. R., 2006, ApJ, 653, 2, 936
  * Schneider & Sluse (2014) Schneider P., Sluse D., 2014, A&A, 564, A103
  * Suyu et al. (2017) Suyu S. H., Bonvin V., Courbin F., et al., 2017, MNRAS, 468, 3, 2590
  * Vegetti et al. (2012) Vegetti S., Lagattuta D. J., McKean J. P., Auger M. W., Fassnacht C. D., Koopmans L. V. E., 2012, Nature, 481, 7381, 341
  * Wagner (2018) Wagner J., 2018, A&A, 620, A86
  * Wagner & Meyer (2019) Wagner J., Meyer S., 2019, MNRAS, 490, 2, 1913
  * Walls & Williams (2018) Walls L. G., Williams L. L. R., 2018, MNRAS, 481, 1, 655
  * Walsh et al. (1979) Walsh D., Carswell R. F., Weymann R. J., 1979, Nature, 279, 381
  * Wojtak et al. (2019) Wojtak R., Hjorth J., Gall C., 2019, MNRAS, 487, 3, 3342
  * Yahalomi et al. (2017) Yahalomi D. A., Schechter P. L., Wambsganss J., 2017, arXiv e-prints, arXiv:1711.07919
  * Zumalacárregui & Seljak (2018) Zumalacárregui M., Seljak U., 2018, Phys. Rev. Lett., 121, 14, 141101