# Measuring Line-of-sight Distances to Haloes with Astrometric Lensing B-mode

Kaiki Taro Inoue1
1Faculty of Science and Engineering, Kindai University, Higashi-Osaka,
577-8502, Japan<EMAIL_ADDRESS>
###### Abstract

Relative astrometric shifts between multiply lensed images provide a valuable
tool to investigate haloes in the intergalactic space. In strong lens systems
in which a single lens plays the primary role in producing multiple images,
the gravitational force exerted by line-of-sight (LOS) haloes can slightly
change the relative positions of multiply lensed images produced by the
dominant lens. In such cases, a LOS halo positioned sufficiently far from the
dominant lens along the LOS can create a pattern in the reduced deflection
angle that corresponds to the B-mode (magnetic or divergence-free mode). By
measuring both the B-mode and E-mode (electric or rotation-free mode), we can
determine the LOS distance ratios, as well as the ’bare’ convergence and shear
perturbations in the absence of the dominant lens. However, scale variations
in the distance ratio lead to mass-sheet transformations in the background
lens plane, introducing some uncertainty in the distance ratio estimation.
This uncertainty can be significantly reduced by measuring the time delays
between the lensed images. Additionally, if we obtain the redshift values of
both the dominant and perturbing haloes, along with the time delays between
the multiply lensed images that are affected by the haloes, the B-mode can
break the degeneracy related to mass-sheet transformations in both the
foreground and background lens planes. Therefore, measuring the astrometric
lensing B-mode has the potential to substantially decrease the uncertainty in
determining the Hubble constant.

###### keywords:

cosmology: theory - gravitational lensing - dark matter

††pagerange: Measuring Line-of-sight Distances to Haloes with Astrometric
Lensing B-mode–LABEL:lastpage††pubyear: 0000

##  1 Introduction

Cold dark matter (CDM) models encounter challenges on scales below $\sim
100\,$kpc, particularly concerning dwarf galaxies. For instance, the observed
number count of dwarf satellite galaxies in nearby galaxies in the Local Group
falls significantly short of the predicted number of subhaloes capable of
hosting dwarf galaxies as projected by CDM models (Kauffmann et al., 1993;
Klypin et al., 1999; Moore et al., 1999). Hydrodynamic simulations that
incorporate baryonic feedback and cosmic reionization have emerged as
potential solutions to this discrepancy for the Milky Way (Wetzel et al.,
2016; Brooks et al., 2017; Fielder et al., 2019). However, it remains
uncertain whether the Milky Way represents a ”typical galaxy,” and the extent
to which this discrepancy can be explained for other galaxies remains an open
question (Nashimoto et al., 2022). It is plausible that dark matter differs
from CDM, and there may be fewer low-mass haloes on scales below $\sim
100\,$kpc capable of hosting dwarf galaxies than CDM predictions suggest.

Gravitational lensing serves as a potent tool for investigating low-mass
haloes with masses $\lesssim 10^{9}\,M_{\odot}$ in the distant universe. In
particular, the study of quadruply lensed quasar-galaxy and galaxy-galaxy
strong lens systems has proven valuable for probing low-mass haloes. This is
because the strong lensing effect induced by a foreground galaxy amplifies the
weak lensing signals of low-mass haloes, which are otherwise challenging to
detect without such enhancement. Typically, the relative positions of lensed
images can be accurately fitted using a smooth gravitational potential within
a few milliarcseconds. However, in some of these systems, the flux ratios of
the lensed images in radio or mid-infrared wavelengths deviate by 10-40% from
the predictions of the theoretical models. It has been suggested that such
anomalies in flux ratios, particularly in radio or mid-infrared emissions, may
be attributed to subhaloes residing within the host galaxy (Mao & Schneider,
1998; Metcalf & Madau, 2001; Chiba, 2002; Dalal & Kochanek, 2002; Keeton et
al., 2003; Inoue & Chiba, 2003, 2005b; Xu et al., 2009, 2010). Nevertheless,
it is worth noting that these anomalies might also find explanations in the
complex gravitational potential of the foreground galaxy (Evans & Witt, 2003;
Oguri, 2005; Gilman et al., 2017; Hsueh et al., 2017; Hsueh et al., 2018). The
interpretation gains support from discrepancies observed in the relative
astrometric shifts of lensed extended images (Treu & Koopmans, 2004; Koopmans,
2005; Vegetti & Koopmans, 2009; Vegetti et al., 2010; Chantry et al., 2010;
Vegetti et al., 2012, 2014; Inoue et al., 2016; Hezaveh et al., 2016;
Chatterjee & Koopmans, 2018; Çagan Şengül et al., 2020).

However, any small-mass haloes along the line of sight (LOS) can also
influence the flux ratios and relative positions of lensed images (Metcalf,
2005; Xu et al., 2012). Based on a semi-analytical approach, Inoue & Takahashi
(2012) argued that the primary cause of anomalies in flux ratios for quadruply
lensed quasars with an Einstein radius of approximately $1^{\prime\prime}$ is
the presence of small-mass haloes within the intergalactic space, rather than
subhaloes within the foreground galaxy. This assertion was subsequently
confirmed by Takahashi & Inoue (2014), which employed $N$-body simulations
capable of resolving dwarf galaxy-sized haloes within a cosmological context.
Building on semi-analytical approach, Inoue (2016) highlighted that for source
redshifts $z_{\textrm{s}}>1$, the cumulative effect of line-of-sight
structures, including haloes and troughs, on altering the flux ratios of
quadruply lensed images ranges from 60 to 80 percent in the CDM models. A
subsequent analysis arrived at a similar conclusion (Despali et al., 2018).

The observational evidence supporting the significant role of LOS haloes in
’substructure-lensing’ is mounting. Based on $N$-body simulations, Takahashi &
Inoue (2014) pointed out that the ’object X’, previously assumed to be a
satellite galaxy in the lensed quasar MG J0414+0534, may in fact be located
within the intergalactic space. Additionally, Sengül et al. (2022) argued that
a presumed dark perturber, assumed to be a satellite galaxy in the lensed
quasar B1938+666, actually constitutes an intergalactic halo with a mass of
approximately $2\times 10^{9}\,M_{\odot}$. Furthermore, a recent assessment of
the lensing power spectra at an angular wave number of $l=1.2\times 10^{6}$ or
roughly $9\,$kpc within the primary lens plane aligns with CDM models, where
LOS haloes exert the predominant influence on the alterations in flux and
relative positions of lensed images (Inoue et al., 2023).

To validate this claim, it becomes imperative to determine line-of-sight (LOS)
distances to perturbing dark haloes. This necessitates the measurement of
astrometric shifts between perturbed multiple images and their unperturbed
counterparts. These astrometric shifts have the potential to break the
degeneracy between subhalo mass and the LOS distance, provided that they are
resolved at the scale of the Einstein radius of the perturbing object.
Additionally, precise determination of sky positions using a dipole structure
is required (Inoue & Chiba, 2005a).

In practical scenarios, the presence of observable dipole structures or a
fifth image due to the strong lensing effects caused by intergalactic small-
mass haloes is rare. Moreover, numerous perturbers with varying masses,
residing at different redshifts and sky positions, can influence astrometric
shifts. Consequently, modeling each dark perturber individually becomes a
formidable challenge. To achieve precise modeling of astrometric shifts, we
must comprehend the coupling effect between the strong lensing exerted by a
dominant primary lens and the weak lensing introduced by a subdominant
secondary lens located outside the primary lens plane. While frameworks for
modeling ’LOS lensing’ by intergalactic haloes have been explored in previous
literature (Erdl & Schneider, 1993; Bar-Kana, 1996; McCully et al., 2014,
2017; Birrer et al., 2017; Fleury et al., 2021), none of these analyses have
fully addressed the coupling effect for general perturbations, encompassing
all multipole moments (beyond flexions). Particularly, the ambiguity of LOS
distances stemming from multi-plane mass-sheet transformations (Schneider,
2014, 2019) remains poorly understood.

In this paper, we investigate the property of B-modes (magnetic modes) in the
two-dimensional vector field of astrometric shifts as a means to determine the
distance ratio to a perturber. The concept is straightforward: if all
perturbers are confined to the primary lens plane, their reduced deflection
angles can be expressed as gradients of a scalar potential, resulting in no
generation of B-modes. However, if some perturbers are residing at a different
lens plane, the reduced deflection angles of those perturbers that are assumed
to be in the primary lens plane cannot be represented as gradients of a scalar
potential due to the coupling effect. Consequently, this generates B-modes
akin to those induced by weak lensing in the cosmic microwave background (CMB)
polarization (Zaldarriaga & Seljak, 1998; Lewis & Challinor, 2006). It should
be noted that, in our scenario, perturbers can exist in the background of the
primary lens, whereas in weak lensing of the CMB, perturbers can only reside
in the foreground of the CMB. In Section 2, we present the theoretical
framework of our method. In Section 3, we explore how multi-plane mass-sheet
transformations impact the LOS distance ratio and time delay in a double lens
system. In Section 4, we examine the astrometric Lensing B-mode and evaluate
the accuracy of distance ratio estimators using simple toy models. Finally, in
Section 5, we offer conclusions and discuss the observational feasibility of
our proposed method.

##  2 Theory

###  2.1 E/B decomposition

Helmholtz’ theorem states that a smooth three dimensional vector field
${\bm{\alpha}}({\bm{r}}),{\bm{r}}=(x_{1},x_{2},x_{3})$ that decays faster than
$|{\bm{r}}|^{-2}$ for $|{\bm{r}}|\rightarrow\infty$ can be uniquely decomposed
as a sum of a rotation free ’electric’ field ${\bm{\alpha}}^{\textrm{E}}$ and
a divergence free ’magnetic’ field ${\bm{\alpha}}^{\textrm{B}}$ as
${\bm{\alpha}}={\bm{\alpha}}^{\textrm{E}}+{\bm{\alpha}}^{\textrm{B}}$. By
constraining $\alpha_{3}$ to be a constant, we can apply the theorem to two
dimensional vectors: In terms of an ’electric’ potential $\psi^{\textrm{E}}$
and a ’magnetic’ potential $\psi^{\textrm{B}}$, a smooth two dimensional
vector field ${\bm{\alpha}}({\bm{x}}),{\bm{x}}=(x_{1},x_{2})$ can be
decomposed as

$\displaystyle\alpha_{i}$ $\displaystyle=$
$\displaystyle\psi^{\textrm{E}}_{,i}+\varepsilon_{i}^{j}\psi^{\textrm{B}}_{,j},$
(1) $\displaystyle=$
$\displaystyle\alpha^{\textrm{E}}_{i}+\alpha^{\textrm{B}}_{i},~{}~{}(i,j=1,2),$
(2)

where $\varepsilon$ is the Levi-Civita tensor and $,i$ is a partial derivative
operator in the direction $x_{i}$. For brevity, we denote the two dimensional
rotation as
$\nabla\times{\bm{\alpha}}^{\textrm{B}}\equiv\alpha^{\textrm{B}}_{2,1}-\alpha^{\textrm{B}}_{1,2}=-\psi^{\textrm{B}}_{,1,1}-\psi^{\textrm{B}}_{,2,2}$
and divergence as
$\nabla\cdot{\bm{\alpha}}^{\textrm{E}}\equiv\alpha^{\textrm{E}}_{1,1}+\alpha^{\textrm{B}}_{2,2}=\psi^{\textrm{E}}_{,1,1}+\psi^{\textrm{E}}_{,2,2}$.

In the following, we consider a lens system in which a source at a redshift
$z_{\textrm{s}}$ is lensed by a primary (dominant) lens at a redshift
$z_{\textrm{d}}$ with a reduced deflection angle $\alpha$ and a secondary lens
(perturber) at a redshift $z_{\textrm{f}}(<z_{\textrm{d}})$ or
$z_{\textrm{b}}(>z_{\textrm{d}})$ with a small reduced deflection angle
$\delta\mbox{\boldmath$a$}\,(|\delta\mbox{\boldmath$\alpha$}|\ll|\mbox{\boldmath$\alpha$}|)$.
If the perturber resides at the lens plane of the dominant lens, i.e.,
$z_{\textrm{P}}=z_{\textrm{d}}$, $\nabla\times\delta\mbox{\boldmath$\alpha$}$
is null. However, if the perturber does not reside at the lens plane, the
rotation of the effective deflection angle111The effective deflection angle is
defined as the astrometric shift due to the secondary lens with respect to the
reduced deflection angle of the dominant lens.
$\nabla\times\delta\mbox{\boldmath$\alpha$}_{\textrm{eff}}$ is not null. We
denote the angular diameter distance between a dominant lens and the source as
$D_{\textrm{sd}}$, between a foreground perturber and a dominant lens as
$D_{\textrm{df}}(z_{\textrm{f}}<z_{\textrm{d}})$, between a dominant lens and
a background perturber as $D_{\textrm{bd}}(z_{\textrm{b}}>z_{\textrm{d}})$,
between a foreground perturber and the source as $D_{\textrm{sf}}$, between a
background perturber and the source as $D_{\textrm{sb}}$.

###  2.2 Rotation and divergence

First, we consider a double lens system in which a perturber resides in front
of a dominant lens, i.e., $z_{\textrm{f}}<z_{\textrm{d}}$. Then, the angular
position $y$ of the source can be written as a function of the angular
position $x$ of the lensed image as

${\bm{y}}={\bm{x}}-\delta\mbox{\boldmath$\alpha$}(\mbox{\boldmath$x$})-\mbox{\boldmath$\alpha$}({\bm{y}}_{\textrm{d}}),$
(3)

where
${\bm{y}}_{\textrm{d}}\equiv\mbox{\boldmath$x$}-\beta\,\delta\mbox{\boldmath$\alpha$}(\mbox{\boldmath$x$}),(0<\beta<1)$
and

$\beta=\frac{D_{\textrm{df}}D_{\textrm{s}}}{D_{\textrm{d}}D_{\textrm{sf}}},$
(4)

is the LOS distance ratio parameter, which encodes the information of the LOS
distance from an observer to the perturber in the foreground with respect to
the dominant lens. If we model the system as the one in which a perturber with
an effective deflection angle
$\delta\mbox{\boldmath$\alpha$}_{\textrm{eff}}^{\textrm{f}}$ resides at the
dominant lens plane, equation (3) can be written as

$y$ $\displaystyle=$
$\displaystyle\mbox{\boldmath$x$}-{\bm{\alpha}}({\bm{x}})-\delta{\bm{\alpha}}_{\textrm{eff}}^{\textrm{f}}({\bm{x}}),$
$\displaystyle\delta{\bm{\alpha}}_{\textrm{eff}}^{\textrm{f}}({\bm{x}})$
$\displaystyle\equiv$
$\displaystyle\delta{\bm{\alpha}}({\bm{x}})+{\bm{\alpha}}({\bm{y}}_{\textrm{d}}({\bm{x}}))-{\bm{\alpha}}(\mbox{\boldmath${\bm{x}}$}).$
(5)

Because of coupling between the strong lensing by the dominant lens and the
weak lensing by the foreground perturber, the effective deflection angle
$\delta{\bm{\alpha}}_{\textrm{eff}}^{\textrm{f}}$ has a magnetic component if
and only if $\beta>0$. The rotation of
$\delta{\bm{\alpha}}_{\textrm{eff}}^{\textrm{f}}$ is

$\nabla\times\delta{\bm{\alpha}}_{\textrm{eff}}^{\textrm{f}}({\bm{x}})=2\beta(\tilde{\gamma}_{1}\delta\gamma_{2}-\tilde{\gamma}_{2}\delta\gamma_{1}),$
(6)

where $\gamma_{i}$ is the $i$-component of the shear tensor of the dominant
lens at $\tilde{{\bm{x}}}$ and $\delta\tilde{\gamma}_{i}$ is the $i$-component
of the shear tensor of the secondary lens at ${\bm{x}}$.

Second, we consider a double lens system, in which a perturber resides in the
background of a dominant lens, i.e., $z_{\textrm{b}}>z_{\textrm{d}}$. Then,
the angular position ${\bm{y}}$ of the source can be written as a function of
the angular position ${\bm{x}}$ of the lensed image as

${\bm{y}}={\bm{x}}-{\bm{\alpha}}({\bm{x}})-\delta{\bm{\alpha}}({\bm{y}}_{\textrm{b}}),$
(7)

where
${\bm{y}}_{\textrm{b}}\equiv{\bm{x}}-\xi\,{\bm{\alpha}}({\bm{x}}),(0<\xi<1)$
and

$\xi=\frac{D_{\textrm{bd}}D_{\textrm{s}}}{D_{\textrm{b}}D_{\textrm{sd}}},$ (8)

is the LOS distance ratio parameter, which encodes the information of the LOS
distance from an observer to perturber in the background of the dominant lens.

As in the foreground case, we can model the system as the one in which a
secondary perturber with an effective deflection angle
$\delta\mbox{\boldmath$\alpha$}_{\textrm{eff}}^{\textrm{b}}$ resides at the
dominant lens plane, equation (7) can be written as

$\displaystyle{\bm{y}}$ $\displaystyle=$
$\displaystyle{\bm{x}}-{\bm{\alpha}}({\bm{x}})-\delta{\bm{\alpha}}_{\textrm{eff}}^{\textrm{b}}({\bm{x}}),$
$\displaystyle\delta{\bm{\alpha}}_{\textrm{eff}}^{\textrm{b}}({\bm{x}})$
$\displaystyle\equiv$
$\displaystyle\delta{\bm{\alpha}}({\bm{y}}_{\textrm{b}}).$ (9)

Because of the difference between the angular position of the lensed image
${\bm{x}}$ and that of the background perturber ${\bm{y}}_{\textrm{b}}$,
$\delta{\bm{\alpha}}_{\textrm{eff}}^{\textrm{b}}({\bm{x}})$ has a magnetic
component if and only if $\xi>0$. The rotation of
$\delta{\bm{\alpha}}_{\textrm{eff}}^{\textrm{b}}$ is

$\nabla\times\delta{\bm{\alpha}}_{\textrm{eff}}^{\textrm{b}}({\bm{x}})=-2\xi(\gamma_{1}\delta\gamma_{2}-\gamma_{2}\delta\gamma_{1}),$
(10)

where $\gamma_{i}$ is the $i$-component of the shear tensor of the dominant
lens at ${\bm{x}}$ and $\delta\gamma_{i}$ is the $i$-component of the shear
tensor of the dominant lens at ${\bm{y}}_{\textrm{b}}$

Ignoring second order terms, equations (6) and (10) can be combined to yield

$\eta=\frac{\nabla\times\delta{\bm{\alpha}}_{\textrm{eff}}({\bm{x}})}{2(\gamma_{1}\delta\gamma_{2}-\gamma_{2}\delta\gamma_{1})},$
(11)

where $-1<\eta<1$ and
$\delta{\bm{\alpha}}_{\textrm{eff}}({\bm{x}})=\delta{\bm{\alpha}}_{\textrm{eff}}^{\textrm{f}}({\bm{x}})$
and $\eta=\beta$ for $\eta>0$ and
$\delta{\bm{\alpha}}_{\textrm{eff}}({\bm{x}})=\delta{\bm{\alpha}}_{\textrm{eff}}^{\textrm{b}}({\bm{x}})$
and $\eta=-\xi$ for $\eta<0$. Thus, from measured rotation of the effective
deflection angle and shears of the dominant lens and the perturber, we can
measure the dimensionless distance ratio parameter $\eta$, that encodes the
information of LOS distance to the perturber. However, it should be noted that
$\eta$ is not a directly measured quantity as $\delta\gamma_{1}$ and
$\delta\gamma_{2}$ are not directly measured. Therefore, we need to estimate
$\eta$ using a certain approximation, which we will discuss in the next
section. If $\eta$ is negative/positive, the perturber resides in the
background/foreground of the dominant lens. In shear-aligned coordinates
$(x^{\prime}_{1},x^{\prime}_{2})$ in which the shear of the dominant lens is
aligned with $x^{\prime}_{1}$ or $x^{\prime}_{2}$ axis, i.e.
$\gamma_{2}^{\prime}=0$, the rotation is proportional to
$\delta\gamma_{2}^{\prime}$. In other words, the amplitude of rotation is
proportional to the non-diagonal shear component of the perturber.

In a similar manner, we can derive the divergence of the effective deflection
$\delta{\bm{\alpha}}_{\textrm{eff}}$. In the foreground and background cases,
the divergences are

$\nabla\cdot\delta{\bm{\alpha}}_{\textrm{eff}}^{\textrm{f}}({\bm{x}})=2\delta\kappa-2\beta(\tilde{\kappa}\delta\kappa+\tilde{\gamma_{1}}\delta\gamma_{1}+\tilde{\gamma_{2}}\delta\gamma_{2}),$
(12)

and

$\nabla\cdot\delta{\bm{\alpha}}_{\textrm{eff}}^{\textrm{b}}({\bm{x}})=2\delta\kappa-2\xi(\kappa\delta\kappa+\gamma_{1}\delta\gamma_{1}+\gamma_{2}\delta\gamma_{2}),$
(13)

respectively. Ignoring second order terms, in terms of $\eta$, equations (12)
and (13) give a modified Poisson equation,

$\displaystyle\nabla\cdot\delta{\bm{\alpha}}_{\textrm{eff}}({\bm{x}})$
$\displaystyle=$ $\displaystyle
2\delta\kappa-2|\eta|(\kappa\delta\kappa+\gamma_{1}\delta\gamma_{1}+\gamma_{2}\delta\gamma_{2}),$
(14) $\displaystyle=$ $\displaystyle 2\delta\kappa_{\textrm{eff}},$

where $\delta\kappa_{\textrm{eff}}$ is the effective convergence perturbation
which encodes the information of coupling between the dominant lens and the
perturber.

The effective deflection angle can be decomposed as magnetic and electric
components,

$\delta{\bm{\alpha}}_{\textrm{eff}}=\delta{\bm{\alpha}}^{\textrm{B}}_{\textrm{eff}}+\delta{\bm{\alpha}}^{\textrm{E}}_{\textrm{eff}},$
(15)

where
$\nabla\cdot\delta{\bm{\alpha}}^{\textrm{B}}_{\textrm{eff}}=\nabla\times\delta{\bm{\alpha}}^{\textrm{E}}_{\textrm{eff}}=0$.

Using equations (11) and (14), we can estimate the ratio of the magnetic
component to the Electric one in the limit $\eta\rightarrow 0$ as

$\displaystyle r^{\textrm{BE}}$ $\displaystyle\equiv$
$\displaystyle\frac{|\delta{\bm{\alpha}}^{\textrm{B}}_{\textrm{eff}}|}{|\delta{\bm{\alpha}}^{\textrm{E}}_{\textrm{{eff}}}|}$
(16) $\displaystyle\sim$
$\displaystyle\frac{|\nabla\times\delta{\bm{\alpha}}^{\textrm{B}}_{\textrm{eff}}|}{|\nabla\cdot\delta{\bm{\alpha}}^{\textrm{E}}_{\textrm{{eff}}}|}$
$\displaystyle\sim$
$\displaystyle\frac{|\eta\,\gamma^{\prime}_{1}|}{\sqrt{2}},$

where prime denotes the coordinates in which the magnification matrix of the
dominant lens is diagonalized, i.e., $\gamma^{\prime}_{2}=0$ and we assumed
that $|\delta\gamma^{\prime}_{2}|\sim|\delta\kappa^{\prime}|/\sqrt{2}$. Let us
suppose a typical lens system in which the dominant lens is modelled by a
singular isothermal sphere (SIS). In the vicinity of an Einstein ring, the
shear is $\gamma^{\prime}_{1}\sim 0.5$. Then the ratio is

$r^{\textrm{BE}}\sim\frac{\sqrt{2}|\eta|}{4}.$ (17)

Let us suppose a lens system with a perturber at $|\eta|=0.4$. Then we have
$r^{\textrm{BE}}\sim 1/7$. Thus the contribution from magnetic component is
expected to be subdominant unless the perturber resides too far from the
dominant lens.

###  2.3 Estimators of LOS distance

In order to estimate the LOS distance ratio $\eta$, we need to assess the
components $\delta\kappa$, $\delta\gamma_{1}$, and $\delta\gamma_{2}$ using
the spatial derivatives of the effective deflection
$\delta{\bm{\alpha}}_{\textrm{eff}}$,

$\delta M_{\textrm{eff}}=\begin{pmatrix}A&B\\\
C&D\end{pmatrix}\equiv\begin{pmatrix}\delta\alpha_{\textrm{eff}1,1}&\delta\alpha_{\textrm{eff}1,2}\\\
\delta\alpha_{\textrm{eff}2,1}&\delta\alpha_{\textrm{eff}2,1}\end{pmatrix}.$
(18)

First, we consider the foreground perturber case. From equations (5, 18),
neglecting second order terms, we have

$\displaystyle A$ $\displaystyle=$
$\displaystyle\delta\kappa+\delta\gamma_{1}-\beta[(\kappa+\gamma_{1})(\delta\kappa+\delta\gamma_{1})+\gamma_{2}\delta\gamma_{2}]$
$\displaystyle+$
$\displaystyle\tilde{\kappa}+\tilde{\gamma_{1}}-\kappa-\gamma_{1},$
$\displaystyle B$ $\displaystyle=$
$\displaystyle\delta\gamma_{2}-\beta(\kappa\delta\gamma_{2}+\gamma_{2}\delta\kappa+\gamma_{1}\delta\gamma_{2}-\gamma_{2}\delta\gamma_{1})$
$\displaystyle+$ $\displaystyle\tilde{\gamma}_{2}-\gamma_{2},$ $\displaystyle
C$ $\displaystyle=$
$\displaystyle\delta\gamma_{2}-\beta(\kappa\delta\gamma_{2}+\gamma_{2}\delta\kappa+\gamma_{2}\delta\gamma_{1}-\gamma_{1}\delta\gamma_{2})$
$\displaystyle+$ $\displaystyle\tilde{\gamma}_{2}-\gamma_{2},$ $\displaystyle
D$ $\displaystyle=$
$\displaystyle\delta\kappa-\delta\gamma_{1}-\beta[(\kappa-\gamma_{1})(\delta\kappa-\delta\gamma_{1})+\gamma_{2}\delta\gamma_{2})]$
(19) $\displaystyle+$
$\displaystyle\tilde{\kappa}-\tilde{\gamma_{1}}-\kappa+\gamma_{1}.$

In first order in $\delta{\bm{\alpha}}$, the last terms that represent small
changes in the magnification matrix of the dominant lens can be linearly
approximated as

$\displaystyle c_{1}$ $\displaystyle\equiv$
$\displaystyle\frac{\partial(\kappa+\gamma_{1})}{\partial{\bm{x}}}\biggr{|}_{\bm{x}}\cdot\delta{\bm{\alpha}}({\bm{x}})\approx-\beta^{-1}(\tilde{\kappa}+\tilde{\gamma_{1}}-\kappa-\gamma_{1}),$
$\displaystyle c_{2}$ $\displaystyle\equiv$
$\displaystyle\frac{\partial(\kappa-\gamma_{1})}{\partial{\bm{x}}}\biggr{|}_{\bm{x}}\cdot\delta{\bm{\alpha}}({\bm{x}})\approx-\beta^{-1}(\tilde{\kappa}-\tilde{\gamma_{1}}+\kappa-\gamma_{1}),$
$\displaystyle c_{3}$ $\displaystyle\equiv$
$\displaystyle\frac{\partial\gamma_{2}}{\partial{\bm{x}}}\biggr{|}_{\bm{x}}\cdot\delta{\bm{\alpha}}({\bm{x}})\approx-\beta^{-1}(\tilde{\gamma_{2}}-\gamma_{2}).$

Let us suppose that a perturber resides in the vicinity of the dominant lens,
i.e., $\beta\ll 1$. Then the magnetic component of
$\delta{\bm{\alpha}}_{\textrm{eff}}$ is much smaller than the electric
component and thus
$\delta{\bm{\alpha}}({\bm{x}})\approx\delta{\bm{\alpha}}_{\textrm{eff}}({\bm{x}})$.
Plugging this relation and equation (LABEL:eq:c1c2c3) into equation (19), we
obtain a quadratic equation in $\beta$ whose positive solution

$\displaystyle\hat{\beta}^{\textrm{A}}$ $\displaystyle\equiv$
$\displaystyle\frac{-\hat{b}-\textrm{sgn}\,\hat{c}\sqrt{\hat{b}^{2}-4\hat{a}\hat{c}}}{2\hat{a}}\approx\beta,$
$\displaystyle\hat{a}$ $\displaystyle\equiv$ $\displaystyle
2c_{3}\gamma_{1}+(c_{2}-c_{1})\gamma_{2},$ $\displaystyle\hat{b}$
$\displaystyle\equiv$
$\displaystyle(B+C)\gamma_{1}+(D-A)\gamma_{2}+(C-B)\kappa,$
$\displaystyle\hat{c}$ $\displaystyle\equiv$ $\displaystyle
B-C=-\nabla\times\delta{\bm{\alpha}}_{\textrm{eff}},$ (21)

is an estimator of $\beta$. If
$|\hat{c}|=2\beta|(\gamma_{2}\delta\gamma_{1}-\gamma_{1}\delta\gamma_{2})|\ll\hat{b}^{2}/|\hat{a}|$
and $\beta>0$, equation (21) can be further simplified as

$\displaystyle\hat{\beta}^{\textrm{B}}$ $\displaystyle\equiv$
$\displaystyle-\frac{\hat{c}}{\hat{b}},$ (22) $\displaystyle=$
$\displaystyle\frac{C-B}{(D-A)\gamma_{2}+(B+C)\gamma_{1}+(C-B)\kappa}$
$\displaystyle\approx$ $\displaystyle\beta.$

One can easily confirm that equation (22) is equivalent to equation (11)
provided that $\beta c_{i}$’s are sufficiently small. If not, equation (21),
which include effects from $\beta c_{i}$’s may give a much better
approximation compared to equation (22). Note that equations (21), (22) are
only valid in some limited regions of ${\bm{x}}$ in which the amplitude of
rotation $|\nabla\times\delta{\bm{\alpha}}_{\textrm{eff}}|$ is not ’very’
small and the strong lensing effect due to the dominant lens or the perturber
is not too large. For instance, if $|C-B|$ is sufficiently smaller than the
typical amplitude of linear perturbations
$\delta\kappa,\delta\gamma_{1},\delta\gamma_{2}$, the equations (21), (22)
give a bad approximation. In that case, second order or higher order
correction is necessary to give an accurate estimate. If the gradient of a
projected gravitational potential of either of the dominant or subdominant
lens is too large, equation (25) gives a bad approximation. In the former
case, $c_{i}$’s become too large. Then the linear approximation in equation
(LABEL:eq:c1c2c3) is no longer valid. In the latter case, the subdominant lens
dominates the lensing effect over that of the dominant lens and thus
$\delta\alpha$ is no longer smaller than $\alpha$.

Next, we consider the background perturber case. The components of the
perturbed magnification matrix $\delta M_{\textrm{eff}}$

$\displaystyle A$ $\displaystyle=$
$\displaystyle\delta\kappa+\delta\gamma_{1}-\xi[(\kappa+\gamma_{1})(\delta\kappa+\delta\gamma_{1})+\gamma_{2}\delta\gamma_{2}],$
$\displaystyle B$ $\displaystyle=$
$\displaystyle\delta\gamma_{2}-\xi(\kappa\delta\gamma_{2}+\gamma_{2}\delta\kappa+\gamma_{1}\delta\gamma_{2}-\gamma_{2}\delta\gamma_{1}),$
$\displaystyle C$ $\displaystyle=$
$\displaystyle\delta\gamma_{2}-\xi(\kappa\delta\gamma_{2}+\gamma_{2}\delta\kappa+\gamma_{2}\delta\gamma_{1}-\gamma_{1}\delta\gamma_{2}),$
$\displaystyle D$ $\displaystyle=$
$\displaystyle\delta\kappa-\delta\gamma_{1}-\xi[(\kappa-\gamma_{1})(\delta\kappa-\delta\gamma_{1})+\gamma_{2}\delta\gamma_{2}],$
(23)

yield an exact solution

$\xi=\frac{B-C}{(D-A)\gamma_{2}+(B+C)\gamma_{1}+(B-C)\kappa}.$ (24)

If the magnitude of rotation $|C-B|=|\hat{c}|$ is sufficiently smaller than
$|(D-A)\gamma_{2}+(B+C)\gamma_{1}|$, equations (22) and (24) yield an
estimator of $\eta$,

$\displaystyle\hat{\eta}$ $\displaystyle\equiv$
$\displaystyle\frac{C-B}{(D-A)\gamma_{2}+(B+C)\gamma_{1}}$ (25)
$\displaystyle\approx$ $\displaystyle\eta.$

If $\hat{\eta}$ is positive(negative), it is likely that a perturber resides
in the foreground (background).

###  2.4 Estimators of LOS perturbations

In a similar manner, we can express the ’bare’222Here ’bare’ means not
influenced by a dominant lens perturbations
$\delta\hat{\kappa},\delta\hat{\gamma}_{1},\delta\hat{\gamma}_{2}$ due to a
perturber in the LOS in terms of the derivatives of
$\delta{\bm{\alpha}}_{\textrm{eff}}$. In the following, we assume that
$\beta\ll 1$ and $\xi\ll 1$.

In the foreground perturber case, if $\beta c_{i}$’s are sufficiently small,
we have

$\displaystyle\delta\hat{\kappa}^{\textrm{B}}$ $\displaystyle=$ $\displaystyle
F[2(AC+BD)\gamma_{1}-(A^{2}+B^{2}-C^{2}-D^{2})\gamma_{2}]/2,$
$\displaystyle\delta\hat{\gamma}_{1}^{\textrm{B}}$ $\displaystyle=$
$\displaystyle F[2(AC-BD)\gamma_{1}$ $\displaystyle-$
$\displaystyle(A+B-C-D)(A-B+C-D)\gamma_{2}]/2,$
$\displaystyle\delta\hat{\gamma}_{2}^{\textrm{B}}$ $\displaystyle=$
$\displaystyle F[2BC\gamma_{1}-(AB-CD)\gamma_{2}],$ $\displaystyle F$
$\displaystyle\equiv$
$\displaystyle(\hat{\beta}^{\textrm{B}})^{-1}(C-B)[4BC\gamma_{1}^{2}-2(B+C)(A-D)\gamma_{1}\gamma_{2}$
(26) $\displaystyle-$ $\displaystyle(A+B-C-D)(-A+B-C+D)\gamma_{2}^{2}]^{-1}.$

If $\beta c_{i}$’s are not sufficiently small, the approximated analytic
solutions have a more complicated form (see Appendix).

Interchanging B and C and substituting $\xi$ into $\beta$ in equation (26),
the exact perturbations for a background perturber are

$\displaystyle\delta\kappa$ $\displaystyle=$ $\displaystyle
G[2(AB+CD)\gamma_{1}-(A^{2}-B^{2}+C^{2}-D^{2})\gamma_{2}]/2,$
$\displaystyle\delta\gamma_{1}$ $\displaystyle=$ $\displaystyle
G[2(AB+CD)\gamma_{1}$ $\displaystyle-$
$\displaystyle(A+B-C-D)(A-B+C-D)\gamma_{2}]/2$ $\displaystyle\delta\gamma_{2}$
$\displaystyle=$ $\displaystyle G[2BC\gamma_{1}-(AC-BD)\gamma_{2}],$
$\displaystyle G$ $\displaystyle\equiv$
$\displaystyle\xi^{-1}(B-C)[4BC\gamma_{1}^{2}-2(B+C)(A-D)\gamma_{1}\gamma_{2}$
(27) $\displaystyle-$ $\displaystyle(A+B-C-D)(-A+B-C+D)\gamma_{2}^{2}]^{-1}.$

##  3 Extended Multi-Plane Mass-sheet Transformation

In this section, we study how multi-plane mass-sheet transformation (MMST)
affects the LOS distance ratio and time delay in a double lens system with a
single source at a certain redshift. A scale transformation in the distance
ratio allow a non-zero mass-sheet transformation in both the foreground and
background lens planes. This implies that we have two degrees of freedom in
scale transformation if the redshift of a perturber is not known. We call such
a transformation (multi-plane MST plus a scale transformation in the distance
ratio of a perturber) an ’extended MMST’ (eMMST).

###  3.1 Distance ratio

Let us recall mass-sheet transformation (MST) in a single lens system with a
single source. The position of the source at a certain redshift is given by

${\bm{y}}={\bm{x}}-{\bm{\alpha}}({\bm{x}}).$ (28)

A scale transformation by a factor of $\lambda$

${\bm{y}}\rightarrow{\bm{y}}^{\prime}=\lambda{\bm{y}}$ (29)

and a scale transformation by a factor of $\lambda$ accompanied by an addition
of a deflection by a constant convergence

${\bm{\alpha}}{({\bm{x}})}\rightarrow{\bm{\alpha}}^{\prime}({\bm{x}})=\lambda{\bm{\alpha}}+(1-\lambda)\mbox{\boldmath$x$}$
(30)

leave a lens equation

${\bm{y}}^{\prime}={\bm{x}}-{\bm{\alpha}}^{\prime}({\bm{x}})$ (31)

invariant. The transformed lens system is equivalent to the original lens
system(28) except for the physical size and the intensity of the source. If we
do not know the true size or the intensity of the source, we cannot
distinguish between the two systems with different deflection angles. The set
of transformations (29) and (30) is called as MST for a single lens system.

In what follows, we consider MMST in a double lens system that consists of a
dominant lens and a subdominant lens which acts as a perturber.

First, we assume that the subdominant lens with a distance ratio parameter
$\beta$ resides in the foreground of the dominant lens. We consider a scale
transformation for the dominant lens

$\displaystyle{\bm{y}}\rightarrow{\bm{y}}^{\prime}$ $\displaystyle=$
$\displaystyle\lambda_{\textrm{d}}{\bm{y}}$ (32) $\displaystyle=$
$\displaystyle{\bm{x}}-\delta\tilde{{\bm{\alpha}}}-(1-\lambda_{\textrm{d}}){\bm{x}}-\lambda_{\textrm{d}}{\bm{\alpha}}({\bm{x}}-\beta\delta{\bm{\alpha}})$

and that for the foreground secondary lens

$\displaystyle{\bm{y}}_{\textrm{d}}\rightarrow{\bm{y}}^{\prime}_{\textrm{d}}$
$\displaystyle=$ $\displaystyle\lambda_{\textrm{f}}{\bm{y}}_{\textrm{d}}$ (33)
$\displaystyle=$
$\displaystyle{\bm{x}}-(1-\lambda_{\textrm{f}}){\bm{x}}-\lambda_{\textrm{f}}\beta\delta{\bm{\alpha}}({\bm{x}}).$

If the distance ratio parameter $\beta$ is not known, we need to consider a
scale transformation for $\beta$,

$\beta\rightarrow\beta^{\prime}=\lambda_{\beta}\beta.$ (34)

Then equation (35) can be written as

$\displaystyle{\bm{y}}_{\textrm{d}}\rightarrow{\bm{y}}^{\prime}_{\textrm{d}}$
$\displaystyle=$
$\displaystyle{\bm{x}}-\beta^{\prime}\delta\mbox{\boldmath$\alpha$}^{\prime}({\bm{x}})$
(35) $\displaystyle=$
$\displaystyle{\bm{x}}-\beta^{\prime}\biggl{[}\frac{1-\lambda_{\textrm{f}}}{\beta^{\prime}}{\bm{x}}+\frac{\lambda_{\textrm{f}}}{\lambda_{\beta}}\delta{\bm{\alpha}}({\bm{x}})\biggr{]},$

and equation (32) as

$\displaystyle{\bm{y}}^{\prime}$ $\displaystyle=$
$\displaystyle{\bm{x}}-\delta\tilde{{\bm{\alpha}}}-(1-\lambda_{\textrm{d}}){\bm{x}}-\lambda_{\textrm{d}}{\bm{\alpha}}({\bm{y}}^{\prime}_{\textrm{d}}/\lambda_{\textrm{f}})$
$\displaystyle=$
$\displaystyle{\bm{x}}-\delta\mbox{\boldmath$\alpha$}^{\prime}({\bm{x}})-\tilde{\mbox{\boldmath$\alpha$}}({\bm{y}}^{\prime}_{\textrm{d}}({\bm{x}},\beta^{\prime},\delta\mbox{\boldmath$\alpha$}^{\prime}))$
$\displaystyle-$
$\displaystyle\mbox{\boldmath$\varepsilon$}({\bm{y}}^{\prime}_{\textrm{d}}({\bm{x}},\beta^{\prime},\delta\mbox{\boldmath$\alpha$}^{\prime})),$
$\displaystyle\tilde{\mbox{\boldmath$\alpha$}}({\bm{y}}^{\prime}_{\textrm{d}})$
$\displaystyle\equiv$
$\displaystyle\lambda_{\textrm{d}}{\bm{\alpha}}({\bm{y}}^{\prime}_{\textrm{d}}/\lambda_{\textrm{f}}),$
(36)

where

$\mbox{\boldmath$\varepsilon$}({\bm{y}}^{\prime}_{\textrm{d}}({\bm{x}},\beta^{\prime},\delta\mbox{\boldmath$\alpha$}^{\prime}))$
(37)

is a residual function for which
$\mbox{\boldmath$\alpha$}^{\prime}\equiv\tilde{\mbox{\boldmath$\alpha$}}+\mbox{\boldmath$\varepsilon$}$
is a transformed deflection angle. Here we assume that
$|\mbox{\boldmath$\varepsilon$}|$ is sufficiently larger or smaller than
$|\mbox{\boldmath$\alpha$}^{\prime}|$. From equation (36), the residual
function satisfies

$\mbox{\boldmath$\varepsilon$}({\bm{y}}^{\prime}_{\textrm{d}})=\biggl{[}\frac{-1+\lambda_{\textrm{f}}}{\beta^{\prime}}-\lambda_{\textrm{d}}+1\biggr{]}{\bm{x}}+\biggl{[}\lambda_{\textrm{d}}-\frac{\lambda_{\textrm{f}}}{\lambda_{\beta}}\biggr{]}\delta{\bm{\alpha}}.$
(38)

In terms of ${\bm{x}}$ and $\delta{\bm{\alpha}}$,
$\mbox{\boldmath$y$}^{\prime}_{\textrm{d}}$ can be represented as

$\mbox{\boldmath$y$}^{\prime}_{\textrm{d}}=\lambda_{\textrm{f}}({\bm{x}}-\beta\delta{\bm{\alpha}}),$
(39)

Since $\varepsilon$ is a function of $\mbox{\boldmath$y$}^{\prime}_{d}$, the
ratios of the coefficients of ${\bm{x}}$ and $\delta{\bm{\alpha}}$ in
equations (38) and (39) should be equivalent to $\beta^{-1}$. Then we have

$\lambda_{\beta}=\frac{1}{(1-\lambda_{\textrm{d}})\beta+\lambda_{\textrm{d}}}.$
(40)

If $\beta$ is known, i.e., $\lambda_{\beta}=1$, then equation (40) implies
$\lambda_{\textrm{d}}=1$, which is the known result for MMST in double (main)
lens systems. However, if $\beta$ is not known, any value of
$\beta^{\prime}=\lambda_{\beta}\beta$ can explain most of lensing phenomena by
an appropriate MST (specified by $\lambda_{\textrm{d}}$) for the dominant
(background) lens. Thus we have two degrees of freedom among the three scale
parameters $\lambda_{\textrm{f}},\lambda_{\textrm{d}},$ and $\lambda_{\beta}$.
If $|\mbox{\boldmath$\varepsilon$}|$ is sufficiently smaller than
$|\delta{\bm{\alpha}}|$ then $\varepsilon$ is not necessarily a function of
$\mbox{\boldmath$y$}^{\prime}_{\textrm{d}}$. However, from equation (38), one
can easily show that
$\lambda_{\beta}\approx\lambda_{\textrm{f}}\approx\lambda_{\textrm{d}}\approx
1$. Therefore, we do not have any solutions in which $\beta^{\prime}$
significantly differs from $\beta$ (i.e., $|1-\lambda_{\textrm{d}}|$ is not so
small with respect to unity).

Second, we assume that the subdominant lens with a distance ratio parameter
$\xi$ resides in the background of the dominant lens. We consider a scale
transformation for the dominant lens

$\displaystyle{\bm{y}}\rightarrow{\bm{y}}^{\prime}$ $\displaystyle=$
$\displaystyle\lambda_{\textrm{d}}{\bm{y}}$ (41) $\displaystyle=$
$\displaystyle{\bm{x}}-\tilde{{\bm{\alpha}}}-(1-\lambda_{\textrm{d}}){\bm{x}}-\lambda_{\textrm{d}}\delta{\bm{\alpha}}({\bm{x}}-\xi\delta{\bm{\alpha}}),$

that for the background secondary lens

$\displaystyle{\bm{y}}_{\textrm{b}}\rightarrow{\bm{y}}^{\prime}_{b}$
$\displaystyle=$ $\displaystyle\lambda_{\textrm{b}}{\bm{y}}_{\textrm{b}}$ (42)
$\displaystyle=$
$\displaystyle{\bm{x}}-(1-\lambda_{\textrm{b}}){\bm{x}}-\lambda_{\textrm{b}}\xi\delta{\bm{\alpha}}({\bm{x}}),$

and that for the distance ratio parameter

$\xi\rightarrow\xi^{\prime}=\lambda_{\xi}\xi.$ (43)

By applying a similar argument above to equations (41) , (42), and (43), we
have

$\lambda_{\xi}=\frac{1}{(1-\lambda_{\textrm{d}})\xi+\lambda_{\textrm{d}}}.$
(44)

Thus the mass-sheet degeneracy (MSD) in the distance ratios $\beta$ and $\xi$,
which is related to the scale transformation in the background lens remain.

###  3.2 Time delay

Figure 1: An example of a light path for a foreground perturber. A light ray
emitted from a point source at S in the source plane is deflected by a
dominant lens at D and a subdominant lens at F and reaches to an obserber at
O.

First we examine a double lens system with a foreground perturber. Then the
time delay $\tau$ of a source at $y$, which is defined as the arrival time
difference between a light path that is lensed by a dominant lens at known
redshift $z_{\textrm{d}}$ and a foreground perturber at unknown redshift
$z_{\textrm{f}}$ and an unlensed path is

$\displaystyle c\tau$ $\displaystyle=$ $\displaystyle
c\tau_{\textrm{f}}+c\tau_{\textrm{d}}$ $\displaystyle=$
$\displaystyle(1+z_{\textrm{f}})\frac{D_{\textrm{f}}D_{\textrm{d}}}{D_{\textrm{df}}}\biggl{[}-\nabla^{-1}_{\bm{x}}\bigl{(}\beta\delta{\bm{\alpha}}({\bm{x}})\bigr{)}+\frac{({\bm{x}}-{\bm{y}}_{\textrm{d}})^{2}}{2}\biggr{]}$
$\displaystyle+$
$\displaystyle(1+z_{\textrm{d}})\frac{D_{\textrm{d}}D_{\textrm{s}}}{D_{\textrm{sd}}}\biggl{[}-\nabla^{-1}_{{\bm{y}}_{\textrm{d}}}\bigl{(}{\bm{\alpha}}({\bm{y}}_{\textrm{d}})\bigr{)}+\frac{({\bm{y}}_{\textrm{d}}-{\bm{y}})^{2}}{2}\biggr{]},$

where $\tau_{\textrm{f}}$ and $\tau_{\textrm{d}}$ are the arrival time
differences between light paths DFO and DO, and SDO and SO, respectively (as
depicted in figure 1). The constant $c$ is the light speed. We investigate how
the time delay $\tau$ behaves under eMMST. The scalings (32), (33), and (34)
yield

$\displaystyle c\tau^{\prime}$ $\displaystyle=$ $\displaystyle
c\tau^{\prime}_{\textrm{f}}+c\tau^{\prime}_{\textrm{d}}$ $\displaystyle=$
$\displaystyle(1+z^{\prime}_{\textrm{f}})\frac{D^{\prime}_{\textrm{f}}D_{\textrm{d}}}{D^{\prime}_{\textrm{df}}}\biggl{[}-\nabla^{-1}_{\bm{x}}\bigl{(}\beta^{\prime}\delta{\bm{\alpha}}^{\prime}({\bm{x}})\bigr{)}+\frac{({\bm{x}}-{\bm{y}}^{\prime}_{\textrm{d}})^{2}}{2}\biggr{]}$
$\displaystyle+$
$\displaystyle(1+z_{\textrm{d}})\frac{D_{\textrm{d}}D_{\textrm{s}}}{D_{\textrm{sd}}}\biggl{[}-\nabla^{-1}_{{\bm{y}}^{\prime}_{\textrm{d}}}\bigl{(}{\bm{\alpha}}^{\prime}({\bm{y}}^{\prime}_{\textrm{d}})\bigr{)}+\frac{({\bm{y}}^{\prime}_{\textrm{d}}-{\bm{y}}^{\prime})^{2}}{2}\biggr{]}$
$\displaystyle=$
$\displaystyle(1+z^{\prime}_{\textrm{f}})\frac{D^{\prime}_{\textrm{f}}D_{\textrm{d}}}{D^{\prime}_{\textrm{df}}}\biggl{[}-\lambda_{\textrm{f}}\nabla^{-1}_{\bm{x}}\bigl{(}\beta\delta{\bm{\alpha}}({\bm{x}})\bigr{)}-\frac{(1-\lambda_{\textrm{f}})x^{2}}{2}$
$\displaystyle+$
$\displaystyle\frac{({\bm{x}}-\lambda_{\textrm{f}}{\bm{y}}_{\textrm{d}})^{2}}{2}\biggr{]}+(1+z_{\textrm{d}})\frac{D_{\textrm{d}}D_{\textrm{s}}}{D_{\textrm{sd}}}\biggl{[}-\lambda_{\textrm{f}}\lambda_{\textrm{d}}\nabla^{-1}_{{\bm{y}}_{\textrm{d}}}\bigl{(}{\bm{\alpha}}({\bm{y}}_{\textrm{d}})\bigr{)}$
$\displaystyle-$
$\displaystyle\frac{\lambda_{\textrm{f}}^{2}h\,y^{2}_{\textrm{d}}}{2}+\frac{(\lambda_{\textrm{d}}{\bm{y}}_{\textrm{d}}-\lambda_{\textrm{f}}{\bm{y}})^{2}}{2}\biggr{]},$
$\displaystyle h$ $\displaystyle\equiv$
$\displaystyle-\frac{\lambda_{\textrm{d}}}{\beta\lambda_{\textrm{f}}}+\frac{1}{\lambda_{\beta}\beta}.$
(46)

Since

$\displaystyle(1+z^{\prime}_{\textrm{f}})\frac{D^{\prime}_{\textrm{f}}D_{\textrm{d}}}{D^{\prime}_{\textrm{df}}}\biggl{[}-\frac{(1-\lambda_{\textrm{f}})x^{2}}{2}+\frac{({\bm{x}}-\lambda_{\textrm{f}}{\bm{y}}_{\textrm{d}})^{2}}{2}\biggr{]}$
(47) $\displaystyle+$
$\displaystyle(1+z_{\textrm{d}})\frac{D_{\textrm{d}}D_{\textrm{s}}}{D_{\textrm{sd}}}\biggl{[}-\frac{\lambda_{\textrm{f}}^{2}h\,y^{2}_{\textrm{d}}}{2}+\frac{(\lambda_{\textrm{d}}{\bm{y}}_{\textrm{d}}-\lambda_{\textrm{f}}{\bm{y}})^{2}}{2}\biggr{]}$
$\displaystyle-$
$\displaystyle(1+z^{\prime}_{\textrm{f}})\frac{D^{\prime}_{\textrm{f}}D_{\textrm{d}}}{D^{\prime}_{\textrm{df}}}\frac{\lambda_{\textrm{f}}({\bm{x}}-{\bm{y}}_{\textrm{d}})^{2}}{2}$
$\displaystyle-$
$\displaystyle(1+z_{\textrm{d}})\frac{D_{\textrm{d}}D_{\textrm{s}}}{D_{\textrm{sd}}}\frac{\lambda_{\textrm{d}}\lambda_{\textrm{f}}({\bm{y}}_{\textrm{d}}-{\bm{y}})^{2}}{2}$
$\displaystyle=$
$\displaystyle(1+z^{\prime}_{\textrm{f}})\frac{D^{\prime}_{\textrm{f}}D_{\textrm{d}}}{D^{\prime}_{\textrm{df}}}\frac{\lambda_{\textrm{d}}(\lambda_{\textrm{d}}-\lambda_{\textrm{f}})y^{2}}{2}$
$\displaystyle+$
$\displaystyle(1+z_{\textrm{d}})\frac{D_{\textrm{d}}D_{\textrm{s}}}{D_{\textrm{sd}}}\frac{\lambda_{\textrm{f}}(\lambda_{\textrm{f}}-1)(\lambda_{\textrm{d}}-1)y^{2}_{\textrm{d}}}{2},$

and

$(1+z^{\prime}_{\textrm{f}})\frac{D^{\prime}_{\textrm{f}}D_{\textrm{d}}}{D^{\prime}_{\textrm{df}}}=-(1+z_{\textrm{d}})\frac{D_{\textrm{d}}D_{\textrm{s}}}{D_{\textrm{sd}}}\lambda_{\textrm{d}}\biggl{[}1-\frac{1}{\beta}\biggr{]},$
(48)

equation (46) gives

$\displaystyle c\tau^{\prime}$ $\displaystyle=$
$\displaystyle\lambda_{\textrm{f}}\lambda_{\textrm{d}}c\tau+\frac{\lambda_{\textrm{d}}(\lambda_{\textrm{d}}-\lambda_{\textrm{f}})F_{\textrm{sd}}y^{2}}{2},$
$\displaystyle F_{\textrm{sd}}$ $\displaystyle\equiv$
$\displaystyle(1+z_{\textrm{d}})\frac{D_{\textrm{d}}D_{\textrm{s}}}{D_{\textrm{sd}}}.$
(49)

Since the last term in equation (49) depend only on the position of a source,
it does not contribute to time delay between lensed images. Thus eMMST admits
a scale transformation in the foreground lens plane and another one in the
background lens plane. Let us suppose that quadruply lensed images consist of
two images that are lensed by only a dominant lens and the other two images
that are perturbed by a foreground halo. Then the time delay between
unperturbed images gives $\lambda_{\textrm{d}}$ and that between perturbed
images gives $\lambda_{\textrm{f}}$ given a Hubble constant $H_{0}$.

Figure 2: An example of a light path for a background perturber. A light ray
emitted from a point source at S in the source plane is deflected by a
subdominant lens at B and a dominant lens at F and reaches to an obserber at
O.

Next we examine a double lens system with a background perturber. In this
case, the time delay of a source at ${\bm{y}}$ is

$\displaystyle c\tau$ $\displaystyle=$ $\displaystyle
c\tau_{\textrm{d}}+c\tau_{\textrm{b}}$ $\displaystyle=$
$\displaystyle(1+z_{\textrm{d}})\frac{D_{\textrm{d}}D_{\textrm{b}}}{D_{\textrm{bd}}}\biggl{[}-\nabla^{-1}_{\bm{x}}\bigl{(}\xi{\bm{\alpha}}({\bm{x}})\bigr{)}+\frac{({\bm{x}}-{\bm{y}}_{\textrm{b}})^{2}}{2}\biggr{]}$
$\displaystyle+$
$\displaystyle(1+z_{\textrm{b}})\frac{D_{\textrm{b}}D_{\textrm{s}}}{D_{\textrm{sb}}}\biggl{[}-\nabla^{-1}_{{\bm{y}}_{\textrm{b}}}\bigl{(}\delta{\bm{\alpha}}({\bm{y}}_{\textrm{d}})\bigr{)}+\frac{({\bm{y}}_{\textrm{b}}-{\bm{y}})^{2}}{2}\biggr{]},$

where $\tau_{\textrm{d}}$ and $\tau_{\textrm{b}}$ are the arrival time
differences between light paths BDO and BO, and SBO and SO, respectively (as
depicted in figure 2). By applying a similar argument above, the scalings
(41), (42), and (43) yield

$\displaystyle c\tau^{\prime}$ $\displaystyle=$ $\displaystyle
c\tau^{\prime}_{\textrm{d}}+c\tau^{\prime}_{\textrm{b}}$ $\displaystyle=$
$\displaystyle(1+z_{\textrm{d}})\frac{D_{\textrm{d}}D^{\prime}_{\textrm{b}}}{D^{\prime}_{\textrm{bd}}}\biggl{[}-\nabla^{-1}_{\bm{x}}\bigl{(}\xi^{\prime}{\bm{\alpha}}^{\prime}({\bm{x}})\bigr{)}+\frac{({\bm{x}}-{\bm{y}}^{\prime}_{\textrm{b}})^{2}}{2}\biggr{]}$
$\displaystyle+$
$\displaystyle(1+z^{\prime}_{\textrm{b}})\frac{D^{\prime}_{\textrm{b}}D_{\textrm{s}}}{D^{\prime}_{\textrm{sb}}}\biggl{[}-\nabla^{-1}_{{\bm{y}}^{\prime}_{b}}\bigl{(}\delta{\bm{\alpha}}^{\prime}({\bm{y}}^{\prime}_{b})\bigr{)}+\frac{({\bm{y}}^{\prime}_{b}-{\bm{y}}^{\prime})^{2}}{2}\biggr{]}$
$\displaystyle=$
$\displaystyle(1+z_{\textrm{d}})\frac{D_{\textrm{d}}D^{\prime}_{\textrm{b}}}{D^{\prime}_{\textrm{bd}}}\biggl{[}-\lambda_{\textrm{d}}\nabla^{-1}_{\bm{x}}\bigl{(}\xi{\bm{\alpha}}({\bm{x}})\bigr{)}-\frac{(1-\lambda_{\textrm{d}})x^{2}}{2}$
$\displaystyle+$
$\displaystyle\frac{({\bm{x}}-\lambda_{\textrm{d}}{\bm{y}}_{\textrm{b}})^{2}}{2}\biggr{]}+(1+z^{\prime}_{\textrm{b}})\frac{D^{\prime}_{\textrm{b}}D_{\textrm{s}}}{D^{\prime}_{\textrm{sb}}}\biggl{[}-\lambda_{\textrm{d}}\lambda_{\textrm{b}}\nabla^{-1}_{{\bm{y}}_{\textrm{b}}}\bigl{(}\delta{\bm{\alpha}}({\bm{y}}_{\textrm{b}})\bigr{)}$
$\displaystyle-$
$\displaystyle\frac{\lambda_{\textrm{d}}^{2}g\,y^{2}_{b}}{2}+\frac{(\lambda_{\textrm{d}}{\bm{y}}_{\textrm{b}}-\lambda_{\textrm{b}}{\bm{y}})^{2}}{2}\biggr{]},$
$\displaystyle g$ $\displaystyle\equiv$
$\displaystyle-\frac{\lambda_{\textrm{b}}}{\xi\lambda_{\textrm{d}}}+\frac{1}{\lambda_{\xi}\xi}$
(51)

Since

$\displaystyle(1+z_{\textrm{d}})\frac{D_{\textrm{d}}D^{\prime}_{\textrm{b}}}{D^{\prime}_{\textrm{bd}}}\biggl{[}-\frac{(1-\lambda_{\textrm{d}})x^{2}}{2}+\frac{({\bm{x}}-\lambda_{\textrm{d}}{\bm{y}}_{\textrm{b}})^{2}}{2}\biggr{]}$
(52) $\displaystyle+$
$\displaystyle(1+z^{\prime}_{\textrm{b}})\frac{D^{\prime}_{\textrm{b}}D_{\textrm{s}}}{D^{\prime}_{\textrm{sd}}}\biggl{[}-\frac{\lambda_{\textrm{d}}^{2}g\,y^{2}_{b}}{2}+\frac{(\lambda_{\textrm{d}}{\bm{y}}_{\textrm{b}}-\lambda_{\textrm{b}}{\bm{y}})^{2}}{2}\biggr{]}$
$\displaystyle-$
$\displaystyle(1+z_{\textrm{d}})\frac{D_{\textrm{d}}D^{\prime}_{\textrm{b}}}{D^{\prime}_{\textrm{bd}}}\frac{\lambda_{\textrm{d}}({\bm{x}}-{\bm{y}}_{\textrm{b}})^{2}}{2}$
$\displaystyle-$
$\displaystyle(1+z^{\prime}_{\textrm{b}})\frac{D^{\prime}_{\textrm{b}}D_{\textrm{s}}}{D^{\prime}_{\textrm{sd}}}\frac{\lambda_{\textrm{b}}\lambda_{\textrm{d}}({\bm{y}}_{\textrm{b}}-{\bm{y}})^{2}}{2}$
$\displaystyle=$
$\displaystyle(1+z^{\prime}_{\textrm{b}})\frac{D^{\prime}_{\textrm{b}}D_{\textrm{s}}}{D^{\prime}_{\textrm{sb}}}\frac{\lambda_{\textrm{b}}(\lambda_{\textrm{b}}-\lambda_{\textrm{d}})y^{2}}{2}$
$\displaystyle+$
$\displaystyle(1+z^{\prime}_{\textrm{b}})\frac{D^{\prime}_{\textrm{b}}D_{\textrm{s}}}{D^{\prime}_{\textrm{sb}}}\frac{\lambda_{\textrm{b}}\lambda_{\textrm{d}}(\lambda_{\textrm{d}}-1)(1-\xi^{-1})y_{\textrm{b}}^{2}}{2}$
$\displaystyle+$
$\displaystyle(1+z_{\textrm{d}})\frac{D_{\textrm{d}}D^{\prime}_{\textrm{b}}}{D^{\prime}_{\textrm{bd}}}\frac{\lambda_{\textrm{d}}(\lambda_{\textrm{d}}-1)y^{2}_{b}}{2}$

and

$\displaystyle(1+z_{\textrm{d}})\frac{D_{\textrm{d}}D^{\prime}_{\textrm{b}}}{D^{\prime}_{\textrm{bd}}}$
$\displaystyle=$
$\displaystyle(1+z_{d})\frac{D_{s}D_{d}}{D_{sd}}\frac{1}{\lambda_{\xi}\xi},$
$\displaystyle(1+z^{\prime}_{\textrm{b}})\frac{D^{\prime}_{\textrm{b}}D_{\textrm{s}}}{D^{\prime}_{\textrm{sb}}}$
$\displaystyle=$
$\displaystyle-(1+z_{d})\frac{D_{s}D_{d}}{D_{sd}}\frac{1}{\lambda_{\xi}\xi-1},$
(53)

equation (51) gives

$\displaystyle c\tau^{\prime}$ $\displaystyle=$
$\displaystyle\frac{\lambda_{\textrm{d}}}{\lambda_{\xi}}c\tau+\frac{\lambda_{\textrm{b}}(\lambda_{\textrm{b}}-\lambda_{\textrm{d}})F^{\prime}_{\textrm{sb}}y^{2}}{2},$
$\displaystyle F^{\prime}_{\textrm{sb}}$ $\displaystyle\equiv$
$\displaystyle(1+z^{\prime}_{\textrm{b}})\frac{D^{\prime}_{\textrm{b}}D_{\textrm{s}}}{D^{\prime}_{\textrm{sb}}},$
(54)

where $\lambda_{\xi}$ is given by equation (44). Since the last term in
equation (54) depends solely on the position of the source, it does not
contribute to the time delay between lensed images. In the case of background
perturbers, eMMST admits a scale transformation in the foreground lens plane
and another one in the background lens plane. However, the transformed time
delay between lensed images is inversely proportional to the scale factor
$\lambda_{\xi}$. If $\lambda_{\textrm{d}}$ is determined by the observation of
time delays between unperturbed images, it is possible to measure
$\lambda_{\xi}$ with astrometric lensing B-mode for an assumed Hubble constant
$H_{0}$.

In scenarios where the redshifts of both foreground and background perturbers,
as well as the dominant lens and the source, are known, observations of
astrometric lensing B-mode can break the MSD. This leads to a reduction in
systematic errors in our estimated value of $H_{0}$. The reason is as follows:
Let’s consider a scenario where a quadruple lens system, with images A, B, C,
and D, is perturbed by a foreground perturber affecting image A and a
background perturber affecting image B. We assume that their gravitational
effects are confined to the vicinity of lensed images A and B, respectively.
Additionally, we assume that the gravitational influence of all other
perturbers, apart from the dominant lens, can be considered negligible.In this
context, a measurement of astrometric lensing B-mode in the de-lensed image A
can provide us with $\lambda_{\textrm{d}}$, as we already have knowledge of
$\beta$. Similarly, a measurement of astrometric lensing B-mode in the de-
lensed image B can furnish us with either $\lambda_{\xi}$ or
$\lambda_{\textrm{b}}$ since $\xi$ is a known quantity. Subsequently, by
measuring the time delay between BC or BD, we can determine $H_{0}$, while
measuring the time delay between AC or AD allows us to ascertain
$\lambda_{\textrm{f}}$.

##  4 Toy models

We investigate the property of astrometric lensing B-mode and the accuracy of
the approximated distance ratios
$\hat{\beta}^{\textrm{A}},\hat{\beta}^{\textrm{B}},$ and $\hat{\eta}$ for a
foreground perturber using simple toy models. As a model of a galaxy halo, we
adopt an singular isothermal sphere (SIS) whose reduced deflection angle is
given by ${\bm{\alpha}}({\bm{x}})=\theta_{\textrm{E}}{\bm{x}}/x$.

One can easily show that any constant shift perturbation in a secondary lens
plane can be explained by a translation in the source plane without a
secondary lens. Since we do not have any information about the original
position of the source, we cannot measure a constant shift. Moreover, a
constant external convergence perturbation cannot be measured due to the
eMMST. Therefore, in order to probe the effects of large-scale perturbation,
we analyse the dominant SIS with a constant external shear at an arbitrary
redshift. In order to probe the effects of small-scale perturbation, we
analyse the dominant SIS with another subdominant SIS with an Einstein radius
of $\theta_{\textrm{{Ep}}}\ll\theta_{\textrm{E}}$ that acts as a perturber at
an arbitrary redshift.

Figure 3: Line-of-sight effect for an SIS model with an external constant
shear in the shear-aligned coordinates. The dominant lens is an SIS with
Einstein radius $\theta_{\textrm{E}}=1\,^{\prime\prime}$ centered at $(0,0)$.
The arrows show the effective deflection angle
$\delta{\bm{\alpha}}_{\textrm{eff}}$ due to an external shear with
$\delta\gamma=0.1$. The colors show the amplitudes of the shift in arcsec.
Figure 4: Astrometric lensing B and E-modes for an external shear of
$\delta\gamma=0.1$. The arrows show the B and E components in the effective
astrometric shift $\delta{\bm{\alpha}}_{\textrm{eff}}$. The model parameters
are the same as in figure 3. The colors show the amplitudes of the shift in
arcsec.

###  4.1 SIS + external shear

Here we study a lens model with a dominant SIS and a constant external shear
$\delta\gamma$. We assume that the center of the effective gravitational
potential of the dominant SIS and that of the constant shear coincide. If a
perturber resides in the foreground of a dominant SIS, the effective
deflection angle of the perturber is

$\delta{\bm{\alpha}}_{\textrm{eff}}^{\textrm{f}}({\bm{x}})=\delta\gamma\mathcal{I}{\bm{x}}+\theta_{\textrm{E}}\frac{{\bm{x}}-\beta\delta\gamma\mathcal{I}{\bm{x}}}{\lVert{\bm{x}}-\beta\delta\gamma\mathcal{I}{\bm{x}}\rVert}-\theta_{\textrm{E}}\frac{{\bm{x}}}{x},$
(55)

where $\mathcal{I}$ is a unit matrix. If a perturber resides in the background
of a dominant SIS, the effective deflection angle of the perturber is

$\delta{\bm{\alpha}}_{\textrm{eff}}^{\textrm{b}}({\bm{x}})=\delta\gamma\mathcal{I}\biggl{[}{\bm{x}}-\xi\theta_{\textrm{E}}\frac{{\bm{x}}}{x}\biggr{]}.$
(56)

The ’magnetic’ potential $\psi^{\textrm{B}}$ can be obtained by solving the
Poisson equation

$-\Delta\psi^{\textrm{B}}=\nabla\times\delta{\bm{\alpha}}_{\textrm{eff}}^{\textrm{f}}({\bm{x}}).$
(57)

To solve equation (57) numerically, we used the Finite Element Method (FEM).
We discretized the inner region of a disk with a radius of
$\theta=100\,\theta_{\textrm{E}}$ into $\sim 2\times 10^{5}$ triangle meshes
in the polar coordinates $(\theta,\phi)$. We used finer meshes in regions with
$\theta<2\,\theta_{\textrm{E}}$ in order to resolve sudden changes in the
perturbed magnetic potential. We imposed a Dirichlet boundary condition at
$\theta=100\,\theta_{\textrm{E}}$. The amplitude of rotation can be estimated
as
$|\nabla\times\delta{\bm{\alpha}}_{\textrm{eff}}^{\textrm{f}}|\sim\beta\delta\gamma/\theta$
or
$|\nabla\times\delta{\bm{\alpha}}_{\textrm{eff}}^{\textrm{b}}|\sim\xi\delta\gamma/\theta$.
Then, the amplitude at $\theta=100\,\theta_{\textrm{E}}$ is expected to be
$\lesssim 10^{-3}$ for $\delta\gamma\leq 0.1$. We numerically checked that the
tiny error at the boundary of the disk does not significantly affect the
accuracy of the obtained solution at $\theta\lesssim\theta_{\textrm{E}}$. In
order to avoid a singularity at the center of the dominant lens, we used a
cored isothermal sphere with very small core radius
$\theta_{\textrm{c}}\ll\theta_{\textrm{E}}$. The effective deflection angle is
given by
${\bm{\alpha}}({\bm{x}})=\theta_{\textrm{E}}(\sqrt{x^{2}+\theta_{\textrm{c}}^{2}}-\theta_{\textrm{c}}){\bm{x}}/x$.
Except for the neighbourhood of singular points (centres of SIS), the absolute
errors in the Poisson equation (57) for numerically obtained solutions of
astrometric lensing B-mode $\delta{\bm{\alpha}}^{\textrm{B}}_{\textrm{eff}}$
were typically $\lesssim 10^{-2}$.

The effective deflection angles $\delta{\bm{\alpha}}_{\textrm{eff}}$ in an SIS
with an external shear of $\gamma\\!=\\!0.1$ are shown in figure 3. If objects
that causes an external shear reside in the foreground(background) of a
dominant lens, the constant shear produces a cross shape (ring shape) pattern
in the field of $\lVert\delta{\bm{\alpha}}_{\textrm{eff}}\rVert$. Except for
the central region in the center of an SIS, the constant shear reduces the
amplitude of $\delta{\bm{\alpha}}_{\textrm{eff}}$. The amplitudes of B-mode
$|\delta{\bm{\alpha}}_{\textrm{eff}}^{\textrm{B}}|$ is as large as $\sim
0\aas@@fstack{\prime\prime}03$ in the vicinity of coordinate axes and
negligibly small in the vicinity of diagonal lines if $\beta=0.5$ or
$\xi=0.5$. The direction of $\delta{\bm{\alpha}}_{\textrm{eff}}^{\textrm{B}}$
in the foreground case is opposite to that in the background case. In
contrast, the amplitudes of E-mode
$|\delta{\bm{\alpha}}_{\textrm{eff}}^{\textrm{E}}|$ are largest in the
diagonal lines and smallest in the coordinate axes except for the central
region (figure 4). For a given distance from the center, the amplitudes of
rotation $\nabla\times\delta{\bm{\alpha}}_{\textrm{eff}}^{\textrm{B}}$ are
largest in the diagonal lines (figure fig:SISES-rot). The ’direction’ of
rotation in the foreground case is opposite to that in the background case.
From the ’direction’ of rotation, we can determine whether the perturber
resides in the foreground or background of the dominant lens.

Figure 5: Rotation of effective deflection angle for an external shear of
$\delta\gamma=0.1$. The colors show the amplitudes of rotation. The model
parameters are the same as in figure 3.  Figure 6: The ratio of B and E-modes
for an external shear of $\delta\gamma=0.1$ and
$\theta_{\textrm{E}}=1^{\prime\prime}$. The orange full curves denote the
ratio of the averaged amplitude of B-mode astrometric shift and that of E-mode
astrometric shift as a function of $\beta$ (top) and $\xi$ (bottom). Both the
modes are averaged over a circle with $\theta=\theta_{\textrm{E}}$. The blue
dashed curves show the ratio of averaged rotation and averaged convergence.
Both the two contributions are averaged over a circle with
$\theta=\theta_{\textrm{E}}$. The blue dot-dashed curves represent the
approximations of the ratios for the foreground case $\sqrt{2}\beta/4$ (top)
and the background case $\sqrt{2}\xi/4$ (bottom).  Figure 7: Relative errors
of approximated distance ratios in the lens plane for $\delta\gamma=0.1$,
$\beta=0.5$, and $\theta_{\textrm{E}}=1^{\prime\prime}$. The colors show the
relative errors $\hat{\beta}^{\textrm{A}}/\beta-1$ (left),
$\hat{\beta}^{\textrm{B}}/\beta-1$ (middle), and $\hat{\eta}/\beta-1$ (right).
Figure 8: Relative errors of approximated distance ratio as a function of
$\beta$ and $\delta\gamma$. The colors show the relative errors averaged over
an Einstein ring with radius $\theta_{\textrm{E}}$ for a given set of
$\delta\gamma$ and $\beta$.

In figure 6, we show the ratio of the B-mode to E-mode
$r^{\textrm{BE}}=|\delta{\bm{\alpha}}_{\textrm{eff}}^{\textrm{B}}|/|\delta{\bm{\alpha}}_{\textrm{eff}}^{\textrm{E}}|$,
averaged over the Einstein ring with radius
$\theta_{\textrm{E}}=1^{\prime\prime}$. The ratio $r^{\textrm{BE}}$ is a
monotonically increasing function of $\beta$ or $\xi$, which equals $\sim 0.5$
for $\beta=0.5$, and $\sim 0.3$ for $\xi=0.5$. We found that the ratio of
rotation to divergence,
$|\nabla\times\delta{\bm{\alpha}}_{\textrm{eff}}^{\textrm{B}}|/|\nabla\cdot\delta{\bm{\alpha}}_{\textrm{eff}}^{\textrm{E}}|$
averaged over the Einstein ring gives a good approximation of
$r^{\textrm{BE}}$. The analytic formula in equation (17) gives a good estimate
for small distance ratio $\beta\lesssim 0.1$ or $\xi\lesssim 0.1$, but the
difference is conspicuous for $\beta\gtrsim 0.1$ or $\xi\gtrsim 0.1$. Next, we
show the relative errors of the distance ratio estimators
$\hat{\beta}^{\textrm{A}}$, $\hat{\beta}^{\textrm{B}}$, and $\hat{\eta}$ in
the dominant lens plane in figure 7. We found that these estimators give a
worse approximation in the central region within the Einstein radius
$\theta_{\textrm{E}}$ and a good approximation in the outer regions
$\theta>\theta_{\textrm{E}}$, especially in the vicinity of the horizontal
axis $x_{1}$. At $\theta=\theta_{\textrm{E}}=1^{\prime\prime}$, for
$\beta\lesssim 0.4$, $\hat{\beta}^{\textrm{A}}$ gives the best approximation
but for $\beta\gtrsim 0.4$, $\hat{\eta}$ gives the best approximation (figure
8). Thus, a hybrid use of $\hat{\beta}^{\textrm{A}}$ and $\hat{\eta}$ may be a
best way to estimate the distance ratio $\beta$ in the SIS + external shear
model.

###  4.2 SIS + SIS

Here we study lens models with a dominant SIS with an Einstein radius of
$\theta_{\textrm{E}}=1^{\prime\prime}$ and a subdominant SIS with an Einstein
radius of $\theta_{\textrm{E0}}\ll\theta_{\textrm{E}}$ apparently centred at
${\bm{x}}_{0}=(x_{1},x_{2})=(1^{\prime\prime},0)$. If a subdominant SIS
resides in the foreground of the dominant SIS, the effective deflection angle
of the perturber is

$\displaystyle\delta{\bm{\alpha}}_{\textrm{eff}}^{\textrm{f}}({\bm{x}})$
$\displaystyle=$
$\displaystyle\theta_{\textrm{E0}}\frac{{\bm{x}}-{\bm{x}}_{0}}{\lVert{\bm{x}}-{\bm{x}}_{0}\rVert}+\theta_{\textrm{E}}\biggl{(}\frac{{\bm{y}}_{\textrm{d}}}{y_{\textrm{d}}}-\frac{{\bm{x}}}{x}\biggr{)},$
$\displaystyle{\bm{y}}_{\textrm{d}}$ $\displaystyle=$
$\displaystyle{\bm{x}}-\beta\theta_{\textrm{E0}}\frac{{\bm{x}}-{\bm{x}}_{0}}{\lVert{\bm{x}}-{\bm{x}}_{0}\rVert}.$
(58)

If an SIS resides in the background of a dominant SIS, the effective
deflection angle of the perturber is

$\displaystyle\delta{\bm{\alpha}}_{\textrm{eff}}^{\textrm{b}}({\bm{x}})$
$\displaystyle=$
$\displaystyle\theta_{\textrm{E0}}\frac{{\bm{y}}_{\textrm{b}}-{\bm{y}}_{\textrm{b0}}}{\lVert{\bm{y}}_{\textrm{b}}-{\bm{y}}_{\textrm{b0}}\rVert},$
$\displaystyle{\bm{y}}_{\textrm{b}}$ $\displaystyle=$
$\displaystyle{\bm{x}}-\xi\theta_{\textrm{E}}\frac{{\bm{x}}}{x},$
$\displaystyle{\bm{y}}_{\textrm{b0}}$ $\displaystyle=$
$\displaystyle{\bm{x}}_{0}-\xi\theta_{\textrm{E}}\frac{{\bm{x}}_{0}}{x_{0}}.$
(59)

The effective deflection angles $\delta{\bm{\alpha}}_{\textrm{eff}}$ in an
SIS+SIS model with $\theta_{E}=1^{\prime\prime}$ and
$\theta_{E0}=0\aas@@fstack{\prime\prime}1$ are shown in figure 9. If a
subdominant SIS resides in the background of a dominant lens, the amplitude of
the effective deflection angle $\delta{\bm{\alpha}}_{\textrm{eff}}$ is
constant but the direction is anisotropic. The directions resemble those of an
electric field of a dipole. If a subdominant SIS resides in the foreground of
a dominant lens, the amplitude of the effective deflection angle
$\delta{\bm{\alpha}}_{\textrm{eff}}$ is not constant and the direction is
anisotropic. The streamlines resemble those in an SIS model with a background
SIS with the same Einstein radius except for those in the central region of
the dominant SIS (figure 10). For $\beta=0.5$ and $\xi=0.5$, the maximum
amplitudes of the astrometric lensing B-mode
$\delta{\bm{\alpha}}_{\textrm{eff}}^{\textrm{B}}$ are $\sim
0\aas@@fstack{\prime\prime}03$. As shown in figure 11, the contours of the
amplitude of lensing B-mode have a dumbbell-like shape with a spindle-shaped
void for both the foreground and background cases. In contrast, the contours
of the amplitude of lensing E-mode
$\delta{\bm{\alpha}}_{\textrm{eff}}^{\textrm{B}}$ have a complex structure
that depends on the position of the subdominant SIS. The amplitude is largest
in the central region of the dominant SIS in the foreground case whereas the
amplitude is largest in the central region of the subdominant SIS. This
implies that any fitting without considering the difference between the
subdominant lens plane and the dominant lens plane may lead to systematic
residual in the positions of quadruple images of a point-like source.

As shown in figure 12, the rotation of the effective deflection angle shows an
octopole pattern that consists of a pair of quadrupoles centred at the centres
of the SISs. The amplitudes of rotations are maximised in the diagonal lines
and minimised in the vertical and horizontal directions of each SIS. The
’directions’ of a rotation in the foreground SIS is opposite to that in the
background SIS.

In figure 13 and figure 14, we show the ratio of the B-mode to E-mode
$r^{\textrm{BE}}=|\delta{\bm{\alpha}}_{\textrm{eff}}^{\textrm{B}}|/|\delta{\bm{\alpha}}_{\textrm{eff}}^{\textrm{E}}|$,
averaged over an arc with $\theta=1.1\,$arcsec and $\phi=0.1\,$arcsec. We did
not use rotations in the Einstein ring of the dominant lens for averaging
because the rotations in such regions are extremely small especially in the
vicinity of the centre of the subdominant SIS. The ratio $r^{\textrm{BE}}$ is
a monotonically increasing function of $\beta$ or $\xi$, which equals $\sim
0.3$ for $\beta=0.5$, and $\sim 0.2$ for $\xi=0.5$. We found that the ratio of
rotation to divergence,
$|\nabla\times\delta{\bm{\alpha}}_{\textrm{eff}}^{\textrm{B}}|/|\nabla\cdot\delta{\bm{\alpha}}_{\textrm{eff}}^{\textrm{E}}|$
averaged over the arc are larger than $r^{\textrm{BE}}$. The analytic formula
in equation (17) gives a good estimate up to moderate values of $\beta\sim
0.5$ and $\xi\sim 0.5$.

Next, we show the values of the distance ratio estimators
$\hat{\beta}^{\textrm{A}}$, $\hat{\beta}^{\textrm{B}}$, and $\hat{\eta}$ in
the dominant lens plane in figure 15. We found that these estimators give a
good approximation in the vicinity of the subdominant SIS.
$\hat{\beta}^{\textrm{B}}$ and $\hat{\eta}$ give a worse result around the
Einstein ring of the dominant lens in which the rotation is almost zero. For
$\beta=0.5$, $\hat{\beta}^{\textrm{A}}$ gives the best result.

Finally, we show the effect of sample region, which was used to estimate the
distant ratio using $\hat{\beta}^{\textrm{A}}$ in figure 16. The center of a
subdominant SIS is fixed to $(x_{1},x_{2})=(1^{\prime\prime},0)$. We consider
two types of sample regions: 1) an arc with a radius of
$1^{\prime\prime}+\theta_{\textrm{Ep}}$ that subtends an azimuthal angle of
$\phi=\theta_{\textrm{Ep}}$, 2) an arc with a radius of
$1\aas@@fstack{\prime\prime}1$ that subtends an azimuthal angle of $\phi=0.1$
rad. Note that we avoided the neighbourhood of an Einstein ring of the
dominant lens as the rotation is very small. The errors in the case 1) is much
smaller than the case 2) as the sample region is much closer to the centre of
the subdominant SIS. For $\beta<0.5$, the relative errors in both 1) and 2)
were found to be $\lesssim 0.1$. Even for large distance ratios $\beta>0.5$,
the relative errors are $<0.1$ for
$\theta_{\textrm{Ep}}<0\aas@@fstack{\prime\prime}04$ in case 1). Thus, the
foreground distance ratio estimator $\hat{\beta}^{\textrm{A}}$ gives a
relatively accurate approximation if it is used for estimating the property of
the central region of a perturbing halo. Even if the lensed arc of an extended
source is not on the centre of a perturbing halo, $\hat{\beta}^{\textrm{A}}$
can be used as a good estimator if $\beta$ is sufficiently small.

Figure 9: Line-of-sight effect for an SIS+SIS model. The dominant lens is an
SIS with Einstein radius $\theta_{\textrm{E}}=1^{\prime\prime}$ centred at
$(0,0)$. The arrows show the effective deflection angle
$\delta{\bm{\alpha}}_{\textrm{eff}}$ due to another subdominant SIS with
$\theta_{\textrm{Ep}}=0\aas@@fstack{\prime\prime}1$ at $(1^{\prime\prime},0)$.
The colors in the middle panel show the amplitudes in arcsec.  Figure 10:
Streamlines of the effective deflection angles for an SIS. The model
parameters are the same as in figure 9.  Figure 11: Astrometric lensing B and
E-modes for an SIS. The arrows show the B and E components in the effective
deflection angle $\delta{\bm{\alpha}}_{\textrm{eff}}$. The model parameters
are the same as in figure 9. The colors show the amplitudes of the shift in
arcsec.  Figure 12: Rotation of effective deflection angle for an SIS with an
Einstein radius of $\theta_{\textrm{E0}}=0\aas@@fstack{\prime\prime}1$ centred
at $(1^{\prime\prime},0)$. The colors show the amplitudes of rotation. The
model parameters are the same as in figure 9.  Figure 13: The ratio of B and
E-modes for an SIS with an Einstein radius of
$\theta_{\textrm{E0}}=0\aas@@fstack{\prime\prime}1$ centred at
$(1^{\prime\prime},0)$. The orange full curves denote the ratios of the
averaged amplitude of B-mode astrometric shift and that of E-mode astrometric
shift as a function of $\beta$ (top) and $\xi$ (bottom). Both the modes are
averaged over an arc with $\theta=1.1\,$arcsec and $\phi=0.1\,$arcsec. The
blue dashed curves show the ratio of averaged rotation and averaged
convergence. Both the two contributions are averaged over the same arc. The
blue dot-dashed curves represent the approximations of the ratios for the
foreground case $\sqrt{2}\beta/4$ (left) and the background case
$\sqrt{2}\xi/4$ (right). The model parameters are the same as in figure 9.
Figure 14: The ratio of B and E-modes for an SIS with an Einstein radius of
$\theta_{\textrm{Ep}}=0\aas@@fstack{\prime\prime}03$ centred at
$(1^{\prime\prime},0)$. The lines are the same in figure 13. The model
parameters except for $\theta_{\textrm{Ep}}$ are the same as in figure 9.
Figure 15: Approximated distance ratios in the lens plane for an SIS with
$\beta=0.5$. The colors show the approximated distance ratios for
$\hat{\beta}^{\textrm{A}}$ (left), $\hat{\beta}^{\textrm{B}}$ (middle), and
$\hat{\eta}$ (right). Green curves denote the Einstein radius
$\theta_{\textrm{E}}$ of the dominant SIS and the blue curves represent the
Einstein radius $\theta_{\textrm{Ep}}$ of the subdominant SIS around the
centre. The model parameters are the same as in figure 9. Figure 16: The
effect of sample region for foreground SIS perturbers. The colors show the
common logarithm of mean relative errors of the distance ratio estimator
$\hat{\beta}^{\textrm{A}}$ as a function of $\beta$ and
$\theta_{\textrm{Ep}}$. The center of a subdominant SIS is fixed to
$(x_{1},x_{2})=(1^{\prime\prime},0)$. The mean is obtained from averaging the
relative errors on an arc with a radius of
$1^{\prime\prime}+\theta_{\textrm{Ep}}$ that subtends an azimuthal angle of
$\phi=\theta_{\textrm{Ep}}$ (left) and those on an arc with a radius of
$1\aas@@fstack{\prime\prime}1$ that subtends an azimuthal angle of $\phi=0.1$
rad (right).

##  5 Conclusion and Discussion

In this study, we investigated the characteristics of astrometric lensing
B-mode in strong lensing systems that consist of a dominant and a subdominant
lens residing at distinct redshifts. The B-mode arises from the coupling
between strong lensing induced by a dominant lens, and weak lensing generated
by a subdominant lens. By measuring both B and E-modes, we can deduce the
distance ratio, and ’bare’ convergence and shear perturbations attributed to
the subdominant lens.

In cases where a subdominant lens is located behind the dominant lens, we can
derive an exact formula if the dominant lens is perfectly modelled. However,
when the situation reverses, with a dominant lens resides behind a subdominant
lens, we cannot obtain an exact formula even with a perfect model of the
dominant lens. In such cases, we employ certain approximations that yield an
exact value when the distance between the dominant and subdominant lenses
approaches zero.

We demonstrated that any scale transformation in the distance ratio of a
subdominant lens corresponds to a mass-sheet transformation in the background
lens plane. Consequently, determining the distance ratio necessitates
assumptions about the values of the mass-sheet within the background lens
plane. Nevertheless, if we measure time delays between perturbed multiple
lensed images and know the redshifts of a subdominant lens, we can break the
mass-sheet degeneracy for a given $H_{0}$, enabling us to determine it without
uncertainty. Moreover, if the redshifts of both foreground and background
perturbers, as well as the dominant lens and the source, are known,
observations of astrometric lensing B-mode can break the mass-sheet
degeneracies. This would lead to a reduction in systematic errors in the
estimated value of $H_{0}$.

Our analysis focuses on systems with a single subdominant lens whose
deflection angle is significantly smaller than that of the dominant lens. In
reality, the gravitational influence of multiple subdominant lenses along each
photon’s path must be considered. When both the second and third dominant
lenses are located either in the foreground or the background of the dominant
lens, the rotation signal may be amplified. Conversely, if they are positioned
separately in the foreground and background, the signal could weaken due to
cancellation. Investigating these effects in systems comprising three or more
lenses falls beyond the scope of this paper.

Assuming that the impact of model degeneracy due to the extended multi-plane
mass-sheet transformation (eMMST) we discussed is negligible, the measurement
of astrometric lensing B-modes holds the potential to constrain the abundance
of intergalactic dark haloes with masses of $\lesssim 10^{8}\,M_{\odot}$ along
the line of sight (LOS) in quasar-galaxy or galaxy-galaxy strong lens systems.
Our toy models suggest that such feasibility can be assessed as follows: As a
reference lens model, we employ a dominant SIS with an Einstein radius of
$1^{\prime\prime}$. Then, a subdominant SIS with an Einstein radius of
$\gtrsim 0\aas@@fstack{\prime\prime}03$ in the LOS would produce B-modes with
a shift of $\gtrsim 0\aas@@fstack{\prime\prime}003$ for $\beta>0.2$ or
$\xi>0.2$. These shifts can be observed with telescopes featuring an angular
resolution of $\lesssim 0\aas@@fstack{\prime\prime}03$, assuming a typical
magnification of $\mu\sim 10$ for lensed image separations of $\lesssim
0\aas@@fstack{\prime\prime}5$. Hence, instruments like the Atacama Large
Millimeter/Submillimeter Array (ALMA) possess the capability to detect
astrometric lensing B-modes resulting from less massive LOS haloes, provided
that the distance between the dominant and subdominant lenses is sufficiently
large and the signal-to-noise ratio of intensity in the lens plane is suitably
high.

In the near future, we plan to investigate the practicality of measuring
astrometric lensing B-modes with ALMA, utilizing more sophisticated models and
taking into account observational capabilities.

## References

  * Bar-Kana (1996) Bar-Kana R., 1996, ApJ, 468, 17
  * Birrer et al. (2017) Birrer S., Welschen C., Amara A., Refregier A., 2017, J. Cosmology Astropart. Phys., 2017, 049
  * Brooks et al. (2017) Brooks A. M., Papastergis E., Christensen C. R., Governato F., Stilp A., Quinn T. R., Wadsley J., 2017, ApJ, 850, 97
  * Chantry et al. (2010) Chantry V., Sluse D., Magain P., 2010, A&A, 522, A95
  * Chatterjee & Koopmans (2018) Chatterjee S., Koopmans L. V. E., 2018, MNRAS, 474, 1762
  * Chiba (2002) Chiba M., 2002, ApJ, 565, 17
  * Dalal & Kochanek (2002) Dalal N., Kochanek C. S., 2002, ApJ, 572, 25
  * Despali et al. (2018) Despali G., Vegetti S., White S. D. M., Giocoli C., van den Bosch F. C., 2018, MNRAS, 475, 5424
  * Erdl & Schneider (1993) Erdl H., Schneider P., 1993, A&A, 268, 453
  * Evans & Witt (2003) Evans N. W., Witt H. J., 2003, MNRAS, 345, 1351
  * Fielder et al. (2019) Fielder C. E., Mao Y.-Y., Newman J. A., Zentner A. R., Licquia T. C., 2019, MNRAS, 486, 4545
  * Fleury et al. (2021) Fleury P., Larena J., Uzan J.-P., 2021, arXiv e-prints, p. arXiv:2104.08883
  * Gilman et al. (2017) Gilman D., Agnello A., Treu T., Keeton C. R., Nierenberg A. M., 2017, MNRAS, 467, 3970
  * Hezaveh et al. (2016) Hezaveh Y. D., et al., 2016, ApJ, 823, 37
  * Hsueh et al. (2017) Hsueh J. W., et al., 2017, MNRAS, 469, 3713
  * Hsueh et al. (2018) Hsueh J.-W., Despali G., Vegetti S., Xu D., Fassnacht C. D., Metcalf R. B., 2018, MNRAS, 475, 2438
  * Inoue (2016) Inoue K. T., 2016, MNRAS, 461, 164
  * Inoue & Chiba (2003) Inoue K. T., Chiba M., 2003, ApJ, 591, L83
  * Inoue & Chiba (2005a) Inoue K. T., Chiba M., 2005a, ApJ, 633, 23
  * Inoue & Chiba (2005b) Inoue K. T., Chiba M., 2005b, ApJ, 634, 77
  * Inoue & Takahashi (2012) Inoue K. T., Takahashi R., 2012, MNRAS, 426, 2978
  * Inoue et al. (2016) Inoue K. T., Minezaki T., Matsushita S., Chiba M., 2016, MNRAS, 457, 2936
  * Inoue et al. (2023) Inoue K. T., Minezaki T., Matsushita S., Nakanishi K., 2023, ApJ, 954, 197
  * Kauffmann et al. (1993) Kauffmann G., White S. D. M., Guiderdoni B., 1993, MNRAS, 264, 201
  * Keeton et al. (2003) Keeton C. R., Gaudi B. S., Petters A. O., 2003, ApJ, 598, 138
  * Klypin et al. (1999) Klypin A., Kravtsov A. V., Valenzuela O., Prada F., 1999, ApJ, 522, 82
  * Koopmans (2005) Koopmans L. V. E., 2005, MNRAS, 363, 1136
  * Lewis & Challinor (2006) Lewis A., Challinor A., 2006, Phys. Rep., 429, 1
  * Mao & Schneider (1998) Mao S., Schneider P., 1998, MNRAS, 295, 587
  * McCully et al. (2014) McCully C., Keeton C. R., Wong K. C., Zabludoff A. I., 2014, MNRAS, 443, 3631
  * McCully et al. (2017) McCully C., Keeton C. R., Wong K. C., Zabludoff A. I., 2017, ApJ, 836, 141
  * Metcalf (2005) Metcalf R. B., 2005, ApJ, 629, 673
  * Metcalf & Madau (2001) Metcalf R. B., Madau P., 2001, ApJ, 563, 9
  * Moore et al. (1999) Moore B., Ghigna S., Governato F., Lake G., Quinn T., Stadel J., Tozzi P., 1999, ApJ, 524, L19
  * Nashimoto et al. (2022) Nashimoto M., Tanaka M., Chiba M., Hayashi K., Komiyama Y., Okamoto T., 2022, ApJ, 936, 38
  * Oguri (2005) Oguri M., 2005, MNRAS, 361, L38
  * Schneider (2014) Schneider P., 2014, A&A, 568, L2
  * Schneider (2019) Schneider P., 2019, A&A, 624, A54
  * Sengül et al. (2022) Sengül A. Ç., Dvorkin C., Ostdiek B., Tsang A., 2022, MNRAS, 515, 4391
  * Takahashi & Inoue (2014) Takahashi R., Inoue K. T., 2014, MNRAS, 440, 870
  * Treu & Koopmans (2004) Treu T., Koopmans L. V. E., 2004, ApJ, 611, 739
  * Vegetti & Koopmans (2009) Vegetti S., Koopmans L. V. E., 2009, MNRAS, 392, 945
  * Vegetti et al. (2010) Vegetti S., Koopmans L. V. E., Bolton A., Treu T., Gavazzi R., 2010, MNRAS, 408, 1969
  * 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, 341
  * Vegetti et al. (2014) Vegetti S., Koopmans L. V. E., Auger M. W., Treu T., Bolton A. S., 2014, MNRAS, 442, 2017
  * Wetzel et al. (2016) Wetzel A. R., Hopkins P. F., Kim J.-h., Faucher-Giguère C.-A., Keres D., Quataert E., 2016, ApJ, 827, L23
  * Xu et al. (2009) Xu D. D., et al., 2009, MNRAS, 398, 1235
  * Xu et al. (2010) Xu D. D., Mao S., Cooper A. P., Wang J., Gao L., Frenk C. S., Springel V., 2010, MNRAS, 408, 1721
  * Xu et al. (2012) Xu D. D., Mao S., Cooper A. P., Gao L., Frenk C. S., Angulo R. E., Helly J., 2012, MNRAS, 421, 2553
  * Zaldarriaga & Seljak (1998) Zaldarriaga M., Seljak U., 1998, Phys. Rev. D, 58, 023003
  * Çagan Şengül et al. (2020) Çagan Şengül A., Tsang A., Diaz Rivero A., Dvorkin C., Zhu H.-M., Seljak U., 2020, Phys. Rev. D, 102, 063502

##  Appendix A Analytic Solutions for convergence and shear perturbations

In the following, we present analytic solutions for approximated ’bare’
convergence perturbation $\delta\kappa$ and shear perturbations
$\delta\gamma_{1},\delta\gamma_{2}$ for a foreground perturber, which are
valid for $\beta\ll 1$. We assume that the magnification matrix of the
dominant lens can be lineally approximated by three constants $c_{1},c_{2},$
and $c_{3}$ at ${\bm{x}}$ as described in equation (LABEL:eq:c1c2c3). Plugging
an approximated solution $\hat{\beta}^{\textrm{A}}$ in equation (21) into
equation (19), we have

$\displaystyle\delta\kappa$ $\displaystyle=$
$\displaystyle-K^{-1}[2(-1+\hat{\beta}^{\textrm{A}}\kappa)(\hat{\beta}^{\textrm{A}}\gamma_{2}(B+\hat{\beta}^{\textrm{A}}c_{3})-(A+\hat{\beta}^{\textrm{A}}c_{1})$
$\displaystyle\times$
$\displaystyle(-1+\hat{\beta}^{\textrm{A}}(\gamma_{1}+\kappa)))-(A+\hat{\beta}^{\textrm{A}}c1-\hat{\beta}^{\textrm{A}}c_{2}-D)$
$\displaystyle\times$
$\displaystyle(-(\hat{\beta}^{\textrm{A}})^{2}\gamma_{2}^{2}-(-1+\hat{\beta}^{\textrm{A}}(\gamma_{1}+\kappa))^{2})],$
$\displaystyle K$ $\displaystyle\equiv$ $\displaystyle
2(-1+\hat{\beta}^{\textrm{A}}(\gamma_{1}+\kappa))(-1+2\hat{\beta}^{\textrm{A}}\kappa+(\hat{\beta}^{\textrm{A}})^{2}(\gamma_{1}^{2}+\gamma_{2}^{2}-\kappa^{2})),$
$\displaystyle\delta\gamma_{1}$ $\displaystyle=$ $\displaystyle
L^{-1}(2(B+\hat{\beta}^{\textrm{A}}c_{3})\gamma_{1}-(A\hat{\beta}^{\textrm{A}}c_{1}-\hat{\beta}^{\textrm{A}}c_{2}-D)\gamma_{2})$
$\displaystyle+$ $\displaystyle
M^{-1}[2\gamma_{1}(-((B+\hat{\beta}^{\textrm{A}}c_{3})(1+\hat{\beta}^{\textrm{A}}\gamma_{1}))-\hat{\beta}^{\textrm{A}}\gamma_{2}(\hat{\beta}^{\textrm{A}}c_{2}+D)$
$\displaystyle+$
$\displaystyle\hat{\beta}^{\textrm{A}}\kappa(B+\hat{\beta}^{\textrm{A}}c_{3}))],$
$\displaystyle L$ $\displaystyle\equiv$ $\displaystyle
2\gamma_{2}(-1+\hat{\beta}^{\textrm{A}}(\gamma_{1}+\kappa)),$ $\displaystyle
M$ $\displaystyle\equiv$ $\displaystyle
2\gamma_{2}(-1+2\hat{\beta}^{\textrm{A}}\kappa+(\hat{\beta}^{\textrm{A}})^{2}(\gamma_{1}^{2}+\gamma_{2}^{2}-\kappa^{2})),$
(61)

and

$\displaystyle\delta\gamma_{2}$ $\displaystyle=$ $\displaystyle
N^{-1}[B+\hat{\beta}^{\textrm{A}}(c_{3}+(\hat{\beta}^{\textrm{A}}c_{2}+D)\gamma_{2}+\hat{\beta}^{\textrm{A}}c_{3}(\gamma_{1}-\kappa))$
$\displaystyle+$ $\displaystyle
B\hat{\beta}^{\textrm{A}}(\gamma_{1}-\kappa)],$ $\displaystyle N$
$\displaystyle\equiv$ $\displaystyle
1-\hat{\beta}^{\textrm{A}}(2\kappa+\hat{\beta}^{\textrm{A}}(\gamma_{1}^{2}+\gamma_{2}^{2}-\kappa^{2})).$
(62)


[a]Sandip Maiti

# Three-dimensional Gross-Neveu model with two flavors of staggered fermions

Debasish Banerjee  Shailesh Chandrasekharan  Marina Krstic Marinkovic

###### Abstract

We introduce a strongly interacting lattice field theory model containing two
flavors of massless staggered fermions with two kind of interactions: (1) a
lattice current-current interaction, and (2) an on-site four-fermion
interaction. At weak couplings, we expect a massless fermion phase since our
interactions become irrelevant at long distances. At strong couplings, based
on previous studies, we argue that our lattice model contains two different
massive fermion phases with different mechanisms of fermion mass generation.
In one phase, fermions become massive through Spontaneous Symmetry Breaking
(SSB) via the formation of a fermion bilinear condensate. In the other phase,
fermion mass arises through a more exotic mechanism without the formation of
any fermion bilinear condensate. Our lattice model is free of sign problems
and can be studied using the fermion bag algorithm. The longer term goal here
is to study both these mass generation phenomena in a single model and
understand how different phases come together.

##  1 Introduction

The traditional approach for massless fermions to acquire mass is through the
formation of fermion bilinear condensates. In the standard model of particle
physics, this occurs through the mechanism of spontaneous symmetry breaking
(SSB), which gives masses to quarks and leptons [1]. In the electro-weak
sector of the Standard Model, the Higgs field plays a central role in
generating masses for fermions through Yukawa couplings [2]. Due to SSB of the
electro-weak symmetry, the Higgs field gets an expectation value, which gives
rise to fermion bilinear mass terms in the low energy effective theory. In the
strong interaction sector, spontaneous breakdown of chiral symmetries is
responsible to create a fermion bilinear quark condensate. On the other hand,
studies of lattice Yukawa models have found exotic phases at strong couplings
[3, 4], where fermions seem to acquire mass without breaking of any symmetries
of the theory. Such phases have been called the paramagnetic strong phase
(PMS-phase) to contrast it with the phase at weak couplings where fermions
remain massless due to lack of any symmetry breaking, which is referred to as
the paramagnetic weak phase (PMW-phase). In the PMS phase, all fermion
bilinear condensates vanish, but fermions are still massive [5, 6]. The
traditional spontaneously broken phase is called the ferromagnetic phase (FM-
phase).

In this work, we introduce a new lattice field theory model with four-fermion
interactions to explore all these three phases under one framework. An
important feature of our model is that it is free from sign problems and can
be studied using the fermion bag algorithms without introducing a bare mass
for the fermions [7]. We focus on $2+1$ dimensions, since the four-fermion
interactions are known to produce many interesting quantum critical points.
For example, in the so-called lattice Thirring model with a current-current
interaction $U$, an interesting second order phase transition between the PMW
and FM phases was observed even with a single flavor of staggered fermions [8,
9]. With two flavors of staggered fermions but with a single site four-fermion
interaction $U^{\prime}$, an exotic second order transition between the PMW
and PMS phases was observed [10, 9]. Our model brings together both of these
transitions into a single phase diagram by introducing a model with both
interactions $U$ and $U^{\prime}$, thus giving us a way to study how the two
mechanisms for fermion mass generation can come together and if there are
interesting new phase transitions between them.

Figure 1: Two possible phase diagrams of our model, containing two types of
interactions $U$ and $U^{\prime}$.

The possible phase diagrams in our model are shown in Fig. 1. The left phase
diagram implies that $U=0$ is a special axis where, due to enhanced
symmetries, a direct second order phase transition between PMW and PMS phase
was possible. Since this generically not possible, introduction of $U$
introduces an intermediate phase. On the other hand, the right phase diagram
implies that the direct phase transition between the PMW and PMS phases is
more robust and is not a consequence of the special symmetries of the $U=0$
axis. Other interesting questions are related to the nature of the phase
transitions between the PMW and FM phases and the PMS and FM phases. If they
are second order, we can compute the critical exponents of these transitions.
The nature of the multi-critical point where all the three phases meet is
another possible subject for investigation.

##  2 Model

We study a $2+1$ dimensional lattice model containing two flavors of massless
staggered fermions interacting with each other via four-fermion interactions.
The Euclidean action of our model is given as a sum of three terms
$S=S_{0}+S^{\text{int}}_{1}+S^{\text{int}}_{2}$, where $S_{0}$ denotes the
free staggered fermion action for each of the flavors $u$ and $d$ and is given
by

$\displaystyle S_{0}$ $\displaystyle=\sum_{\langle
x,y\rangle}\left(\bar{u}_{x}M_{x,y}u_{y}+\bar{d}_{x}M_{x,y}d_{y}\right),$ (1)

$S^{\text{int}}_{1}$ is a nearest neighbor current-current interaction term
within each flavor

$\displaystyle S^{\text{int}}_{1}$ $\displaystyle=-U\sum_{\langle
x,y\rangle}\left(\bar{u}_{x}u_{x}\bar{u}_{y}u_{y}+\bar{d}_{x}d_{x}\bar{d}_{y}d_{y}\right),$
(2)

and $S^{\text{int}}_{2}$ is a single site interaction between the two flavors
given by

$\displaystyle S^{\text{int}}_{2}$
$\displaystyle=-U^{\prime}\sum_{x}\left(\bar{u}_{x}u_{x}\bar{d}_{x}d_{x}\right).$
(3)

The matrix $M$ is the massless staggered fermion matrix defined by

$\displaystyle M_{x,y}$
$\displaystyle=\sum\limits_{\hat{\mu}}\frac{\eta_{x,\hat{\mu}}}{2}(\delta_{x+\hat{\mu},y}-\delta_{x-\hat{\mu},y}),$
(4)

where $x\equiv(x_{0},x_{1},x_{2})$ denotes a lattice site on a 3 dimensional
cubic lattice and $\hat{\mu}$ represent unit vectors in the three directions.
The $\eta_{x,\hat{\mu}}$ are the usual staggered phase factors defined as:
$\eta_{x,\hat{0}}=1$, $\eta_{x,\hat{1}}=(-1)^{x_{0}}$ and
$\eta_{x,\hat{2}}=(-1)^{x_{0}+x_{1}}$. When $U^{\prime}=0$, the action $S$
represents two decoupled single flavored lattice Thirring models that has been
studied earlier [8]. The coupling $U^{\prime}$ couples the two flavors with
each other. When $U=0$, the lattice model is the same as the one that was
studied in Ref. [10].

Our lattice model has a rich internal symmetry structure. The free fermion
action $S_{0}$ is symmetric under $SU(4)\times U(1)$ transformations, while
$S^{\text{int}}_{1}$ is invariant under $SU(2)\times U(1)\times SU(2)\times
U(1)$ and $S^{\text{int}}_{2}$ is symmetric under $SU(4)$ transformation.
Hence, the action $S$ of our model is invariant under $SU(2)\times SU(2)\times
U(1)$ transformations. The observables we wish to measure include local four-
fermion condensates given by

$\displaystyle N_{u}=\frac{2}{V}\sum_{\langle
x,y\rangle}\Big{\langle}\bar{u}_{x}u_{x}\bar{u}_{y}u_{y}\Big{\rangle},\quad
N_{d}=\frac{2}{V}\sum_{\langle
x,y\rangle}\left\langle\bar{d}_{x}d_{x}\bar{d}_{y}d_{y}\right\rangle,\quad
N_{i}=\frac{1}{V}\sum_{x}\left\langle\bar{u}_{x}u_{x}\bar{d}_{x}d_{x}\right\rangle,$
(5)

where the sum $\langle x,y\rangle$ is over nearest neighbor bonds. Similarly,
fermion bilinear susceptibilities given by

$\displaystyle\chi_{uu}=\frac{1}{V}\sum_{x,y}\Big{\langle}\bar{u}_{x}u_{x}\bar{u}_{y}u_{y}\Big{\rangle},\quad\chi_{dd}=\frac{1}{V}\sum_{x,y}\Big{\langle}\bar{d}_{x}d_{x}\bar{d}_{y}d_{y}\Big{\rangle},\quad\chi_{ud}=\frac{1}{V}\sum_{x,y}\Big{\langle}\bar{u}_{x}u_{x}\bar{d}_{y}d_{y}\Big{\rangle}.$
(6)

The symmetries of our model imply that $N_{u}=N_{d}$ and
$\chi_{uu}=\chi_{dd}$. In the above expressions the expectation values are
defined as

$\displaystyle\left\langle\mathscr{O}\right\rangle$
$\displaystyle=\frac{1}{Z}\int[\mathcal{D}\bar{u}\mathcal{D}u\mathcal{D}\bar{d}\mathcal{D}d]\hskip
2.84544pt\mathscr{O}\hskip 2.84544pte^{-S[\bar{u},u,\bar{d},d]}$ (7)

with $Z$ being the partition function. The susceptibilities $\chi_{uu}$ and
$\chi_{ud}$ diverge as $V$ in the FM-phase with non-zero fermion bilinear
condensate but saturate in both the PMW phase and the PMS phase.

##  3 Fermion Bag Approach

The traditional approach to solve four-fermion models is by introducing an
auxiliary field, where we can use the Hubbard-Stratonovich transformation to
convert the four-fermion coupling into a fermion bilinear by introducing an
auxiliary bosonic field [11]. Methods of gradient flow, generalizing the
Lefschetz thimble method to the lattice Thirring model at finite density, have
also been recently explored [12, 13]. In this work, we use an alternate
approach called the Fermion Bag approach [7]. The basic idea of this alternate
method is to regroup fermion worldlines inside regions denoted as fermion
bags. The partition function can then be expressed as a sum over
configurations of fermion bags. However, since there are several ways to
regroup fermion worldlines, the approach is not unique and depending on the
physics the regrouping must be done thoughtfully. When this is possible, this
approach can be more efficient than the auxiliary field approach.

For our model, several possible ways to identify fermion bags. The first step
involves expanding the partition function

$\displaystyle
Z=\int[\mathcal{D}\bar{u}\mathcal{D}u\mathcal{D}\bar{d}\mathcal{D}d]\
e^{-\sum\limits_{x,y}(\bar{u}_{x}M_{x,y}u_{y}+\bar{d}_{x}M_{x,y}d_{y})+U\sum\limits_{x,y}(\bar{u}_{x}u_{x}\bar{u}_{y}u_{y}+\bar{d}_{x}d_{x}\bar{d}_{y}d_{y})+U^{\prime}\sum\limits_{x}\bar{u}_{x}u_{x}\bar{d}_{x}d_{x}}$
(8)

as a sum over powers of interactions, as is done in bare perturbation theory.
This is done by expanding

$\displaystyle
e^{U\bar{u}_{x}u_{x}\bar{u}_{y}u_{y}}=(1+U\bar{u}_{x}u_{x}\bar{u}_{y}u_{y})=\sum\limits_{b_{x,y}=0,1}(U\bar{u}_{x}u_{x}\bar{u}_{y}u_{y})^{b_{x,y}}$
(9) $\displaystyle
e^{U^{\prime}\bar{u}_{x}u_{x}\bar{d}_{x}d_{x}}=(1+U^{\prime}\bar{u}_{x}u_{x}\bar{d}_{x}d_{x})=\sum\limits_{i_{x}=0,1}(U^{\prime}\bar{u}_{x}u_{x}\bar{d}_{x}d_{x})^{i_{x}}.$
(10)

We can introduce a new set of variables representing the interactions, which
in our case would be _dimers_ or _bonds_ $b_{x,y}$ for nearest neighbor
interactions and _instantons_ $i_{x}$ for single site interactions. Thus,
every configuration $[b,i]$ divides all lattice sites into either dimers $[b]$
and $[i]$. Sites that do not belong to a dimer or an instanton will be
referred to as free sites. The latter does not have any instantons or bonds.
Figure 2 illustrates a configuration on a $2^{3}$ lattice.

Figure 2:  An illustration of a possible configuration $[b,i]$ on a $2^{3}$
lattice. The red circles denote the instanton sites and red links denote
dimers. Dimers on $u$ and $d$ layers can occur on different sites, while
instanton positions on both the layers are same and can be visualized as bonds
connecting the two layers. The black circles denote free sites and do not
contain any instantons or dimers.

One way to identify fermion bags in our model would be to group Grassmann
variables that appear at every dimer or instanton and integrate them out. This
gives us factors of $U$ or $U^{\prime}$ for every bond and instanton. The
remaining free sites then form their own group of sites and are integrated out
separately. The partition function can then be written as a sum over all
configurations of $[b,i]$ as [14],

$\displaystyle Z$
$\displaystyle=\sum_{[b,i]}{U^{\prime}}^{N_{i}}U^{N_{u}+N_{d}}{\rm
det}(W_{u})~{}~{}{\rm det}(W_{d})$ (11)

where $N_{i}$ represents the number of instantons, while $N_{u}$ and $N_{d}$
represent the number of $u$-dimers and $d$-dimers in the configuration. The
matrices $W_{u}$ and $W_{d}$ are free staggered fermion matrix restricted to
the free sites. We call this approach strong coupling fermion bag since in
this approach the sizes of the matrices $W_{u}$ and $W_{d}$ will be small.

Figure 3: The variation of the average bond density ($N_{u}$) and the average
instanton density ($N_{i}$) with couplings $U$ and $U^{\prime}$ for a lattice
of linear dimension 2.

An alternate way of constructing the fermion bag is to consider each
configuration $[b,i]$ as a term in the perturbative expansion and write the
partition function as

$\displaystyle Z=\sum_{[b,i]}U^{{}^{\prime}N_{i}}U^{N_{u}+N_{d}}$
$\displaystyle\Bigg{\\{}\int[\mathcal{D}\bar{u}\mathcal{D}u]e^{-\sum\limits_{x,y}(\bar{u}_{x}M_{x,y}u_{y})}\bar{u}_{z_{1}}u_{z_{1}}...\bar{u}_{z_{k}}u_{z_{k}}\Bigg{\\}}$
$\displaystyle\times\Bigg{\\{}\int[\mathcal{D}\bar{d}\mathcal{D}d]e^{-\sum\limits_{x,y}(\bar{d}_{x}M_{x,y}d_{y})}\bar{d}_{w_{1}}d_{w_{1}}...\bar{d}_{w_{\ell}}d_{w_{\ell}}\Bigg{\\}}.$
(12)

In this expression $k=2N_{u}+N_{i}$ and $\ell=2N_{d}+N_{i}$. The terms in the
brackets in Eq. (12) can be computed using Wick’s theorem and summed over all
contractions. This yields a simple result

$\displaystyle Z$ $\displaystyle=[\mathrm{det}(M)]^{2}\
\sum_{[b,i]}{U^{\prime}}^{N_{i}}U^{N_{u}+N_{d}}{\rm det}(G_{u})~{}~{}{\rm
det}(G_{d})$ (13)

where $G_{u}$ and $G_{d}$ are $k\times k$ and $\ell\times\ell$ propagator
matrix between the sites $z_{1},...z_{k}$ and $w_{1},...,w_{\ell}$
respectively. We call this approach the weak coupling fermion bag, since the
sizes of the matrices $G_{u}$ and $G_{d}$ will be small at weak couplings. For
the efficiency of the fermion bag algorithm, it is often useful to switch
between the two fermion bag formulations.

Figure 4: The behavior of $\chi_{uu}$ and $\chi_{ud}$ as a function of the
couplings $U$ and $U^{\prime}$.

##  4 Results

We are currently developing the fermion bag Monte Carlo algorithm to explore
the phase diagram of our model. In order to confirm the correctness of the
algorithm, it would be useful to compare the results with exact calculations.
Here we show the results from these exact calculations on a $2^{3}$ lattice
size with anti-periodic boundary conditions in all three directions. The
behavior of the observables $N_{u}$ and $N_{i}$ as a function of the couplings
$U$ and $U^{\prime}$ are shown in Fig. 3. These four-point condensates are
smooth functions, increasing from 0 for small couplings, and approaching 1 for
large couplings.

To explore the phase diagram of our model further, it is more natural to study
the susceptibilities $\chi_{uu}$ and $\chi_{dd}$. These are plotted as
functions of the couplings $U$ and $U^{\prime}$ in Fig. 4. In the context of
phase transitions, a peak in the plot of the susceptibilities as a function of
the coupling for a fixed system size, can indicate a phase transition in the
system. However, we must also show that the peak actually diverges in the
thermodynamic limit to confirm the existence of the phase transition. As can
be seen in Fig. 4, both susceptibilities show peaks. It is important to point
out the fact that for very large couplings both susceptibilities as we have
defined it becomes small and is related to the fact that the lattice is being
filled with bonds and instantons. This however does not imply that the
condensate vanishes. To measure the true condensate, we will need to study the
dependence on the system size. We expect this to be dramatically different, at
large U and small $U^{\prime}$ as compared to large $U^{\prime}$ and small
$U$. The former will show that $\chi_{uu}$ will grow as a function of the
system size while the latter will saturate. These can be inferred from
previous studies.

##  5 Conclusions and Future work

We have introduced a three-dimensional Gross-Neveu model with two couplings
$U$ and $U^{\prime}$ such that two mechanisms of fermion mass generation are
possible and can compete. We have shown some exact results for $2^{3}$ lattice
and already see some rudimentary features of the phase diagram. In order to
understand the physics of our model in more detail, we plan to extend our
results for larger lattice size using the fermion bag algorithm. We hope to
explore the phase diagram of our model, and study the nature of the phase
transitions. If some of these are second order, we plan to calculate the
critical exponents.

## References

  * [1] Heather E. Logan.  TASI 2013 lectures on Higgs physics within and beyond the Standard Model.  6 2014.
  * [2] Antonio Pich.  The Standard Model of Electroweak Interactions.  In 2010 European School of High Energy Physics, pages 1–50, 1 2012.
  * [3] Sinya Aoki, I-Hsiu Lee, and Robert E Shrock.  Charged fermion correlation functions in a lattice scalar-fermion model with global chiral U(1) symmetry.  Nuclear Physics B, 388(1):229–242, 1992.
  * [4] I-Hsiu Lee, Junko Shigemitsu, and Robert E. Shrock.  Study of different lattice formulations of a Yukawa model with a real scalar field.  Nuclear Physics B, 334(1):265–278, 1990.
  * [5] Simon Catterall.  Fermion mass without symmetry breaking.  JHEP, 01:121, 2016.
  * [6] Ana Hasenfratz, Wei-qiang Liu, and Thomas Neuhaus.  Phase Structure and Critical Points in a Scalar Fermion Model.  Phys. Lett. B, 236:339–343, 1990.
  * [7] Shailesh Chandrasekharan.  The Fermion bag approach to lattice field theories.  Phys. Rev. D, 82:025007, 2010.
  * [8] Shailesh Chandrasekharan and Anyi Li.  Fermion bags, duality and the three dimensional massless lattice Thirring model.  Phys. Rev. Lett., 108:140404, 2012.
  * [9] Shailesh Chandrasekharan and Anyi Li.  Quantum critical behavior in three dimensional lattice Gross-Neveu models.  Phys. Rev. D, 88:021701, 2013.
  * [10] Venkitesh Ayyar and Shailesh Chandrasekharan.  Massive fermions without fermion bilinear condensates.  Phys. Rev. D, 91(6):065035, 2015.
  * [11] F. F. Assaad, M. Bercx, F. Goth, A. Götz, J. S. Hofmann, E. Huffman, Z. Liu, F. Parisen Toldin, J. S. E. Portela, and J. Schwab.  The ALF (Algorithms for Lattice Fermions) project release 2.0. Documentation for the auxiliary-field quantum Monte Carlo code.  12 2020.
  * [12] Andrei Alexandru, Gokce Basar, Paulo F. Bedaque, Gregory W. Ridgway, and Neill C. Warrington.  Monte Carlo calculations of the finite density Thirring model.  Phys. Rev. D, 95(1):014502, 2017.
  * [13] Andrei Alexandru, Paulo F. Bedaque, Henry Lamm, Scott Lawrence, and Neill C. Warrington.  fermions at finite density in $2+1$ dimensions with sign-optimized manifolds.  Phys. Rev. Lett., 121:191602, Nov 2018.
  * [14] Venkitesh Ayyar, Shailesh Chandrasekharan, and Jarno Rantaharju.  Benchmark results in the 2D lattice Thirring model with a chemical potential.  Phys. Rev. D, 97:054501, Mar 2018.