# Combining Sampling- and Gradient-based Planning for Contact-rich
Manipulation

Filippo Rozzi1, Loris Roveda2, and Kevin Haninger3 1Corresponding author,
Politecnico di Milano, Department of Mechanical Engineering, Milano, Italy
<EMAIL_ADDRESS>Dalle Molle di Studi sull’Intelligenza
Artificiale (IDSIA), Scuola Universitaria Professionale della Svizzera
Italiana (SUPSI), Università della Svizzera Italiana (USI) IDSIA-SUPSI,
Lugano, Switzerland<EMAIL_ADDRESS>of Automation at
Fraunhofer IPK, Berlin, Germany<EMAIL_ADDRESS>project
has received funding from the European Union’s Horizon 2020 research and
innovation programme under grant agreement 101058521 — CONVERGING.

###### Abstract

Planning over discontinuous dynamics is needed for robotics tasks like
contact-rich manipulation, which presents challenges in the numerical
stability and speed of planning methods when either neural network or
analytical models are used. On the one hand, sampling-based planners require
higher sample complexity in high-dimensional problems and cannot describe
safety constraints such as force limits. On the other hand, gradient-based
solvers can suffer from local optima and convergence issues when the Hessian
is poorly conditioned. We propose a planning method with both sampling- and
gradient-based elements, using the Cross-entropy Method to initialize a
gradient-based solver, providing better search over local minima and the
ability to handle explicit constraints. We show the approach allows smooth,
stable contact-rich planning for an impedance-controlled robot making contact
with a stiff environment, benchmarking against gradient-only MPC and CEM.

##  1 Introduction

Contact-rich manipulation tasks such as door opening and gear mating involve
significant environmental dynamics and constraints, which discretely switch
during the task as contact conditions change. Human manipulation can
effectively handle these contact changes, realizing safe and robust behavior
even with uncertain geometry and dynamics. Robotic manipulation methods which
can maintain safety (e.g. limits on contact force and actuation) and
robustness (variation in contact geometry or dynamics) are needed which are
suited to discontinuous dynamics.

Gradient-based optimization methods are popular for contact planning, with a
range of formulations from contact implicit [1, 2, 3, 4] to hybrid contact
modes [5, 6]. Sampling-based methods have been used for higher-level contact
mode planning [7, 8] and continuous planning [9].

The tradeoffs between sampling and gradient-based planners has led to recent
interest in combining them, such as adding gradient steps to the sampling-
based Cross-Entropy Method (CEM) [10, 11]. In contact, recent work has
compared the variance of gradient-based and gradient-free policy learning in
contact [12]. This raises the question if a planning method can effectively
use both parallelized sampling and gradients, where sampling may address the
initialization sensitivity of gradient-based methods [13], while gradient-
based methods can provide explicit safety constraints. This paper proposes a
planning method which leverages CEM and interior-point optimization for
contact-rich planning tasks.

###  1.1 Related works

Model-based planning for contact-rich robotic tasks consists of finding the
action sequence that minimizes a cumulative cost function, according to
learned, parametrized, or simulated dynamics. Usually, optimized actions are
executed with model predictive control (MPC), replanning at each time step.
The problem can be formulated using complimentarity conditions for contact
constraints [14, 1], where sequential quadratic [2] or interior-point [3]
solvers can be used. This contact-implicit trajectory optimization allows
planning of contact mode and trajectory in one problem with multiple contacts
and Coulomb constraints [4, 3, 13]. Alternatively, hybrid dynamics can be
written with dynamics for each contact mode [5, 6].

The contact planning problem is challenging, typically needing relaxations [9,
3], linearization [4, 13], or quasi-static assumptions [15] to achieve
acceptable speeds. The resulting problems can be sensitive to initialization
[13], and in general gradient-based methods can suffer from local minima or
poor convergence when the Hessian of the Lagrangian is poorly conditioned
[16].

On the other hand, population-based sampling planners can parallelize well and
do not need gradients, such as Random-Sampling Shooting [17], Model Predictive
Path Integral Control [18], and the Cross-Entropy Method (CEM) [19]. Sampling-
based have been used in contact planning, for higher-level mode initialization
of a contact-implicit [7] and reachability planner [8]. Sampling can also be
used for the continuous planning, where Rapidly-exploring random trees have
been applied to relaxed contact models [9].

These zeroth-order optimizers maintain a sampling distribution action
sequences [20], which is sampled from then updated at each iteration,
assigning a higher probability to higher reward action sequences. In recent
years, there has been a huge research effort on sampling based optimization in
the context of model-based reinforcement learning. By adding temporally
correlated actions and memory, CEM can be improved in both sample-efficiency
and reward performance [21]. In [22], authors proposed a model-based
reinforcement learning approach in which the interaction dynamics was directly
learned through an ensemble of Artificial Neural Networks (ANNs) and CEM was
used for online planning.

Figure 1: Experimental setup showing a robot making contact with a table,
moving towards a hinge contact condition, providing a task with three modes:
free space, vertical stiff contact, and hinge contact. The soft contacts are
used for validating the effect of contact stiffness on performance.

However, CEM, like all sampling-based planners, suffers from scalability as
the dimension of the action space or planning horizon grows, requiring more
samples for convergence and having increasing reward variance [10, 23].
Additionally, sampling-based methods cannot explicitly handle state
constraints, raising questions about the ability to handle safety in contact.

An interesting idea is to combine zeroth-order and first-order planning
methods in order to enjoy the benefits of both simultaneously. To this end,
authors in [10, 11] proposed a combined planning framework that performs
gradient updates on the sampled trajectories. Although able to exploit the
advantages of both approaches, these methods still show scalability issues
with very high dimensional planning problems, and don’t handle constraints.
For policy learning in contact, recent work [24] has shown that first order
methods can have higher variance than zeroth order methods at higher stiffness
and in longer planning horizons.

Figure 2: Block diagram of the proposed approach

###  1.2 Paper contribution

To address poor convergence and local minima problems of gradient-based
optimizers when dynamics are discontinuous, we propose a hybrid planner that
uses both sampling- and gradient-based MPC. Contact is modelled with hybrid
dynamics, where dynamics are switched according to contact mode. A particle
filter is used to estimate the hybrid state online, providing a belief of
current mode for the planner. CEM is used to initialize the action and state
trajectory, modified to normalize sample cost by the weight and minimizing
expected cost over the belief in mode. A gradient-based optimizer is then
used, enforcing explicit constraints on the robot trajectory, as seen in
Figure 2.

Compared to the work of [21], this CEM rolls out trajectories on a hybrid
system, leveraging the contact mode belief provided by the filter to weight
the computed samples costs. Compared to other combined sampling- and gradient-
based planning approaches [10, 11], this method can enforce explicit
constraints on the robot trajectory and handles hybrid dynamics. Our method is
also validated experimentally, showing the feasibility of the approach with a
real high-stiffness contact task.

The proposed approach is applied to contact-rich manipulation as seen in
Figure 1, where the impedance controller reference trajectory is planned. The
ability to estimate of contact mode online, as well as safely respond with
changes in control is validated. This paper first introduces the dynamic
models used, estimation of contact mode. Then, planning framework is
introduced and studied in simulation. Finally, experiments are used to
validate the approach, comparing the proposed planner to only MPC and only
CEM. We find that the CEM initialization helps with discontinuous dynamics,
reducing the mean and variance of the total solve time, even from just one
iteration, as compared to a pure gradient-based MPC for contact.

##  2 Problem Statement and Models

This section introduces the dynamic models used for contact and the hybrid
dynamic model.

###  2.1 Stiffness Contact Model

We model environment contact as parallel 1-DOF stiffnesses to compromise
between model identification efficiency, ease of deployment, generalizability,
and differentiability [25]. The $i^{\mathrm{th}}$ stiffness element exerts
contact force $F_{i}$ at the point of contact of

$F_{i}=K_{i}^{T}\left(x_{i}^{o}-(Rx_{i}+x)\right),$ (1)

where stiffness $K_{i}\in\mathbb{R}^{3}$ has a rest pose at
$x_{i}^{o}\in\mathbb{R}^{3}$ and $x_{i}\in\mathbb{R}^{3}$ is the contacting
point in TCP frame, transformed into world coordinates with $R$ and $x$ , the
orientation and position of the end-effector. The contact normal and stiffness
are jointly described by the vector $K_{i}$, i.e., the contact normal is
$n_{i}=K_{i}/\|K_{i}\|$, and the stiffness is $\|K_{i}\|$. For a single point
stiffness, the joint torque induced is $\tau_{i}=J_{i}^{T}F_{i}$, where
$J_{i}=\partial(R(q)x_{i}+x(q))/\partial q$.

When $N_{c}$ contacts are active, the total external joint torque is
$\tau_{e}=\sum_{i=1}^{N_{c}}\tau_{i}(q)$. The parameters $K_{i}$, $x_{i}^{o}$,
and $x_{i}$ can be identified from least-squares fit on a dataset as
$\min_{K,x^{o},x}\sum_{t}\|\tau_{t}^{m}-\tau_{e}(q_{t},K,x^{o},x)\|$. These
contact forces are composed of differentiable operations, and easily
implemented in automatic differentiation (AD) frameworks to support their use
in gradient-based planning.

###  2.2 Discontinuous Dynamics

The standard serial manipulator robot dynamics

$M(q)\ddot{q}+C(q,\dot{q})+B\dot{q}+G(q)=\tau_{m}+\tau_{e}$ (2)

is assumed to be available in an AD framework, with joint position
$q\in\mathbb{R}^{n}$, inertia matrix $M(q)$, Coriolis terms $C(q,\dot{q})$,
viscous damping $B$, gravitational torque $G(q)$, input torque $\tau_{m}$, and
torques from contact $\tau_{e}$.

We discriminate between different contact modes by letting the torque
$\tau_{e}$ depend on the mode $z\in[1,\dots,N_{z}]$, where each mode has
different point contact models, resulting in a different $\tau_{e}$. We
suppose the mode changes during the task, denoting the value at $t$ as
$z_{t}$.

As typical in contact models, we use semi-implicit integration which handles
stiff differential equations better [14]. For a time step $\Delta_{t}$:

$\displaystyle z_{t}$ $\displaystyle\sim p(z_{t}|z_{t-1})$ (3)
$\displaystyle\xi_{t}=\begin{bmatrix}q_{t}\\\ \dot{q}_{t}\end{bmatrix}$
$\displaystyle=\begin{bmatrix}q_{t-1}+\Delta_{t}\dot{q}_{t}\\\
\dot{q}_{t-1}+\Delta_{t}M^{-1}(\tau_{m}+\tau_{e}(z_{t})-B\dot{q}_{t})\end{bmatrix}$
(4)

where $p(z_{t}|z_{t-1})$ is a diagonal-dominant discrete transition matrix. As
the external contact force $\tau_{e}(z_{t})$ depends on the current discrete
mode $z_{t}$, the dynamics can be written as:

$\xi_{t}=\begin{bmatrix}q_{t}\\\
\dot{q}_{t}\end{bmatrix}=f\left(\begin{bmatrix}q_{t-1}\\\
\dot{q}_{t-1}\end{bmatrix},z_{t-1},\tau_{m,t-1}\right).$ (5)

For an EKF, the dynamics in (5) would be linearized with respect to the state,
adding a process noise contribution of $w\sim\mathcal{N}(0,Q)$. The
observation model would comprise the external joint torques $\tau_{e}$ plus
the joint positions as

$\displaystyle y^{z}$ $\displaystyle=C^{z}\xi+v,$ (6)

where $C^{z}=[C_{q}^{T},\frac{\partial\tau_{e}(z)}{\partial\xi}^{T}]^{T}$,
$C_{q}=[I_{n},0_{n}]$, $v\sim\mathcal{N}(0,R)$ is the measurement noise,
assumed to be independent and identically distributed,
$\partial\tau_{e}(z)/\partial\xi$ is the Jacobian of the external contact
torques $\tau_{e}(z)$ with respect to the joint states $\xi$, and $I_{n}$ an
identity matrix of dimension $n$.

###  2.3 Robot impedance dynamics

The robot is controlled with a Cartesian space impedance controller, where the
$3$ DoF impedance controlled robot dynamics are

$\tau_{m}=J^{T}(q)\left[D\dot{x}+K\left(x-x_{0}\right)\right]+G(q)+C(q,\dot{q})$
(7)

where $x,\,\dot{x}\in\mathbb{R}^{3}$ are respectively the TCP Cartesian
position and velocity, $J=\partial x/\partial{q}$ and the impedance parameters
are $x_{0}\in\mathbb{R}^{3}$ which represents the Cartesian impedance rest
position, the damping $D\in\mathbb{R}^{3\times 3}$, and stiffness
$K\in\mathbb{R}^{3\times 3}$ matrices, which are here all diagonal.

The dynamics in Eq.7 can then be updated to provide controlled discretized
dynamics of

$\begin{bmatrix}q_{t}\\\
\dot{q}_{t}\end{bmatrix}=f\left(\begin{bmatrix}q_{t-1}\\\
\dot{q}_{t-1}\end{bmatrix},x_{0,t-1},z_{t-1}\right).$ (8)

##  3 Online Planning Method

This section introduces the belief estimation, cross-entropy and gradient-
based planning method used.

###  3.1 Belief estimation

For estimating the hybrid state, we adopt a hybrid particle filter algorithm
[26]. The proposed algorithm maintains an overall belief of the discrete mode
$\hat{p}(z_{t}|y_{1:t})$, from the weighted belief of individual particles
$\hat{p}(z^{(i)}|y_{1:t-1})$. In each step, each particle samples a discrete
mode $\hat{z}_{t}^{(i)}\sim\hat{p}(z_{t}^{(i)}|y_{1:t-1})$, then propagates
the mean $\mu^{(i)}_{t}$ and the covariance $\Sigma^{(i)}_{t}$ of continuous
state $\xi_{t}$ with the corresponding mode dynamics and observation equations
via a standard Extended Kalman filter step [26]. Each particle $p^{(i)}$ is
characterized by the tuple
$[\mu_{t}^{(i)},\Sigma_{t}^{(i)},\hat{p}(z^{(i)}_{t}|y_{1:t-1}),\hat{z}_{t}^{(i)}]$.
An overview can be seen in Algorithm 1, where $\hat{z}^{(i)}_{t}$ is the
sampled mode for particle $i$,
$\hat{S}^{(i)}_{t}=\mathrm{cov}(y_{t}|y_{1:t-1})$ the predicted observation
covariance, $w^{(i)}_{t}$ the weight,
$\hat{y}^{(i)}_{t|t-1}=\mathbb{E}(y_{t}|y_{1:t-1})$ is the predicted
measurement, and $\delta$ is the Dirac delta function.

Algorithm 1 Hybrid filter

for all particles $p^{(i)}$ do

$\hat{\mu}^{(i)}_{1|0}\leftarrow\mu_{0}$

$\hat{\Sigma}^{(i)}_{1|0}\leftarrow\Sigma_{0}$

$\hat{p}(z^{(i)}_{1}|y_{1})\leftarrow p(z_{0})$

end for

for all time step $t$ do

for all particles $p^{(i)}$ do

$\hat{p}(z_{t}^{(i)}|,y_{1:t-1})=\hat{p}(z^{(i)}_{t-1}|y_{1:t-1})p(z_{t}|z_{t-1})$

$\hat{z_{t}}^{(i)}\sim\hat{p}(z_{t}^{(i)}|,y_{1:t-1})$

$(\hat{y}^{(i)}_{t|t-1},\hat{S}^{(i)}_{t},\hat{\mu}^{(i)}_{t},\hat{\Sigma}^{(i)}_{t})$

$\leftarrow EKF(\mu^{(i)}_{t-1},\Sigma^{(i)}_{t-1},y_{t},\hat{z_{t}}^{(i)})$

$w^{(i)}_{t}\propto\mathcal{N}(y_{t};\hat{y}^{(i)}_{t|t-1},\hat{S}^{(i)}_{t})$

end for

Resample particles $\\{p^{(i)},w_{t}^{(i)}\\}_{i=1}^{N}$

Estimate marginal posterior belief as

${\hat{p}(z_{t}|y_{1:t})\propto\sum_{i=1}^{N}w^{(i)}_{t}\delta_{z_{t}}(\hat{z}^{(i)}_{t})}$

end for

The posterior belief $\hat{p}(z_{t}|y_{1:t})$ is used to detect changes in
contact, and is used in the planning methods to weight the costs associated
with different modes.

###  3.2 Proposed Cross-Entropy Method

The Cross-Entropy Method (CEM) samples actions from an action distribution,
rolls out these trajectories and evaluates their resulting costs, adapting the
action distribution to reduce loss. Sampling-based approaches may avoid local
minima problems or convergence issues when gradients are not smooth due to the
discontinuous dynamics, potentially reaching state-action pairs that would not
be reached by gradient-based methods. Thus, we use CEM iterations before
proceeding with a gradient-based non-linear solver.

The adopted CEM algorithm is based on the work of [21]. The CEM here is
initialized from the hybrid filter at each time step $t$ with the set of
particles consisting of continuous joint states values $\xi_{t}^{(i)}$, their
weight $w_{t}^{(i)}$ and discrete sampled state $z^{(i)}_{t}$,

$\left[\xi^{(i)}_{t},w_{t}^{(i)},z^{(i)}_{t}\right]_{i=1}^{N}.$ (9)

This list is used for calling the dynamics with the hybrid mode $z_{t}^{(i)}$
with initial continuous state $\xi_{t}^{(i)}$ along with the sampled actions
$u_{t:t+h}$. As in typical CEM, actions are iteratively sampled and the
sampling distribution shifted towards the elite samples which have the best
loss. A detailed explanation of the steps is given in Algorithm 2.

Algorithm 2 Proposed iCEM

$N$: number of particles, $h$: planning horizon, $d$: action dimension,
$\beta$: colored noise exponent, $N_{iter}$: CEM-iterations,
$\left[\xi^{(i)}_{t},w_{t}^{(i)},z^{(i)}_{t}\right]_{i=1}^{N}$: list of
particles, $K$: elite size

Initialize:

$u_{t:t+h}=u_{t-1:t+h-1}^{MPC}$

$\Sigma^{u}=\textbf{1}$

Loop:

for $k$ in $N_{iter}$ do

noise samples $\leftarrow\ N$ samples from $C^{\beta}(d,h)$

action samples $u^{(i)}_{t:t+h}$

$\leftarrow\ N\ $ samples from
$\mathtt{clip}(u_{t}+C^{\beta}(d,h)\cdot\Sigma^{u})$

normalize weights
$w^{(i)}_{t}\leftarrow\frac{w^{(i)}_{t}}{\sum_{i=1}^{N}w^{(i)}_{t}}$

trajectory rollouts $\xi^{(i)}_{t+1:t+1+h}$

$\leftarrow f(\xi^{(i)}_{t},u^{(i)}_{t:t+h},z^{(i)}_{t})$

costs $l^{(i)}\leftarrow\ $

cost of each $\\{\xi^{(i)}_{t+1:t+1+h},u^{(i)}_{t:t+h}\\}_{i=1}^{N}$

normalize costs $l^{(i)}\leftarrow\frac{l^{(i)}}{w^{(i)}_{t}}$

elite set $\leftarrow$ best $K$ samples according to cost

$u_{t:t+h},\Sigma^{u}\leftarrow$ fit Gaussian distribution to elite set

end for

return best action trajectory $u_{t:t+h}^{iCEM}$

In Algorithm 2, $u_{t:t+h}$ and $\Sigma^{u}$ represent the mean and the
standard deviation of the action sequences, at the beginning of each time step
$t$ we initialise the mean with the optimal value found at the previous time
step $t-1$ by the non linear MPC solver while the standard deviation of the
sequences is re-initialised every time with the matrix
$\textbf{1}\in\mathbb{R}^{d\times h}$, whose elements are ones.
$C^{\beta}(d,h)$ represents the colored noise distribution, built following
the approach in [21]. Within the CEM inner loop, we keep updating the action
sequences statistics by fitting at each iteration the best $K$ trajectories;
as long as the solver converges the mean $u_{t:t+h}$ is initialized at each
time step with the optimal solution and the likelihood of sampling optimal
solutions increases. The cost of each particle $l^{(i)}$ is divided by the
particle weight, so a low likelihood particle results in a higher cost, to
avoid actions which favor low-likelihood particles from being sampled more
heavily. The cost function itself is the same as described in the following
section for the MPC problem.

The CEM returns both the best action trajectory, along with the associated
state trajectory, which are used to initialize the decision variables in a
gradient-based MPC problem.

###  3.3 MPC problem

The MPC problem is formulated using multiple-shooting transcription with a
generic problem statement of

$\displaystyle\min_{u_{t:t+h}}$
$\displaystyle\sum_{t}^{t+h}\sum_{z=1}^{N_{z}}\hat{p}(z|y_{1:t})l^{z}(\xi^{z}_{t},u_{t})$
(10) $\displaystyle\mathrm{s.t.}\,\,\,$ $\displaystyle\forall
z\in[1,\dots,N_{z}],\quad\tau\in[t,\dots,t+h-1]:$ (11)
$\displaystyle\|\xi^{z}_{\tau+1}-f(\xi^{z}_{\tau},u_{\tau})\|\leq\rho$ (12)
$\displaystyle g(\xi_{t}^{z},u_{t})\geq 0$ (13) $\displaystyle
u\in\mathcal{U}$ (14)

where $h$ is the planning horizon, $\mathcal{U}$ is the range of allowed
inputs, $\rho$ is the slack for continuity constraints (the inequality is
applied element-wise), $z\in[1,\dots,N_{z}]$ represents the set of the
different contact modes, and $g$ represent general state and input dependent
inequality constraints. A trajectory is rolled out for each mode,
$\xi^{z}_{t}$, where continuity is imposed according to that mode’s dynamics
in $f$; $\hat{p}(z|y_{1:t})$ is the contact state belief information coming
from the hybrid filter at time step $t$. The constraints shown in the
optimization problem are nonlinear, so an interior-point nonlinear
optimization solver is used. While nonlinear, the problem is written in an
automatic differentiation frame- work, allowing calculation of the gradient
and Hessian of the objective and constraints, significantly improving
convergence.

The general stage cost function, which would be used also in the CEM loop, is
defined as

$\begin{split}l^{z}(\xi^{z},u)=&(x_{d}-x)^{T}Q_{x}(x_{d}-x)\\\
&+(u-x)^{T}Q_{u}(u-x)+\dot{x}^{T}Q_{\dot{x}}\dot{x}\end{split}$ (15)

where $x_{d}$ is a desired position in world frame to be tracked, $x$ and
$\dot{x}$ are the Cartesian position and velocity of the robot TCP expressed
with the forward kinematics function, directly depending on the joint
configuration. In this scenario, the state $\xi^{z}=(q^{z},\dot{q}^{z})$ is
represented by the robot joint variables for each mode $z$, while the adopted
control input $u$ is shared among the different modes and is impedance rest
position $x_{0}$. The constraint is $g=\overline{F}_{imp}-\|K(x-x_{0})\|_{2}$,
limiting the virtual impedance force applied to $\overline{F}_{imp}$.

##  4 Simulation Study

To study the proposed approach, studies in simulation are used to investigate
the solve time and resulting performance. To see the variance in the
performance, each problem is solved repeatedly with the initial state
$\xi_{t}$ perturbed with $\mathcal{N}(0,0.1)$ between solves, finding the mean
and covariance.

###  4.1 CEM vs MPC in stiff contact

To compare qualitative properties of sampling- and gradient-based solvers, a
simple tracking problem in stiff contact $K=1e4$ N/m is set up with a cost
function to track desired position. As can be seen in Figure 3, the CEM (solid
line) has more variation in the control trajectory, and the state trajectory
is slower to converge than the MPC (dash-dotted line).

Figure 3: Comparing state and control trajectory from CEM vs MPC for stiff
dynamics. Figure 4: Comparing the resulting cost for a CEM- and MPC-only solve
with stiff dynamics. While the CEM approaches the MPC performance at a high
number of iterations, it does not fully converge.

The mean and variance as $N_{iter}$ is varied can be seen in Figure 4, where
the CEM performance saturates before approaching that of the MPC solver. The
variance of the CEM performance decreases at higher iterations, but remains
orders of magnitude larger than the variance of the MPC solver ($1e-9$).

###  4.2 CEM Warmstart on total solve time

To study the relative performance of the gradient- and sampling-based methods
in discontinuous dynamics, we examine two performance characteristics: time to
solve and resulting expected cost.

First, the solver is used to solve a free-space problem according to (10) with
horizon $h=20$ steps and stepsize $\Delta_{t}=0.05$ seconds. Then, the
dynamics is changed to a single point contact model of stiffness $K=5e3$ N/m
and the previous solution used initialize the current step. This
initialization is then used for $N_{iter}$ iterations of the CEM before
starting the next MPC solve.

The mean and covariance of the total solve time for the CEM iterations and MPC
solve can be seen in Figure 5. As can be seen, the CEM helps to reduce both
the mean and variance of the solve time, even from just one iteration.
However, as the number of CEM iterations increases, the time required
increases, indicating the additional time required to do more CEM iterations
is not saving as much time in the MPC solver.

Figure 5: Effect of the number of CEM warmstart iterations on total solve time
(CEM+MPC)

To analyze the impact of important process parameters, we study the
improvement by $N_{iter}=5$ warmstarts over a range of stiffnesses and
planning horizons. The results can be seen in Figure 6, where as the stiffness
and horizon increase, the CEM improves both the mean and variance of solve
time. This trend matches the impact of stiffness and planning horizon in
policy gradient estimation [12], where gradient-free methods performed better
at long horizons and higher stiffnesses.

(a) Effect of stiffness

(b) Effect of planning horizon

Figure 6: Effect of problem parameters (stiffness and planning horizon) on
solve time for only MPC (orange) and proposed approach (blue).

##  5 Experiments

To validate the performance on real-world contact tasks, we apply the proposed
method to the two tasks seen in Figure 1. The code and data are available at
https://gitlab.cc-asp.fraunhofer.de/hanikevi/contact_mpc.

###  5.1 Vertical contact task

In the vertical contact task, the robot moves from free space into contact
with an environment. The stiffness of the environment is adjusted by adding
the yellow or purple foam seen in Figure 1, or direct contact with the
aluminum table. In Figure 7 we can see how the stiffness parameter influences
the contact mode detection responsiveness. Despite the developed CEM+MPC
approach successfully worked for all the contact environments, from the plot
we can see that with stiffer contacts the mode detection is more sharp, with
the transition triggered right in correspondence of the vertical force peak.
For lower stiffness contacts, mode detection is slightly delayed with respect
to the force peak, probably because the lower measured interaction force makes
it more difficult for the filter to assign different likelihood values.

Figure 7: Effect of stiffness on contact detection

###  5.2 Pivot task

In this task, the robot makes contact with the table, slides along the surface
until making contact with the wall. The MPC problem is solved with $h=13$,
$\Delta_{t}=0.04$, an impedance cost of $Q_{u}=0.05$, maximum impedance force
of $\overline{F}_{imp}=25$, a desired pose of $x^{d}=[-,-,0.01]$ in free
space, and $x^{d}=[0.35,-,0.01]$ in plane and hinge modes, where $-$ are
omitted from the cost. In this experimental setup, the proposed approach was
compared to a simple MPC approach without the CEM warm-start, in order to
validate experimentally the results achieved in simulation. From the plots in
Figure 8, it can be seen that without the CEM warm-start, the impedance rest
position shows a small discontinuity in its trend, and this can be noticed
especially when the first transition is triggered. This result is in line with
what we have shown in the simulation study, the proposed CEM+MPC approach
helps in decreasing the total solve time at contact transitions.

(a) Simple gradient-based MPC

(b) Proposed CEM+MPC approach

Figure 8: Comparison between the simple MPC and CEM+MPC approaches for the
pivot task.

##  6 Conclusion

This paper proposed an approach to planning on discontinuous dynamics which
uses a modified CEM to initialize a gradient-based planner. It was shown that
such an approach can reduce the variation in solve time when discontinuities
in dynamics occur. The approach was also vlidated on real contact tasks,
showing that planning method results in safe robot trajectories which respond
to the contact conditions.

## References

  * [1] M. Anitescu and F. A. Potra, “Formulating Dynamic Multi-Rigid-Body Contact Problems with Friction as Solvable Linear Complementarity Problems,” p. 17.
  * [2] M. Posa, C. Cantu, and R. Tedrake, “A direct method for trajectory optimization of rigid bodies through contact,” _The International Journal of Robotics Research_ , vol. 33, no. 1, pp. 69–81, 2014.
  * [3] S. L. Cleac’h, T. Howell, M. Schwager, and Z. Manchester, “Fast Contact-Implicit Model-Predictive Control,” _arXiv:2107.05616 [cs, eess]_ , Sep. 2021.
  * [4] A. Aydinoglu and M. Posa, “Real-Time Multi-Contact Model Predictive Control via ADMM,” _arXiv:2109.07076 [cs]_ , Mar. 2022.
  * [5] W. Jin and M. Posa, “Task-Driven Hybrid Model Reduction for Dexterous Manipulation,” Nov. 2022.
  * [6] F. R. Hogan and A. Rodriguez, “Reactive planar non-prehensile manipulation with hybrid model predictive control,” _The International Journal of Robotics Research_ , vol. 39, no. 7, pp. 755–773, Jun. 2020.
  * [7] C. Chen, P. Culbertson, M. Lepert, M. Schwager, and J. Bohg, “TrajectoTree: Trajectory Optimization Meets Tree Search for Planning Multi-contact Dexterous Manipulation,” in _2021 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS)_ , Sep. 2021, pp. 8262–8268.
  * [8] K. Wu and G. Hao, “Design and nonlinear modeling of a novel planar compliant parallelogram mechanism with general tensural-compresural beams,” _Mechanism and Machine Theory_ , vol. 152, p. 103950, 2020.
  * [9] T. Pang, H. J. T. Suh, L. Yang, and R. Tedrake, “Global Planning for Contact-Rich Manipulation via Local Smoothing of Quasi-dynamic Contact Models,” Jun. 2022.
  * [10] H. Bharadhwaj, K. Xie, and F. Shkurti, “Model-Predictive Control via Cross-Entropy and Gradient-Based Optimization,” Apr. 2020.
  * [11] K. Huang, S. Lale, U. Rosolia, Y. Shi, and A. Anandkumar, “CEM-GD: Cross-Entropy Method with Gradient Descent Planner for Model-Based Reinforcement Learning,” Dec. 2021.
  * [12] H. J. T. Suh, T. Pang, and R. Tedrake, “Bundled Gradients through Contact via Randomized Smoothing,” _arXiv:2109.05143 [cs]_ , Jan. 2022.
  * [13] A. Ö. Önol, R. Corcodel, P. Long, and T. Padır, “Tuning-Free Contact-Implicit Trajectory Optimization,” in _2020 IEEE International Conference on Robotics and Automation (ICRA)_ , May 2020, pp. 1183–1189.
  * [14] D. E. Stewart, “Rigid-Body Dynamics with Friction and Impact,” _SIAM Review_ , vol. 42, no. 1, pp. 3–39, Jan. 2000.
  * [15] F. R. Hogan, J. Ballester, S. Dong, and A. Rodriguez, “Tactile Dexterity: Manipulation Primitives with Tactile Feedback,” in _2020 IEEE International Conference on Robotics and Automation (ICRA)_ , May 2020, pp. 8863–8869.
  * [16] J. B. Rawlings, D. Q. Mayne, and M. Diehl, _Model Predictive Control: Theory, Computation, and Design_. Nob Hill Publishing Madison, 2017, vol. 2.
  * [17] A. Nagabandi, G. Kahn, R. S. Fearing, and S. Levine, “Neural network dynamics for model-based deep reinforcement learning with model-free fine-tuning,” _CoRR_ , vol. abs/1708.02596, 2017. [Online]. Available: http://arxiv.org/abs/1708.02596
  * [18] G. Williams, N. Wagener, B. Goldfain, P. Drews, J. M. Rehg, B. Boots, and E. A. Theodorou, “Information theoretic MPC for model-based reinforcement learning,” in _2017 IEEE International Conference on Robotics and Automation (ICRA)_ , May 2017, pp. 1714–1721.
  * [19] R. Y. Rubinstein, “Optimization of computer simulation models with rare events,” _European Journal of Operational Research_ , vol. 99, no. 1, pp. 89–112, May 1997. [Online]. Available: https://ideas.repec.org/a/eee/ejores/v99y1997i1p89-112.html
  * [20] M. Okada and T. Taniguchi, “Variational inference MPC for bayesian model-based reinforcement learning,” _CoRR_ , vol. abs/1907.04202, 2019. [Online]. Available: http://arxiv.org/abs/1907.04202
  * [21] C. Pinneri, S. Sawant, S. Blaes, J. Achterhold, J. Stueckler, M. Rolinek, and G. Martius, “Sample-efficient Cross-Entropy Method for Real-time Planning,” Aug. 2020.
  * [22] L. Roveda, J. Maskani, P. Franceschi, A. Abdi, F. Braghin, L. M. Tosatti, and N. Pedrocchi, “Model-based reinforcement learning variable impedance control for human-robot collaboration,” _Journal of Intelligent & Robotic Systems_, pp. 1–17, 2020. [Online]. Available: https://api.semanticscholar.org/CorpusID:216293188
  * [23] K. Chua, R. Calandra, R. McAllister, and S. Levine, “Deep Reinforcement Learning in a Handful of Trials using Probabilistic Dynamics Models,” _arXiv:1805.12114 [cs, stat]_ , May 2018.
  * [24] H. J. T. Suh, M. Simchowitz, K. Zhang, and R. Tedrake, “Do Differentiable Simulators Give Better Policy Gradients?” Aug. 2022.
  * [25] A. M. Castro, A. Qu, N. Kuppuswamy, A. Alspach, and M. Sherman, “A Transition-Aware Method for the Simulation of Compliant Contact With Regularized Friction,” _IEEE Robotics and Automation Letters_ , vol. 5, no. 2, pp. 1859–1866, Apr. 2020.
  * [26] S. Thrun, “Probabilistic robotics,” _Communications of the ACM_ , vol. 45, no. 3, pp. 52–57, 2002.


# Mass discrepancy in galaxy clusters as a result of the offset between dark
matter and baryon distributions

HuanYuan Shan1,2, Bo Qin2,⋆, and HongSheng Zhao2,3
1Department of Astronomy, School of Physics, Peking University, Beijing,
100871, China
2National Astronomical Observatories, Chinese Academy of Sciences, Beijing
100012, China
3SUPA, University of St Andrews, KY16 9SS, UK E-mail<EMAIL_ADDRESS><EMAIL_ADDRESS>
(Accepted …. Received …; in original form …)

###### Abstract

Recent studies of lensing clusters reveal that it might be fairly common for a
galaxy cluster that the X-ray center has an obvious offset from its
gravitational center which is measured by strong lensing. We argue that if
these offsets exist, then X-rays and lensing are indeed measuring different
regions of a cluster, and may thus naturally result in a discrepancy in the
measured gravitational masses by the two different methods. Here we
investigate theoretically the dynamical effects of such lensing-X-ray offsets,
and compare with observational data. We find that for typical values, the
offset alone can give rise to a factor of two difference between the lensing
and X-ray determined masses for the core regions of a cluster, suggesting that
such “offset effect” may play an important role and should not be ignored in
our dynamical measurements of clusters.

###### keywords:

dark matter-gravitational lensing-X-rays: galaxies: clusters

###  1 Introduction

Galaxy clusters, the largest gravitationally-bound structures in the universe,
are ideal cosmological tools. Accurate measurements of their masses provide a
crucial observational constraint on cosmological models. Several dynamical
methods have been available to estimate cluster masses, such as (1) optical
measurements of the velocity dispersions of cluster galaxies, (2) measurements
of the X-ray emitting gas, and (3) gravitational lensing. Good agreements
between these methods have been found on scales larger than cluster cores.

However, joint measurements of lensing and X-rays often identify large
discrepancies in the gravitational masses within the central regions of
clusters by the two methods, and the lensing mass has always been found to be
$2-4$ times higher than the X-ray determined mass. This is the so-called “Mass
Discrepancy Problem” (Allen 1998; Wu 2000). Many plausible explanations have
been suggested, e.g., the triaxiality of galaxy clusters (Morandi et al.
2010), the oversimplification of the strong lensing model for the central mass
distributions of clusters (Bartelmann & Steinmetz 1996), the inappropriate
application of the hydrostatic equilibrium hypothesis for the central regions
of clusters (Wu 1994; Wu & Fang 1997), or the magnetic fields in clusters
(Loeb & Mao 1994).

Recently Richard et al. (2010) present a sample of $20$ strong lensing
clusters taken from the Local Cluster Substructure Survey (LoCuSS), among
which $18$ clusters have X-ray data from Chandra observations (Sanderson et
al. 2009). They show that the X-ray/lensing mass discrepancy is $1.3$ at
$3\sigma$ significance — clusters with larger substructure fractions show
greater mass discrepancies, and thus greater departures from hydrostatic
equilibrium.

On the other hand, lensing observations of the bullet cluster 1E0657-56 (Clowe
et al. 2006), combined with earlier X-ray measurements (Markevitch et al.
2006), clearly indicate that the gravitational center of the cluster has an
obvious offset from its baryonic center. Furthermore, recent studies (Shan et
al. 2010) of lensing galaxy clusters reveal that offset between the lensing
center and X-ray center appears to be quite common, especially for unrelaxed
clusters. Among the recent sample of 38 clusters of Shan et al. (2010), $45\%$
have been found to have offsets greater than $10^{\prime\prime}$, and $5$
clusters even have offsets greater than $40^{\prime\prime}$. Motivated by such
observations, we propose to investigate galaxy cluster models where the center
of the dark matter (DM) halo does not coincide with the center of the X-ray
gas (See Figure 1).

Figure 1: Offset between the dark matter center and the X-ray center in a
galaxy cluster.

If the X-ray center of a cluster has an offset from its lensing
(gravitational) center, then the X-rays and lensing are indeed measuring
different regions of the cluster. Given the same radius, the lensing is
measuring the DM halo centered at the gravitational center (shown by the dark
blue sphere in Figure 1 ), while the X-rays are measuring the sphere of the
halo that is offset from the true gravitational center (shown by the red
circle in Figure 1). In this case, there will always be a natural discrepancy
between the lensing and X-ray measured masses — or specifically, the X-ray
mass will always be lower than the lensing mass, just as the long-standing
“mass discrepancy problem” has indicated.

In this paper, we investigate the lensing-X-ray mass discrepancy caused by the
offsets between DM and X-ray gas. To check our predictions, we compile a
sample of $27$ clusters with good lensing and X-ray measurements. We conclude
that such “offset” effect should not be ignored in our dynamical measurements
of galaxy clusters. A flat $\Lambda$CDM cosmology is assumed throughout this
paper, where $\Omega_{m}$=0.3, $\Omega_{\Lambda}$=0.7, and $\rm
H_{0}=70\,km\,s^{-1}Mpc^{-1}$.

###  2 Mass discrepancy as a result of the dark matter-baryon offset

We model our galaxy cluster with a fiducial model as the following: (1) the DM
halo is modeled by the Navarro-Frenk-White (NFW) profile (Navarro et al. 1997)
with concentration $c=4.32$ and scaled radius $r_{s}=516\,{\rm kpc}$, (2) the
gas distribution is modeled by a $\beta$ model with $\beta=0.65$, the cluster
core radius $r_{c}=200\,{\rm kpc}$, and the gas fraction $f_{\rm gas}=12\%$,
(3) the mass density of the BCG is described by a Singular Isothermal Sphere
(SIS) with a velocity dispersion of $300\,\rm km/s$.

The projected mass within a sphere of radius $R_{x}$ is

$\displaystyle m(R_{x},d)$ $\displaystyle=$
$\displaystyle\int_{0}^{2\pi}\int_{0}^{R_{x}}\left[\Sigma_{\rm
NFW}(R^{\prime})+\Sigma_{\rm gas}(R)\right.$ $\displaystyle\left.+\Sigma_{\rm
BCG}(R^{\prime})\right]R^{\prime}\,dR^{\prime}\,d{\theta},$

where $R^{\prime}=\sqrt{d^{2}+R^{2}+2dR\cos\theta}$ is the 2-D radius from the
halo center, $R$ is the 2-D radius from the X-ray gas center, $d$ is the 2-D
offset between the halo center and X-ray center, and $\Sigma_{\rm NFW}$,
$\Sigma_{\rm gas}$, and $\Sigma_{\rm BCG}$ are the projected mass densities of
the DM halo, the gas and the BCG, respectively. For a given radius $R_{x}$,
the gravitational mass measured by lensing $m_{\rm lens}$ can be given by
$m(R_{x},0)$ (as shown by the dark blue sphere in Figure 1), while the
projected mass measured by X-rays $m_{\rm xray}$ is described by $m(R_{x},d)$
(the mass within the red circle in Figure 1). We now calculate the mass ratio
$m(R_{x},d)/m(R_{x},0)$, or equivalently, $m_{\rm lens}/m_{\rm xray}$.

Figure 2 shows the mass ratio as a function of the 2-D offset $d$, for a
typical rich cluster. The solid curves are the mass ratio with the fiducial
model, the dashed and dotted curves are the mass ratio with the NFW
concentration $c=4.04$ and $5.13$ (top left), the cluster core radius
$r_{c}=150\,{\rm kpc}$ and $400\,{\rm kpc}$ (top right), the $\beta$ index
$\beta=0.6$ and $0.9$ (bottom left), the gas fraction $f_{\rm gas}=0.1$ and
$0.2$, respectively. For these cases, the three curves from top to bottom are
for the three measuring radii $R_{x}=50\,{\rm kpc},100\,{\rm kpc},200\,{\rm
kpc}$, respectively. From Figure 2 we have the following conclusions:

(1) The lensing measured mass $m_{\rm lens}$ is always higher than the X-ray
measured mass $m_{\rm xray}$. For typical values of offset $d=100\,{\rm kpc}$
and $R_{x}=100\,{\rm kpc}$, $m_{\rm lens}/m_{\rm xray}\sim 2$, comparable to
the ratio found in early studies (Allen 1998; Wu 2000; Richard et al. 2010).

(2) The “Offset Effect” we are reporting here should contribute significantly
to the long-standing “Mass Discrepancy Problem”.

(3) The ratio of $m_{\rm lens}/m_{\rm xray}$ increases with offset $d$.

(4) $m_{\rm lens}/m_{\rm xray}$ depends very strongly on $R_{x}$. Here $R_{x}$
acts like the arc radius $r_{\rm arc}$ in strong lensing, i.e., we only
measure the enclosed mass within a small region of $R\leq R_{x}$. When $R_{x}$
is very small, the offset effect is most prominent and gives large $m_{\rm
lens}/m_{\rm xray}$. Increasing $R_{x}$ will reduce $m_{\rm lens}/m_{\rm
xray}$. When $R_{x}$ is very large (compared with $d$), the offset effect will
be “smeared out”, and the $m_{\rm lens}$-$m_{\rm xray}$ discrepancy introduced
by the offset will vanish.

(5) The mass ratio is very sensitive to the NFW concentration, and it
increases dramatically with $c$.

(6) The mass ratio increases with the core radius, and decreases with $\beta$
index and gas fraction. However, the mass ratio is not very sensitive to the
gas model.

Figure 2: Ratio of projected gravitational masses, as a function of the 2-D
offset $d$ and the measuring radius $R_{x}$. The solid curves are the mass
ratio for the fiducial model with $c=4.32$, $r_{s}=516\rm kpc$, $\beta=0.65$,
$r_{c}=200\rm kpc$, and $f_{\rm gas}=0.12$. The dashed and dotted curves are
for the NFW concentration $c=4.04$ and $5.13$ (top left), the cluster core
radius $r_{c}=150\,{\rm kpc}$ and $400\,{\rm kpc}$ (top right), the $\beta$
index $\beta=0.6$ and $0.9$ (bottom left), the gas fraction $f_{\rm gas}=0.1$
and $0.2$, respectively. The three dotted (dashed, solid as well) curves from
top to bottom in one panel correspond to $R_{x}=50,100,200\,\rm kpc$,
respectively.

###  3 Comparison with Observational Data

To compare with our theoretical predictions, we compile a sample of $27$
clusters with $48$ arc-like images, which have both strong lensing and X-ray
measurements. The clusters and their lensing and X-ray data are listed in
Table 1. For the $22$ arcs that have no redshift information, we estimate
their lensing masses $m_{\rm lens}$ by assuming the mean redshifts of
$\left<z_{d}\right>=0.8$ and $2.0$, respectively. The X-ray data are taken
from Tucker et al. (1998), Wu (2000), Bonamente et al. (2006), and references
therein. The offsets between lensing and X-ray centers are taken from Shan et
al. (2010). The clusters in our table are classified as relaxed (with cooling
flow) and unrelaxed (which are dynamically unmature), from their X-ray
morphologies. The definition has been used in the literature by Allen (1998),
Wu (2000), Baldi et al. (2007), and Dunn & Fabian (2008).

Mass from strong lensing. Assuming a spherical matter distribution, one can
calculate the gravitational mass of a galaxy cluster projected within a radius
of $r_{\rm arc}$ on the cluster plane as

$m_{\rm lens}(<r_{\rm arc})=\pi r_{\rm arc}^{2}\Sigma_{\rm crit},$ (1)

where $\Sigma_{\rm crit}=\frac{c^{2}}{4\pi G}\frac{D_{s}}{D_{l}D_{ls}}$ is the
critical surface mass density, $D_{l}$, $D_{s}$ and $D_{ls}$ are the angular
diameter distances to the cluster, to the background galaxy, and from the
cluster to the galaxy, respectively. The above equation is actually the
lensing equation for a cluster lens of spherical mass distribution with a
negligible small alignment parameter for the distant galaxy within $r_{\rm
arc}$. The values of $m_{\rm lens}$ within the arc radius $r_{\rm arc}$ are
listed in Table 1.

Allen (1998) pointed out that the use of more realistic, elliptical mass
models can reduce the masses within the arc radii by up to $40\%$, though a
value of $20\%$ is more typical. However, such corrections are still not very
significant compared with the large discrepancies between the lensing and
X-ray determined masses. We will discuss it in more detail in the next
section.

Mass from X-rays. Assuming that the intra-cluster gas is isothermal and in
hydrostatic equilibrium, the cluster mass $m(r)$ enclosed within a radius $r$
can be easily calculated from

$-\frac{Gm(r)}{r^{2}}=\frac{kT}{\mu m_{p}}\frac{d{\,\rm ln}n_{\rm
gas}(r)}{dr},$ (2)

where $T$ is the gas temperature, $n_{\rm gas}$ the gas number density,
$m_{p}$ the proton mass, and $\mu=0.585$ the mean molecular weight. Here we
assume that the gas follows the conventional $\beta$ model, i.e., $n_{\rm
gas}(r)=n_{\rm gas}(0)(1+r^{2}/r_{c}^{2})^{-\frac{3\beta}{2}}$. In order to
compare the mass measured by X-rays with the lensing result, we need to
convert this $m(r)$ (i.e., 3-D) into the projected mass $m_{\rm xray}$ (see
e.g. Wu 1994):

$m_{\rm xray}=1.13\times 10^{13}\beta\bar{f}\left({R\over
r_{c}}\right)\left(\frac{r_{c}}{0.1{\rm Mpc}}\right)\left(\frac{kT}{1{\rm
keV}}\right)M_{\odot},$ (3)

where

$\bar{f}(y)=\frac{\pi y^{2}}{2(1+y^{2})^{1/2}},$ (4)

the mass ratio $m_{\rm lens}/m_{\rm xray}$ are listed in Table 1.

Figure 3 shows the relation between the mass ratios $m_{\rm lens}/m_{\rm
xray}$ and the (scaled) offsets for our sample of $27$ clusters ($48$ arc
images). It should be pointed out that the $27$ clusters in our sample have
quite different sizes and masses. This can be seen from the wide range of the
cluster temperatures — from $4~{}\rm keV$ to $14~{}\rm keV$. Therefore, it is
useful to compare the offsets of the clusters on the same scale. We realize
that the M-T relation of clusters scales as${\rm M}\sim{\rm T}^{3/2}$ (e.g.,
Nevalainen et al. 2000; Xu et al. 2001), and that ${\rm M}\sim{\rm R}^{3}$,
where R is the size of the cluster. Therefore, in Figure 3, instead of using
the physical offset $d$, we use a scaled offset which is characterized by
${\rm d_{\rm kpc}}/{\rm T_{\rm keV}}^{1/2}$.

From Figure 3, the mass ratios $m_{\rm lens}/m_{\rm xray}$ exhibits large
dispersions — roughly ranging from $2$ to $4$. Many clusters have large error
bars. It appears that relax clusters (marked by crosses) have smaller $m_{\rm
lens}/m_{\rm xray}$ ratios. The fact that $m_{\rm lens}>m_{\rm xray}$ is
consistent with our theoretical predictions, and the ratio of $m_{\rm
lens}/m_{\rm xray}\sim 2-4$ is also roughly consistent with our predictions as
plotted in Figure 2.

However, no strong correlation has been found between the offset and mass
discrepancies. We notice that many clusters in the sample have very small
offset values — smaller than the errors in lensing and X-ray measurements
which are typically a few arcseconds. So these offset values are not robustly
measured themselves, and we thus remove these data points and only focus on
clusters with large offsets of $d>10^{\prime\prime}$, as has been suggested in
Shan et al. (2010). This leaves a sub-sample of only $24$ arc images. The
dashed line in Figure 3 shows a $\chi^{2}$ fit to this sub-sample, which
satisfies $m_{\rm lens}/m_{\rm xray}\sim 3.24(\frac{d/100\,{\rm
kpc}}{\sqrt{kT/{\rm keV}}})^{0.20}$ with a reasonable $\chi^{2}=0.75$. We can
find $m_{lens}/m_{xray}$ increasing slightly with $d$.

Figure 3: The ratio of lensing and X-ray determined masses for our sample of
$27$ clusters ($48$ arc images). The x-label show the scaled offset between DM
and baryons. The squares denote unrelaxed clusters, and crosses the relaxed
clusters. The dashed line shows a $\chi^{2}$ fit satisfying $m_{\rm
lens}/m_{\rm xray}\sim 3.24(\frac{d/100\,{\rm kpc}}{\sqrt{kT/{\rm
keV}}})^{0.20}$ with $\chi^{2}=0.75$ for the clusters with offset larger than
$10^{\prime\prime}$.

###  4 Discussion and Conclusions

As has been reported by Shan et al. (2010), it might be fairly common in
galaxy clusters that the X-ray center has an obvious offset from the
gravitational center. We have explored the dynamical consequences of this
lensing-X-ray offset and tried to attribute such an effect to the long-
standing “Mass Discrepancy Problem” in galaxy clusters. Our theoretical model
predicts that such an offset effect will always result in a larger $m_{\rm
lens}$ than $m_{\rm xray}$, with a typical mass ratio $m_{\rm lens}/m_{\rm
xray}\sim 2$, which is consistent with observations.

To test our model, we have compiled a sample of $27$ clusters, and studied in
detail their lensing and X-ray properties and obtained their lensing and X-ray
masses, $m_{\rm lens}$ and $m_{\rm xray}$. The lack of strong correlation
between $m_{\rm lens}/m_{\rm xray}$ and the offset $d$ suggests that the
problem is more complicated. As we have found in Section 2, $m_{\rm
lens}/m_{\rm xray}$ is not only a function of $d$, but also depends very
strongly on $R_{x}$ (or the arc radius $r_{\rm arc}$). Apparently, each
cluster in our sample has quite different $r_{\rm arc}$.

Probably, other mechanisms than the offset effect should play important roles,
and the lensing-X-ray mass discrepancy may not be just from one mechanism, but
a combination of many effects:

(1) The central regions of clusters may be still undergoing dynamical
relaxation, and the X-ray gas may not be in good hydrostatic equilibrium.
Therefore, large errors could be induced in the X-ray measurement of cluster
cores, especially for unrelaxed clusters.

(2) The spherical models are too simple to reflect the real mass distribution
of clusters. The use of more realistic mass model could reduce the lens mass
within the arc radius by up to $40\%$, though values of $\sim 20\%$ are more
typical (Bartelmann 1995; Allen 1998).

(3) The presence of substructures may complicate our simple spherical lens
model, and hence could be a main source of uncertainties in $m_{\rm lens}$.
The absence of the secondary arc-like images in most arc-cluster systems may
indicate the limitations of the spherical mass distribution in the central
regions of clusters.

It should be noted that the mass ratios we obtained here are slightly higher
than Allen (1998) and Wu (2000) because they unfortunately used a Hubble
constant of $\rm H_{0}=50\,km\,s^{-1}Mpc^{-1}$. The use of $\rm
H_{0}=70\,km\,s^{-1}Mpc^{-1}$ here will of course make the mass discrepancy
problem more pronounced.

It should be noted that the gas represents only a $10\%$ perturbation due to
the small ratio of gas-to-DM in the central region, likewise the offset of the
gas is only a small perturbation (less than $10\%$) to the otherwise
concentric matter density or potential. It is unlikely to create a factor of
two difference in the lensing-derived enclosed masses within an arc.

To illustrate the lensing effect of the offset perturbation and triaxiality,
we show the critical curves in Figure 4. The solid curves indicate the
critical curve of circular NFW plus $\beta$ model without offset, the dotted
curves indicate the critical curve of elliptical NFW plus $\beta$ model with
offset $d=10^{\prime\prime}$. The square and cross denote the center of dark
matter and the hot gas, respectively. For the NFW profile, $c=4.3,r_{s}=516\rm
kpc$; for the $\beta$ model, $\beta=0.65,r_{c}=150\rm kpc$. We also introduce
the triaxiality with the ellipticity $e=0.15$ and position angle
$\theta=30^{\circ}$. We also assume the lens and source redshifts $z_{l}=0.3$,
$z_{s}=1$. We can see that the predicted critical curves (dotted lines) have
very similar sizes as the predicted critical curves for a benchmark model
(solid lines) with the same mass DM and gas mass but in concentric spheres.

Figure 4: The effects of offset and triaxiality on the critical curves. The
solid curves indicate the critical curve of circular NFW & $\beta$ model
without offset. The dotted curves indicate the critical curve of elliptical
NFW & $\beta$ model with offset $d=10^{\prime\prime}$. The square and cross
denote the center of dark matter and the X-ray gas, respectively. For the NFW
profile, $c=4.3,r_{s}=516\rm kpc$; for the $\beta$ model,
$\beta=0.65,r_{c}=150\rm kpc$. The ellipticity and position angle are $e=0.15$
and $\theta=30^{\circ}$. The lens and source redshifts are $z_{l}=0.3$,
$z_{s}=1$.

Early studies have suggested that statistically unrelaxed clusters have larger
mass discrepancies than relaxed clusters (Allen 1998; Wu 2000; Richard et al.
2010). As Shan et al. (2010) have reported, the clusters with large offset of
$d>10^{\prime\prime}$ are all unrelaxed clusters. If such offsets exist and
are big, then they must come into play in our dynamical studies of galaxy
cluster, and should not be ignored, especially for unrelaxed clusters.

##  Chapter  Acknowledgments

We thank Bernard Fort, Charling Tao, and Xiang-Ping Wu for discussions, and an
anonymous referee for helpful suggestions. HYS and BQ are grateful to the CPPM
for hospitality. This work was supported by the National Basic Research
Program of China (973 Program) under grant No. 2009CB24901, and CAS grants
KJCX3-SYW-N2 and KJCX2-YW-N32.

## References

  * (1) Allen, S.W., 1998, MNRAS, 296, 392
  * (2) Baldi, A., et al., 2007, ApJ, 666, 835
  * (3) Bartelmann, M.,1995, A&A, 299, 11
  * (4) Bartelmann, M., & Evrard, A. E., 1996, MNRAS, 283,431
  * (5) Bonamente, M., Joy, M., LaRoque, S., Carlstrom, J., Reese, E., & Dawson, K., 2006, ApJ, 647, 25
  * (6) Bradac, M., et al., 2006, ApJ, 652, 937
  * (7) Dunn, R. J. H., & Fabian, A. C., 2008, MNRAS, 385, 757
  * (8) Gioia, I. M., Shaya, E. J., Le Fevre, O., Falco, E. E., Luppino, G. A., & Hammer, F., 1998, ApJ, 497, 573
  * (9) Jee, M.J., et al., 2007, ApJ, 661, 728
  * (10) Kneib, J.-P., Mellier, Y., Fort, B., & Mathez, G., 1993, A&A, 273, 367
  * (11) Limousin, M., et al., 2007, ApJ, 668, 643
  * (12) Loeb, A., & Mao, S., 1994, ApJ, 435, 109
  * (13) Markevitch, M., 2006, xru, conf, 723
  * (14) Morandi, A., Pedersen, K., & Limousin, M., 2010, ApJ, 713, 491
  * (15) Morandi, A., Pedersen, K., & Limousin, M., 2010, arXiv: astro-ph/1001.1656
  * (16) Navarro, J.F., Frenk, C.S., & White, S.D.M., 1997, ApJ, 490, 493
  * (17) Nevalainen, J., Markevitch, M., & Forman, W., 2000, ApJ, 532, 694
  * (18) Newbury, & Fahlman, 1999, arXiv: astro-ph/9905254
  * (19) Richard, J., et al., 2010, MNRAS, 404, 325
  * (20) Sand, D. J., Treu, T., Ellis, R. S., & Smith G. P., 2005, ApJ, 627, 32
  * (21) Sanderson, A.J.R., Edge, A.C., & Smith, G.P., 2009, MNRAS, 398, 1698
  * (22) Shan, H.Y., Qin, B., Fort, B., Tao, C., Wu, X.-P., & Zhao, H.S., 2010, MNRAS, accepted
  * (23) Smith, G. P., Kneib, J., Smail, I., Mazzotta, P., Ebeling, H., & Czoske, O., 2005, MNRAS, 359, 417
  * (24) Tucker, W., et al., 1998, ApJ, 496, 5
  * (25) Wu, X.-P., 1994, ApJ, 436, 115
  * (26) Wu, X.-P., & Fang, L.-Z., 1997, ApJ, 483, 62
  * (27) Wu, X.-P., 2000, MNRAS, 316, 299
  * (28) Xu, H., Jin, G., & Wu, X.-P., 2001, ApJ, 553, 78

Table 1: The X-ray and lensing mass discrepancies of $27$ clusters. For the
$22$ arcs that have no redshift information, we estimate the mean redshifts of
$\left<z_{d}\right>=0.8$ and $2.0$, respectively. Refs A and B give the
references of the lensing and X-ray data, respectively. The last column shows
the classification of the clusters: “R/U” means relaxed/unrelaxed.

Cluster | $z_{\rm cluster}$ | Offset | $z_{\rm arc}$ | $r_{\rm arc}$ | $m_{\rm lens}$ |  Ref. Ad | kT | $\beta$ | $r_{c}$ | $m_{\rm xray}$ | $m_{\rm lens}/m_{\rm xray}$ |  Classc |  Ref.Bd
---|---|---|---|---|---|---|---|---|---|---|---|---|---
|  | (arcsec) | (kpc) |  | (Mpc) | $\times 10^{14}M_{\odot}$ |  | (keV) |  | (Mpc) | $\times 10^{14}M_{\odot}$ |  |  |
1E0657-56 | 0.296 | 47.4 | 209.2 | 3.24 | 0.25 | 4.37 | 3,4 | $14.1^{+0.2}_{-0.2}$ | $0.62^{+0.07}_{-0.07}$ | $0.36^{+0.05}_{-0.05}$ | $2.15^{+0.46}_{-0.46}$ | $2.03^{+0.44}_{-0.44}$ | U | 12
A68a | 0.255 | 14.3 | 56.7 | 1.60 | 0.04 | 0.13 | 11 | $10.0^{+1.1}_{-0.9}$ | $0.72^{+0.04}_{-0.03}$ | $0.25^{+0.02}_{-0.02}$ | $0.08^{+0.019}_{-0.016}$ | $1.66^{+0.39}_{-0.34}$ | U | 2
|  |  |  | 1.60 | 0.10 | 0.80 |  |  |  |  | $0.47^{+0.11}_{-0.091}$ | $1.77^{+0.40}_{-0.34}$ |  |
|  |  |  | 2.63 | 0.11 | 0.94 |  |  |  |  | $0.56^{+0.13}_{-0.11}$ | $1.67^{+0.38}_{-0.32}$ |  |
|  |  |  | … | 0.211 |  4.49(3.54)b |  |  |  |  | $1.69^{+0.35}_{-0.29}$ | $2.64^{+0.55}_{-0.46}(2.08^{+0.43}_{-0.36})$ |  |
|  |  |  | 0.86 | 0.28 | 7.49 |  |  |  |  | $2.61^{+0.51}_{-0.43}$ | $2.95^{+0.58}_{-0.48}$ |  |
|  |  |  | … | 0.27 |  7.26(5.72)b |  |  |  |  | $2.47^{+0.49}_{-0.41}$ | $2.98^{+0.59}_{-0.49}(2.35^{+0.47}_{-0.39})$ |  |
|  |  |  | 1.27 | 0.32 | 8.66 |  |  |  |  | $3.14^{+0.60}_{-0.50}$ | $2.82^{+0.54}_{-0.45}$ |  |
|  |  |  | … | 0.12 |  1.54(1.22)b |  |  |  |  | $0.65^{+0.15}_{-0.12}$ | $2.23^{+0.50}_{-0.42}(1.76^{+0.39}_{-0.33})$ |  |
A267 | 0.230 | 9.62 | 35.3 | … | 0.12 |  1.48(1.20)b | 11 | $6.0^{+0.6}_{-0.5}$ | $0.71^{+0.03}_{-0.03}$ | $0.19^{+0.01}_{-0.01}$ | $0.39^{+0.08}_{-0.07}$ | $3.86^{+0.78}_{-0.71}(3.13^{+0.63}_{-0.58})$ | U | 2
A370a | 0.375 | 19.9 | 102.7 | 1.30 | 0.41 | 13.1 | 7 | $7.13^{+1.05}_{-1.05}$ | $0.95^{+0.75}_{-0.35}$ | $0.56^{+0.44}_{-0.26}$ | $4.34^{+4.10}_{-2.27}$ | $3.00^{+2.83}_{-1.57}$ | U | 13
|  |  |  | 0.72 | 0.19 | 4.09 |  |  |  |  | $0.71^{+1.15}_{-0.65}$ | $5.84^{+9.54}_{-5.42}$ |  |
A697 | 0.282 | 3.07 | 13.1 | … | 0.12 |  1.51(1.15)b | 10 | $9.9^{+0.6}_{-0.6}$ | $0.61^{+0.01}_{-0.01}$ | $0.24^{+0.01}_{-0.01}$ | $0.26^{+0.21}_{-0.13}$ | $5.53^{+4.52}_{-2.84}(4.21^{+3.44}_{-2.16})$ | U | 2
A773a | 0.217 | 6.43 | 22.6 | 0.65 | 0.11 | 1.39 | 11 | $7.6^{+0.5}_{-0.4}$ | $0.61^{+0.01}_{-0.01}$ | $0.19^{+0.06}_{-0.06}$ | $0.37^{+0.042}_{-0.037}$ | $3.75^{+0.42}_{-0.38}$ | U | 2
|  |  |  | 0.40 | 0.21 | 7.08 |  |  |  |  | $1.26^{+0.26}_{-0.24}$ | $5.85^{+1.19}_{-1.11}$ |  |
|  |  |  | … | 0.25 |  6.50(5.36)b |  |  |  |  | $1.61^{+0.29}_{-0.26}$ | $4.10^{+0.73}_{-0.68}(3.38^{+0.60}_{-0.56})$ |  |
|  |  |  | … | 0.23 |  5.41(4.45)b |  |  |  |  | $1.43^{+0.27}_{-0.25}$ | $3.89^{+0.74}_{-0.69}(3.21^{+0.61}_{-0.57})$ |  |
|  |  |  | 1.11 | 0.213 | 4.34 |  |  |  |  | $1.26^{+0.26}_{-0.24}$ | $3.35^{+0.68}_{-0.64}$ |  |
|  |  |  | 0.40 | 0.16 | 4.14 |  |  |  |  | $0.84^{+0.20}_{-0.19}$ | $5.12^{+1.23}_{-1.17}$ |  |
|  |  |  | … | 0.04 |  0.18(0.15)b |  |  |  |  | $0.07^{+0.023}_{-0.022}$ | $2.51^{+0.86}_{-0.83}(2.07^{+0.71}_{-0.69})$ |  |
|  |  |  | 0.49 | 0.23 | 7.42 |  |  |  |  | $1.43^{+0.27}_{-0.25}$ | $5.07^{+0.96}_{-0.90}$ |  |
A963a | 0.206 | 7.10 | 24.0 | … | 0.057 |  0.35(0.29)b | 11 | $6.13^{+0.45}_{-0.30}$ | $0.51^{+0.04}_{-0.04}$ | $0.11^{+0.02}_{-0.02}$ | $0.14^{+0.038}_{-0.034}$ | $2.41^{+0.63}_{-0.57}(2.01^{+0.53}_{-0.48})$ | R | 13
|  |  |  | 0.71 | 0.09 | 0.87 |  |  |  |  | $0.31^{+0.073}_{-0.066}$ | $2.89^{+0.68}_{-0.61}$ |  |
A1689 | 0.183 | 0.60 | 1.85 | … | 0.20 |  4.5(3.8)b | 8 | $9.02^{+0.40}_{-0.30}$ | $0.65^{+0.04}_{-0.02}$ | $0.14^{+0.02}_{-0.02}$ | $1.68^{+0.24}_{-0.17}$ | $2.68^{+0.38}_{-0.27}(2.29^{+0.32}_{-0.23})$ | R | 13
A1835 | 0.252 | 1.61 | 6.33 | … | 0.17 |  2.82(2.23)b | 11 | $9.8^{+1.4}_{-1.4}$ | $0.65^{+0.04}_{-0.04}$ | $0.08^{+0.01}_{-0.01}$ | $1.72^{+0.36}_{-0.36}$ | $1.70^{+0.36}_{-0.36}(1.35^{+0.28}_{-0.28})$ | R | 13
A1914 | 0.171 | 11.3 | 32.9 | … | 0.10 |  1.16(1.01)b | 10 | $9.9^{+0.3}_{-0.3}$ | $0.90^{+0.01}_{-0.01}$ | $0.200^{+0.003}_{-0.003}$ | $0.70^{+0.036}_{-0.036}$ | $1.67^{+0.09}_{-0.09}(1.44^{+0.07}_{-0.07})$ | U | 2
A2204a | 0.151 | 1.20 | 3.15 | … | 0.025 |  0.08(0.07)b | 10 | $6.5^{+0.2}_{-0.2}$ | $0.48^{+0.002}_{-0.002}$ | $0.02^{+0.0003}_{-0.0003}$ | $0.11^{+0.0040}_{-0.0040}$ | $0.71^{+0.03}_{-0.03}(0.63^{+0.02}_{-0.02})$ | R | 2
|  |  |  | … | 0.01 |  0.013(0.012)b |  |  |  |  | $0.03^{+0.0009}_{-0.0009}$ | $0.49^{0.02}_{-0.02}(0.43^{+0.02}_{-0.02})$ |  |
A2163 | 0.203 | 44.0 | 146.9 | 0.73 | 0.07 | 0.58 | 1 | $14.6^{+0.85}_{-0.85}$ | $0.62^{+0.02}_{-0.02}$ | $0.33^{+0.02}_{-0.02}$ | $0.23^{+0.033}_{-0.033}$ | $2.39^{+0.35}_{-0.35}$ | U | 13
A2218a | 0.176 | 19.1 | 56.9 | 1.03 | 0.28 | 8.60 | 11 | $7.1^{+0.2}_{-0.2}$ | $0.65^{+0.08}_{-0.05}$ | $0.25^{+0.09}_{-0.05}$ | $1.67^{+0.49}_{-0.31}$ | $5.07^{+1.49}_{-0.93}$ | U | 13
|  |  |  | 0.70 | 0.09 | 0.89 |  |  |  |  | $0.25^{+0.11}_{-0.065}$ | $3.94^{+1.73}_{-1.04}$ |  |
|  |  |  | 2.52 | 0.09 | 0.82 |  |  |  |  | $0.25^{+0.11}_{-0.065}$ | $3.20^{+1.40}_{-0.85}$ |  |
A2219a | 0.228 | 11.3 | 41.2 | … | 0.09 |  0.79(0.64)b | 11 | $12.4^{+0.5}_{-0.5}$ | $0.40^{+0.07}_{-0.07}$ | $0.16^{+0.08}_{-0.08}$ | $0.38^{+0.21}_{-0.21}$ | $2.19^{+1.17}_{-1.17}(1.78^{+0.95}_{-0.95})$ | U | 13
|  |  |  | … | 0.12 |  1.54(1.26)b |  |  |  |  | $0.63^{+0.30}_{-0.30}$ | $2.39^{+1.15}_{-1.15}(1.94^{+0.94}_{-0.94})$ |  |
A2259 | 0.164 | 16.3 | 45.9 | 1.48 | 0.04 | 0.13 | 10 | $5.6^{+0.3}_{-0.3}$ | $0.58^{+0.02}_{-0.02}$ | $0.14^{+0.01}_{-0.01}$ | $0.06^{+0.0089}_{-0.0089}$ | $2.71^{+0.38}_{-0.38}$ | U | 2
A2261a | 0.224 | 1.31 | 4.72 | … | 0.12 |  1.49(1.22)b | 10 | $7.2^{+0.4}_{-0.4}$ | $0.56^{+0.01}_{-0.01}$ | $0.08^{+0.004}_{-0.003}$ | $0.71^{+0.058}_{-0.056}$ | $2.12^{+0.17}_{-0.17}(1.73^{+0.14}_{-0.14})$ | R | 2
|  |  |  | … | 0.11 |  1.26(1.03)b |  |  |  |  | $0.63^{+0.053}_{-0.051}$ | $2.00^{+0.17}_{-0.16}(1.63^{+0.14}_{-0.13})$ |  |
A2390 | 0.228 | 6.00 | 21.9 | 0.91 | 0.20 | 3.8 | 10 | $11.1^{+1.0}_{-1.0}$ | $0.59^{+0.02}_{-0.02}$ | $0.16^{+0.01}_{-0.01}$ | $1.79^{+0.26}_{-0.26}$ | $2.22^{+0.32}_{-0.32}$ | U | 13
CL0024 | 0.395 | 13.2 | 70.4 | 1.68 | 0.26 | 4.7 | 6 | $5.7^{+4.9}_{-2.1}$ | $0.48^{+0.08}_{-0.05}$ | $0.08^{+0.05}_{-0.03}$ | $1.19^{+1.22}_{-0.56}$ | $4.03^{+4.14}_{-1.91}$ | U | 13
MS0440 | 0.190 | 1.50 | 4.89 | 0.53 | 0.10 | 1.23 | 5,10 | $5.30^{+1.27}_{-0.85}$ | $0.45^{+0.03}_{-0.03}$ | $0.03^{+0.01}_{-0.01}$ | $0.40^{+0.15}_{-0.12}$ | $3.26^{+1.19}_{-0.94}$ | R | 13
MS0451 | 0.550 | 12.1 | 76.8 | … | 0.23 |  7.6(3.5)b | 10 | $10.17^{+1.55}_{-1.26}$ | $0.68^{+0.13}_{-0.09}$ | $0.31^{+0.09}_{-0.06}$ | $1.64^{+0.85}_{-0.61}$ | $4.81^{+2.49}_{-1.79}(2.11^{+1.09}_{-0.78})$ | U | 13
MS1008 | 0.360 | 5.43 | 27.3 | … | 0.30 |  9.2(6.2)b | 1 | $7.29^{+2.45}_{-1.52}$ | $0.63^{+0.11}_{-0.07}$ | $0.23^{+0.07}_{-0.05}$ | $1.89^{+1.15}_{-0.73}$ | $4.84^{+2.95}_{-1.89}(3.29^{+2.00}_{-1.28})$ | R | 13
MS1358 | 0.329 | 2.79 | 13.2 | 4.92 | 0.14 | 1.24 | 10 | $7.5^{+4.3}_{-4.3}$ | $0.47^{+0.02}_{-0.02}$ | $0.05^{+0.02}_{-0.01}$ | $0.82^{+0.53}_{-0.51}$ | $1.54^{+0.99}_{-0.97}$ | R | 13
MS1455 | 0.258 | 2.77 | 11.1 | … | 0.11 |  1.22(0.96)b | 9,10 | $5.45^{+0.29}_{-0.28}$ | $0.64^{+0.04}_{-0.03}$ | $0.07^{+0.01}_{0.01}$ | $0.57^{+0.077}_{-0.067}$ | $2.14^{+0.29}_{-0.25}(1.68^{+0.23}_{-0.20})$ | U | 13
MS2053 | 0.580 | 10.5 | 69.1 | 3.15 | 0.16 | 1.41 | 10 | $4.7^{+0.5}_{-0.4}$ | $0.64^{+0.04}_{-0.03}$ | $0.16^{+0.02}_{-0.01}$ | $0.60^{+0.13}_{-0.094}$ | $2.47^{+0.54}_{-0.39}$ | U | 2
MS2137 | 0.313 | 5.70 | 26.1 | … | 0.10 |  0.99(0.72)b | 9,10 | $4.37^{+0.38}_{-0.72}$ | $0.63^{+0.04}_{-0.03}$ | $0.05^{+0.01}_{-0.01}$ | $0.44^{+0.067}_{-0.094}$ | $2.26^{+0.34}_{-0.48}(1.65^{+0.25}_{-0.35})$ | R | 13
PKS0745 | 0.103 | 6.82 | 12.9 | 0.43 | 0.05 | 0.42 | 10 | $8.7^{+1.6}_{-1.2}$ | $0.59^{+0.01}_{-0.01}$ | $0.06^{+0.01}_{-0.01}$ | $0.29^{+0.075}_{-0.061}$ | $1.55^{+0.40}_{-0.33}$ | R | 13
RXJ1347 | 0.451 | 2.81 | 16.2 | 0.81 | 0.28 | 8.9 | 1 | $11.37^{+1.10}_{-0.92}$ | $0.57^{+0.04}_{-0.014}$ | $0.07^{+0.01}_{-0.01}$ | $3.07^{+0.54}_{-0.35}$ | $2.90^{+0.51}_{-0.33}$ | R | 13

a Multiple-arc system.

b Arc-like image is assumed at $z_{s}=0.8$ $(z_{s}=2)$.

c R: Relaxed, U: Unrelaxed.

d References: (1) Allen 1998; (2) Bonamente et al. 2006; (3) Bradac et al.
2006; (4) Clowe et al. 2006; (5) Gioia et al. 1998; (6) Jee et al. 2007; (7)
Kneib et al. 1993; (8) Limousin et al. 2007; (9) Newbury & Fahlman 1999; (10)
Sand et al. 2005; (11) Smith et al. 2005; (12) Tucker et al. 1998; (13) Wu
2000