depend on the details of the models we have considered. In order to preserve
generality, we shall present this calculation in $d$ spatial dimensions. When
$d=2$, this is the approximation that is used to analyze the interacting
Hamiltonians in Eq. 15 and Eq. 68 in the main body of the text. When $d=3$,
this is the approximation that is used to analyze the interacting Hamiltonians
in Eq. 40 and Eq. 60 in the main body of the text.

To start, we will rewrite the Hubbard-Stratonovich fields as

$\begin{split}&\lambda_{1x}=\lambda_{0}\cos(\nu)\cos(\theta_{x})e^{i\phi_{1x}},\\\
&\lambda_{2x}=\lambda_{0}\cos(\nu)\sin(\theta_{x})e^{i\phi_{2x}},\\\
&\lambda_{1y}=\lambda_{0}\sin(\nu)\cos(\theta_{y})e^{i\phi_{1y}},\\\
&\lambda_{2y}=\lambda_{0}\sin(\nu)\sin(\theta_{y})e^{i\phi_{2y}},\\\
\end{split}$ (74)

where $\lambda_{0}$ is real and positive, and we have suppressed the
dependence on the lattice coordinates for simplicity. In terms of these
fields, the Lagrangian can be written as

$\begin{split}\mathcal{L}=&\bm{c}^{\dagger}G_{f}^{-1}\bm{c}\\\
&+\frac{2}{V}\lambda^{2}_{0}\cos(\phi_{1x}+\phi_{2x}-A_{xy})\cos^{2}(\nu)\sin(2\theta_{x})\\\
&+\frac{2}{V}\lambda^{2}_{0}\cos(\phi_{1y}+\phi_{2y}-A_{xy})\sin^{2}(\nu)\sin(2\theta_{y})],\end{split}$
(75)

where $G_{f}^{-1}$ is the matrix which contains all terms that are quadratic
in the fermion operators. Due to the Hubbard-Stratonovich decomposition of the
ring exchange interaction, $G_{f}^{-1}$ depends on the Hubbard Stratonovich
fields in Eq. 74.

To proceed, we will make use of the fact that the self-consistent values of
the Hubbard-Stratonovich fields are also those that minimize the effective
potential for the Hubbard-Stratonovich fields, once the fermions are
integrated out. To find the effective potential for the Hubbard-Stratonovich
fields, we will assume that they acquire a uniform expectation value. In this
approximation, the Lagrangian is quadratic in the fermion operators, and
$G^{-1}_{f}$ can be written diagonally in momentum and frequency space,
$G^{-1}_{f}\rightarrow G^{-1}_{f}(\omega,k)$. Using this, the resulting
effective potential for the Hubbard Stratonovich fields can be written as

$\begin{split}\mathcal{H}_{eff}=-&\log[\text{Det}(G_{f}^{-1})]\\\
&-\frac{2}{V}\lambda^{2}_{0}\cos(\phi_{1x}+\phi_{2x}-A_{xy})\cos^{2}(\nu)\sin(2\theta_{x})\\\
&-\frac{2}{V}\lambda\cos(\phi_{1y}+\phi_{2y}-A_{xy})\sin^{2}(\nu)\sin(2\theta_{y}),\end{split}$
(76)

where

$\text{Det}(G_{f}^{-1})\equiv\int\frac{d\omega
d^{d}k}{(2\pi)^{d+1}}\text{Det}(G_{f}^{-1}(\omega,k)).$ (77)

In terms of $\mathcal{H}_{eff}$, the self-consistent values of the Hubbard
Stratonovich fields satisfy

$\begin{split}&\partial_{\phi_{1x}}\mathcal{H}_{eff}=0,\\\
&\partial_{\phi_{2x}}\mathcal{H}_{eff}=0,\\\
&\partial_{\phi_{1y}}\mathcal{H}_{eff}=0,\\\
&\partial_{\phi_{2y}}\mathcal{H}_{eff}=0,\\\ \end{split}$ (78)
$\partial_{\theta_{x}}\mathcal{H}_{eff}=0,$ (79)
$\partial_{\theta_{y}}\mathcal{H}_{eff}=0,$ (80)
$\partial_{\nu}\mathcal{H}_{eff}=0,$ (81)
$\partial_{\lambda_{0}}\mathcal{H}_{eff}=0.$ (82)

For all the models we consider in this paper, the constraints in Eq. 78 are
solved by setting

$\begin{split}&\phi_{1x}=\phi_{x}+A_{xy},\\\ &\phi_{2x}=-\phi_{x},\\\
&\phi_{1y}=\phi_{y}+A_{xy},\\\ &\phi_{2y}=-\phi_{y},\\\ \end{split}$ (83)

where $\phi_{x}$ and $\phi_{y}$ are two new phase fields we have introduced
here. Similarly, for all the models we consider in this paper, the constraints
in Eq. 79 and 80 are solved by setting

$\theta_{x}=\theta_{y}=\pi/4.$ (84)

The solution to Eq. 81 depends on the details of the model we are considering.
For the Hamiltonians in Eq. 15 and 40, Eq. 81 is solved by

$\nu=\pi/4.$ (85)

For the Hamiltonians in Eq. 60 and 68, Eq. 81 is solved by

$\nu=0,$ (86)

provided that $\gamma>\gamma^{\prime}$. For all the lattice models, Eq. 82
must be solved numerically as a function of the lattice parameters. In order
to avoid any unnecessary numeric constants in our definitions, it will also be
useful to define the new variables

$\begin{split}&\lambda\equiv\lambda_{0}/2,\\\
&\lambda_{x}\equiv\lambda_{0}/\sqrt{2}.\end{split}$ (87)

In addition to the effective potential, we can also consider the effective
kinetic term for the Hubbard Stratonovich fields. These terms can be found by
evaluating the polarization diagram of the fermions. At one loop order, we
find that when $d=2$, the fermions generate a kinetic term of the form

$\begin{split}\mathcal{H}_{kin-
eff}=&-g_{x}\lambda_{2x}(\bm{r}+\hat{y})\lambda_{1x}(\bm{r})\\\
&-g_{y}\lambda_{2y}(\bm{r}+\hat{x})\lambda_{1y}(\bm{r})+h.c.,\end{split}$ (88)

where we have used the original definitions of the Hubbard-Stratonovich
fields, and added back in the dependence on the lattice coordinates. When
$d=3$, we also find a similar term, which is equivalent to Eq. 88 upon
exchanging the $2$D lattice coordinate $\bm{r}$ with the $3$D lattice
coordinate $\bm{R}$. Combining Eq. 88, with our results from analyzing Eq.
78-82, we find the the kinetic energy is minimized when

$\begin{split}&\Delta_{y}\phi_{x}=A_{xy}\\\ &\Delta_{x}\phi_{y}=A_{xy},\\\
\end{split}$ (89)

where $\Delta_{x}$ and $\Delta_{y}$ are the lattice derivatives.

Using Eq. 83-87 and 89, along with the numeric solutions for the self-
consistent values of $\lambda_{0},$ we are able to determine the self-
consistent values of the Hubbard-Stratonovich fields $\lambda_{ai}$.

##  Appendix B Quadrupole Moment of the Mean Field Hamiltonian

In this appendix we will discuss how to determine the quadrupole moment for
the $2$D mean field models we discussed in the main text. For a Hamiltonian
coupled to a spatially varying gauge field $A_{xy}(\bm{r})$, the change in the
quadrupole moment during an adiabatic process is given by

$\begin{split}\frac{\partial}{\partial\theta}Q_{xy}=\lim_{\epsilon\rightarrow
0}\frac{i}{\epsilon L_{x}L_{y}}\sum_{n\neq
0}&\left[\frac{\bra{0}\frac{\partial
H}{\partial\theta}\ket{n}\bra{n}\sum_{\bm{r}}\frac{\partial
H}{\partial{A_{xy}(\bm{r})}}\ket{0}}{\epsilon+E_{0}-E_{n}}\right.\\\
&-\left.\frac{\bra{0}\sum_{\bm{r}}\frac{\partial
H}{\partial{A_{xy}(\bm{r})}}\ket{n}\bra{n}\frac{\partial
H}{\partial\theta}\ket{0}}{\epsilon+E_{n}-E_{0}}\right].\end{split}$ (90)

In the limit that the gauge field is flat ($A_{xy}(\bm{r})=$const.), Eq. 90
reduces to Eq. 11. The expectation values in Eq. 90 are calculated with
vanishing gauge field $A_{0}(\bm{r})=A_{xy}(\bm{r})=0$.

Here, we shall consider the mean field limit of the $2$D interacting model

$\begin{split}H^{Q}&=\sum_{\bm{r}}\left\\{\bm{c}^{\dagger}(\bm{r})h^{o}\bm{c}(\bm{r})-A_{0}(\bm{r})\bm{c}^{\dagger}(\bm{r})\bm{c}(\bm{r})\right.\\\
&\phantom{=}-Vc_{1}^{\dagger}(\bm{r})c^{\dagger}_{2}(\bm{r}+\hat{x}+\hat{y})c_{3}(\bm{r}+\hat{x})c_{4}(\bm{r}+\hat{y})e^{iA_{xy}(\bm{r})}\\\
&\phantom{=}+\left.h.c.\right\\},\\\
h^{o}&=\mu\Gamma_{0}+t(\Gamma_{2}+\Gamma_{4}).\end{split}$ (91)

For this Hamiltonian

$\begin{split}\frac{\partial
H^{Q}}{\partial{A_{xy}(\bm{r})}}=&-iVc_{1}^{\dagger}(\bm{r})c^{\dagger}_{2}(\bm{r}+\hat{x}+\hat{y})c_{3}(\bm{r}+\hat{x})c_{4}(\bm{r}+\hat{y})\\\
&+h.c.,\end{split}$ (92)

where we have evaluated the derivative at $A_{xy}(\bm{r})=0$. To use mean
field theory, we shall decompose the ring exchange term using the following
Hubbard-Stratonovich transformation

$\begin{split}&-Vc_{1}^{\dagger}(\bm{r})c_{2}^{\dagger}(\bm{r}+\hat{x}+\hat{y})c_{3}(\bm{r}+\hat{x})c_{4}(\bm{r}+\hat{y})e^{iA_{xy}(\bm{r})}\\\
&\rightarrow\lambda_{1x}(\bm{r})c_{2}^{\dagger}(\bm{r}+\hat{x}+\hat{y})c_{4}(\bm{r}+\hat{y})e^{iA_{xy}(\bm{r})}\\\
&\phantom{\rightarrow}+\lambda_{2x}(\bm{r})c_{1}^{\dagger}(\bm{r})c_{3}(\bm{r}+\hat{x})\\\
&\phantom{\rightarrow}+\lambda_{1y}(\bm{r})c_{2}^{\dagger}(\bm{r}+\hat{x}+\hat{y})c_{3}(\bm{r}+\hat{x})e^{iA_{xy}(\bm{r})}\\\
&\phantom{\rightarrow}-\lambda_{2y}(\bm{r})c_{1}^{\dagger}(\bm{r})c_{4}(\bm{r}+\hat{y})\\\
&\phantom{\rightarrow}-\frac{2}{V}\lambda_{1x}(\bm{r})\lambda_{2x}(\bm{r})-\frac{2}{V}\lambda_{1y}(\bm{r})\lambda_{2y}(\bm{r}).\end{split}$
(93)

This is equivalent to the Hubbard-Stratonovich transformation used in the main
text upon shifting $\lambda_{1x}\rightarrow\lambda_{1x}e^{iA_{xy}}$ and
$\lambda_{1y}\rightarrow\lambda_{1y}e^{iA_{xy}}$. After this transformation,

$\begin{split}\frac{\partial
H^{Q}}{\partial{A_{xy}(\bm{r})}}=&[i\lambda_{1x}(\bm{r})c_{2}^{\dagger}(\bm{r}+\hat{x}+\hat{y})c_{4}(\bm{r}+\hat{y})\\\
&+i\lambda_{1y}(\bm{r})c_{2}^{\dagger}(\bm{r}+\hat{x}+\hat{y})c_{3}(\bm{r}+\hat{x})h.c.].\end{split}$
(94)

In the mean field limit, the Hubbard-Stratonovich fields acquire expectation
values of the form

$\begin{split}&\lambda_{1x}(\bm{r})=\lambda e^{i\phi_{x}(\bm{r})},\\\
&\lambda_{2x}(\bm{r})=\lambda e^{-i\phi_{x}(\bm{r})},\\\
&\lambda_{1y}(\bm{r})=\lambda e^{i\phi_{y}(\bm{r})},\\\
&\lambda_{2y}(\bm{r})=\lambda e^{-i\phi_{y}(\bm{r})},\end{split}$ (95)

where $\lambda$ is a constant. Using this, the mean field Hamiltonian is

$\begin{split}H^{Q}_{MF}=&\sum_{\bm{r}}\left\\{\bm{c}^{\dagger}(\bm{r})h^{o}\bm{c}(\bm{r})-A_{0}(\bm{r})\bm{c}^{\dagger}(\bm{r})\bm{c}(\bm{r})\right.\\\
&+[\lambda
e^{i\phi_{x}(\bm{r})+iA_{xy}(\bm{r})}c_{2}^{\dagger}(\bm{r}+\hat{x}+\hat{y})c_{4}(\bm{r}+\hat{y})\\\
&+\lambda
e^{-i\phi_{x}(\bm{r})}c_{1}^{\dagger}(\bm{r})c_{3}(\bm{r}+\hat{x})\\\
&+\lambda
e^{i\phi_{y}(\bm{r})+iA_{xy}(\bm{r})}c_{2}^{\dagger}(\bm{r}+\hat{x}+\hat{y})c_{3}(\bm{r}+\hat{x})\\\
&-\left.\lambda
e^{-i\phi_{y}(\bm{r})}c_{1}^{\dagger}(\bm{r})c_{4}(\bm{r}+\hat{y})+h.c.]\right\\},\end{split}$
(96)

and

$\begin{split}\frac{\partial H^{Q}_{MF}}{\partial{A_{xy}(\bm{r})}}=&[i\lambda
e^{i\phi_{x}(\bm{r})}c_{2}^{\dagger}(\bm{r}+\hat{x}+\hat{y})c_{4}(\bm{r}+\hat{y})\\\
&+i\lambda_{1y}e^{i\phi_{y}(\bm{r})}c_{2}^{\dagger}(\bm{r}+\hat{x}+\hat{y})c_{3}(\bm{r}+\hat{x})h.c.].\end{split}$
(97)

As a consistency check, we note that combining Eq. 94 and 95 also leads to Eq.
97.

To proceed, we will introduce a new non-interacting Hamiltonian

$\begin{split}H^{Q}_{q}=&\sum_{\bm{r}}\left\\{\bm{c}^{\dagger}(\bm{r})h^{o}\bm{c}(\bm{r})\right.\\\
&+[\lambda
e^{i\phi_{x}(\bm{r})+iq}c_{2}^{\dagger}(\bm{r}+\hat{x}+\hat{y})c_{4}(\bm{r}+\hat{y})\\\
&+\lambda
e^{-i\phi_{x}(\bm{r})}c_{1}^{\dagger}(\bm{r})c_{3}(\bm{r}+\hat{x})\\\
&+\lambda
e^{i\phi_{y}(\bm{r})+iq}c_{2}^{\dagger}(\bm{r}+\hat{x}+\hat{y})c_{3}(\bm{r}+\hat{x})\\\
&-\left.\lambda
e^{-i\phi_{y}(\bm{r})}c_{1}^{\dagger}(\bm{r})c_{4}(\bm{r}+\hat{y})+h.c.]\right\\},\end{split}$
(98)

which is equivalent to Eq. 96 at $A_{0}(\bm{r})=0$ and $A_{xy}(\bm{r})=q$
(const). In particular, at $A_{xy}(\bm{r})=q=0,$ Eq. 96 and 98 are equal to
each other, and by extension have the same ground states and excited states.
Additionally, we also have that

$\sum_{\bm{r}}\frac{\partial
H^{Q}_{MF}}{\partial{A_{xy}(\bm{r})}}=\frac{\partial H^{Q}_{q}}{\partial q},$
(99)

where the derivatives are evaluated at $A_{xy}(\bm{r})=q=0$. Combining this
with Eq. 90 we find that the change in the quadrupole moment in the mean field
limit is

$\begin{split}\frac{\partial}{\partial\theta}Q_{xy}=\lim_{\epsilon\rightarrow
0}\frac{i}{\epsilon L_{x}L_{y}}\sum_{n\neq
0}&\left[\frac{\bra{0}\frac{\partial
H^{Q}_{q}}{\partial\theta}\ket{n}\bra{n}\frac{\partial H^{Q}_{q}}{\partial
q}\ket{0}}{\epsilon+E_{0}-E_{n}}\right.\\\ &-\left.\frac{\bra{0}\frac{\partial
H^{Q}_{q}}{\partial q}\ket{n}\bra{n}\frac{\partial
H^{Q}_{q}}{\partial\theta}\ket{0}}{\epsilon+E_{n}-E_{0}}\right],\end{split}$
(100)

where $\ket{n}$ are the energy eigenstates of $H^{Q}_{q}$ at $q=0$ (which are
equal to energy eigenstates of $H^{Q}_{MF}$ at $A_{xy}(\bm{r})=0$). After some
algebra, we can express the total change in the quadrupole moment in terms of
the eigenfunctions of $H^{Q}_{q}$ as

$\begin{split}\Delta Q_{xy}=\Delta\left[(-i)\sum_{\alpha\in
occ}\int\frac{d^{2}\vec{k}}{4\pi^{2}}\bra{\vec{k},\alpha}\partial_{q}\ket{\vec{k},\alpha}\right],\end{split}$
(101)

where$\ket{\vec{k},\alpha}$ is a single particle eigenfunction of $H^{Q}_{q}$
with momentum $\vec{k}$ and band index $\alpha$, and the sum is over the
occupied bands. We have also integrated over $\theta$ to derive Eq. 101. This
formula can be used to calculate the change in the quadrupole moment during a
given adiabatic process, such as the one in Sec. III.2.

We can also use this procedure to find the change in the quadrupole moment for
the mean field limit of Eq. 68. To do this, we will introduce the Hamiltonian

$\begin{split}H^{M}_{q^{\prime}}=&\sum_{\bm{r}}\left[\bm{c}^{\dagger}(\bm{r})h^{oM}\bm{c}(\bm{r})-A_{0}(\bm{r})\bm{c}^{\dagger}(\bm{r})\bm{c}(\bm{r})\right.\\\
&+[\lambda_{x}e^{i\phi_{x}+iq^{\prime}}c_{2}^{\dagger}(\bm{r}+\hat{x}+\hat{y})c_{4}(\bm{r}+\hat{y})\\\
&+\left.\lambda_{x}e^{-i\phi_{x}}c_{1}^{\dagger}(\bm{r})c_{3}(\bm{r}+\hat{x})h.c.]\right],\end{split}$
(102)

where $h^{oM}$ is defined as in Eq. 68. Here, we have included a parameter
$q^{\prime}$, which we will use to find the quadrupole moment. Following the
same steps as before, we find that the change in the quadrupole moment is
given by

$\begin{split}\Delta Q_{xy}=\Delta\left[(-i)\sum_{\alpha\in
occ}\int\frac{d^{2}\vec{k}}{4\pi^{2}}\bra{\vec{k},\alpha}\partial_{q^{\prime}}\ket{\vec{k},\alpha}\right],\end{split}$
(103)

where $\ket{\vec{k},\alpha}$ are the energy eigenstates of the Hamiltonian in
Eq. 102.

##  Appendix C Role of the Phase Fields in the Mean Field Limit

In decomposing the ring exchange interaction into the sum of terms that are
quadratic in the fermion creation/annihilation operators, we found it was
necessary to introduce phase fields $\phi_{x}$ and $\phi_{y}$. Under phase
shifts that depend linearly on position, $c({\bf{r}})\rightarrow
c({\bf{r}})e^{i\beta_{1}x+i\beta_{2}y}$, the phase fields transform as
$\phi_{x}\rightarrow\phi_{x}+\beta_{1}$ and
$\phi_{y}\rightarrow\phi_{y}+\beta_{2}$. Similarly, under a gauge transform
$\Lambda$, the phase fields transform as
$\phi_{i}\rightarrow\phi_{i}+\Delta_{i}\Lambda$ for $i=x,y$. These fields
enter the mean field Hamiltonians $h_{MF}$, by shifting the lattice momentum
$k_{i}\rightarrow k_{i}+\phi_{i}$. In this appendix, we shall show that
integrating over the phase fields projects out states with non-vanishing many-
body electric polarization. Here, we shall only consider the $\phi_{x}$ field.
The analysis for the $\phi_{y}$ field can be done analogously.

To show that the integration projects out states with non-vanishing many-body
polarization, let us consider the $h_{MF}$. If we set $\phi_{x}=0$, then the
Hamiltonian $h_{MF}(k_{x})$ is diagonalized by the single-particle states

$\gamma^{\dagger}_{n}(k_{x})=u^{a\dagger}_{n}(k_{x})c^{\dagger}_{a}(k_{x}),$
(104)

where $n$ is the band index. Here we are suppressing all the dependence on all
momenta except $k_{x}$ for simplicity. In terms of these single-particle
eigenstates, the ground state of the Hamiltonian with $\phi_{x}=0$ is given by

$\begin{split}\ket{0,0}&=\prod_{\\{n,k_{x}\\}\in\text{occ}}u^{a\dagger}_{n}(k_{x})c^{\dagger}_{a}(k_{x})\ket{vac}\\\
&=\prod_{\\{n,k_{x}\\}\in\text{occ}}u^{a\dagger}_{n}(k_{x})\sum_{x}\frac{e^{ik_{x}x}}{\sqrt{L_{x}}}c^{\dagger}_{a}(x)\ket{vac}.\end{split}$
(105)

For $\phi_{x}\neq 0$, the mean field Hamiltonian is given by
$h_{MF}(k_{x}+\phi_{x})$. Hence for $\phi_{x}\neq 0$, the state in Eq. 105
becomes

$\begin{split}\ket{0,\phi_{x}}&=\prod_{\\{n,k_{x}\\}\in\text{occ}}u^{a\dagger}_{n}(k_{x}+\phi_{x})c^{\dagger}_{a}(k_{x})\ket{vac}\\\
&=\prod_{\\{n,k_{x}\\}\in\text{occ}}u^{a\dagger}_{n}(k_{x}+\phi_{x})\sum_{x}\frac{e^{ik_{x}x}}{\sqrt{N_{x}}}c^{\dagger}_{a}(x)\ket{vac},\end{split}$
(106)

where $N_{x}$ is the length of the lattice system in the $x$-direction.

To proceed we will make two assumptions. First, we will assume that the single
particle bands defined via Eq. 104 are completely filled in the ground state
(as they are in non-interacting band insulators). Second, we will assume that
either $\phi_{x}$ can be written in the form $2\pi n/N_{x}$ for
$n\in\mathbb{Z}$ (which means that the dipole gauge transformations must obey
periodic boundary conditions), or that we are in the thermodynamic limit
($N_{x}\rightarrow\infty$). When these conditions are true, we can shift
$k_{x}\rightarrow k_{x}-\phi_{x}$, and rewrite the state $\ket{0,\phi_{x}}$ as

$\begin{split}\ket{0,\phi_{x}}&=\prod_{\\{n,k_{x}\\}\in\text{occ}}u^{a\dagger}_{n}(k_{x})c^{\dagger}_{a}(k_{x}-\phi_{x})\ket{vac}\\\
&=\prod_{\\{n,k_{x}\\}\in\text{occ}}u^{a\dagger}_{n}(k_{x})\sum_{x}\frac{e^{ik_{x}x-i\phi_{x}x}}{\sqrt{L_{x}}}c^{\dagger}_{a}(x)\ket{vac}\\\
&=e^{-i\phi_{x}xn(x)}\prod_{\\{n,k_{x}\\}\in\text{occ}}u^{a\dagger}_{n}(k_{x})\sum_{x}\frac{e^{ik_{x}x}}{\sqrt{L_{x}}}c^{\dagger}_{a}(x)\ket{vac}\\\
&=e^{-i\phi_{x}xn(x)}\ket{0,0},\end{split}$ (107)

where $n(x)=\sum_{a}c^{\dagger}_{a}(x)c_{a}(x)$. The operator
$e^{-i\phi_{x}xn(x)}$ from Eq. 107 is related to the expectation value of the
polarization operator for periodic systems via

$P_{x}=-\frac{1}{\phi_{x}N_{x}}\text{Im}\log\langle
e^{-i\phi_{x}xn(x)}\rangle.$ (108)

where $P_{x}$ is electronic contribution to the polarization of the system.

We will now rewrite the $\phi_{x}=0$ ground state from Eq. 105 as a sum of
eigenstates of the operator $e^{-i\phi_{x}xn(x)}$

$\begin{split}\ket{0,0}&=\sum_{n,X}a_{X,n}\ket{X,n}\\\ \end{split}$ (109)

where $\ket{X,n}$ satisfies
$e^{-i\phi_{x}xn(x)}\ket{X,n}=e^{-i\phi_{x}X}\ket{X,n}$. From this we can
conclude that the state $\ket{X,n}$ has total polarization $X/N_{x}+P_{ion}$,
where $P_{ion}$ is the contribution to the polarization from the ions. The
ionic contribution must be included in order for the system is charge neutral,
and for the polarization to be invariant with respect to a change of
coordinates. For the models we are considering, a simple calculation shows
that $P_{ion}=0$ mod$(1)$.

Combining Eq. 107 and Eq. 109, the $\phi_{x}\neq 0$ state is given by

$\begin{split}\ket{0,\phi_{x}}&=e^{-i\phi_{x}xn(x)}\sum_{n,X}a_{X,n}\ket{X,n}\\\
&=\sum_{n,X}a_{X,n}e^{i\phi_{x}X}\ket{X,n}.\end{split}$ (110)

From this we can see that integrating over $\phi_{x}$ will project out any
states with $X/N_{x}\neq 0\text{ mod(1)}$. These are exactly the unpolarized
many-body states.

##  Appendix D Continuum Analysis of $M\mathcal{T}$ model

In this appendix, we will analyze the $3$D $M\mathcal{T}$ symmetric model
given in Eq. 60 in the continuum, near $V=0$. At $V=0$, the system is
quasi-1d, i.e., a decoupled $3$D array of $1$D wires oriented in the
$z$-direction. If we pass to the continuum along the $z$-direction, each of
these wires consists of 4 massless fermions (which correspond to the
fluctuations of the lattice fermions with momentum near $k_{z}=\pi$), and 4
massive fermions (which correspond to the fluctuations of the lattice fermions
with momentum near $k_{z}=0$). At $V=0$, the continuum Lagrangian for these
fields is

$\begin{split}&\mathcal{L}=\sum_{\bm{r}}\bm{\psi}^{\dagger}(\bm{r})G^{-1}_{0}(\omega+A_{0}(\bm{r}),p_{z}+A_{z}(\bm{r}))\bm{\psi}(\bm{r})\\\
&\phantom{====}+\bm{\Psi}^{\dagger}(\bm{r})G^{-1}_{\Lambda}(\omega+A_{0}(\bm{r}),p_{z}+A_{z}(\bm{r}))\bm{\Psi}(\bm{r})\\\
&G^{-1}_{0}(\omega,k_{z})=\omega\text{I}+p_{z}\Gamma^{0},\\\
&G^{-1}_{\Lambda}(\omega,k_{z})=\omega\text{I}-p_{z}\Gamma^{0}-M(\cos(\kappa)\Gamma^{2}+\sin(\kappa)\Gamma^{4}),\end{split}$
(111)

where $\bm{\psi}=(\psi_{1},\psi_{2},\psi_{3},\psi_{4})$ are the light fermion
operators, $\bm{\Psi}=(\Psi_{1},\Psi_{2},\Psi_{3},\Psi_{4})$ are the heavy
fermions operators, and $\kappa$ is defined such that
$\tan(\kappa)=\gamma^{\prime}/\gamma$. We will take $M$ to be the UV cutoff
for this theory.

We can now consider the interactions in the continuum. As before, we will use
a Hubbard-Stratonovich transformation to decompose the lattice ring exchange
interaction. Using Eq. 42 we see that in the continuum the Hubbard-
Stratonovich fields $\lambda_{1/2x}$ and $\lambda_{1/2y}$ couple to both the
light fermions $\psi$ and the heavy fermions $\Psi$. Since they are gapped,
the heavy fermions can be integrated out, leaving a Lagrangian in terms of the
light fermions and the Hubbard-Stratonovich fields. After integrating out the
heavy fermions, the Lagrangian can be written as

$\mathcal{L}^{dCS}=\mathcal{L}_{\psi}+\mathcal{L}_{\lambda\lambda},$ (112)

where $\mathcal{L}_{\psi}$ contains all terms involving the light fermions
$\psi$ and $\mathcal{L}_{\lambda\lambda}$ contains all couplings between the
Hubbard-Stratonovich fields. If we ignore any terms that osculate like
$(-1)^{z}$, $\mathcal{L}_{\psi}$ is given by

$\begin{split}\mathcal{L}_{\psi}=&\sum_{\bm{r}}\bm{\psi}^{\dagger}(\bm{r})G^{-1}_{0}(\omega+A_{0}(\bm{r}),p_{z}+A_{z}(\bm{r}))\bm{\psi}(\bm{r})\\\
&-[\lambda_{2x}(\bm{r})\psi_{1}^{\dagger}(\bm{r})\psi_{3}(\bm{r}+\hat{x})\\\
&\phantom{=}+\lambda_{1x}(\bm{r})\psi_{2}^{\dagger}(\bm{r}+\hat{x}+\hat{y})\psi_{4}(\bm{r}+\hat{y})\\\
&\phantom{=}+\lambda_{2y}(\bm{r})\psi_{1}^{\dagger}(\bm{r})\psi_{4}(\bm{r}+\hat{y})\\\
&\phantom{=}+\lambda_{1y}(\bm{r})\psi_{2}^{\dagger}(\bm{r}+\hat{x}+\hat{y})\psi_{3}(\bm{r}+\hat{x})+h.c.].\end{split}$
(113)

To determine $\mathcal{L}_{\lambda\lambda}$ we must integrate out the heavy
fermions $\Psi$. At one loop order, $\mathcal{L}_{\lambda\lambda}$ is given by

$\begin{split}\mathcal{L}_{\lambda\lambda}=&\sum_{\bm{r}}-u_{1}[\lambda_{1x}(\bm{r})\lambda^{*}_{1x}(\bm{r})+\lambda_{2x}(\bm{r})\lambda^{*}_{2x}(\bm{r})\\\
&+\lambda_{1y}(\bm{r})\lambda^{*}_{1y}(\bm{r})+\lambda_{2y}(\bm{r})\lambda^{*}_{2y}(\bm{r})]\\\
&+[u^{x}_{2}\lambda_{1x}(\bm{r})\lambda_{2x}(\bm{r}+\hat{y})+u^{y}_{2}\lambda_{1y}(\bm{r})\lambda_{2y}(\bm{r}+\hat{x})\\\
&\phantom{=}+\frac{2}{V}\lambda_{1x}(\bm{r})\lambda_{2x}(\bm{r})e^{-iA_{xy}(\bm{r})}\\\
&\phantom{=}+\frac{2}{V}\lambda_{1y}(\bm{r})\lambda_{2y}(\bm{r})e^{-iA_{xy}(\bm{r})}+h.c.],\end{split}$
(114)

where $u_{1}=\frac{\log(4)-1}{16\pi},$
$u^{x}_{2}=\frac{\cos^{2}(\kappa)}{16\pi}$, and
$u^{y}_{2}=\frac{\sin^{2}(\kappa)}{16\pi}$.

As discussed in the main text, we shall employ the self-consistent mean field
theory approximation. Similar to before, the Hubbard-Stratonovich fields
acquire expectation values of the form

$\begin{split}\lambda_{1x}(\bm{r})&=\lambda_{x}e^{i\phi_{x}(\bm{r})+iA_{xy}(\bm{r})},\\\
\lambda_{2x}(\bm{r})&=\lambda_{x}e^{-i\phi_{x}(\bm{r})},\\\
\lambda_{1y}(\bm{r})&=0,\\\ \lambda_{2y}(\bm{r})&=0.\\\ \end{split}$ (115)

Here, $\lambda_{1y}$ and $\lambda_{2y}$ have vanishing expectation values,
while $\lambda_{1y}$ and $\lambda_{2y}$ do not, since $u^{x}_{2}>u^{y}_{2}$ in
Eq. 114. The value of $\lambda_{x}$ is determined by the effective potential

$\begin{split}H_{\lambda\lambda}=&\frac{4}{V_{x}^{\prime}}\lambda_{x}^{2}+\frac{1}{2\pi}\lambda_{x}^{2}\Big{(}\log(\frac{\lambda_{x}^{2}}{M^{2}})-1\Big{)},\end{split}$
(116)

where $V_{x}^{\prime}=(V^{-1}+\frac{u_{1}+u^{x}_{2}}{2})^{-1}$. The effective
potential is minimized by $\lambda_{x}=Me^{-\frac{4\pi}{V_{x}^{\prime}}}$. In
agreement with the numeric results, we find that $\lambda_{x}$ vanishes when
$V\rightarrow 0,$ and increases monotonically with increasing $V$.
Additionally, due to the $u^{x}_{2}$ term in Eq. 114, at low energies

$\begin{split}&\phi_{x}(\bm{r}+\hat{y})-\phi_{x}(\bm{r})=\Delta_{y}\phi_{x}(\bm{r})=A_{xy}(\bm{r}).\end{split}$
(117)

As noted before, under a gauge transformation $\Lambda$, the phase fields
$\phi_{x}$ transform as
$\phi_{x}(\bm{r})\rightarrow\phi_{x}(\bm{r})+\Delta_{x}\Lambda(\bm{r})$.
Because of this, all terms in Eq. 117 have the same gauge transformation, as
desired.

To find the effective response theory for the continuum model in the mean
field limit, the effective response action can be found by evaluating the
current-current correlation functions. Here we are interested in the
correlation functions between the currents associated with $A_{0}(\bm{r})$,
$A_{z}(\bm{r})$, $A_{xy}(\bm{r}),$ and $\phi_{x}(\bm{r})$ . These currents are
found by taking the functional derivative of the mean field fermionic action
with respect to one of the aforementioned fields. In the bulk, we find that in
the low frequency and momentum limit, the effective response theory response
is

$\begin{split}\mathcal{L}_{eff,\text{bulk}}=\frac{1}{4\pi}&(A_{xy}-\partial_{y}\phi_{x})\partial_{t}A_{z}\\\
&-(A_{xy}-\partial_{y}\phi_{x})\partial_{z}A_{0},\end{split}$ (118)

where we have passed to the continuum in the $x$ and $y$ directions as well.
Using Eq. 117, we find that the bulk action vanishes. However, at the
boundaries normal to the $x$ and $y$-directions, as well as the hinges
separating these surfaces, we find that the response Lagrangians are

$\begin{split}L_{eff,\pm x}=0\\\ L_{eff,\pm
y}=\pm&\frac{1}{2\pi}\big{[}A_{0}\partial_{x}A_{z}+\phi_{x}\partial_{z}A_{0}-\phi_{x}\partial_{t}A_{y}\big{]},\\\
L_{eff,\pm x,\pm y}=-&\frac{1}{4\pi}A_{0}A_{z}.\\\ L_{eff,\pm x,\mp
y}=+&\frac{1}{4\pi}A_{0}A_{z}.\end{split}$ (119)

This is consistent with what we found using the mean field analysis in Sec.
VI.1.


# Hubble constant from the cluster-lensed quasar system SDSS J1004+4112:
investigation of the lens model dependence

Yuting<EMAIL_ADDRESS>Institute for Frontiers in Astronomy and
Astrophysics, Beijing Normal University, Beijing 102206, China  Department of
Astronomy, Beijing Normal University, Beijing 100875, China  Department of
Physics, University of Tokyo, Tokyo 113-0033, Japan  Masamune
<EMAIL_ADDRESS>Center for Frontier Science, Chiba
University, Chiba 263-8522, Japan  Department of Physics, Graduate School of
Science, Chiba University, Chiba 263-8522, Japan  Kavli Institute for the
Physics and Mathematics of the Universe, University of Tokyo, Kashiwa, Chiba
277-8583, Japan  Shuo<EMAIL_ADDRESS>Institute for Frontiers in
Astronomy and Astrophysics, Beijing Normal University, Beijing 102206, China
Department of Astronomy, Beijing Normal University, Beijing 100875, China

###### Abstract

As a fundamental parameter for modern cosmology, the Hubble constant $H_{0}$
is experiencing a serious crisis. In this paper, we explore an independent
approach to measure $H_{0}$ based on the time-delay cosmography with strong
gravitational lensing of a quasar by a galaxy cluster. Specifically we focus
on the strong lensing system SDSS J1004+4112 with the maximum image separation
of 14.62′′, the first system of a quasar lensed by a galaxy cluster with five
multiple images. Incorporating the latest time-delay measurements, we
investigate the lens model dependence from the combination of 16 different
lens mass models. We find that the lens model dependence is indeed large, with
the combined measurement of the Hubble constant of
$H_{0}=67.5^{+14.5}_{-8.9}\mathrm{~{}km~{}s^{-1}~{}Mpc^{-1}}$ that is obtained
by summing posteriors of the Hubble constant from the 16 models with equal
weighting. Interestingly, our results show that the value of Hubble constant
decreases as the complexity of the perturbation around the lens increases,
although weighting based on positional errors of quasar images does not
significantly improve the $H_{0}$ constraint. We find that the 16 different
mass models predict largely different shapes of the lensed quasar host galaxy
as well as other lensed galaxies behind the cluster. By selecting two mass
models that best reproduces those shapes, the constraint on the Hubble
constant is significantly tightened to
$H_{0}=59.1^{+3.6}_{-3.5}\mathrm{~{}km~{}s^{-1}~{}Mpc^{-1}}$. While we caution
that our analysis still does not fully explore all the possible mass model
uncertainty, our results highlight the importance of including as many
constraints as possible such as extended shapes of lensed galaxies for
obtaining tight constraints on the Hubble constant from cluster-lensed quasar
lens systems.

cosmological parameters, gravitational lensing, quasars, galaxy clusters

##  I Introduction

The Hubble constant $H_{0}$, which quantifies the expansion rate of the
Universe today, is one of the most important cosmological parameters that is
also related to the age and the total energy density of the Universe Weinberg
et al. (2013); Freedman et al. (2001). Although $H_{0}$ has been measured with
high precision, the determination of its accurate value is extremely
challenging, followed with an increasing tension. In particular, the
discrepancy of measured $H_{0}$ values between different measurements is
becoming more pronounced. The cosmic microwave background radiation of the
early Universe from the Planck satellite provides a tight constraint on the
Hubble constant as $H_{0}=67.4\pm 0.5\mathrm{~{}km~{}s^{-1}~{}Mpc^{-1}}$,
which is based on the so-called standard cosmology, $\Lambda$-dominated cold
dark matter ($\Lambda$CDM) model Planck Collaboration et al. (2020). Such
concordance $\Lambda$CDM model has withstood most popular observational
evidences Riess et al. (1998); Perlmutter et al. (1999); Cao et al. (2015);
Troxel et al. (2018); Scolnic et al. (2018); Bianchini et al. (2020); Han et
al. (2021); Asgari et al. (2021); Mossa et al. (2020). An independent
measurement of $H_{0}$ has been obtained by type Ia supernovae (SNe Ia)
calibrated via the distance ladder in the local Universe, with the value of
$H_{0}=73.04\pm 1.04\mathrm{~{}km~{}s^{-1}~{}Mpc^{-1}}$ released by SH0ES
(SNe, H0, for the Equation of State of dark energy) collaboration Riess et al.
(2022). It is obvious that the $H_{0}$ constraints derived by these two
methods are inconsistent. Such significant 5$\sigma$ tension has sparked
debates about the validity of the $\Lambda$CDM model.

As a consequence, any other methodology to measure the Hubble constant is very
necessary to clarify the origin of the discrepancy. The tip of the red giant
branch (TRGB) as a local standard candle provides another way to determine the
Hubble constant. Ref. Freedman et al. (2019) presented $H_{0}=69.8\pm
0.8\mathrm{~{}km~{}s^{-1}~{}Mpc^{-1}}$ by measuring TRGB in nine SNe Ia hosts
and calibrating TRGB in the Large Magellanic Cloud. Several recent
calibrations of the TRGB method were further combined, which are internally
self-consistent at the 1% level Freedman (2021). Making geometric distance
measurements to megamaser galaxies in the Hubble flow is another particular
approach to constrain the Hubble constant Braatz and Gugliucci (2008); Reid et
al. (2013); Kuo et al. (2013, 2015); Gao et al. (2017). Recently, Ref. Pesce
et al. (2020) employed such approach to derive the value of $H_{0}=73.9\pm
3.0\mathrm{~{}km~{}s^{-1}~{}Mpc^{-1}}$. In addition, gravitational waves from
compact binary mergers provide a useful means of directly measuring luminosity
distances, and therefore provide an independent approach to constrain $H_{0}$,
which is sometimes referred to as the standard siren technique. For instance,
from the neutron star merger event GW170817 with its electromagnetic
counterparts AT2017gfo and GRB170817A, Ref. Dietrich et al. (2020) obtained
the Hubble constant value of
$H_{0}=66.2^{+4.4}_{-4.2}\mathrm{~{}km~{}s^{-1}~{}Mpc^{-1}}$. The LIGO-Virgo-
KAGRA Scientific Collaboration obtained
$H_{0}=68^{+12}_{-8}\mathrm{~{}km~{}s^{-1}~{}Mpc^{-1}}$ with 47 gravitational-
wave sources from the Third LIGO-Virgo-KAGRA Gravitational-Wave Transient
Catalog (GWTC3) The LIGO Scientific Collaboration et al. (2021).

It is extremely important to further derive measurements of $H_{0}$ using
independent techniques such as strong gravitational lensing. Ref. Refsdal
(1964) proposed that observations of time delays between multiple images of
gravitationally lensed supernovae could be used to determinate $H_{0}$
independently of the distance ladder. Very recently, this method finally
applied to the gravitationally lensed supernova SN Refsdal to obtain the
constraint on the Hubble constant of
$64.8^{+4.4}_{-4.3}\mathrm{~{}km~{}s^{-1}~{}Mpc^{-1}}$ Kelly et al. (2023).

Given the rarity of lensed supernovae, lensed quasars have so far been mostly
used to constrain the Hubble constant. For instance, the H0LiCOW ($H_{0}$
Lenses in COSMOGRAIL’s Wellspring) collaboration presented the estimation of
the Hubble constant from the analysis of time delays for six lensed quasars as
$H_{0}=73.3^{+1.7}_{-1.8}\mathrm{~{}km~{}s^{-1}~{}Mpc^{-1}}$ Wong et al.
(2020). While the H0LiCOW result was obtained by assuming rather simplistic
radial mass distributions for lensing galaxies, the TDCOSMO (Time-Delay
COSMOgraphy) team addressed the effect of lens model degeneracies by
considering more flexible radial mass distributions and breaking the
degeneracies by stellar kinematics, yielding a somewhat weaker constraint of
$H_{0}=74.5^{+5.6}_{-6.1}\mathrm{~{}km~{}s^{-1}~{}Mpc^{-1}}$ Birrer et al.
(2020). In the same paper, a slightly different value of
$H_{0}=67.4^{+4.1}_{-3.2}\mathrm{~{}km~{}s^{-1}~{}Mpc^{-1}}$ was also obtained
by adding a prior from Sloan ACS survey galaxy-galaxy strong lens systems.
Recently, Ref. Shajib et al. (2023) obtained
$H_{0}=77.1^{+7.3}_{-7.1}\mathrm{~{}km~{}s^{-1}~{}Mpc^{-1}}$ from the lens
RXJ1131-1231 using spatially resolved stellar kinematics of the lens galaxy.
Clearly a challenge lies in how to break degeneracies between different mass
models for obtaining tight constraints on the Hubble constant.

We note that galaxy-scale lensed quasars are used in the H0LiCOW and TDCOSMO
analysis. Since lens modeling of galaxy- and cluster-scale lenses involve
quite different systematic errors, it is worth exploring the possibility of
obtaining the complementary constraint on the Hubble constant from cluster-
lensed quasar systems, which is a relatively new area of research partly
because of the small sample of such strong lens systems. One notable
difference between galaxy- and cluster-scale lensed quasars is that
constraints on the lens mass model from additional multiple images of galaxies
behind the lensing cluster are commonly available for the latter Sharon et al.
(2005, 2017); Acebron et al. (2022a, b). While it is known that cluster strong
lens mass modeling exhibits relatively large positional offsets between
observed and model-predicted multiple images, a careful analysis with mock
strong lens data from ray-tracing simulations Meneghetti et al. (2017) has
demonstrated that it is possible to recover input Hubble constant accurately
even with such large positional offsets Kelly et al. (2023). Ref. Denzel et
al. (2021) used the cluster-scale lensed quasar SDSS J1004+4112 Inada et al.
(2003); Oguri et al. (2004); Inada et al. (2005); Fohlmeister et al. (2008)
together with 7 galaxy-scale lensed quasars to obtain
$H_{0}=71.8^{+3.9}_{-3.3}\mathrm{~{}km~{}s^{-1}~{}Mpc^{-1}}$. Ref. Napier et
al. (2023a) focused more specifically on cluster-scale lensed quasars by
studying three cluster-lensed quasars SDSS J1004+4112, SDSS J1029+2623 Inada
et al. (2006); Oguri et al. (2008); Fohlmeister et al. (2013), and SDSS
J2222+2745 Dahle et al. (2013, 2015), to produce a combined measurement of
Hubble constant as $H_{0}=71.5\pm 6.1\mathrm{~{}km~{}s^{-1}~{}Mpc^{-1}}$.
However, in this analysis only one specific mass model is considered for each
lens system, and as a result the degeneracy among different mass models has
not been fully explored. Very recently, Ref. Martínez-Arrizabalaga et al.
(2023) presented a new constraint on the Hubble constant from the SDSS J1004 +
4112, $H_{0}=74^{+9}_{-13}\mathrm{~{}km~{}s^{-1}~{}Mpc^{-1}}$, using a free-
form lens model.

In this paper, we focus on the first cluster-lensed quasar system SDSS
J1004+4112 Inada et al. (2003); Oguri et al. (2004); Inada et al. (2005) to
derive an independent $H_{0}$ measurement, paying a careful attention to the
dependence of the $H_{0}$ constraint on the assumed lens mass model.
Specially, the latest time-delay measurements between different images Muñoz
et al. (2022) and the combination of 16 different lens mass models are
considered in order to address the lens mass dependence on derived values of
$H_{0}$. Our work is complementary to Ref. Forés-Toribio et al. (2022) in
which the latest time-delay measurements are used to place tight constraints
on the lens mass distribution. We also discuss possibilities to break the lens
model degeneracy by adding other observables to obtain better constraints on
$H_{0}$. Our work will be useful for future cosmological application of
cluster-lensed quasar systems that are actively being searched and discovered
recently Shu et al. (2018, 2019); Martinez et al. (2023); Napier et al.
(2023b).

This paper is arranged as follows. In Sec. II, we briefly describe the
cluster-lensed quasar system SDSS J1004+4112. We introduce the methodology of
lens modeling in Sec. III. Results of our analysis are presented in Sec. IV.
The possible improvement of the Hubble constant with additional constraints is
discussed in Sec. V. We summarize the main conclusion in Sec. VI.

Figure 1: The HST/ACS F814W image of the cluster-lensed quasar system SDSS J1004+4112. Red dots labeled with A-E are five multiple images of a lensed quasar. Red squares are images of three lensed background galaxies at redshifts 3.33, 2.74, and 3.28, respectively. Table 1: Observational constraints from multiple images of lensed quasar and galaxies. Source redshift, $x$\- and $y$-positions with respect to the quasar image A, the magnitude difference, and the relative time delay are listed in columns 2–5, respectively.  Name | $z_{\rm s}$ |  $\Delta x$ [′′] |  $\Delta y$ [′′] | $\Delta m$ |  $\Delta t$ [days]
---|---|---|---|---|---
A | 1.734 | 0.000 | 0.000 | $\equiv 0$ | $825.99\pm 2.10$
B |  | $-1.317$ | 3.532 | $0.35\pm 0.3$ | $781.92\pm 2.20$
C |  | 11.039 | $-4.492$ | $0.87\pm 0.3$ | $\equiv 0$
D |  | 8.399 | 9.707 | $1.50\pm 0.3$ | $2456.99\pm 5.55$
E |  | 7.197 | 4.603 | $6.30\pm 0.8$ | $\cdots$
A1.1 | 3.33 | 3.93 | $-2.78$ |  |
A1.2 |  | 1.33 | 19.37 |  |
A1.3 |  | 19.23 | 14.67 |  |
A1.4 |  | 18.83 | 15.87 |  |
A1.5 |  | 6.83 | 3.22 |  |
A2.1 | 3.33 | 4.13 | $-2.68$ |  |
A2.2 |  | 1.93 | 19.87 |  |
A2.3 |  | 19.43 | 14.02 |  |
A2.4 |  | 18.33 | 15.72 |  |
A2.5 |  | 6.83 | 3.12 |  |
A3.1 | 3.33 | 4.33 | $-1.98$ |  |
A3.2 |  | 2.73 | 20.37 |  |
A3.3 |  | 19.95 | 13.04 |  |
A3.4 |  | 18.03 | 15.87 |  |
A3.5 |  | 6.83 | 3.02 |  |
B1.1 | 2.74 | 8.88 | $-2.16$ |  |
B1.2 |  | $-5.45$ | 15.84 |  |
B1.3 |  | 8.33 | 2.57 |  |
B2.1 | 2.74 | 8.45 | $-2.26$ |  |
B2.2 |  | $-5.07$ | 16.04 |  |
B2.3 |  | 8.33 | 2.57 |  |
C1.1 | 3.28 | 10.25 | $-3.06$ |  |
C1.2 |  | $-7.55$ | 15.39 |  |
C1.3 |  | 8.49 | 2.72 |  |
C2.1 | 3.28 | 9.95 | $-3.36$ |  |
C2.2 |  | $-7.30$ | 15.44 |  |
C2.3 |  | 8.49 | 2.72 |  |

##  II Observations of SDSS J1004+4112

Table 2: Assumed positional errors of quasar images $\sigma_{\rm q}$ and
galaxy images $\sigma_{\rm g}$ for the 16 different lens mass models.

Label | Model * |  $\sigma_{\rm q}$[′′] |  $\sigma_{\rm g}$[′′]
---|---|---|---
m111 |  anfw+gals+pert | 0.12 | 1.2
m112 |  anfw+gals+pert+mpole(m=3) | 0.07 | 0.7
m113 |  anfw+gals+pert+mpole(m=3,4) | 0.05 | 0.5
m114 |  anfw+gals+pert+mpole(m=3,4,5) | 0.04 | 0.4
m121 |  anfw+ahern+gals+pert | 0.07 | 0.7
m122 |  anfw+ahern+gals+pert+mpole(m=3) | 0.05 | 0.5
m123 |  anfw+ahern+gals+pert+mpole(m=3,4) | 0.04 | 0.4
m124 |  anfw+ahern+gals+pert+mpole(m=3,4,5) | 0.03 | 0.3
m211 |  jaffe+gals+pert | 0.12 | 1.2
m212 |  jaffe+gals+pert+mpole(m=3) | 0.09 | 0.9
m213 |  jaffe+gals+pert+mpole(m=3,4) | 0.07 | 0.7
m214 |  jaffe+gals+pert+mpole(m=3,4,5) | 0.04 | 0.4
m221 |  jaffe+ahern+gals+pert | 0.10 | 1.0
m222 |  jaffe+ahern+gals+pert+mpole(m=3) | 0.08 | 0.8
m223 |  jaffe+ahern+gals+pert+mpole(m=3,4) | 0.06 | 0.6
m224 |  jaffe+ahern+gals+pert+mpole(m=3,4,5) | 0.05 | 0.5

  * *

For Model components, anfw refers to the NFW profile for the halo, jaffe
refers to the pseudo-Jaffe profile for the halo, gals refers to being the BCG
treated same as the other member galaxies, ahern refers to the Hernquist
profile for the BCG, pert refers to an external shear, mpole refers to
multipole perturbations of order m. See the text for a detailed explanation on
the label.

As the first system of a quasar lensed by a cluster, the large-separation
gravitationally lensed quasar SDSS J1004+4112 Inada et al. (2003); Oguri et
al. (2004); Inada et al. (2005) was discovered from the Sloan Digital Sky
Survey (SDSS) York et al. (2000); Abazajian et al. (2004) data, which has
generated much interest due to its rarity and particularity. This system is
composed of five multiple images of a quasar at $z_{\rm s}=1.734$ Inada et al.
(2008) with the maximum separation angle between the multiple images of
14.62′′. A cluster of galaxies as the lensing object has been identified at
$z_{\rm l}=0.68$ centered among the five images, the brightest galaxy and more
galaxy cluster members have also been observed. In Fig. 1, we present the five
multiple images of the lensed quasar as well as seven lensed background
galaxies at three different redshifts, which have also been spectroscopically
confirmed Sharon et al. (2005); Liesenborgs et al. (2009); Oguri (2010). Such
a wealth of observational constraints suggest that it may be possible to
tightly constrain the Hubble constant from this lens system.

In our investigation, we employ positions of the five quasar images in Ref.
Oguri (2010) which was adapted according to the Hubble Space Telescope
Advanced Camera for Surveys (HST/ACS) F814W images Inada et al. (2005). For
images of lensed background galaxies, we include the constraints identified in
Refs. Sharon et al. (2005); Liesenborgs et al. (2009); Oguri (2010). There are
7 multiply-imaged knots originating from three galaxies (A, B, and C) at
different redshifts. The galaxy A contains three knots, each of which is
lensed into five images. Galaxies B and C each have two knots in a galaxy with
three lensed images. In Table 1, we summarize positions of the multiple
images, which are distributed in a very wide range of radius. This remarkable
feature helps us model this cluster strong lensing system in great detail and
to constrain the Hubble constant. Moreover, the newly measured time delay
between lensed quasar images D and C ($t_{\rm DC}=2456.99\pm 5.55$ days) is
included in our analysis in addition to the previously measured time delays
between images A-C and B-C Muñoz et al. (2022); Fohlmeister et al. (2007,
2008) as well as magnitude differences between the quasar multiple images. For
the brightest cluster galaxy (BCG), the position is fixed to the observed
position at ($\Delta x$, $\Delta y$)=(7.104′′, 4.362′′) with respect to the
quasar image A. We also include the positions, ellipticities, position angles
and luminosity ratios (with respect to the BCG) of 90 cluster galaxy members,
which are selected from the HST/ACS F435W and F814W images with the red-
sequence feature and whose properties are measured also in the HST/ACS F814W
image using the SExtractor software Bertin and Arnouts (1996).

##  III Lens Modeling

We compute the lens model with the publicly available software glafic Oguri
(2010, 2021), which is designed to model strong lens systems and can compute
efficiently lensed image positions of point sources with adaptive-meshing.

###  III.1 Dark matter halo

The Navarro-Frenk-White (NFW) profile is employed to model the mass
distribution of a dark matter halo Navarro et al. (1997). The radial density
profile of the NFW profile is described by

$\rho=\frac{\rho_{\rm s}}{(r/r_{\rm s})(1+r/r_{\rm s})^{2}},$ (1)

where $r_{\rm s}$ indicates the scale radius and $\rho_{\rm s}$ is the
characteristic density. The projected NFW profile that is used for lensing
analysis is parameterized by mass $M$, the position $x$ and $y$, the
ellipticity $e$, the position angle $\theta_{\rm e}$, and the concentration
parameter $c\equiv r_{\rm vir}/r_{\rm s}$. In our analysis, we optimize all
these parameters with the range for the mass of
$10^{14}-10^{15}h^{-1}M_{\odot}$ and concentration parameter from $1.0$ to
$10.0$. We adopt the fast approximation of the lensing calculation of the NFW
profile proposed in Ref. Oguri (2021), which is named anfw in glafic.

In addition to the NFW profile, the Pseudo-Jaffe Ellipsoid (jaffe in glafic)
profile Jaffe (1983); Keeton (2001) is also used to model the dark matter
halo. In this model, the radial density profile is described as

$\rho=\frac{\rho_{\mathrm{s}}}{\\{1+(r/r_{\rm core})^{2}\\}\\{1+(r/r_{\rm
trun})^{2}\\}}.$ (2)

The projected density profile is parameterized by the velosity dispersion
$\sigma$, the position, the ellipicity, the position angle, the core radius,
and the truncation radius. In our analysis, we optimize all the parameters,
but we restrict the range of $\sigma$ from $500\mathrm{~{}km~{}s^{-1}}$ to
$800\mathrm{~{}km~{}s^{-1}}$, $r_{\rm trun}$ from $5.0^{{}^{\prime\prime}}$ to
$50.0^{{}^{\prime\prime}}$ and $r_{\rm core}$ is smaller than
$15.0^{{}^{\prime\prime}}$.

###  III.2 The brightest cluster galaxy and other cluster galaxy members

In some of our mass models, the brightest cluster galaxy (BCG) is
parameterized by the Hernquist profile Hernquist (1990), which takes the
following radial form

$\rho=\frac{M_{\rm tot}}{2\pi(r/r_{\rm b})(1+r/r_{\rm b})^{3}},$ (3)

where the scale radius $r_{\rm b}$ is related to the effective radius $R_{\rm
e}$ of a project surface mass density as $r_{\rm b}=0.551R_{\rm e}$. The
projected density profile can be parametrized by the total mass $M_{\rm tot}$,
the position, the ellipticity, position angle and the scale radius. In our
analysis, we fix the position $x$ and $y$ to the observed BCG centroid as
described in Sec. II and we optimize the total mass $M_{\rm tot}$ and the
scale radius $r_{\rm b}$ with the range of $10^{12}-10^{13}h^{-1}M_{\odot}$
and $1.0^{{}^{\prime\prime}}-10.0^{{}^{\prime\prime}}$, respectively. We adopt
the fast approximated calculation of lensing properties of the Hernquist
profile (ahern in glafic) proposed in Ref. Oguri (2021). In addition, we add
Gaussian priors to the ellipticity and the position angle of the BCG,
$e=0.355\pm 0.05$ and $\theta_{\rm e}=-20.47^{\circ}\pm 5^{\circ}$,
respectively, based on the observed shape of the light profile of the BCG.

For the other cluster galaxy members, the scaled Pseudo-Jaffe ellipsoids (gals
in glafic) are included to model. The velocity dispersion $\sigma$ and
truncation radius $r_{\rm trun}$ for each galaxy are scaled based on its
luminosity relative to that of the BCG as

$\displaystyle\frac{\sigma}{\sigma_{*}}$
$\displaystyle=\left(\frac{L}{L_{*}}\right)^{1/4},$ (4)
$\displaystyle\frac{r_{\rm trun}}{r_{\rm trun,*}}$
$\displaystyle=\left(\frac{L}{L_{*}}\right)^{\eta}.$ (5)

During our optimization, $\sigma_{*}$, $r_{\rm trun,*}$, and $\eta$ are all
allowed to vary with the range for the $\sigma_{*}$ of
$100.0-400.0\mathrm{~{}km~{}s^{-1}}$, $r_{\rm trun,*}$ from
$1.0^{{}^{\prime\prime}}$ to $100.0^{{}^{\prime\prime}}$ and $\eta$ from $0.2$
to $1.5$.

To check the dependence of our results on the treatment of the BCG, in some
lens mass models we do not explicitly include the BCG as a separate component,
but include it as one of Pseudo-Jaffe ellipsoids of cluster member galaxies
(gals) with its velocity dispersion and the truncation radius determined
following the same scaling relations for cluster member galaxies.

###  III.3 External perturbations and multipole perturbations

In order to achieve a better fit, multipole perturbations are sometimes
included in mass modeling. In this paper, we consider multipole perturbations
of order up to 5 at most (e.g., $m=2$, $3$, $4$, and $5$) with the potential
(Evans and Witt, 2003; Kawano et al., 2004; Congdon and Keeton, 2005; Yoo et
al., 2006)

$\phi=-\frac{\epsilon}{m}r^{2}\cos m(\theta-\theta_{\epsilon}-\pi/2).$ (6)

For each order, $\epsilon$ and $\theta_{\epsilon}$ are model parameters, which
we optimize without any prior ranges. The perturbation with $m=2$ is referred
to as an external shear (Keeton et al., 1997) and named pert in glafic, while
the terms with $m\geq 3$ are named mpole.

###  III.4 Model optimization

In this paper, we consider 16 different lens mass models to fit this cluster-
lensed quasar system in order to check how the constraint on $H_{0}$ depends
on the different assumption on the lens mass model. Specifically, we consider
the following three variation of assumptions on the mass model, (1) which
model is used to describe the dark matter halo; (2) whether to remove the BCG
from all cluster member galaxies and include a separate Hernquist component to
model it; (3) how many orders of perturbations are considered.

More specifically, we use the NFW profile (anfw in GLAFIC) and the Pseudo-
Jaffe Ellipsoid profile (jaffe in GLAFIC) to describe the mass distribution of
dark matter halos. For convenience, we label modeling of the dark matter halo
with anfw or jaffe as 1 or 2. In addition, if we include the brightest cluster
galaxy (BCG) as one of the cluster members, we apply the scaled Pseudo-Jaffe
Ellipsoids profile (gals in GLAFIC) to model mass distribution; we label this
situation as 1. On the contrary, BCG as a separate Hernquist profile (ahern in
GLAFIC) component from the other cluster member galaxies is labeled as 2. As
for perturbations, we consider not only external shear (pert in GLAFIC) but
also multipole perturbations (mpole in GLAFIC) of order up to 5 at most.
Conveniently, labels 1-4 are used to represent which corresponds to pert only,
pert plus mpole with $m=3$, pert plus mpole with $m=3$, $4$, and pert plus
mpole with $m=3$, $4$, $5$, respectively.

We summarize all our mass models and corresponding labels, which describe the
difference of the models more concisely, in Table 2. The label m$ijk$ means
the label $i$ for the dark matter halo, the label $j$ for the cluster member
galaxies (BCG and other members), and the label $k$ for perturbations. For
instance, the label ‘m112’ means that we employ the NFW profile (anfw) to
model the dark matter halo, include BCG as one of the cluster members that are
described by the scaled Pseudo-Jaffe Ellipsoid profile (gals), and the
external shear (pert) and third-pole perturbation (mpole(m=3)) are included.
Therefore we also note it as model m112 (anfw+gals+pert+mpole(m=3)). As
another example, the difference of the model m122
(anfw+ahern+gals+pert+mpole(m=3)) with the model m112 is that we describe the
BCG separately with the Hernquist profile (ahern) instead of treating the BCG
similarly with other cluster member galaxies by the scaled Pseudo-Jaffe
Ellipsoid profile (gals). As yet another example, the difference of the model
m222 (jaffe+ahern+gals+pert+mpole(m=3)) with the model m122 is that the
Pseudo-Jaffe Ellipsoid profile (jaffe) is used to parameterised the dark
matter halo instead of the NFW (anfw) profile. If we include an additional
order of multipole perturbations to the model m222, the new model is referred
to as m223 (jaffe+ahern+gals+pert+mpole(m=3,4)).

In cluster strong lens modeling, it is customary to assume positional errors
that are much larger than measurement errors of positions of multiple images,
in order to take account of the complexity of cluster mass distributions due
to e.g., substructures in clusters and the line-of-sight matter fluctuations.
For massive clusters, a typical scatter of multiple image positions due to
such complexity is $\sim 0.5^{\prime\prime}-1^{\prime\prime}$ (e.g., Remolina
González et al. (2018)). Because of the lack of the guidance on what
positional errors should be assumed for SDSS J1004+4112, we choose positional
errors for each of our mass models such that the value of reduced-$\chi^{2}$
for the best-fitting model becomes roughly equal to one. Following Oguri
(2010), we also choose smaller positional errors for quasar multiple images
than for galaxy multiple images, because accurate quasar image positions are
important for accurate predictions of time delays between multiple images.
Table 2 summarizes specific positional errors that we adopt for the 16
different lens mass models.

In our mass modeling, we have 71 observational constraints. Naturally,
different mass models use different model parameters, and as a result the
degree of freedom (dof) is different for different mass models. For example,
for the model m111 (anfw+gals+pert) in which we use the NFW profile to
describe the dark matter halo, do not employ the Hernquist profile to model
the BCG separately from all the cluster member galaxies, and just consider the
external shear (i.e., pert), there are 28 model parameters and the degree of
freedom is 43. For model m224 (jaffe+ahern+gals+pert+mpole(m=3,4,5)) in which
the dark matter halo is parameterised by the Pseudo-Jaffe Ellipsoid profile,
the Hernquist profile is employed to model the BCG, and we not only consider
external shear but also third-, fourth-, and fifth-order multipole
perturbations, we have 39 model parameters leaving with 32 degree of freedom.

In our all calculations and model parameter optimizations, we use a standard
$\chi^{2}$ minimization method implemented in the glafic software and estimate
$\chi^{2}$ in the source plane by taking the full account of the magnification
tensor at each image position (see Oguri (2010)) in order to speed up the
calculations.

Table 3:  The Hubble constant with 68.3$\%$ confidence interval and reduced-$\chi^{2}$ ($\chi^{2}$/dof) values for the 16 different lens mass models.  Model |  $H_{0}$ |  $\chi^{2}$/dof
---|---|---
|  [$\mathrm{~{}km~{}s^{-1}~{}Mpc^{-1}}$] |
m111 |  $71.4^{+5.0}_{-5.0}$ |  46.7/43
m112 |  $59.2^{+1.5}_{-3.5}$ |  36.0/41
m113 |  $57.5^{+2.0}_{-3.0}$ |  41.6/39
m114 |  $58.0^{+2.0}_{-3.0}$ |  31.1/37
m121 |  $71.1^{+2.0}_{-3.0}$ |  38.2/39
m122 |  $62.4^{+3.5}_{-3.5}$ |  37.5/37
m123 |  $60.9^{+2.5}_{-2.0}$ |  38.7/35
m124 |  $60.8^{+2.0}_{-2.0}$ |  31.0/33
m211 |  $88.9^{+4.0}_{-3.0}$ |  53.8/42
m212 |  $88.7^{+3.5}_{-3.0}$ |  39.1/40
m213 |  $79.5^{+4.0}_{-2.5}$ |  38.6/38
m214 |  $79.6^{+3.0}_{-3.0}$ |  38.1/36
m221 |  $76.7^{+3.5}_{-3.5}$ |  45.8/38
m222 |  $74.8^{+3.5}_{-3.5}$ |  42.8/36
m223 |  $64.2^{+3.5}_{-3.5}$ |  39.6/34
m224 |  $62.1^{+3.0}_{-2.5}$ |  29.7/32

Figure 2:  Top: The summary of the constraint on the Hubble constant for all
the 16 different lens mass models. Bottom: Similar to the top panel, but for
assumed position errors of quasar images $\sigma_{\rm q}$ in 16 different lens
mass models with measurements of the Hubble constant.

Figure 3: Top: The probability distribution functions (PDFs) of the Hubble
constant for all the 16 different lens mass models (dashed), as well as the
PDF after combining all the 16 mass models with equal weighting (solid). The
vertical gray dotted lines and shaded region indicate the median and $68.3\%$
confidence interval for the PDF of the combined result. Middle: Similar to the
top panel, but combining all the mass models with different weighting. Bottom:
Also similar to the top panel, but for 2 best mass models that best reproduce
the observed shapes of the lensed quasar host galaxy and the lensed galaxies
behind the cluster.

##  IV Hubble constant from the time delay cosmology

Figure 4: Radial convergence (left) and tangential shear (right) profilers for
the 16 different lens mass models, which are computed assuming the source
redshift of $z_{\rm s}=1$. Grey dots with errors represent tangential shear
measurements from previous weak gravitational leasing analysis Oguri et al.
(2012).

Thanks to the progress of observations as well as theoretical understanding of
the lens model dependence, time-delay cosmography has emerged as an important
independent technique to constrain the Hubble constant. For a strong lensing
system with a background quasar as a source and a cluster as a lens,
variability patterns of a quasar for different multiple images have some time
lag because of different light paths for the different multiple images. The
time delay between image $i$ and $j$ is calculated as (e.g., Schneider (1985))

$\Delta t_{ij}=\frac{D_{\Delta{\rm t}}}{c}\Delta\phi_{ij},$ (7)

where the time delay distance $D_{\Delta{\rm t}}$ is defined as (e.g.,
(Refsdal, 1964; Schneider et al., 1992))

$D_{\Delta{\rm t}}\equiv(1+z_{\rm l})\frac{D_{\rm l}^{\rm A}D_{\rm s}^{\rm
A}}{D_{\rm ls}^{\rm A}},$ (8)

with $z_{\rm l}$ being the redshift of the lens, $D_{\rm l}^{\rm A}$, $D_{\rm
s}^{\rm A}$, and $D_{\rm ls}^{\rm A}$ being angular diameter distances from
the observer to the lens, from the observer to the source, and between the
lens and the source, respectively. Representing the angular image position in
the image plane and the angular source position in the source plane by
$\vec{\theta}$ and $\vec{\beta}$, respectively, the Fermat potential
difference $\Delta\phi_{ij}$ is calculated as (e.g., Blandford and Narayan
(1986))

$\Delta\phi_{ij}=\phi(\vec{\theta}_{i};\vec{\beta})-\phi(\vec{\theta}_{j};\vec{\beta}),$
(9)
$\phi(\vec{\theta};\vec{\beta})=\frac{1}{2}(\vec{\theta}-\vec{\beta})^{2}-\psi(\vec{\theta}),$
(10)

where $\psi(\vec{\theta})$ indicates the gravitational lens potential.

Given the definition of the angular diameter distance in the Friedmann-
Lemaitre-Robertson-Walker (FLRW) metric

$D^{\rm A}(z_{\rm 1},z_{2})=\frac{1}{1+z_{2}}f_{\rm K}[X(z_{1},z_{2})],$ (11)

where the spatial curvature is described by $K$,

$f_{\rm
K}(X)=\left\\{\begin{array}[]{lll}K^{-1/2}\sin(K^{1/2}X),~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}K>0\\\
X,~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}K=0\\\
(-K)^{-1/2}\sinh[(-K)^{1/2}X],~{}K<0,\\\ \end{array}\right.$ (12)

and

$X(z_{1},z_{2})=\frac{c}{H_{0}}\int_{z_{1}}^{z_{2}}E(z^{\prime})dz^{\prime},$
(13)

with the dimensionless Friedman equation being $E(z)$, the Hubble constant is
associated with time delays in the flat Universe ($K=0$) as follows

$H_{0}=\frac{1}{\Delta t_{ij}}\frac{\int_{0}^{z_{\rm
l}}E(z^{\prime})dz^{\prime}\int_{0}^{z_{\rm
s}}E(z^{\prime})dz^{\prime}}{\int_{z_{\rm l}}^{z_{\rm
s}}E(z^{\prime})dz^{\prime}}\Delta\phi_{ij}.$ (14)

It is clearly seen that $H_{0}$ can be constrained from the measurements of
time delays, as long as the Fermat potential $\Delta\phi_{ij}$ is well
constrained from lens mass modeling. Throughout our analysis, we assume a flat
Universe with $\Omega_{\rm m}=0.3$ and $\Omega_{\Lambda}=0.7.$

In practical analysis, we include $H_{0}$ as a model parameter and vary it
together with all the other model parameters to simultaneously fit the image
positions, magnitude differences, and time delays that are summarized in Table
1. The error on $H_{0}$ is estimated by changing the $H_{0}$ value around the
best-fit value with a step size of $0.5\mathrm{~{}km~{}s^{-1}~{}Mpc^{-1}}$,
and for each value of $H_{0}$ we optimize the other model parameters to obtain
the best-fit $\chi^{2}$, and derive the $1\sigma$ error by the range where the
difference of the best-fit $\chi^{2}$, $\Delta\chi^{2}$, is smaller than $1$,
assuming Gaussian errors.

We summarize our results in Table 3 and Fig. 2. We find that the best-fit
$H_{0}$ values are significantly different for different lens mass models and
are much larger than the statistical errors derived for individual lens mass
models. For instance, when the NFW profile is used to model the dark matter
halo and using the external shear to describe the perturbation (m111), we
obtain $H_{0}=71.4^{+5.0}_{-5.0}\mathrm{~{}km~{}s^{-1}~{}Mpc^{-1}}$. When more
multipole perturbations are added, it becomes
$H_{0}=59.2^{+1.5}_{-3.5}\mathrm{~{}km~{}s^{-1}~{}Mpc^{-1}}$ for m112, and
$H_{0}=57.5^{+2.0}_{-3.0}\mathrm{~{}km~{}s^{-1}~{}Mpc^{-1}}$ for m113. We find
that the best-fit Hubble constant becomes smaller when more multipole
perturbations are added. Interestingly, this trend appears to be common for
all the mass models we consider in this paper, and is also seen in the bottom
panel of Fig. 2.

We find that changing the halo density profile from the NFW profile to the
Pseudo-Jaffe Ellipsoid profile has a larger impact on the resulting $H_{0}$
value. For instance, we obtain
$H_{0}=79.6^{+3.0}_{-3.0}\mathrm{~{}km~{}s^{-1}~{}Mpc^{-1}}$ in lens modeling
combination m214, which is significantly larger than $H_{0}$ for m114. Fig. 2
clearly indicates that $H_{0}$ values tend to be much higher in general when
the Pseudo-Jaffe Ellipsoid profile is used to model the dark matter halo.
Modelling the BCG separately as the Hernquist profile also has some impact on
the result, e.g., $H_{0}=62.1^{+3.0}_{-2.5}\mathrm{~{}km~{}s^{-1}~{}Mpc^{-1}}$
is derived from m224, which is much smaller than the $H_{0}$ value obtained
for m214.

Combining results from these different lens mass models require the
determination of how the weight is assigned to each lens mass model. For
instance, in the analysis of SN Refsdal Kelly et al. (2023), a weight for each
mass model is assigned based on the reproducibility of flux ratios and time
delay ratios of multiple supernova images that are not included in the mass
modeling. In our case, we cannot use the flux ratios and time delay ratios for
weighting as they are explicitly included as model constraints. Given the lack
of the clear guidance on how to assign weights to individual lens mass models,
first we simply sum up posteriors $\exp(-\Delta\chi^{2}/2)$ of the Hubble
constant from individual mass models with equal weighting. The combined
posterior probability distribution function (PDF) of $H_{0}$, as well as PDFs
for individual lens mass models, are shown in Fig. 3. The constraint on
$H_{0}$ from this combined PDF is
$H_{0}=67.5^{+14.5}_{-8.9}\mathrm{~{}km~{}s^{-1}~{}Mpc^{-1}}$, which
represents a weak constraint on $H_{0}$ due to the large scatter among
different lens mass models.

To see how the result depends on weighting, we next combine all the lens mass
models with weighting based on the assumed positional errors of quasars images
$\sigma_{\rm q}$. Specifically, we adopt the weighting value of
$w_{i}=1/{\sigma^{2}_{\rm q}}_{,i}$, where ${\sigma_{\rm q}}_{,i}$ are
positional errors of individual mass models quoted in Table 2. With this
weighting, we obtain
$H_{0}=62.0^{+15.7}_{-4.2}\mathrm{~{}km~{}s^{-1}~{}Mpc^{-1}}$ (see Fig. 3),
which represents a slightly tighter constraint than for equal weighting. While
the slight improvement comes from the trend of $H_{0}$ with the positional
errors as shown in Fig. 2, we conclude that the error on the Hubble constant
is still large even with weighting based on positional errors of quasar
images.

We note that the 16 lens mass models explored in this paper do not cover the
full lens model uncertainty, and combining these 16 lens mass model results
with equal weighting or weighting with positional errors is also not well
motivated, indicating that this constraint on $H_{0}$ should be taken with
caution. The multiple peaks in the PDF indeed imply that we do not yet fully
marginalize over the full lens model uncertainty. Nevertheless, our analysis
at least suggests that the lens mass model uncertainty is quite large for the
current observational constraints and the treatment of positional errors.
These findings motivate us to further discuss other possibilities to
distinguish or exclude some of the lens mass models to further improve
constraints on the Hubble constant.

Figure 5: The predicated shapes of the lensed host galaxy of the quasar and
lensed background galaxies for all the 16 different lens mass models. Shapes
enclosed in red ellipses indicate arcs whose shapes are relatively more
sensitive to the difference of lens mass models.

##  V Possibility of improving Hubble constant with additional constraints

First, it is possible that different lens mass models predict quite different
outer density profiles of the lensing cluster, which can be constrained by
weak gravitational lensing. To explore this possibility, we calculate radial
convergence and tangential shear profiles of all the 16 different lens mass
models. Since the mass distribution is not circular symmetric, we compute the
azimuthal average of the convergence at each radius to derive radial
convergence and tangential shear profiles. The results shown in Fig. 4
indicate that the convergence and tangential shear profiles are quite similar
between different lens mass models, implying that they cannot be distinguished
by weak gravitational lensing. For comparison, we plot results from the weak
gravitational lensing analysis in Ref. Oguri et al. (2012) using Subaru
Suprime-Cam images, which indicates that measurements errors are much larger
compared with the difference of tangential shear profiles among different
models. In order to distinguish these different models with weak lensing, we
need much improved weak lensing measurements with deep imaging from space
telescopes such as James Webb Space Telescope, and hence is challenging.

In addition, we predict the shapes of lensed quasar host galaxies and lensed
background galaxies, which are shown in Fig. 5. We find that there are obvious
differences of those shapes between different lens mass models. For instance,
as multipole perturbations increase, the arcs become more extended in most
cases. There are also notable differences between different models to describe
the dark matter halo such that the arcs are smaller when using the Pseudo-
Jaffe Ellipsoid profile than using the NFW profile. In addition, some slight
differences are presented when we model the BCG with the Hernquist profile
separately. By comparing our models prediations with observations, we find
that the lens mass models m113 and m123 best match the observed shapes of
lensed quasar host galaxy and lensed galaxies.

When we restrict our mass models to m113 and m123 and exclude the other
models, the combined constraint becomes much tighter,
$H_{0}=59.1^{+3.6}_{-3.5}\mathrm{~{}km~{}s^{-1}~{}Mpc^{-1}}$. The PDF of
$H_{0}$ for this case is also shown in Fig. 3. Again we caution that these two
mass models perhaps do not fully cover the lens model uncertainty. However we
can at least conclude that adding constraints from shapes of lensed galaxies
can greatly help breaking degeneracies between different lens mass models.

##  VI Discussions and Conclusion

In this paper, we explore the possibility of constraining the Hubble constant
using time-delay cosmography with a cluster-lensed quasar system. As a
specific example, we focus on the first discovered cluster-lensed quasar lens
system SDSS J1004+4112, in which a background quasar is lensed into five
multiple images. In addition, several multiple images of galaxies behind the
lensing cluster have also been identified.

In order to check how the constraint on the Hubble constant depends on the
choice of the lens mass model, we employ 16 different lens mass models in
which different assumptions on profiles of the dark matter halo, the BCG, and
multipole perturbations are made. We find that the variation of the best-fit
values of the Hubble constant is very large among the 16 lens mass models.
Interestingly, we find that the value of the best-fit Hubble constant
decreases as the complexity of the multipole perturbation increases, whose
origin is to be understood. By summing posteriors of the Hubble constant with
equal weighting, the combined constraint is
$H_{0}=67.5^{+14.5}_{-8.9}\mathrm{~{}km~{}s^{-1}~{}Mpc^{-1}}$. We also
consider different weighting values based on positional errors of quasar
images, $w_{i}=1/{\sigma^{2}_{\rm q}}_{,i}$, to obtain a slightly tighter
constraint on $H_{0}=62.0^{+15.7}_{-4.2}\mathrm{~{}km~{}s^{-1}~{}Mpc^{-1}}$,
which does not represent a significant improvement.

We discuss possible ways to distinguish different lens mass models by adding
other observational constraints. We find that there are no obvious differences
in radial density profiles and tangential shear profiles in the outer region
that can be constrained by weak gravitational lensing. In contrast, we find
that there are distinct differences in critical curves, and also predicted
shapes of quasar host galaxies as well as other lensed galaxies. By comparing
those shapes with observations, we find that there are two models (m113 and
m123) that best match the observed shapes. By combining the constraints from
these two models, we obtain much tighter constraint of
$H_{0}=59.1^{+3.6}_{-3.5}\mathrm{~{}km~{}s^{-1}~{}Mpc^{-1}}$. We emphasize
that our analysis with 16 mass models does not fully cover the whole mass
model uncertainty and our results are also dependent on how to assign weight
to each mass model. Therefore our result should be taken with caution. We
expect that our results at least indicate that incorporating as many lensing
constraints as possible, most notably those from shapes of lensed extended
objects, is important for obtaining tight constraints on the Hubble constant
from cluster-lensed quasar lens systems.

The importance of including constraints from extended arcs has also been
emphasized for galaxy-scale quasar lens systems Kochanek et al. (2001); Suyu
et al. (2010, 2017). Our analysis indicates that this may also be the case for
the cluster-scale lens system. However, one technical difficulty in the
analysis of cluster-scale lens systems lies in the fact that we cannot
perfectly reproduced observed image positions in cluster lens modeling, which
mainly originates from small-scale matter distributions that cause random
scatter in multiple image positions. In presence of significant mismatches of
image positions between model predictions and observations, it is not obvious
how to rigorously incorporating such shape information in lens modeling,
because we cannot simply conduct pixel-by-pixel fitting of lensed arc shapes.
How to efficiently and rigorously include such extended arc information in
cluster strong lens modeling, which requires understanding of how such small-
scale matter distributions not only affect their positions but also their
shapes, is an open question, which needs to be overcome for turning cluster-
lensed quasar systems into useful and accurate probe of the Hubble constant.
In order to obtain accurate constraints on the Hubble constant from time
delays, accurate and realistic error bars on time delays are also important.

A caveat is that the situation could be different for more massive clusters
(e.g., MACS J1149+2223) where there are many more multiple images, typically
more than 100, to tightly constrain lens mass distributions. The large mass
model dependence for SDSS J1004+4112 may originate from the small number of
multiple images that we can use. Deep imaging of cluster-lensed quasar systems
with e.g., James Webb Space Telescope may increase the number of multiple
images for SDSS J1004+4112 and other cluster-lensed quasar systems, which can
help reduce the lens mass model dependence on the Hubble constant measurement.

###### Acknowledgements.

We would like to thank the referee for useful comments and suggestion. This
work was supported by the National Natural Science Foundation of China under
Grants Nos. 12021003, 11690023, and 11920101003; the Strategic Priority
Research Program of the Chinese Academy of Sciences, Grant No. XDB23000000;
the Interdiscipline Research Funds of Beijing Normal University. YL was
supported by the Interdiscipline Research Funds of Beijing Normal University
(Grant No. BNUXKJC2017) and China Scholarship Council (Grant No.
202106040084). MO was supported by JSPS KAKENHI Grant Numbers JP22H01260,
JP20H05856, and JP22K21349.

## References

  * Weinberg et al. (2013) D. H. Weinberg, M. J. Mortonson, D. J. Eisenstein, C. Hirata, A. G. Riess, and E. Rozo, Phys. Rep. 530, 87 (2013), eprint 1201.2434.
  * Freedman et al. (2001) W. L. Freedman, B. F. Madore, B. K. Gibson, L. Ferrarese, D. D. Kelson, S. Sakai, J. R. Mould, J. Kennicutt, Robert C., H. C. Ford, J. A. Graham, et al., ApJ 553, 47 (2001), eprint astro-ph/0012376.
  * Planck Collaboration et al. (2020) Planck Collaboration, N. Aghanim, Y. Akrami, M. Ashdown, J. Aumont, C. Baccigalupi, M. Ballardini, A. J. Banday, R. B. Barreiro, N. Bartolo, et al., A&A 641, A6 (2020), eprint 1807.06209.
  * Riess et al. (1998) A. G. Riess, A. V. Filippenko, P. Challis, A. Clocchiatti, A. Diercks, P. M. Garnavich, R. L. Gilliland, C. J. Hogan, S. Jha, R. P. Kirshner, et al., AJ 116, 1009 (1998).
  * Perlmutter et al. (1999) S. Perlmutter, G. Aldering, G. Goldhaber, R. A. Knop, P. Nugent, P. G. Castro, S. Deustua, S. Fabbro, A. Goobar, D. E. Groom, et al., ApJ 517, 565 (1999).
  * Cao et al. (2015) S. Cao, M. Biesiada, R. Gavazzi, A. Piórkowska, and Z.-H. Zhu, ApJ 806, 185 (2015), eprint 1509.07649.
  * Troxel et al. (2018) M. A. Troxel, N. MacCrann, J. Zuntz, T. F. Eifler, E. Krause, S. Dodelson, D. Gruen, J. Blazek, O. Friedrich, S. Samuroff, et al., Phys. Rev. D 98, 043528 (2018).
  * Scolnic et al. (2018) D. M. Scolnic, D. O. Jones, A. Rest, Y. C. Pan, R. Chornock, R. J. Foley, M. E. Huber, R. Kessler, G. Narayan, A. G. Riess, et al., ApJ 859, 101 (2018).
  * Bianchini et al. (2020) F. Bianchini, W. L. K. Wu, P. A. R. Ade, A. J. Anderson, J. E. Austermann, J. S. Avva, L. Balkenhol, E. Baxter, J. A. Beall, A. N. Bender, et al., Phys. Rev. D 102, 083504 (2020).
  * Han et al. (2021) D. Han, N. Sehgal, A. MacInnis, A. van Engelen, B. D. Sherwin, M. S. Madhavacheril, S. Aiola, N. Battaglia, J. A. Beall, D. T. Becker, et al., J. Cosmology Astropart. Phys. 2021, 031 (2021).
  * Asgari et al. (2021) M. Asgari, C.-A. Lin, B. Joachimi, B. Giblin, C. Heymans, H. Hildebrandt, A. Kannawadi, B. Stölzner, T. Tröster, J. L. van den Busch, et al., A&A 645, A104 (2021).
  * Mossa et al. (2020) V. Mossa, K. Stöckel, F. Cavanna, F. Ferraro, M. Aliotta, F. Barile, D. Bemmerer, A. Best, A. Boeltzig, C. Broggini, et al., Nature 587, 210 (2020).
  * Riess et al. (2022) A. G. Riess, W. Yuan, L. M. Macri, D. Scolnic, D. Brout, S. Casertano, D. O. Jones, Y. Murakami, G. S. Anand, L. Breuval, et al., ApJ 934, L7 (2022), eprint 2112.04510.
  * Freedman et al. (2019) W. L. Freedman, B. F. Madore, D. Hatt, T. J. Hoyt, I. S. Jang, R. L. Beaton, C. R. Burns, M. G. Lee, A. J. Monson, J. R. Neeley, et al., ApJ 882, 34 (2019), eprint 1907.05922.
  * Freedman (2021) W. L. Freedman, ApJ 919, 16 (2021), eprint 2106.15656.
  * Braatz and Gugliucci (2008) J. A. Braatz and N. E. Gugliucci, ApJ 678, 96 (2008), eprint 0804.4136.
  * Reid et al. (2013) M. J. Reid, J. A. Braatz, J. J. Condon, K. Y. Lo, C. Y. Kuo, C. M. V. Impellizzeri, and C. Henkel, ApJ 767, 154 (2013), eprint 1207.7292.
  * Kuo et al. (2013) C. Y. Kuo, J. A. Braatz, M. J. Reid, K. Y. Lo, J. J. Condon, C. M. V. Impellizzeri, and C. Henkel, ApJ 767, 155 (2013), eprint 1207.7273.
  * Kuo et al. (2015) C. Y. Kuo, J. A. Braatz, K. Y. Lo, M. J. Reid, S. H. Suyu, D. W. Pesce, J. J. Condon, C. Henkel, and C. M. V. Impellizzeri, ApJ 800, 26 (2015), eprint 1411.5106.
  * Gao et al. (2017) F. Gao, J. A. Braatz, M. J. Reid, J. J. Condon, J. E. Greene, C. Henkel, C. M. V. Impellizzeri, K. Y. Lo, C. Y. Kuo, D. W. Pesce, et al., ApJ 834, 52 (2017), eprint 1610.06802.
  * Pesce et al. (2020) D. W. Pesce, J. A. Braatz, M. J. Reid, A. G. Riess, D. Scolnic, J. J. Condon, F. Gao, C. Henkel, C. M. V. Impellizzeri, C. Y. Kuo, et al., ApJ 891, L1 (2020), eprint 2001.09213.
  * Dietrich et al. (2020) T. Dietrich, M. W. Coughlin, P. T. H. Pang, M. Bulla, J. Heinzel, L. Issa, I. Tews, and S. Antier, Science 370, 1450 (2020).
  * The LIGO Scientific Collaboration et al. (2021) The LIGO Scientific Collaboration, the Virgo Collaboration, the KAGRA Collaboration, R. Abbott, H. Abe, F. Acernese, K. Ackley, N. Adhikari, R. X. Adhikari, V. K. Adkins, et al., arXiv e-prints arXiv:2111.03604 (2021), eprint 2111.03604.
  * Refsdal (1964) S. Refsdal, MNRAS 128, 307 (1964).
  * Kelly et al. (2023) P. L. Kelly, S. Rodney, T. Treu, M. Oguri, W. Chen, A. Zitrin, S. Birrer, V. Bonvin, L. Dessart, J. M. Diego, et al., Science 380, abh1322 (2023), eprint 2305.06367.
  * Wong et al. (2020) K. C. Wong, S. H. Suyu, G. C. F. Chen, C. E. Rusu, M. Millon, D. Sluse, V. Bonvin, C. D. Fassnacht, S. Taubenberger, M. W. Auger, et al., MNRAS 498, 1420 (2020), eprint 1907.04869.
  * Birrer et al. (2020) S. Birrer, A. J. Shajib, A. Galan, M. Millon, T. Treu, A. Agnello, M. Auger, G. C. F. Chen, L. Christensen, T. Collett, et al., A&A 643, A165 (2020), eprint 2007.02941.
  * Shajib et al. (2023) A. J. Shajib, P. Mozumdar, G. C. F. Chen, T. Treu, M. Cappellari, S. Knabel, S. H. Suyu, V. N. Bennert, J. A. Frieman, D. Sluse, et al., A&A 673, A9 (2023), eprint 2301.02656.
  * Sharon et al. (2005) K. Sharon, E. O. Ofek, G. P. Smith, T. Broadhurst, D. Maoz, C. S. Kochanek, M. Oguri, Y. Suto, N. Inada, and E. E. Falco, ApJ 629, L73 (2005), eprint astro-ph/0507360.
  * Sharon et al. (2017) K. Sharon, M. B. Bayliss, H. Dahle, M. K. Florian, M. D. Gladders, T. L. Johnson, R. Paterno-Mahler, J. R. Rigby, K. E. Whitaker, and E. Wuyts, ApJ 835, 5 (2017), eprint 1609.08848.
  * Acebron et al. (2022a) A. Acebron, C. Grillo, P. Bergamini, A. Mercurio, P. Rosati, G. B. Caminha, P. Tozzi, G. B. Brammer, M. Meneghetti, A. Morelli, et al., ApJ 926, 86 (2022a), eprint 2111.05871.
  * Acebron et al. (2022b) A. Acebron, C. Grillo, P. Bergamini, G. B. Caminha, P. Tozzi, A. Mercurio, P. Rosati, G. Brammer, M. Meneghetti, M. Nonino, et al., A&A 668, A142 (2022b), eprint 2208.13788.
  * Meneghetti et al. (2017) M. Meneghetti, P. Natarajan, D. Coe, E. Contini, G. De Lucia, C. Giocoli, A. Acebron, S. Borgani, M. Bradac, J. M. Diego, et al., MNRAS 472, 3177 (2017), eprint 1606.04548.
  * Denzel et al. (2021) P. Denzel, J. P. Coles, P. Saha, and L. L. R. Williams, MNRAS 501, 784 (2021), eprint 2007.14398.
  * Inada et al. (2003) N. Inada, M. Oguri, B. Pindor, J. F. Hennawi, K. Chiu, W. Zheng, S.-I. Ichikawa, M. D. Gregg, R. H. Becker, Y. Suto, et al., Nature 426, 810 (2003), eprint astro-ph/0312427.
  * Oguri et al. (2004) M. Oguri, N. Inada, C. R. Keeton, B. Pindor, J. F. Hennawi, M. D. Gregg, R. H. Becker, K. Chiu, W. Zheng, S.-I. Ichikawa, et al., ApJ 605, 78 (2004), eprint astro-ph/0312429.
  * Inada et al. (2005) N. Inada, M. Oguri, C. R. Keeton, D. J. Eisenstein, F. J. Castander, K. Chiu, P. B. Hall, J. F. Hennawi, D. E. Johnston, B. Pindor, et al., PASJ 57, L7 (2005), eprint astro-ph/0503310.
  * Fohlmeister et al. (2008) J. Fohlmeister, C. S. Kochanek, E. E. Falco, C. W. Morgan, and J. Wambsganss, ApJ 676, 761 (2008), eprint 0710.1634.
  * Napier et al. (2023a) K. Napier, K. Sharon, H. Dahle, M. Bayliss, M. D. Gladders, G. Mahler, J. R. Rigby, and M. Florian, arXiv e-prints arXiv:2301.11240 (2023a), eprint 2301.11240.
  * Inada et al. (2006) N. Inada, M. Oguri, T. Morokuma, M. Doi, N. Yasuda, R. H. Becker, G. T. Richards, C. S. Kochanek, I. Kayo, K. Konishi, et al., ApJ 653, L97 (2006), eprint astro-ph/0611275.
  * Oguri et al. (2008) M. Oguri, E. O. Ofek, N. Inada, T. Morokuma, E. E. Falco, C. S. Kochanek, I. Kayo, T. Broadhurst, and G. T. Richards, ApJ 676, L1 (2008), eprint 0802.0002.
  * Fohlmeister et al. (2013) J. Fohlmeister, C. S. Kochanek, E. E. Falco, J. Wambsganss, M. Oguri, and X. Dai, ApJ 764, 186 (2013), eprint 1207.5776.
  * Dahle et al. (2013) H. Dahle, M. D. Gladders, K. Sharon, M. B. Bayliss, E. Wuyts, L. E. Abramson, B. P. Koester, N. Groeneboom, T. E. Brinckmann, M. T. Kristensen, et al., ApJ 773, 146 (2013), eprint 1211.1091.
  * Dahle et al. (2015) H. Dahle, M. D. Gladders, K. Sharon, M. B. Bayliss, and J. R. Rigby, ApJ 813, 67 (2015), eprint 1505.06187.
  * Martínez-Arrizabalaga et al. (2023) J. Martínez-Arrizabalaga, J. M. Diego, and L. J. Goicoechea, arXiv e-prints arXiv:2309.14776 (2023), eprint 2309.14776.
  * Muñoz et al. (2022) J. A. Muñoz, C. S. Kochanek, J. Fohlmeister, J. Wambsganss, E. Falco, and R. Forés-Toribio, ApJ 937, 34 (2022), eprint 2206.08597.
  * Forés-Toribio et al. (2022) R. Forés-Toribio, J. A. Muñoz, C. S. Kochanek, and E. Mediavilla, ApJ 937, 35 (2022), eprint 2206.09856.
  * Shu et al. (2018) Y. Shu, R. Marques-Chaves, N. W. Evans, and I. Pérez-Fournon, MNRAS 481, L136 (2018), eprint 1809.07337.
  * Shu et al. (2019) Y. Shu, S. E. Koposov, N. W. Evans, V. Belokurov, R. G. McMahon, M. W. Auger, and C. A. Lemon, MNRAS 489, 4741 (2019), eprint 1909.02010.
  * Martinez et al. (2023) M. N. Martinez, K. A. Napier, A. P. Cloonan, E. Sukay, K. Gozman, K. Merz, G. Khullar, J. J. Lin, O. S. Matthews Acuña, E. Medina, et al., ApJ 946, 63 (2023), eprint 2209.03972.
  * Napier et al. (2023b) K. Napier, M. Gladders, K. Sharon, H. Dahle, A. P. Cloonan, G. Mahler, I. Escapa, J. Garza, A. Kisare, N. Malagon, et al., arXiv e-prints arXiv:2305.14317 (2023b), eprint 2305.14317.
  * York et al. (2000) D. G. York, J. Adelman, J. Anderson, John E., and S. F. e. a. Anderson, AJ 120, 1579 (2000), eprint astro-ph/0006396.
  * Abazajian et al. (2004) K. Abazajian, J. K. Adelman-McCarthy, M. A. Agüeros, S. S. Allam, K. Anderson, S. F. Anderson, J. Annis, N. A. Bahcall, I. K. Baldry, S. Bastian, et al., AJ 128, 502 (2004), eprint astro-ph/0403325.
  * Inada et al. (2008) N. Inada, M. Oguri, E. E. Falco, T. J. Broadhurst, E. O. Ofek, C. S. Kochanek, K. Sharon, and G. P. Smith, PASJ 60, 27 (2008), eprint 0808.1769.
  * Liesenborgs et al. (2009) J. Liesenborgs, S. de Rijcke, H. Dejonghe, and P. Bekaert, MNRAS 397, 341 (2009), eprint 0904.2382.
  * Oguri (2010) M. Oguri, PASJ 62, 1017 (2010), eprint 1005.3103.
  * Fohlmeister et al. (2007) J. Fohlmeister, C. S. Kochanek, E. E. Falco, J. Wambsganss, N. Morgan, C. W. Morgan, E. O. Ofek, D. Maoz, C. R. Keeton, J. C. Barentine, et al., ApJ 662, 62 (2007).
  * Bertin and Arnouts (1996) E. Bertin and S. Arnouts, A&AS 117, 393 (1996).
  * Oguri (2021) M. Oguri, PASP 133, 074504 (2021), eprint 2106.11464.
  * Navarro et al. (1997) J. F. Navarro, C. S. Frenk, and S. D. M. White, ApJ 490, 493 (1997), eprint astro-ph/9611107.
  * Jaffe (1983) W. Jaffe, MNRAS 202, 995 (1983).
  * Keeton (2001) C. R. Keeton, arXiv e-prints astro-ph/0102341 (2001), eprint astro-ph/0102341.
  * Hernquist (1990) L. Hernquist, ApJ 356, 359 (1990).
  * Evans and Witt (2003) N. W. Evans and H. J. Witt, MNRAS 345, 1351 (2003), eprint astro-ph/0212013.
  * Kawano et al. (2004) Y. Kawano, M. Oguri, T. Matsubara, and S. Ikeuchi, PASJ 56, 253 (2004), eprint astro-ph/0404013.
  * Congdon and Keeton (2005) A. B. Congdon and C. R. Keeton, MNRAS 364, 1459 (2005), eprint astro-ph/0510232.
  * Yoo et al. (2006) J. Yoo, C. S. Kochanek, E. E. Falco, and B. A. McLeod, ApJ 642, 22 (2006), eprint astro-ph/0511001.
  * Keeton et al. (1997) C. R. Keeton, C. S. Kochanek, and U. Seljak, ApJ 482, 604 (1997), eprint astro-ph/9610163.
  * Remolina González et al. (2018) J. D. Remolina González, K. Sharon, and G. Mahler, ApJ 863, 60 (2018), eprint 1807.03291.
  * Oguri et al. (2012) M. Oguri, M. B. Bayliss, H. Dahle, K. Sharon, M. D. Gladders, P. Natarajan, J. F. Hennawi, and B. P. Koester, MNRAS 420, 3213 (2012), eprint 1109.2594.
  * Schneider (1985) P. Schneider, A&A 143, 413 (1985).
  * Schneider et al. (1992) P. Schneider, J. Ehlers, and E. E. Falco, _Gravitational Lenses_ (1992).
  * Blandford and Narayan (1986) R. Blandford and R. Narayan, ApJ 310, 568 (1986).
  * Kochanek et al. (2001) C. S. Kochanek, C. R. Keeton, and B. A. McLeod, ApJ 547, 50 (2001), eprint astro-ph/0006116.
  * Suyu et al. (2010) S. H. Suyu, P. J. Marshall, M. W. Auger, S. Hilbert, R. D. Blandford, L. V. E. Koopmans, C. D. Fassnacht, and T. Treu, ApJ 711, 201 (2010), eprint 0910.2773.
  * Suyu et al. (2017) S. H. Suyu, V. Bonvin, F. Courbin, C. D. Fassnacht, C. E. Rusu, D. Sluse, T. Treu, K. C. Wong, M. W. Auger, X. Ding, et al., MNRAS 468, 2590 (2017), eprint 1607.00017.


A simple model for ultradiscrete Hopf bifurcation

Shousuke Ohmori∗) and Yoshihiro Yamazaki

Department of Physics, Waseda University, Shinjuku, Tokyo 169-8555, Japan

*corresponding author<EMAIL_ADDRESS>

Abstract

Dynamical properties of ultradiscrete Hopf bifurcation, similar to those of
the standard Hopf bifurcation, are discussed by proposing a simple model of
ultradiscrete equations with max-plus algebra. In ultradiscrete Hopf
bifurcation, limit cycles emerge depending on the value of a bifurcation
parameter in the model. The limit cycles are composed of a finite number of
discrete states. Furthermore, the model exhibits excitability. The model is
derived from two different dynamical models with Hopf bifurcation by means of
ultradiscretization; it is a candidate for a normal form for ultradiscrete
Hopf bifurcation.

For analysis of dynamical properties of nonlinear equations, piecewise
linearization of them has often been carried out. Especially,
ultradiscretization, one of piecewise linearization methods, successfully
retains and elucidates the essence of dynamical structures in nonlinear
integrable systems [1, 2]. The ultradiscretization method is summarized as
follows. Difference equations are derived from given continuous nonlinear
equations. By some limiting procedure, the difference equations are converted
into other type of difference equations with max-plus algebra, which is called
ultradiscrete equation. This limiting procedure produces piecewise
linearization of the original equations. Note that the ultradiscretization can
be applied to another type of dynamical systems such as non-integrable non-
equilibrium dissipative systems and reaction-diffusion systems [3, 4, 5, 6, 7,
8, 9, 10].

Recently, we have applied ultradiscretization to bifurcation phenomena in one-
dimensional dynamical systems[10]. Bifurcation phenomena have been hugely
studied from viewpoint of continuous[11, 12, 13] and discrete[14, 15, 16]
dynamical systems. In our study, focusing on the one-dimensional normal forms
of saddle-node, transcritical, and pitchfork bifurcations, their ultradiscrete
equations were derived and the dynamical properties for the obtained
ultradiscrete equations were investigated. In particular, we found that they
possess “ultradiscrete bifurcations”, which has similar properties to the
original bifurcations. These ultradiscrete bifurcations can be visually
understood by piecewise linear graphs of the ultradiscrete equations.

In this letter, we focus on Hopf bifurcation in two-dimensional dynamical
systems. First, we propose a model of ultradiscrete equations which exhibits a
bifurcation similar to Hopf bifurcation. Actually, our proposed model shows a
dynamical transition between a monostable state and a state with limit cycles.
Next, we show that our model can be a normal form of ultradiscrete Hopf
bifurcation by deriving it from two different nonlinear dynamical models.

Let us focus on the following model consisting of max-plus equations with a
bifurcation parameter $B$:

$\displaystyle X_{n+1}$ $\displaystyle=$ $\displaystyle Y_{n}+\max(0,2X_{n}),$
(1) $\displaystyle Y_{n+1}$ $\displaystyle=$ $\displaystyle B-\max(0,2X_{n}).$
(2)

Eqs. (1)-(2) are considered as a discrete dynamical system
$\bm{x}_{n+1}=\bm{F}(\bm{x}_{n})$ for the state variable
$\bm{x}_{n}=(X_{n},Y_{n})$, equipping the evolution operator

$\displaystyle\bm{F}:R^{2}\to R^{2},(x,y)\mapsto(y+\max(0,2x),B-\max(0,2x)).$
(3)

A trajectory
$\\{\bm{x}_{0},\bm{x}_{1},\bm{x}_{2},\dots\\}\left(\equiv\\{\bm{x}_{n}\\}\right)$
from the initial point $\bm{x}_{0}=(X_{0},Y_{0})$ is given by
$\bm{x}_{n}=\underbrace{\bm{F}\circ\cdots\circ\bm{F}}_{n}(\bm{x}_{0})=\bm{F}^{n}(\bm{x}_{0})(n=1,2,\dots)$.
Here we set the following two regions I and II in $(X_{n},Y_{n})$ plane as
shown in Fig.1 (a): (I) $X_{n}>0$ and (II) $X_{n}\leq 0$. In each region, Eqs.
(1)-(2) can be represented as the following matrix form.

(Region I)

When $X_{n}>0$, Eqs. (1)-(2) can be rewritten as

$\displaystyle\left(\begin{array}[]{ccc}X_{n+1}\\\
Y_{n+1}\end{array}\right)=\left(\begin{array}[]{ccc}2&1\\\ -2&0\\\
\end{array}\right)\left(\begin{array}[]{ccc}X_{n}\\\
Y_{n}\end{array}\right)+\left(\begin{array}[]{ccc}0\\\ B\end{array}\right),$
(12)

where Eq. (12) has the fixed point $\bm{\bar{x}}_{I}=(B,-B)$. The matrix
$\bm{A}_{I}=\left(\begin{array}[]{ccc}2&1\\\ -2&0\\\ \end{array}\right)$
satisfies Tr$\bm{A}_{I}=$ det$\bm{A}_{I}=2$, where Tr and det stand for trace
and determinant of a matrix, respectively. Therefore, the trajectory given by
Eq. (12) is characterized as a clockwise spiral source[16], whose center is
the unstable fixed point $\bm{\bar{x}}_{I}$. Fig. 1(b) shows the trajectory of
$\\{\bm{x}_{n}\\}$ with $\bm{x}_{0}=(1,0)$ and $B=0$.

(Region II)

When $X_{n}\leq 0$, the matrix form of Eqs. (1)-(2) is

$\displaystyle\left(\begin{array}[]{ccc}X_{n+1}\\\
Y_{n+1}\end{array}\right)=\left(\begin{array}[]{ccc}0&1\\\ 0&0\\\
\end{array}\right)\left(\begin{array}[]{ccc}X_{n}\\\
Y_{n}\end{array}\right)+\left(\begin{array}[]{ccc}0\\\ B\end{array}\right).$
(21)

Equation (21) has the fixed point $\bm{\bar{x}}_{II}=(B,B)$. For the matrix
$\bm{A}_{II}=\left(\begin{array}[]{ccc}0&1\\\ 0&0\\\ \end{array}\right)$,
Tr$\bm{A}_{II}=$ det$\bm{A}_{II}=0$. The fixed point $\bm{\bar{x}}_{II}$
becomes a stable node. Actually for any $\bm{x}_{0}=(X_{0},Y_{0})$, it is
readily found that $\bm{x}_{1}=(Y_{0},B)$ and
$\bm{x}_{2}=(B,B)=\bm{\bar{x}}_{II}$; Fig. 1 (c) shows an example.

(a) (b) (c)

Figure 1:  (a) Definition of regions I and II. (b) A trajectory in the
vicinity of the unstable focus $\bm{\bar{x}}_{I}$. (c) Trajectories obtained
from Eq. (21).

Taking these dynamical properties in regions I and II into account,
bifurcation for Eqs. (1)-(2), can be grasped as follows. (i) When $B\leq 0$,
both $\bm{\bar{x}}_{I}$ and $\bm{\bar{x}}_{II}$ are in region II;
$\bm{\bar{x}}_{II}$ becomes a fixed point, but $\bm{\bar{x}}_{I}$ does not.
Then, Eqs. (1)-(2) have a unique fixed point $\bm{\bar{x}}_{II}=(B,B)$. (i-a)
When $\bm{x}_{0}=(X_{0},Y_{0})$ belongs to “region II-1”, which means
$X_{0}\leq 0$ and $Y_{0}\leq 0$, $\bm{x}_{1}=(Y_{0},B)$ and
$\bm{x}_{2}=(B,B)=\bm{\bar{x}}_{II}$. Then any initial point in region II-1
reaches $\bm{\bar{x}}_{II}$ with two iteration steps. (i-b) $\bm{x}_{0}$ in
region I moves into region II-1 within four iteration steps as shown in Fig.2.
Then, any initial point in region II-1 reaches $\bm{\bar{x}}_{II}$ within six
steps. (i-c) When $\bm{x}_{0}$ is in “region II-2”, which shows $X_{0}\leq 0$
and $Y_{0}>0$, $\bm{x}_{1}=(Y_{0},B)$ belongs to region I. Therefore,
$\bm{x}_{0}$ reaches $\bm{\bar{x}}_{II}$ within seven steps from the property
(i-b). From (i-a)-(i-c), it is found that the fixed point
$\bm{\bar{x}}_{II}=(B,B)$ is stable. The trajectories of $\\{\bm{x}_{n}\\}$ in
this case are shown in Fig. 3. Figure 3 also shows that excitability occurs
when $\bm{x}_{0}$ is in region II-2. Figure 4 shows time evolutions of
$(X_{n},Y_{n})$ from two different initial conditions. By comparing with
Fig.4(a), it is found that excitability occurs in Fig.4(b).

Figure 2: Flowchart for time evolutions of $\bm{x}_{n}$ from the initial state
$(X_{0},Y_{0})$ in region I. $(X_{0},Y_{0})$ moves into region II-1 within
four iteration steps.

(a) (b)

Figure 3:  Trajectories obtained from Eqs. (1)-(2) with $B\leq 0$. (a)
Schematic explanation. (b) Numerical results from four different initial
states described by filled circles and squares. Trajectories proceed in the
direction of the arrows and finally reach the stable fixed point
$\bm{\bar{x}}_{II}$ shown by the black star.

(a) (b)

Figure 4:  $(X_{n},Y_{n})$ as a function of $n$ for Eqs. (1)-(2) with $B=-1$
from two different initial conditions.

(ii) When $B>0$, $\bm{\bar{x}}_{I}=(B,-B)$ becomes a unique unstable fixed
point. We find that there exist only two different clockwise periodic
solutions, $\mathcal{C}$ and $\mathcal{C}_{s}$ as shown in Fig.5(a), which are
composed of the following seven points, respectively:
$\mathcal{C}=\\{(B,B)\rightarrow(3B,-B)\rightarrow(5B,-5B)\rightarrow(5B,-9B)\rightarrow(B,-9B)\rightarrow(-7B,-B)\rightarrow(-B,B)\left[\rightarrow(B,B)\right]\\}$,
$\mathcal{C}_{s}=\\{\left(\frac{B}{15},B\right)\rightarrow\left(\frac{17B}{15},\frac{13B}{15}\right)\rightarrow\left(\frac{47B}{15},-\frac{19B}{15}\right)\rightarrow\left(5B,-\frac{79B}{15}\right)\rightarrow\left(\frac{71B}{15},-9B\right)\rightarrow\left(\frac{7B}{15},-\frac{127B}{15}\right)\rightarrow\left(-\frac{113B}{15},\frac{1B}{15}\right)\left[\rightarrow\left(\frac{B}{15},B\right)\right]\\}$.
It is also found that a trajectory with any initial condition except for
$\bm{\bar{x}}_{I}$ is finally absorbed in $\mathcal{C}$ or $\mathcal{C}_{s}$.
Therefore, $\mathcal{C}$ and $\mathcal{C}_{s}$ are considered as limit cycles.
Figure 5(b) shows trajectories from four different initial conditions; they
finally converge into $\mathcal{C}$.

(a) (b)

Figure 5:  (a) The two limit cycles $\mathcal{C}$ (open circles) and
$\mathcal{C}_{s}$ (filled circles). (b) Examples of trajectories starting from
four different filled squares. The trajectories finally converge into
$\mathcal{C}$. The star in each figure shows $\bm{\bar{x}}_{I}$.

Here we discuss basins for $\mathcal{C}$ and $\mathcal{C}_{s}$. It is
confirmed that any trajectory has a point with $Y_{n}=B$ at a certain
iteration step $n$. In other words, a point on the line $l_{B}$,
$\bm{x}_{n}=(X_{n},B)$ where $X_{n}\in(-\infty,\infty)$, exists in all
trajectories. It is found that all trajectories with a point on $l_{B}$ except
for $\left(\frac{B}{15},B\right)$ and $\left(\frac{17B}{15},B\right)$ are
finally absorbed into $\mathcal{C}$. And any trajectory with
$\left(\frac{B}{15},B\right)$ or $\left(\frac{17B}{15},B\right)$ finally goes
into $\mathcal{C}_{s}$. Therefore, there exist only two limit cycles
$\mathcal{C}$ and $\mathcal{C}_{s}$ in this model.

Now we show that Eqs. (1)-(2) become a candidate of a normal form of
ultradiscrete Hopf bifurcation. Through ultradiscretization, Eqs. (1)-(2) can
be derived from the following two different dynamical systems described by
partial differential equations.
(i) Sel′kov model[12, 17]

$\displaystyle{}\displaystyle\frac{dx}{dt}$ $\displaystyle=$
$\displaystyle-x+ay+x^{2}y,$ (22) $\displaystyle\displaystyle\frac{dy}{dt}$
$\displaystyle=$ $\displaystyle b-ay-x^{2}y,$ (23)

where $a$ and $b$ are positive. They can be considered as bifurcation
parameters for Hopf bifurcation. Actually, eqs. (22)-(23) exhibit Hopf
bifurcation when $a$ and $b$ satisfy $b^{2}=\frac{1}{2}(1-2a\pm\sqrt{1-8a})$.
Ultradiscretiazation of Sel′kov model can be performed in the following way.
By tropical discretization[4], the following difference equations are adopted
for eqs. (22)-(23),

$\displaystyle{}x_{n+1}=\frac{x_{n}+\Delta t(ay_{n}+x^{2}_{n}y_{n})}{1+\Delta
t},$ (24) $\displaystyle y_{n+1}=\frac{y_{n}+\Delta tb}{1+\Delta
t(a+x_{n}^{2})},$ (25)

where $\Delta t$ is the discretized time interval. $x_{n}=x(n\Delta t)$,
$y_{n}=y(n\Delta t)$, where $n$ is positive integer. The variable
transformations,

$\displaystyle\begin{cases}\Delta
t=e^{T/\varepsilon},\;\;\;x_{n}=e^{X_{n}/\varepsilon},\;\;\;y_{n}=e^{Y_{n}/\varepsilon},\\\
a=e^{A/\varepsilon},\;\;\;b=e^{B/\varepsilon},\;\;\;\end{cases}$ (26)

are applied to eqs.(24)-(25), and the ultradiscrete limits

$\displaystyle\begin{cases}\displaystyle\lim_{\varepsilon\to+0}\varepsilon\log(e^{A/\varepsilon}+e^{B/\varepsilon}+\cdot\cdot\cdot)~{}=~{}\max(A,B,\dots),\\\
\displaystyle\lim_{\varepsilon\to+0}\varepsilon\log(e^{A/\varepsilon}\cdot
e^{B/\varepsilon}\cdot\dots)=A+B+\dots.\end{cases}$ (27)

are performed. Then, the ultradiscrete equations for Sel′kov model are
obtained as

$\displaystyle X_{n+1}$ $\displaystyle=$
$\displaystyle\max(X_{n},T+\max(A+Y_{n},2X_{n}+Y_{n}))-\max(0,T),$ (28)
$\displaystyle Y_{n+1}$ $\displaystyle=$
$\displaystyle\max(Y_{n},T+B)-\max(0,T+\max(A,2X_{n})).$ (29)

Assuming that $T\geq\max\\{0,-A,Y_{n}-B,-(X_{n}+Y_{n})\\}$ for all $n$ and
$A=0$, eqs. (28)-(29) are identical to eqs.(1)-(2).
(ii) Lengyel model [18, 19],

$\displaystyle{}\displaystyle\frac{dx}{dt}$ $\displaystyle=$ $\displaystyle
r-x-\frac{4xy}{1+x^{2}},$ (30) $\displaystyle\displaystyle\frac{dy}{dt}$
$\displaystyle=$ $\displaystyle x(1-\frac{sy}{1+x^{2}}),$ (31)

where $r$ and $s$ are positive bifurcation parameters for Hopf bifurcation.
Ultradiscrete equations of Eqs. (30)-(31) are obtained in the similar way to
the case of Eqs. (22)-(23). By tropical discretization, the following
difference equations are derived from Eqs. (30)-(31),

$\displaystyle{}x_{n+1}=\frac{x_{n}+\Delta t(r-x_{n})}{1+\Delta
t(\frac{4y_{n}}{1+x^{2}_{n}})},$ (32) $\displaystyle
y_{n+1}=\frac{y_{n}+\Delta tx_{n}}{1+\Delta t(\frac{sx_{n}}{1+x_{n}^{2}})}.$
(33)

The variable transformations,

$\displaystyle\begin{cases}\Delta t=e^{T/\varepsilon},\;\;\;1-\Delta
t=e^{M/\varepsilon},\;\;\;x_{n}=e^{X_{n}/\varepsilon},\;\;\;y_{n}=e^{Y_{n}/\varepsilon},\\\
r=e^{R/\varepsilon},\;\;\;s=e^{S/\varepsilon},\;\;\;\end{cases}$ (34)

and ultradiscrete limits (27) are carried out. Then, the ultradiscrete
equations for Lengyel model are obtained as

$\displaystyle X_{n+1}$ $\displaystyle=$ $\displaystyle
R-Y_{n}+\max(0,2X_{n}),$ (35) $\displaystyle Y_{n+1}$ $\displaystyle=$
$\displaystyle-S+\max(0,2X_{n}),$ (36)

on the assumption $T\geq\max\\{0,X_{n}+M-R,Y_{n}-X_{n},|X_{n}|-S\\}$ for all
$n$. It is noted that Eqs.(35)-(36) have counterclockwise limit cycles
reflected on the signs in the equations. Therefore taking the variable
transformations $Y_{n}\to-Y_{n}$ and setting $R=0,S=B$ in Eqs.(35)-(36),
Eqs.(1)-(2) can be reproduced.

In conclusion, we have proposed the simple model for ultradiscrete Hopf
bifurcation. Depending on the value of the bifurcation parameter, the model
exhibits excitability and possesses the limit cycles. The model is derived
from the two different nonlinear dynamical models, Sel′kov model and Lengyel
model; it can be a normal form of ultradiscrete Hopf bifurcation. Further
investigation for basins of $\mathcal{C}$ and $\mathcal{C}_{s}$ and
characterization for their periodicity will be future problems.

Acknowledgement

The authors are grateful to Prof. M. Murata, at Tokyo University of
Agriculture and Technology, Prof. K. Matsuya at Musashino University, Prof. D.
Takahashi, Prof. T. Yamamoto, and Prof. Emeritus A. Kitada at Waseda
University for useful comments and encouragements. This work was supported by
Sumitomo Foundation, Grant Number 200146.

## References

  * [1] T. Tokihiro, D. Takahashi, J. Matsukidaira, and J. Satsuma, Phys. Rev. Lett. 76, 3247 (1996).
  * [2] B. Grammaticos, Y. Ohta, A. Ramani, D. Takahashi, and K. M. Tamizhmani, Phys. Lett. A 226, 53 (1997).
  * [3] T. Nagatani, Phys. Rev. E 58 700 (1998).
  * [4] M. Murata, J. Differ. Equations Appl. 19 1008 (2013).
  * [5] S. Ohmori and Y. Yamazaki, Prog. Theor. Exp. Phys. 08A01 (2014).
  * [6] K. Matsuya and M. Murata, Discrete Contin. Dyn. Syst. B 20 173 (2015).
  * [7] M. Murata, J. Phys. A Math, Theor. 48 255202 (2015).
  * [8] S. Gibo and H. Ito, J. Theor. Biol. 378 89 (2015).
  * [9] S. Ohmori and Y. Yamazaki, J. Phys. Soc. Jpn. 85 045001 (2016).
  * [10] S. Ohmori and Y. Yamazaki, J. Math. Phys. 61 122702 (2020)
  * [11] J. Guckenheimer and P. Holmes, _Nonlinear Oscillations, Dynamical Systems, and Bifurcations of Vector Fields_ (Springer, New York, 1983).
  * [12] Steven. H. Strogatz, _Nonlinear Dynamics and Chaos_ (Westview Press, U.S. 1994).
  * [13] G. Nicolis, _Introduction to Nonlinear Science_ (Cambridge Univ. Press 1995).
  * [14] C. Robinson, _Dynamical systems -Stability, Symbolic Dynamics, and Chaos-, 2ed edition_ (CRC Press, Florida 1999).
  * [15] Yuri A. Kuznetsov, _Elements of Applied Bifurcation Theory_ (Springer-Verlag, New York, 2010).
  * [16] O. Galor, _Discrete Dynamical Systems_ (Springer, New York 2010).
  * [17] E. E. Sel′kov, Eur. J. Biochem. 4 79 (1968).
  * [18] I. Lengyel, G. Rabai, and I. R. Epstein, J. Am. Chem. Soc, 112 9104 (1990).
  * [19] I. Lengyel and I. R. Epstein, Science, 251 650 (1991).