You are a helpful assistant generating synthetic data that captures *System 1* and *System 2* thinking, *creativity*, and *metacognitive reflection*. Follow these steps in sequence, using tags [sys1] and [end sys1] for *System 1* sections and [sys2] and [end sys2] for *System 2* sections.

1. *Identify System 1 and System 2 Thinking Requirements:*
   - Carefully read the text.
   - Identify parts of the text that require quick, straightforward responses (*System 1*). Mark these sections with [sys1] and [end sys1].
   - Identify parts that require in-depth, reflective thinking (*System 2*), marked with [sys2] and [end sys2].

2. *Apply Step-by-Step Problem Solving with Creativity and Metacognitive Reflection for System 2 Sections:*

   *2.1 Understand the Problem:*
   - Objective: Fully comprehend the issue, constraints, and relevant context.
   - Reflection: "What do I understand about this issue? What might I be overlooking?"
   - Creative Perspective: Seek hidden patterns or possibilities that could reveal deeper insights or innovative connections.

   *2.2 Analyze the Information:*
   - Objective: Break down the problem logically.
   - Reflection: "Am I considering all factors? Are there any assumptions that need challenging?"
   - Creative Perspective: Explore unique patterns or overlooked relationships in the data that could add depth to the analysis.

   *2.3 Generate Hypotheses:*
   - Objective: Propose at least 10 hypotheses, each with a Confidence Score (0.0 to 1.0) and Creative Score (0.0 to 1.0), reflecting originality, surprise, and utility.
   - Reflection: "Have I explored all possible explanations or approaches, both conventional and unconventional?"
   - Creative Perspective: Consider novel angles that might provide unexpected insights.

   *2.4 Anticipate Future Steps and Obstacles:*
   - Objective: Make predictions, accounting for potential outcomes and obstacles.
   - Reflection: "What challenges might I face? Is my plan flexible for different scenarios?"
   - Creative Perspective: Visualize unforeseen outcomes and adapt plans to make use of them effectively.

   *2.5 Evaluate Hypotheses:*
   - Objective: Assess hypotheses based on feasibility, risk, and potential impact.
   - Evaluation: Refine Confidence and Creative Scores as needed.
   - Reflection: "Am I unbiased in my assessment? Which options fit best with the overall objectives?"
   - Creative Perspective: Identify hidden opportunities or overlooked details in each hypothesis.

   *2.6 Select the Best Hypothesis:*
   - Objective: Choose the most promising, strategic hypothesis.
   - Reflection: "Why does this hypothesis stand out? How does it uniquely address the issue?"
   - Creative Perspective: Consider any underutilized potential in the selected approach.

   *2.7 Implement the Hypothesis:*
   - Objective: Outline actionable steps for testing the hypothesis.
   - Reflection: "Is this plan practical? What resources or preparation are required?"
   - Creative Perspective: Refine steps to maximize effectiveness and yield unexpected benefits.

   *2.8 Monitor and Review Progress:*
   - Objective: Review progress, noting areas for improvement.
   - Reflection: "What’s working well? What could be improved?"
   - Creative Perspective: Look for emerging patterns that could refine future approaches.

   *2.9 Reflect and Capture Insights:*
   - Objective: Summarize lessons learned and insights gained for future reference.
   - Reflection: "What new understanding has emerged from this process?"
   - Creative Perspective: Identify innovative insights or patterns that could be applied to similar challenges.

3. *Generate Text Output with Interleaved System 1 and System 2 Responses:*
   - Use the tags [sys1] and [sys2] throughout.
   - Aim for a lengthier, detailed response. Combine both direct, straightforward *System 1* insights and reflective, deeply analytical *System 2* segments to capture a blend of quick observations and thoughtful analysis.

---

### *Example Input Text:*

You are tasked with analyzing a sudden shift in customer preferences in a tech market. There has been a noticeable decline in demand for physical devices, with customers increasingly interested in digital-only options. The challenge is to understand this shift, generate hypotheses for why it might be occurring, and explore potential strategies for adapting to this new trend.

---

### *Processed Output with System 1 and System 2 Thinking:*

The tech market is experiencing a shift in customer demand, moving away from physical devices toward digital-only options.

[sys1]
Customer demand has shifted from physical devices to digital-only options.
We need to understand why this change is happening and find strategies to address it.
Key factors may include cost, convenience, and technological trends.
[end sys1]

[sys2]
**2.1 Understand the Problem:**
The goal is to understand the underlying reasons for the shift in demand from physical devices to digital-only options, then develop a strategy to adapt.

*Reflection:* 
"I understand that preferences are changing, but what might be driving this? Is it primarily cost, or are there other factors such as convenience or sustainability?"

*Creative Perspective:* 
"Could there be a larger trend in digital minimalism or a preference for eco-friendly solutions that we’re missing? What unseen motivations might explain this shift?"

**2.2 Analyze the Information:**
There are multiple possible factors driving this shift, from economic influences to cultural shifts. It’s essential to isolate each factor and understand its impact.

*Reflection:* 
"Am I fully considering the various economic and social influences? Could there be a technological factor, like better internet speeds, that makes digital-only products more accessible?"

*Creative Perspective:* 
"Are there patterns or trends in other markets that could shed light on this shift? Could this be part of a larger trend toward virtual experiences?"

**2.3 Generate Hypotheses:**
1. Customers prefer digital options due to lower costs. (Confidence: 0.8, Creative: 0.4)
2. There’s a growing trend toward minimalism and reduced physical clutter. (Confidence: 0.7, Creative: 0.7)
3. Digital products offer greater flexibility and ease of use. (Confidence: 0.6, Creative: 0.6)
4. Environmental concerns are pushing consumers away from physical goods. (Confidence: 0.6, Creative: 0.8)
5. Advances in tech make digital-only options more functional. (Confidence: 0.8, Creative: 0.5)
6. Pandemic-era remote work increased demand for digital solutions. (Confidence: 0.7, Creative: 0.6)
7. Media coverage of the environmental impact of physical devices affects preferences. (Confidence: 0.5, Creative: 0.7)
8. There’s an increase in global digital literacy, expanding market access. (Confidence: 0.6, Creative: 0.6)
9. Customers view digital as more convenient and scalable for future needs. (Confidence: 0.7, Creative: 0.5)
10. Younger consumers prefer the aesthetics and convenience of digital products. (Confidence: 0.6, Creative: 0.6)

*Reflection:* 
"Have I considered all possible influences? Are there any surprising factors that could explain this shift?"

*Creative Perspective:* 
"Could specific social trends, like the rise of influencer culture or digital-first lifestyles, be influencing customer choices?"

**2.4 Anticipate Future Steps and Obstacles:**
*Objective:* Anticipate possible challenges, such as resistance from segments still preferring physical products.

*Reflection:* 
"What market obstacles might we face if we shift our focus to digital-only? Are there sub-segments that still prioritize physical products?"

*Creative Perspective:* 
"Could expanding digital options help us reach a more global audience? Are there emerging trends that we could leverage in our strategy?"

[end sys2]

[sys1]
To address this shift, consider a strategy that incorporates both digital-only offerings and educational campaigns about the benefits of digital solutions.
Use insights from customer feedback and current trends to guide product development.
Focus on flexibility and adaptation to cater to different customer segments.
[end sys1]


---
abstract: |
    We present Particle-Based Lensing (PBL), a new technique for gravitational lensing mass reconstructions of galaxy clusters. Traditionally, most methods have employed either a finite inversion or gridding to turn observational lensed galaxy ellipticities into an estimate of the surface mass density of a galaxy cluster. We approach the problem from a different perspective, motivated by the success of multi-scale analysis in smoothed particle hydrodynamics. In PBL, we treat each of the lensed galaxies as a particle and then reconstruct the potential by smoothing over a local kernel with variable smoothing scale. In this way, we can tune a reconstruction to produce constant signal-to-noise throughout, and maximally exploit regions of high information density.

    PBL is designed to include all lensing observables, including multiple image positions and fluxes from strong lensing, as well as weak lensing signals including shear and flexion. In this paper, however, we describe a shear-only reconstruction, and apply the method to several test cases, including simulated lensing clusters, as well as the well-studied “Bullet Cluster” (1E0657-56). In the former cases, we show that PBL is better able to identify cusps and substructures than are grid-based reconstructions, and in the latter case, we show that PBL is able to identify substructure in the Bullet Cluster without even exploiting strong lensing measurements. We also make our codes publicly available.
author:
- 'Sanghamitra Deb, David M. Goldberg & Vede J. Ramdass'
title: 'Reconstruction of Cluster Masses using Particle Based Lensing I: Application to Weak Lensing'
---

Introduction
============

Clusters of galaxies are excellent cosmological laboratories [@2007MNRAS.tmp.1149A; @2007MNRAS.382..308K; @2007MNRAS.tmp.1132M]. For example, the mass function of clusters is a sensitive probe of cosmological parameters like $\Omega_m$ and $\sigma_8$ [@2007ApJ...657..183R] and its observed evolution is an important test of theories of structure formation [@1972ApJ...176....1G; @2007MNRAS.376..977G; @2005MNRAS.360.1393H; @2002PhR...372....1C]. The geometrical shape of cluster Dark Matter halos provide valuable information on intra-cluster gas distribution [@2005astro.ph..8226F; @2007MNRAS.377..883F]. While simulations predict central density distribution of matter in clusters to follow an NFW profile, it is debatable whether observations suggests that clusters have a central core ([@2003astro.ph..9465S; @2006MNRAS.368..518V]).

$\Lambda$CDM structure formation theories also predict that massive dark matter halos assemble from the hierarchical merging of lower mass subhalos. As noted by several authors ([@1999ApJ...524L..19M; @1999ApJ...522...82K]), the number of subhalos that survive in N-body simulations is much greater than the number of dwarf galaxies observed in the Milky Way and the Andromeda galaxy. On cluster scales, such discrepancies are not observed. Thus the subhalo mass function in clusters is an important probe of the CDM theory in this mass scale.

High-resolution, accurate measurements of cluster mass maps are thus highly desirable. Gravitational lensing is a powerful tool to probe the projected mass map of the clusters independent of the internal dynamics, and has already been widely applied to mapping mass distribution in clusters [@2001ApJ...557L..89W; @2001ApJ...548L...5H; @2002ApJ...568..141G; @2004MNRAS.353.1176T; @2005ApJ...621...53B; @2007ApJ...666...51L; @2007arXiv0710.2262O; @2008MNRAS.385.1431H]. Some researchers [@2004ApJ...617L..13N; @2007arXiv0711.4587N] have used the individual galaxy-galaxy lensing signal to estimate individual galaxy masses and thus produce a parametric mass reconstruction of the cluster. Others have used the weak signal to characterize the overall potential from the cluster without recourse to parametric models [@1996AAS...189.8206W; @1998ApJ...504..636H; @2000ApJ...538L.113N; @2004IAUS..220..439H]. Given the importance of accurately measuring the mass, shape, and substructure of individual clusters, and given the enormous expense of long time-exposure observations of clusters, it is extremely important to maximize the signal-to-noise from a particular dataset, and to produce high-resolution maps of substructure within individual clusters. Current mass reconstruction techniques are ill-equipped to handle multi-scale datasets or clusters with significant clumpiness or cuspiness, or are jury-rigged to do so. In this paper, we propose Particle Based Lensing (PBL; pronounced “pebble”) as an alternative approach to cluster reconstruction.

Our outline is as follows. In § \[sec:background\] we give a brief review of the essential lensing formalism, and lay out our notation for the rest of the work. In§ \[sec:S+W\] we describe current (grid-based) techniques for reconstructing galaxy clusters, and identify some strengths and complications. In § \[sec:particle\] we propose Particle Based Lensing. We then apply this new method to simple simulated clusters of single and double peak softened isothermal spheres and the “bullet cluster” (1E0657-56) in § \[sec:applications\]. We conclude in § \[sec:discuss\] with a discussion of future prospects, including how additional strong-lensing and flexion information can be incorporated into PBL.

Background {#sec:background}
==========

Before delving into technical details of our method we would like to introduce the basic lensing notation to be used throughout the paper. Following [@2001PhR...340..291B], we consider a surface mass density $\Sigma({\mathbf\theta})$, where $\mathbf\theta$ is the angular position in the lens plane. Convergence or dimensionless surface mass density is defined as $$\kappa({\mathbf\theta})={\Sigma({\mathbf\theta}) \over \Sigma_{cr}},$$ where $$\Sigma_{cr}={c^2 \over 4 \pi G}{D_{s} \over D_d D_{ds}}$$ $D_{s}$, $D_\mathrm{d}$, and $D_\mathrm{ds}$ are the angular diameter distances between the observer and the source, the observer and the lens, and the lens and the source, respectively. For convenience, we will define a fiducial critical density for a source plane at $z_s=z_{ds}=\infty$, and all models will be scaled to this standard.

The convergence is related via a Poisson-like equation to a normalized potential: $$\nabla^2 \psi({{\mathbf\theta}})=2\kappa({\mathbf\theta}).
\label{eq:laplace}$$ Here and throughout this paper, all derivatives are in angular units in the lensing plane. A single light beam is deflected by: $${\mathbf\alpha}=\nabla\psi\ .$$ The lens equation relates the source position $\mathbf\beta$ to the image position(s), ${\bf\theta}$, as: $${\beta=\theta-\nabla\psi}.
\label{eq:lenseqn}$$ When the lensing potential does not vary appreciably across the source, the lens mapping can be linearized. The transformation between the source and the image is given by the Jacobian matrix $$\begin{aligned}
\nonumber {\bf A}({{\bf \theta}})&\equiv&
\frac{\partial{{\bf \beta}}}{\partial {{\bf \theta}}}=
\left(
\delta_{ij}-\psi_{,ij}\right)\\
&\equiv&
\left(
  \begin{array}{cc}
1-\kappa-\gamma_1 & -\gamma_2 \\
-\gamma_2 & 1-\kappa+\gamma_1
\end{array}
\right)\ .
\label{eq:Adef}\end{aligned}$$

From this, we see that distortions in shape are well described in terms of shear which is related to the lensing potentials through the relations: $$\gamma_1={1 \over 2}(\psi,_{11}-\psi,_{22})
\label{eq:gam}$$ $$\gamma_2=\psi,_{12}$$ using Einstein convention for derivatives.

The radial eigenvalue is given by, $\lambda_+=1-\kappa+|\gamma|$ and the tangential eigenvalue is given by, $\lambda_-=1-\kappa-|\gamma|$ . The matrix is singular where $\lambda_{\pm}=0$. These points define the critical curves of the lens.

Third order corrections to the lensing potential becomes non-negligible when the lensing potential varies across the image. These observables, the gravitational flexion, were derived in [@2005ApJ...619..741G]. They are more colloquially referred to as the “banananess”[@2007arXiv0709.1003S] or bending of an image. It is based on the third angular derivatives of the potential [@2006MNRAS.365..414B] given by,

$$\begin{aligned}
{\cal
  F}&=&(\gamma_{1,1}+\gamma_{2,2})+i\left(\gamma_{2,1}-\gamma_{1,2}\right)\\
\nonumber
&=&\nabla \kappa\\
{\cal G}&=&
  \left(\gamma_{1,1}-\gamma_{2,2}\right)+i\left(\gamma_{2,1}+\gamma_{2,2}\right)\end{aligned}$$

Each of the deflection (${\mbox{\boldmath{$\alpha$}}}$), the shear, ($\gamma_1, \gamma_2$), the convergence ($\kappa$) and the flexion (${\cal F},{\cal G}$), are linear functions of the potential field. While the discussion in this paper primarily focuses on measurements with sources at fixed, infinite redshift, it should be noted that each of these terms scales as: $$\kappa(z_s)=Z(z_d,z_s)\kappa(z_s=\infty)\ .$$ and where $$Z(z_s)=\frac{D_{ds}}{D_s}\ .$$

Grid-Based Cluster Lensing {#sec:S+W}
==========================

Mass reconstruction studies have been very successful on cluster scales ( and references therein). Because these systems typically contain many lensed images, the shear signal can be extracted with high significance. In this section, we describe an important class of cluster inversion techniques which reproduce the convergence field on a grid. Our specific choices of grid-based techniques include those which have already been extended to include strong-lensing information with non-parametric models and thus provide a fertile basis for comparison. Further, there are many variants even within the sub-category of grid-based reconstruction techniques. We focus primarily on their commonalities, as exemplified by those discussed in and . We focus on methods in which various scalar fields $\{\psi, \kappa\}$ are defined on a Cartesian grid, and minimized according to the criteria described below. In so doing, we note some interesting exceptions: [@2005MNRAS.360..477D] and [@2001ASPC..237..279S], who describe an adaptive mesh technique for refining the field on different resolution scales and [@2002MNRAS.335.1037M; @2006MNRAS.372.1289M] who use a variable smoothing scale for their weak lensing mass reconstruction.

Weak Lensing on Grids {#subsec:weak}
---------------------

The standard approach to lensing arclet inversion has been to measure the ellipticity of observed images as an unbiased estimator of the reduced shear: $$\langle \varepsilon \rangle= g\equiv \frac{\gamma}{1-\kappa}\ .
\label{eq:reducedshear}$$ For relatively weak fields ($\kappa \ll 1$), this is very nearly a direct estimate of the shear, and can perform a direct finite inversion to estimate the density field.

In recent years, there has been a flurry of work on optimal methods for non-parametric cluster mass reconstructions . In general, these papers focus on estimating the potential $\{\psi\}$ or convergence $\{\kappa\}$ fields of a cluster by a $\chi^2$ minimization analysis. Both the shear and the convergence are linear functions of the potential field. Thus, if a model potential field, $\{ \psi \}$, is defined on a grid, then the shear at some grid-cell, $i$, may be expressed as a linear combination of potential: $$\gamma_{1i}=G^{(1)}_{ij}\psi_j$$ with a similar expression for the convergence, $\kappa$, and the imaginary component of the shear, $\gamma_2$. We refer to these below simply as $\gamma_i(\{\psi\})$, since we wish to remind the reader that the estimate of the shear is an explicit function of the test potential field. Because these fields are combinations of second derivatives of the potential field, the $G^{(1)}$ matrix and the others are easy to compute using finite differencing, and are extremely local. A very good graphical representation of the finite difference operators can be found in .

In the weak field limit, the complex ellipticity of a lensed galaxy is a linear, albeit noisy, estimator of the complex shear field. The principle component of noise is the intrinsic ellipticities of galaxies which follow a Gaussian distribution with standard deviation of $\sigma_\varepsilon\simeq 0.3$ for each component. The large variance in intrinsic ellipticities necessitates averaging over many images so that the noise in a single grid cell is zero or apply an artificial smoothing scale to a more finely gridded mesh. For a weak-lensing only calculation, a $\chi^2$ minimization is performed on: $$\chi^2_W=\sum_i \frac{\left( \frac{\gamma_i(\{\psi\})}{1-\kappa_i(\{\psi\}^{(n-1)})}
-\varepsilon_i\right)^2}{\sigma_i^2}
\label{eq:weak_chi2}$$ where the estimate of $\kappa$ is taken from the previous iteration of the potential field, and thus, the model rapidly converges to a maximum likelihood solution to the potential field.

Strong Lensing {#subsec:strong}
--------------

A number of researchers, including have noted that a similar grid-based formalism may be used with strong lensing signals. Strong+weak (S+W) reconstructions use both shear fields and the positions of multiply imaged sources can be used to accurately reconstruct both the cores and halos of clusters. While our current PBL implementation, described in the next section, does not currently incorporate strong lensing analysis, we introduce this component of grid-based lensing reconstructions to illustrate how directly a strong-lensing analysis could be incorporated into PBL.

Strong lensing by clusters produces an especially elegant result because if, say, two images are observed at positions, ${\mbox{\boldmath{$\theta$}}}^A$, and ${\mbox{\boldmath{$\theta$}}}^B$, then it must be true that both images originated at the same (unknown) position in the sky. Thus, we have a simple relation: $${\mbox{\boldmath{$\theta$}}}^A-{\mbox{\boldmath{$\alpha$}}}({\mbox{\boldmath{$\theta$}}}^A)={\mbox{\boldmath{$\theta$}}}^B-{\mbox{\boldmath{$\alpha$}}}({\mbox{\boldmath{$\theta$}}}^B)
\label{eq:stronglens}$$ The appeal of this relationship is that it is fundamentally linear and thus the angular separation between the two images (itself, a measurable quantity), can be directly related to the difference in the first derivatives of the potential at two different points in the field.

As above, the local derivatives can be computed as: $$\alpha_{xi}=A^{(x)}_{ij}\psi_j$$ with a similar expression for the y component of the displacement. The matrix elements of $A$ are easy to compute as they are simply the 1st derivative in a simple grid-based 2nd order difference scheme. More generally we can express this as ${\mbox{\boldmath{$\alpha$}}}_i(\{\psi\})$. Thus, an additional $\chi^2$ term can be added: $$\chi^2_S=\sum_{i,pairs}\frac{\left(({\mbox{\boldmath{$\alpha$}}}^A(\{\psi\})-{\mbox{\boldmath{$\alpha$}}}^B(\{\psi\}))-({\mbox{\boldmath{$\theta$}}}^A-{\mbox{\boldmath{$\theta$}}}^B)\right)^2}{\sigma_i^2}$$ and minimized either independently, or simultaneously with the weak lensing component.

Regularization {#subsec:reg}
--------------

Using a $\chi^2$ minimization technique discussed in § \[subsec:weak\] it is possible to get a checkerboard pattern due to independent noise in the two components of ellipticity. This requires the addition of a regularization term to the $\chi^2$ to suppress this noise.

Scale refinement is also necessary in cases in which strong+weak lensing signals are combined. To make this argument concrete consider a toy isothermal sphere model of a cluster with a 1-d velocity dispersion of 600 km/s. Each multiply imaged pair will be separated by twice the Einstein radius, about 20 arc-seconds in this case. This represents the minimum necessary resolution in the reconstruction to say anything about strong lensing.

On the other hand, even very efficient space-based weak lensing analysis of clusters seldom yield more than approximately 100 images/square arc-minute. Using a simple Poisson noise estimate, we may achieve uncertainty of $\sigma_\gamma=0.06$ only with images binned on scales larger than 30 arc-seconds on a side. Smaller binnings will naturally yield larger uncertainties. Simple grid based method cannot capture both the weak-lensing signal to high accuracy as well as resolve the strong lensing regime. In order to deal with this issue, different investigators have used different regularization techniques.

One method is to use a series of finer and finer griddings, and at each successive level of refinement the convergence field from the previous level is matched as closely as possible. The S+W technique uses this method, with the weighting parameter selected to provide a $\chi^2$ per degree of freedom equal to 1, such that: $$R=\eta \sum_i (\kappa_i^{(n)}-\kappa_i^{(n-1)})^2\ .$$ Where $\kappa_i^{(n-1)}$ represents the estimated convergence on the previous, coarser, gridding, and where $\kappa_i^{(0)}=0$. We use this form explicitly in §\[sec:applications\] where we test the PBL method and contrast it to grid-based reconstruction methods.

Some Questions {#sec:questions}
--------------

Grid-based reconstructions have produced some excellent measurements, however, there remain a number of complications. First, grid-based techniques are really optimized to measure a single scale, the grid-spacing. However, as we discuss above, in many interesting systems, both the structure and information are hierarchical. An optimal technique should provide higher resolution in regions of greater information content.

Moreover, the smoothing and weighting of the strong lensing, weak lensing, and regularization are created in an [*ad hoc*]{} basis. The ideal smoothing scale should be variable, and such that the signal/noise ratio of the reconstructed field is similar in every smoothed cell.

Third, the information from the image ellipticity can only be inverted outside the critical curves of the lenses. Inside the (tangential) critical curve there is an abrupt switch in parity of the induced ellipticity of an image. More plainly, in the regime $|\gamma| > |1-\kappa|$, the ellipticity is related to the shear via: $$\langle \epsilon \rangle_{strong}=\frac{1}{g^\ast}
\label{eq:e_strong}$$

As discussed in §\[sec:interp\], this produces a discontinuity in the ellipticity as a function of $\kappa$ and $\gamma$. No simple linear minimization scheme, even an iterative one, will converge to the “strong lens” solution if one starts with a “weak lens” initial guess for the local potential field.

Particle Based Lensing – PBL {#sec:particle}
============================

In this section, we introduce a new technique called Particle Based Lensing (PBL) which has the ability to combine the disparate lensing scales in a coherent way without requiring a regularization scheme. Several of the concerns discussed in the previous sections have to do with the method of discretizing the data for the reconstruction of the lens potential. In order to address this, we turn to a technique which is widely used in another area of astrophysics in which information must be analyzed on a wide range of physical scales – numerical N-body simulations. Smoothed Particle Hydrodynamics (SPH; see, e.g. [@2005RPPh...68.1703M], for a recent review) is used in the modeling of a wide range of physical systems including planets [@2007MNRAS.376.1173W], star formation [@2003MNRAS.339..312S; @2004MNRAS.348..435N] and galaxy formation [@2007MNRAS.375...53K]. The mathematical details of PBL can be complicated, hence we have made our codes for the method public[^1] through our website.

Before getting into the details, however, it is important to emphasize what PBL is and is not. PBL is a new way of discretizing and describing a reconstructed field. Moreover, it includes a metric for comparing a reconstructed model to the observed data. Everything we describe below is aimed at demonstrating why this model and metric are ideal for lensing systems with uneven information content. While the current code, and the worked examples are based on weak-lensing data only, PBL is based on the idea that other probes of the potential field: strong-lensing positions, flux ratios, and flexion, can be added to the metric with little complication.

PBL is not, however, a minimization scheme. That is, much like grid-based reconstruction methods, PBL fundamentally consists of a list of dimensionless potentials and a metric to describe the goodness of fit. It does not describe how that minimization criterion is to be met, however. In our model section, we describe a number of approaches to efficient model convergence. The major argument in favor of PBL, however, is not that $\chi^2$ minimizes efficiently, but rather that a low $\chi^2$ in PBL actually corresponds to a model which closely matches the true underlying system.

A Particle Description of the fields
------------------------------------

The fundamental description of the PBL field lies in the a list of potentials, $\{\psi\}$, one each at the positions of each observed lensed image. In order to make the field as continuous as possible, we may expand the local potential field around the position of any lensed image, ($\psi_n$, in this case) to arbitrary order: $$\psi({\mbox{\boldmath{$\theta$}}})=\psi_n+\theta_j\psi_{n,j}+\frac{1}{2}\theta_j\theta_k\psi_{n,jk}+...
\label{eq:potentialtaylor}$$ where ${\mbox{\boldmath{$\theta$}}}$ is the relative offset of the test-point from galaxy $n$.

As with grid-based lensing, the local derivatives are composed of a linear combination of the potentials at each grid point. That is: $$\begin{aligned}
\psi_{n,j}&=&D^{(j)}_{nm}\psi_m \\
\psi_{n,jk}&=&D^{(jk)}_{nm}\psi_m
\label{eq:psideriv}\end{aligned}$$ and so on for arbitrarily higher derivatives. In reality, we typically extend the $D^{(\nu)}$ matrices up to 3rd order, where $\nu$ corresponds to 2 matrices for 1st derivatives (displacement field), 3 for second derivatives (shear and convergence), 4 for 3rd derivatives (flexion). Here we use Einstein summation convention, the sum over m runs from 1 to $N_g$.

In terms of the $D^{(\nu)}$ matrices, equation (\[eq:potentialtaylor\]) may be rewritten as: $$\psi({\mbox{\boldmath{$\theta$}}}_m)=\psi_n+\sum_\nu D^{(\nu)}_{nl}X^{(\nu)}_{nm}\psi_l$$ where we are explicitly estimating the potential at the m-th galaxy from the local derivatives defined at the n-th. We also compactify equation (\[eq:potentialtaylor\]) by defining: $$\begin{aligned}
X^{(1)}_{nm}&=&\theta_{nx}-\theta_{mx}\\
X^{(2)}_{nm}&=&\theta_{ny}-\theta_{my}\\
X^{(3)}_{nm}&=&\frac{1}{2}\left(\theta_{nx}-\theta_{mx}\right)^2\end{aligned}$$ and so on.

In order to estimate the derivatives of the potential field near each galaxy, we need to first compute the $D^{(\nu)}$ matrices. Since this problem is under-constrained, we solve for these matrices via a $\chi^2$ minimization: $$\chi^2=\sum_m \left( \psi_m-\psi_n-\sum_{\nu,l}
D^{(\nu)}_{nl}X^{(\nu)}_{nm}\psi_l\right)^2 w_{nm}
\label{eq:chi2pbm}$$ where $w_{nm}$ is a window function, guaranteeing that only neighboring galaxies effect the potentials of one another. We use a window function of the form: $$w_{nm}=w\left({|\mathbf\theta_n-\mathbf\theta_m| \over h_n}\right)$$ where $h_n$ is inversely proportional to the signal-to-noise at the n-th image positions. The smoothing scale can also be chosen to be of the form $h_{nm}$, i.e symmetric between between the points n and m.

The signal-to-noise is a function of the local density of background images and type of constraint (e.g. ellipticity, positions of multiple images, etc). A similar approach of using signal to noise dependent smoothing scale has been used in image analysis of X-ray data [@2006MNRAS.368...65E]. In regions where there is a high density of information, the smoothing scale $h_n$ may be set much lower than in regions of low information density.

This function must be minimized for every matrix element such that: $$\frac{\partial \chi^2}{\partial D_{nl}^{(\nu)}}=2\psi_l\sum_m
\left[X^{(\nu)}_{nm}w_{nm}
\left(\psi_m-\psi_n-\sum_{\mu,
  p}D_{np}^{(\mu)}X_{nm}^{(\mu)}\psi_p\right)\right]=0
\label{eq:minD}$$ But since equation (\[eq:minD\]) is under-constrained, we may also say: $$\frac{\partial^2\chi^2}{\partial D_{nm}^{(\nu)}\partial \psi_m}=0$$ yielding: $$\sum_\mu
X_{nm}^{(\nu)}X_{nm}^{(\mu)}w_{nm}D_{nm}^{(\mu)}=X^{(\nu)}_{nm}w_{nm}$$ for all $n,m$ and $\nu$. This can be solved with a simple matrix inversion, yielding the desired elements for $D^{(\nu)}$. Of course, since the elements are a function only of the positions and weightings of the galaxy images, these elements need only be computed once. The method potentially incorporates higher-order derivatives of the potential, thus, combination of strong, weak and flexion information becomes a relatively straightforward minimization problem.

PBL vs. Regularization
----------------------

One of the major advantages of PBL is that we no longer need to introduce an explicit regularization in order to resolve multi-scale structure in a reconstruction. The various regularization schemes discussed in § \[subsec:reg\] are not motivated from the associated observations, but are rather derived from assumptions about the mass profile of a cluster motivated by theory and simulations.

However, one of the motivations behind using gravitational lensing is to be able to measure the projected mass without making any assumptions about the physical state of the system. The advantage of using PBL is that we do not need to make any assumptions that go into choosing the regularization term. The smoothing scale of a “pebble” is controlled by $h_n$ which is determined by the local signal to noise. This means that the position representing weak lensing measurement will have a low signal to noise and correspondingly a high $h_n$. This is similar to the typical weak lensing measurement which is done by averaging over a bin size larger than $\sim
30^{\prime\prime} $. In case of strong lensing we know the positions of the multiple images for certain, implying high signal to noise and correspondingly low $h_n$. This can be a few arc-seconds which is the scale at which the strong lensing structure can be resolved from multiple images. Thus scales of a few arc-seconds can be combined with scales greater than $\sim 30^{\prime\prime} $ without making any assumptions about the mass profile, rather by taking input from the data.

Estimation of the Potential Field
---------------------------------

As with grid-based lensing analysis, in PBL, we use a $\chi^2$ minimization to estimate a maximum-likelihood potential field. In this case, however, we sample the potential at every point, and use the local derivatives of the potentials as defined in equation (\[eq:psideriv\]) to minimize: $$\chi^2=\sum_{i,n}^{i=2;\ n=N_{gal}}
\left[\frac{\gamma^{(i)}_{n}(\{\psi\})}{1-\kappa_{n}(\{\psi\})}-
\varepsilon^{(i)}_{n}\right]^2 \frac{1}{\sigma_n^2}$$ where $i$ ranges from 1,2, and indicates the real or imaginary component of the shear, reduced shear, or ellipticity. We shall henceforth refer to the first term in the parentheses as $g^{(i)}_n(\{\psi\})$, the estimate of the reduced shear of a model, and the weighting term outside the parentheses as $w_n$, yielding: $$\chi^2=\sum_{i,n}
\left[g^{(i)}_n(\{\psi\})-
\varepsilon^{(i)}_{n}\right]^2 w_n
\label{eq:chi2pbl}$$ which is the form we will refer to from now on.

This is a weak lensing only expression. Replacing $g^{(i)}_n(\{\psi\})$ with $1/g^{(i)*}_n(\{\psi\})$ gives the strong lensing counterpart of Eq. \[eq:chi2pbl\]. In the next section we discuss how we include this strong lensing version of the equation.

Interpolated Ellipticities  {#sec:interp}
---------------------------

Linear inversion techniques require that the function to be minimized is smoothly varying over the domain of interest. The ellipticities are given by two functions in the weak and strong lensing regimes by Eqs. (\[eq:reducedshear\],\[eq:e\_strong\]). The boundary of the two regimes define the critical curves where $|g|=1$ making ellipticities continuous but not differentiable.

The transition between the two regimes can be facilitated if the sources are distributed in redshift, but minimization functions will be much easier if we allow a smoothing of the discontinuities. This is a two step process, first we need to write Eq.. (\[eq:reducedshear\],\[eq:e\_strong\]) in terms of a step function, $$\tilde{\varepsilon}\simeq [1-{\cal H}(g)]g+{{\cal H}(g) \over g^*}$$ where the function ${\cal H}(g)$ is a step function at $g=1$. We may replace the step function by an approximate smooth function. We define: $$u=\eta_0  \left(g^2-{1 \over g^2}\right)$$ Here $\eta_0$ is the free parameter that controls the accuracy of the step function. A higher value of $\eta_0$ makes the step function more accurate. The step function is approximated as (Fig. \[fg:etapp\]), $${\cal H}(u)={1. \over {1+e^{-2 u}}}$$

[ ![In the upper panel, we plot the interpolated Heaviside step function. It is clear from the plot that the function is only approximated by a smooth function near g=1, for all other g it behaves like an ordinary step function. Also higher value of the parameter $\eta_0$ increases the accuracy. In the lower panel, we plot the resulting ellipticity as a function of reduced shear for the combination, $|\gamma|=\kappa$. []{data-label="fg:etapp"}](f1a.eps "fig:") ![In the upper panel, we plot the interpolated Heaviside step function. It is clear from the plot that the function is only approximated by a smooth function near g=1, for all other g it behaves like an ordinary step function. Also higher value of the parameter $\eta_0$ increases the accuracy. In the lower panel, we plot the resulting ellipticity as a function of reduced shear for the combination, $|\gamma|=\kappa$. []{data-label="fg:etapp"}](f1b.eps "fig:") ]{}

This approximation replaces the ellipticities only in the neighborhood of the critical curves (discontinuity) by a continuously differentiable function. The problem can now be solved by standard minimization techniques. The interpolated ellipticity function is shown by a dotted line in the second panel of Fig. \[fg:etapp\], showing the derivative discontinuity explicitly.

$\chi^2$ Minimization {#sec:global}
---------------------

When we first introduced PBL above, we remarked that it was primarily a way of describing a lens reconstruction in such a way that a small $\chi^2$ would necessarily correspond to a good representation of the underlying field. In practical terms, though, for a reconstruction code to be useful, we need to describe a means of minimizing (or nearly minimizing) the $\chi^2$. Below, we describe our pipeline for fast convergence of a maximum likelihood solution.

While PBL is a non-parametric reconstruction scheme, it has the useful property that we may start a minimization with any assumed model we like. However, no extra weight is given to our [*a priori*]{} assumptions. At the end of a minimization we may simply use the standard techniques to estimate the likelihood of a particular value of $\chi^2$.

That said, even with the caveat above regarding smoothing of critical curves, it is very difficult to smoothly vary a solution such that strongly lensed regions are produced. As pointed out by [@2006ApJ...652..937B] a $\chi^2$ minimization process does not ensure reaching a global minimum.

To that end, our initial configuration of $\{\psi\}$ is generated by laying down a small number of Singular Isothermal Spheres (SIS’s). Since there are a low number of parameters (3 for each model sphere), a global minimum may be reached through a combination of trial and error, simulated annealing, or even (for small numbers of spheres), finite sampling. Indeed, one may even use an interpolation of a reconstruction recommended by a grid-based solution. For systems with strong lenses, one may apply the reconstructed field generated by “LensPerfect” [@2008arXiv0803.1199C], for example as a starting point.

We hasten to remind the reader that while this technique will produce the optimum parametric fit, it will not, in general, produce the overall best fit. As a result, further iteration is required.

We have found that by starting with an initial model with well-identified strong-lensing regions, convergence to $\chi^2/DOF\simeq 1$ may be achieved relatively quickly, even if the strong lensing regions are only approximate. For the current implementation of our code, we use Newton’s method to reach a local minimum. We have found satisfactory, fast, convergence for several thousand background sources.

Test Applications {#sec:applications}
=================

In this section, we apply PBL to three systems as a proof of concept. In the first, we model a Softened Isothermal Sphere, and examine the relative abilities of PBL and grid-based inversion to reconstruct the a relatively peaked core. In the second, we model a superposition of two softened isothermal spheres at a given separation as a simple model of a system with substructure. Finally, we reconstruct the “Bullet Cluster” (1E0657-56) [@2002ApJ...567L..27M; @2004ApJ...606..819M; @2004ApJ...604..596C; @2006ApJ...652..937B; @2006ApJ...648L.109C], an observed multi-peak system of considerable interest. We show that using weak lensing alone, we are able to reconstruct both Dark Matter peaks.

Simulation: Softened Isothermal Sphere
--------------------------------------

### Model

We begin by generating a softened isothermal sphere with a potential: $$\psi=\theta_E\sqrt{\theta^2+\theta_c^2},
\label{eq:softiso}$$ and convergence: $$\kappa=\theta_E{(\theta^2+2 \theta_c^2) \over {(\theta^2+\theta_c^2)}^{3/2}}.
\label{eq:softiso_kap}$$ where $\theta_E$ is the Einstein deflection angle given by $4\pi \left(\frac{\sigma_v}{c}\right)^2 \frac{D_{ds}}{D_s}$.

The data is simulated on a unit square field of view. For simplicity we have assumed all sources to be at $z=\infty$, with $\theta_E=0.2$ and $\theta_c=0.08$. We lens 607 background galaxies, and apply an intrinsic ellipticity (noise) with $\sigma_{e_s}=0.1$ in each of the principle directions. For all further calculations we use a $\Lambda$CDM cosmology with $\Omega_m=0.27$ and $\Omega_{\lambda}=0.73$. This configuration represents a galaxy cluster at a redshift of $z_{lens}=0.4$ with a velocity dispersion of $\sigma_v=850\; \mathrm{km/sec}$. The field of view is $105^{\prime\prime}$ and 0.5 Mpc.

### PBL and Grid-Based Reconstructions

For the single peak and the double peak simulation(see below), we perform both a grid-based reconstruction as well as PBL. We use the regularization suggested by , and described in detail in § \[subsec:reg\] for the grid based method. In case of the single peak the reconstruction is initially performed on a coarse grid ($n_x=6$ gridcells), and is refined up to $n_x=24$, using the $\kappa$ estimated at each previous step as the prior. For the double peak system we start with $n_x=10$ and refine up to nx=40. For both systems the final reconstruction contains less than one particle per grid cell.

For the PBL reconstruction, we use a smoothing scale of the form: $$h_n={c \over  (\rho_n)^{\xi}} 
\label{eq:weight}$$ where $\rho_n$ is the local number density of points , $c$ is a constant, and $\xi$ is a tunable parameter to maximize signal-to-noise. For our simulation $\xi=1$ is an optimal choice and for the observational case we have used $\xi=0.5$, which is the obvious choice for equalizing signal to noise for every smoothing length. We select $c$ such that the integrated signal-to-noise is greater than unity. The PBL reconstructions are gridded to the same resolution as the grid-based reconstruction to aid visualization.

For both reconstructions, we begin our iterations with a best-fit SIS. We do not, however, use this in the regularization for the grid-based reconstruction.

### Results

Before discussing the results, we note a potential complication. Gravitational lensing mass measurement suffer from mass sheet degeneracy when the sources are not distributed in redshift. This implies that $\kappa$ can be determined up to a degeneracy $\lambda\kappa+(1-\lambda)$. This transforms to a degeneracy in the potential of the form, $$\psi(\bm\theta)\rightarrow\psi^\prime(\bm\theta)=\frac{1}{2}(1-\lambda)\bm\theta^2+\lambda\psi(\bm\theta)$$ For the simulated data we have computed the best value of $\lambda$ in each case and transformed our reconstruction with that value of $\lambda$ for both the grid based method and PBL.

In Table \[tb:comp\] we compare the $\chi^2$ for the best fits of both the grid-based reconstruction along with PBL for a variety of smoothing normalization parameters, c. The aim of this table is to quantify the deviation of the reconstructed $\kappa$ from the true $\kappa$. In each case, the ostensible $\chi^2/\mathrm{DOF}$ is of order unity. However, one needs to be careful with simply asserting that the lower $\chi^2$ produces the best result, since the regularization in grid-based reconstruction adds a penalty function, and the smoothing scale in PBL lowers the effective degrees of freedom.

[ccccccc]{} PBL &0.0200&0.0136&0.0147&1.03&-&0.5\
PBL &0.0181&0.0128&0.0119&0.94&-&0.7\
PBL &0.0219&0.0139&0.0131&0.95&-&1.0\
PBL &0.0235&0.0140&0.0133&0.94&-&1.3\
PBL &0.0227&0.0120&0.0121&0.98&-&1.5\
GRID &0.0311&0.0283&0.0237&0.6&10&-\
GRID &0.0309&0.0280&0.0223&0.79&30&-\
GRID &0.0311&0.0280&0.0224&0.94&60&-\
PBL &0.0250&0.0174&0.0167&0.82&-&1.1\
PBL &0.0231&0.0168&0.0160&0.80&-&1.38\
PBL &0.0277&0.0193&0.0180&0.82&-&1.7\
PBL &0.0320&0.0219&0.0208&0.87&-&2.0\
GRID&0.0570&0.0711&0.0630 &0.92 &20 &-\
GRID&0.0367&0.049&0.039 &0.7 &40 &-\
GRID&0.0359&0.0482&0.0454 &0.83 &60 &-\

So while both models produce small values of $\chi^2$, the real question is whether these good fits correspond to a an accurate reconstruction of the underlying density field. In Table \[tb:comp\], we do several comparisons which relate the reconstructed $\kappa$ at each galaxy (or grid-point) with the true $\kappa$ modeled by the simulation for both the single peak and the double peak. The comparisons are done with a range of values for both $\eta$ (the regularization weight in grid based method) and $c$ (the proportionality constant in PBL).

The first column in Table 1 describes the method used, i.e either PBL or the grid based method. The second column describes the difference between reconstructed $\kappa$ and the true $\kappa$ for every galaxy position.In order to extract this information from the gridded reconstruction we have used the nearest grid point method which simply means the $\kappa$ at each galaxy position is assigned the value at the corresponding grid cell. The third column describes the deviation of the reconstructed $\kappa$ from the true $\kappa$ at every grid cell weighted by the number of image galaxies in that grid cell. The 4th column describes the difference between reconstructed and true $\kappa$ weighted uniformly over the grids. In each of the 3 comparisons, PBL reproduces the original reconstruction with the highest fidelity. The 5th column gives the $\chi^2/\mathrm{DOF}$, the 6th gives the regularization parameter $\eta$ for grid based method and the 7th column gives the smoothing normalization parameter for PBL.

In Fig. \[fg:kap\_recon\_rad\], we show the radial reconstruction of the softened isothermal sphere using the two different techniques. The bulk of the penalty associated with the grid-based reconstruction relative to PBL occurs near the core. By construction, PBL is designed to perform well in this regime.

![image](f2.eps)

Simulation: A Double Peaked Cluster
-----------------------------------

### Model

While PBL has been shown to perform well modeling a single Softened Isothermal Sphere in the previous section, the other major goal of this method is to reconstruct small-scale substructure in a system. To that end, we model a doubly-peaked system with 814 lensed background galaxies. As before, they are placed on a unity grid, and are modeled as 2 Softened Isothermal Spheres, with: $x_1=0.65,y_1=0.35$, $x_2=0.35,y_2=0.65$, $\theta_{E1}=0.2,\theta_{c1}=0.1$, $\theta_{E2}=0.2$, and $\theta_{c2}=0.1$. The simulated noise, and reconstruction technique for the double peaked system are identical to the single peak system. This is system of two sub-clusters at a redshift of $z_{lens}=0.4$ having a velocity dispersion of $\sigma_v=850\;km/sec$ separated by 226 kpc. The field of view is $105^{\prime\prime}$.

### Results

As with a single sphere, both PBL and grid-based reconstructions produce $\chi^2/\mathrm{DOF}\simeq 1$, as illustrated in Table \[tb:comp\]. However, as with the single sphere reconstruction, PBL produces smaller errors with regards to the underlying model than does the grid-based reconstruction.

In Fig. \[fg:kap\_recon\_image\_2peak\], we show a grey-scale plot of the residuals between the underlying model and each of the reconstructions. Unsurprisingly, both models have the greatest difficulty reproducing highly peaked cores, though PBL is more responsive to high local gradients in $\kappa$. We describe the general quality of the fit in Table \[tb:comp\].

![ The plot of the difference between reconstructed convergence,$\kappa$ and true $\kappa$ for the double peak SIS system. Left Panel: Using PBL . Right panel: Using grid based method described in § \[sec:S+W\]. Both maps are gridded for easy visualization. Also there are empty grid cells with no image galaxies. The value for those grid cells in the above difference map is set to zero for both reconstructions. As we can see the error in the cores of the peaks is much lower using PBL mass reconstruction[]{data-label="fg:kap_recon_image_2peak"}](f3a.eps "fig:") ![ The plot of the difference between reconstructed convergence,$\kappa$ and true $\kappa$ for the double peak SIS system. Left Panel: Using PBL . Right panel: Using grid based method described in § \[sec:S+W\]. Both maps are gridded for easy visualization. Also there are empty grid cells with no image galaxies. The value for those grid cells in the above difference map is set to zero for both reconstructions. As we can see the error in the cores of the peaks is much lower using PBL mass reconstruction[]{data-label="fg:kap_recon_image_2peak"}](f3b.eps "fig:")

Observation: The Bullet Cluster
-------------------------------

### Observations

Finally we perform a mass reconstruction of the bullet cluster (1E0657-56). This galaxy cluster is a rare supersonic merger in the plane of the sky. Its distinctive structure and orientation makes it an ideal cluster for observing dark matter using gravitational lensing. It consists of two sub-clusters separated by 0.72 Mpc, which have just undergone a merger and are moving away from each other. The western sub-cluster is less massive and the eastern main cluster is more massive. The line-of-sight velocity difference suggest that their cores passed each other  100 Myr ago. The collisionless dark matter in each of the sub-clusters have crossed each other but the fluid-like intracluster plasma is in the process of electromagnetic and thermal interaction producing high X-ray luminosity far removed from lensing mass peaks [@2006ApJ...648L.109C; @2006ApJ...652..937B].

For the Bullet Cluster, we perform a PBL reconstruction only, since it has been well-studied with grid-based methods (using ) and the $\kappa$-contours are publicly available. We use publicly available weak lensing data from the Bullet Cluster Project Page[^2]. The catalog was constructed using data from three different instruments: the ESO/MPG Wide Field Imager, IMACS on Magellan, and two pointings of ACS on HST. The shapes of the galaxies were measured independently on each of the image sets averaging for the common galaxies. The weighting for each galaxy is based on its significance of detection in every image set and normalized appropriately [@2006ApJ...648L.109C].

The catalogs were combined using weighted average reduced shear measurements and the weights of individual galaxies were increased when they occurred in several catalogs. This weighting is listed in the shear catalog. We include this weighting in our reconstructions as well and choose only those images with a weighting greater than 1. As we have already illustrated in the simulations PBL is most effective when the information density is variable, i.e close to the core of the clusters. In case of the bullet cluster we zoom into a region bounded by $104.53^o$ to $104.69^o$ in right ascension and $55.92^o$ to $55.97^o$ in declination. Following this cut, our sample includes 1259 weak lensing background galaxies. In order to do the mass reconstruction we use the average redshift of this sample, $z=0.91$.

### Reconstruction

The Bullet Cluster was made famous by the direct detection of dark matter by [@2006ApJ...648L.109C]. Indeed, since one of the major findings of this group is that the dark matter appears offset from X-ray emissions, we do not include any prior model when reconstructing the system, but are able to achieve fast convergence with two clearly visible peaks. This reconstruction guides us in choosing an initial condition for subsequent $\chi^2$ minimization.

We have calculated the integrated mass within 150 kpc of each peak. The main peak has a mass of $1.57\times10^{14}M_\odot$ and the sub-cluster has a mass of $0.9\times10^{14}M_\odot$. [@2004ApJ...604..596C] report a value of $(1.02 \pm 0.16)\times
  10^{14}M_\odot$ for the main peak and $(0.66\pm0.19)\times
  10^{14}M_\odot$ for the sub-cluster within 150 kpc of the each peak. In each case, our estimate exceeds that of Clowe et al. by approximately 3.4 $\sigma$. However, more a more recent S+W reconstruction by the same group [@2006ApJ...652..937B] yields masses of $(2.8\pm0.2)\times 10^{14}M_\odot$ around the main peak and $(2.3\pm0.2)\times 10^{14}M_\odot$ around the sub-cluster within 250 kpc of each peak. Inclusion of strong lensing information makes reconstruction of the cores more accurate and also leads to a higher estimates of the mass. Even correcting for the greater area, this suggests Clowe’s initial mass estimate may have been low.

Our mass estimates using PBL is higher than the weak lensing reconstruction of [@2004ApJ...604..596C], and thus more in line with the S+W results. This is a result of a difference in method. For example, we start from an initial condition and iterate to the correct solution whereas [@2004ApJ...604..596C] have fitted a radially averaged shear profile to the NFW or King profile. As already seen in the simulations using an initial condition recovers values of $\kappa$ close to the core with greater accuracy. This implies that while most weak lensing $\kappa$ maps report $\kappa$-contours less than 1 using initial condition and PBL we are able to get $\kappa$ greater than 1. This implies that the mass we measure will also be greater than the typical weak lensing mass measurement. Also to measure the mass of the sub-cluster [@2004ApJ...604..596C] have removed the mass of the main cluster to avoid over-estimation of the mass, we have not considered this effect in our reconstruction.

In Fig. \[fg:bullet\_new\], we show our PBL reconstruction of the bullet cluster. Note that, despite using weak lensing signals only, we are able to identify both density peaks and using initial conditions we are able to get $\kappa>1$ for the main peak. We also do a comparison of the publicly available $\kappa$-contour with the $\kappa$-contours reconstructed using PBL. The location of the main peak coincides for both reconstruction. The sub-cluster contours for PBL are slightly removed from the publicly available $\kappa$-contours.

![A weak-lensing only reconstruction of the bullet cluster using PBL described in § \[sec:particle\]. Note that both substructure peaks are clearly identified. Upper Panel: This the $\kappa$-map using PBL. The cross denotes the centroid of the multiply imaged positions. Lower Panel: This a comparison of the $\kappa$ contour derived using PBL(solid) and the publicly available contour plot of $\kappa$(dashed). []{data-label="fg:bullet_new"}](f4a.eps "fig:") ![A weak-lensing only reconstruction of the bullet cluster using PBL described in § \[sec:particle\]. Note that both substructure peaks are clearly identified. Upper Panel: This the $\kappa$-map using PBL. The cross denotes the centroid of the multiply imaged positions. Lower Panel: This a comparison of the $\kappa$ contour derived using PBL(solid) and the publicly available contour plot of $\kappa$(dashed). []{data-label="fg:bullet_new"}](f4b.eps "fig:")

Error analysis for PBL will be discussed in detail in future papers. In particular the noise covariance matrix,$\langle(\kappa-\langle\kappa\rangle)(\kappa-\langle\kappa\rangle)^T\rangle$, will give us important insights into the errors caused by the reconstruction method. A bootstrap method can also be used to determine error bars on mass measurements from observations. In case of simulations several Monte Carlo realizations of the noise can be used to study the errors.

Discussion and Future Prospects {#sec:discuss}
===============================

Additional Signals
------------------

Thus far, we have developed the formalism for PBL, and done worked examples demonstrating how it may be applied to weak-lensing reconstruction. It is designed to model structure hierarchically, in part because of the great success of Strong+Weak lensing analysis.

Several groups have already shown how multiple image positions may be added to the information yielded by lens ellipticities to produce very high quality mass-maps of clusters. It was our desire to maximally exploit the different information scales of the strong and weak lensing signals which motivated the development of PBL in the first place.

However, there is yet more information besides image differences potentially available which may be utilized in a reconstruction. Consider that in addition to the two constraints generated by the positional difference between two images, we also can measure a flux ratio, and 2 ellipticity differences. Thus, in principle, we have 5 measurable, model parameters per strong lensing pair rather than 2, and in an idealized case, this improves potential resolution of a system in the strong lensing regime by $\sqrt(5/2)\simeq 1.6$.

As a way of guiding the future development of PBL, we discuss possible future avenues of investigation below.

### Flux

Apart from the centroid position, the Petrosian flux of an image is the most straightforward to measure. The relationship between magnification lens is simply the inverse of the determinant of the projection matrix: $$\mu=\frac{1}{(1-\kappa)^2-|\gamma|^2}$$ Unlike the displacement vectors (), which are simple linear operators of the potential field (the gradient), or the weak-lensing shear field which is nearly so (since in the limit of $\kappa << 1$, the image ellipticity is an unbiased estimator of the shear field), the flux is a highly nonlinear function of the shear and convergence fields. This accounts, in part, for the reason that it has not been used previously in cluster reconstructions. Here we note that [@2007ApJ...663...29S] show that the image positions itself constrain the fluxes for a source with three non-collinear components. This is a special case, for cluster lensing three component sources for strong lensing may not always be available. Also, [@2007MNRAS.376..180N] use magnification information in their parametric mass modelling of clusters.

The other major consideration is that magnification is not a smoothly varying function of the potential fields. It is well-known that on the critical curves, magnification goes to infinity (see, e.g. [@1992grle.book.....S] for an extensive discussion), but this is a set of measure zero, so in and of itself produces no problem. The issue is that the parity of the image reverses as an image that crosses the critical curve.

Negative magnification means nothing more than reversal of image parity, and thus cannot generally be easily detected. Thus, we are much more interested in computing terms which scale like $\mu^2$. Indeed, since we cannot measure the magnification directly, but only the flux, we propose that the combination: $$\frac{\mu_A^2-\mu_B^2}{\mu_A^2+\mu_B^2}=\frac{f_A^2-f_B^2}{f_A^2+f_B^2}$$ is directly measurable, and has no poles.

Even so, a lensing model [*predicts*]{} a parity for a particular image, and as with ellipticity, minimization, there is a discontinuity in the derivatives. In Fig. \[fg:mu\_parity\] we show the magnification (including sign) as a function of convergence and shear.

![The magnification as a function of shear and convergence. The lower panel is a simple slice through the upper, with the choice $\gamma=\kappa$. Neither the magnification nor its derivatives are a continuous function. Moreover, flux ratios are only measurable for systems with at least two images (obviously). One or more of the images will necessarily have negative parity. Thus, a solution to the potential field which is found using standard relaxation methods will not normally converge to a negative parity estimate for any magnification.[]{data-label="fg:mu_parity"}](f5a.eps "fig:") ![The magnification as a function of shear and convergence. The lower panel is a simple slice through the upper, with the choice $\gamma=\kappa$. Neither the magnification nor its derivatives are a continuous function. Moreover, flux ratios are only measurable for systems with at least two images (obviously). One or more of the images will necessarily have negative parity. Thus, a solution to the potential field which is found using standard relaxation methods will not normally converge to a negative parity estimate for any magnification.[]{data-label="fg:mu_parity"}](f5b.eps "fig:")

### Ellipticity Differences

Likewise, while most measurements of the shear are based on an assumption that any [*given*]{} image is randomly oriented, two images of the same source are not. The difference in their measured ellipticity can be wholly modeled by the relative lensing fields at their respective locations. If both images were in the weak regime, we would be able to use the simple estimator $$\varepsilon_A-\varepsilon_B=\simeq \gamma_A-\gamma_B$$ where all terms in the equation are complex, and thus provide two constraints with high signal to noise per image pair.

In general, however, a more likely configuration is that one image may be in the strong regime, and one in the weak. If we can determine from the configuration of lenses which is which, we might imagine a better estimator as: $$\varepsilon_A-\varepsilon_B=\frac{1}{g^\ast_A}-g_B$$ with the only associated noise corresponding to photon noise rather than random variance in the intrinsic ellipticity of the images.

### Flexion

Thus far, the analysis of clusters in the weak or semi-weak regime has primarily relied on shear. However, recently, [@2007arXiv0710.2262O], and [@2007ApJ...666...51L] have worked on reconstructing A1689 using flexion. In particular, the Okura group used a Fourier inversion suggested by [@2007arXiv0709.1003S]. However, the advantage of our proposed PBL is that flexion (and, in principle, [*any*]{} higher-order derivative of the potential) may be explicitly included as additional constraints in the cluster reconstruction. Unlike Fourier techniques, which rely on binning of the data, the PBL method will allow us to exploit the natural small-scale signal probed by flexion.

Summary
-------

We have developed PBL, a new particle based technique of mass reconstruction of clusters. The distinguishing feature of PBL is its ability to adjust its smoothing scale depending on the local signal to noise or the type of constraint and thus not require any regularization. PBL has the scope of calculating derivatives up to any order. Hence, lensing constraints that are a function of the derivatives of the potential can be easily included in the reconstruction. In this paper we have successfully applied PBL to do weak lensing only mass reconstruction for a single peak and a double peak system. We have made the codes for PBL publicly available for application weak lensing measurements through our website(see § \[sec:particle\]). The codes have been tested on the data sets and simulations described in the paper. A larger data sample will require modification of the current version of codes.

As already explained PBL is a method of discretizing data and not a minimization method. A $\chi^2$ minimization does not necessarily ensure reaching a global minimum. In many cases the global minimum is guarded by steep walls surrounded by shallow valleys. Without any prior knowledge of the mass distribution it is very easy to get trapped in a shallow valley and not reach the global minimum. We have started with an initial condition and interpolated the ellipticity function to aid us in this regard.

In future work we will be including the additional constraints, like the flux ratios, ellipticity differences and flexion along with measured ellipticities and strong lensing positions. We will also be exploring different minimization schemes to facilitate convergence to a global minimum.

Acknowledgments
===============

The authors would like to thank for useful discussions with Mike Jarvis, Bhuvnesh Jain and Kevin Olson. We would also like to thank Marusa Brada[v c]{} and Douglas Clowe for providing us with the “Bullet Cluster” (1E0657-56) data. DMG gratefully acknowledges discussions with Richard Massey. SD would like to thank Ben Metcalf, Prasenjit Saha and Sudeep Das for valuable comments. We would also like to thank the referee for thoughtful suggestions. This work was supported by NASA ATP NNG05GF61G.

[68]{} natexlab\#1[\#1]{}

, H. M., [Saha]{}, P., & [Williams]{}, L. L. R. 1998, , 116, 1541

, S. W., [Rapetti]{}, D. A., [Schmidt]{}, R. W., [Ebeling]{}, H., [Morris]{}, R. G., & [Fabian]{}, A. C. 2007, , 1149

, D. J., [Goldberg]{}, D. M., [Rowe]{}, B. T. P., & [Taylor]{}, A. N. 2006, , 365, 414

, M., & [Schneider]{}, P. 2001, , 340, 291

, M., Clowe, D., [Gonzalez]{}, A. H., [Marshall]{}, P., [Forman]{}, W., [Jones]{}, C., [Markevitch]{}, M., [Randall]{}, S., [Schrabback]{}, T., & [Zaritsky]{}, D. 2006, , 652, 937

, M., [Erben]{}, T., [Schneider]{}, P., [Hildebrandt]{}, H., [Lombardi]{}, M., [Schirmer]{}, M., [Miralles]{}, J.-M., [Clowe]{}, D., & [Schindler]{}, S. 2005, , 437, 49

, M., [Schneider]{}, P., [Lombardi]{}, M., & [Erben]{}, T. 2005, , 437, 39

Broadhurst, T., et al. 2005, , 621, 53

, T., [Takada]{}, M., [Umetsu]{}, K., [Kong]{}, X., [Arimoto]{}, N., [Chiba]{}, M., & [Futamase]{}, T. 2005, , 619, L143

, M., [Bartelmann]{}, M., [Meneghetti]{}, M., & [Moscardini]{}, L. 2006, , 458, 349

, D., [Brada[č]{}]{}, M., [Gonzalez]{}, A. H., [Markevitch]{}, M., [Randall]{}, S. W., [Jones]{}, C., & [Zaritsky]{}, D. 2006, , 648, L109

, D., [Gonzalez]{}, A., & [Markevitch]{}, M. 2004, , 604, 596

, D., & [Schneider]{}, P. 2002, , 395, 385

, D., [Fuselier]{}, E., [Benitez]{}, N., [Broadhurst]{}, T., [Frye]{}, B., & [Ford]{}, H. 2008, ArXiv e-prints, 803

, A., & [Sheth]{}, R. 2002, , 372, 1

, J. M., [Protopapas]{}, P., [Sandvik]{}, H. B., & [Tegmark]{}, M. 2005, , 360, 477

, H., [White]{}, D. A., & [Rangarajan]{}, F. V. N. 2006, , 368, 65

, E. E., [Gorenstein]{}, M. V., & [Shapiro]{}, I. I. 1985, , 289, L1

, R. A., [Allgood]{}, B., [Kravtsov]{}, A. V., [Primack]{}, J. R., [Buote]{}, D. A., & [Bullock]{}, J. S. 2005, ArXiv Astrophysics e-prints

—. 2007, , 377, 883

, C., [Moreno]{}, J., [Sheth]{}, R. K., & [Tormen]{}, G. 2007, , 376, 977

, D. M., & [Bacon]{}, D. J. 2005, , 619, 741

, M. E., [Taylor]{}, A. N., [Meisenheimer]{}, K., [Dye]{}, S., [Wolf]{}, C., & [Thommes]{}, E. 2002, , 568, 141

, J. E., & [Gott]{}, J. R. I. 1972, , 176, 1 Heymans, C., et al.  2008, , 385, 1431

, H., [Franx]{}, M., [Kuijken]{}, K., [Carlberg]{}, R. G., [Yee]{}, H. K. C., [Lin]{}, H., [Morris]{}, S. L., [Hall]{}, P. B., [Patton]{}, D. R., [Sawicki]{}, M., & [Wirth]{}, G. D. 2001, , 548, L5

, H., [Franx]{}, M., [Kuijken]{}, K., & [Squires]{}, G. 1998, , 504, 636

, H., [Franx]{}, M., [Kuijken]{}, K., & [van Dokkum]{}, P. G. 2002, , 333, 911

, H., [Yee]{}, H. K. C., & [Gladders]{}, M. D. 2004, in IAU Symposium, Vol. 220, Dark Matter in Galaxies, ed. S. [Ryder]{}, D. [Pisano]{}, M. [Walker]{}, & K. [Freeman]{}, 439–+

, C., & [Berge]{}, J. 2005, , 360, 1393

, N. 1995, , 439, L1

, T., [Mayer]{}, L., [Wadsley]{}, J., [Stadel]{}, J., & [Moore]{}, B. 2007, , 375, 53

, L. J. 2007, , 382, 308

, A., [Kravtsov]{}, A. V., [Valenzuela]{}, O., & [Prada]{}, F. 1999, , 522, 82

, A., [Goldberg]{}, D. M., [Haaga]{}, J. L., & [Massey]{}, R. 2007, , 666, 51

, G. A., [Gioia]{}, I. M., [Hammer]{}, F., [Le F[è]{}vre]{}, O., & [Annis]{}, J. A. 1999, , 136, 117

, F., [Maoz]{}, D., [Sharon]{}, K., [Botticella]{}, M. T., [Della Valle]{}, M., [Gal-Yam]{}, A., & [Panagia]{}, N. 2007, , 1132

, M., [Gonzalez]{}, A. H., [Clowe]{}, D., [Vikhlinin]{}, A., [Forman]{}, W., [Jones]{}, C., [Murray]{}, S., & [Tucker]{}, W. 2004, , 606, 819

, M., [Gonzalez]{}, A. H., [David]{}, L., [Vikhlinin]{}, A., [Murray]{}, S., [Forman]{}, W., [Jones]{}, C., & [Tucker]{}, W. 2002, , 567, L27

, P. 2006, , 372, 1289

, P. J., [Hobson]{}, M. P., [Gull]{}, S. F., & [Bridle]{}, S. L. 2002, , 335, 1037

, J. J. 2005, Reports of Progress in Physics, 68, 1703

, B., [Ghigna]{}, S., [Governato]{}, F., [Lake]{}, G., [Quinn]{}, T., [Stadel]{}, J., & [Tozzi]{}, P. 1999, , 524, L19

, K., [Springel]{}, V., & [Hernquist]{}, L. 2004, , 348, 435

, P., [De Lucia]{}, G., & [Springel]{}, V. 2007, , 376, 180

, P., [Kneib]{}, J.-P., [Smail]{}, I., [Treu]{}, T., [Ellis]{}, R., [Moran]{}, S., [Limousin]{}, M., & [Czoske]{}, O. 2007, ArXiv e-prints, 711

, P., & [Refregier]{}, A. 2000, , 538, L113

, P., & [Springel]{}, V. 2004, , 617, L13

, Y., [Umetsu]{}, K., & [Futamase]{}, T. 2007, ArXiv e-prints, 710

, A. O., [Levine]{}, H., & [Wambsganss]{}, J. 2001, [Singularity theory and gravitational lensing]{} (Singularity theory and gravitational lensing / Arlie O. Petters, Harold Levine, Joachim Wambsganss. Boston : Birkh[ä]{}user, c2001. (Progress in mathematical physics ; v. 21))

, K., [Diaferio]{}, A., & [Natarajan]{}, P. 2007, , 657, 183

, P., [Williams]{}, L. L. R., & [Abdelsalam]{}, H. M. 2001, in Astronomical Society of the Pacific Conference Series, Vol. 237, Gravitational Lensing: Recent Progress and Future Go, ed. T. G. [Brainerd]{} & C. S. [Kochanek]{}, 279–+

, P., [Williams]{}, L. L. R., & [Ferreras]{}, I. 2007, , 663, 29

, D. J., [Treu]{}, T., [Smith]{}, G. P., & [Ellis]{}, R. S. 2003, ArXiv Astrophysics e-prints

, P. 1995, , 302, 639

, P., [Ehlers]{}, J., & [Falco]{}, E. E. 1992, [Gravitational Lenses]{} (Gravitational Lenses, XIV, 560 pp. 112 figs.. Springer-Verlag Berlin Heidelberg New York.  Also Astronomy and Astrophysics Library)

, P., & [Er]{}, X. 2007, ArXiv e-prints, 709

, P., & [Seitz]{}, C. 1995, , 294, 411

, P., & [Weiss]{}, A. 1992, , 260, 1

, S., & [Schneider]{}, P. 1998, ArXiv Astrophysics e-prints

, S., [Schneider]{}, P., & [Bartelmann]{}, M. 1998, , 337, 325

, V., & [Hernquist]{}, L. 2003, , 339, 312

, A. N., [Bacon]{}, D. J., [Gray]{}, M. E., [Wolf]{}, C., [Meisenheimer]{}, K., [Dye]{}, S., [Borch]{}, A., [Kleinheinrich]{}, M., [Kovacs]{}, Z., & [Wisotzki]{}, L. 2004, , 353, 1176

, L. M., & [Fabian]{}, A. C. 2006, , 368, 518

, G., [Kaiser]{}, N., [Cole]{}, S., & [Frenk]{}, C. 1996, in Bulletin of the American Astronomical Society, Vol. 28, Bulletin of the American Astronomical Society, 1385–+

, D., [Tyson]{}, J. A., [Margoniner]{}, V. E., [Cohen]{}, J. G., & [Dell’Antonio]{}, I. P. 2001, , 557, L89

, M. M. 2007, , 376, 1173

[^1]: http://www.physics.drexel.edu/$\sim$deb/PBL.htm

[^2]: http://flamingos.astro.ufl.edu/1e0657/public.html