text:
# Energetically stable particle-like Skyrmions in a trapped Bose-Einstein condensate

## Abstract

We numerically show that a topologically nontrivial 3D Skyrmion can be energetically stable in a trapped two-component atomic Bose-Einstein condensate, for the parameters of 87Rb condensate experiments. The separate conservation of the two atomic species can stabilize the Skyrmion against shrinking to zero size, while drift of the Skyrmion due to the trap-induced density gradient can be prevented by rotation or by a laser potential.

Localized topological excitations that do not perturb the order parameter field at large distances from the particle, and which are characterized by a topologically invariant winding number, are well-known in nuclear and elementary particle physics SKY61 ; WIT83 ; BAT01 ; FAD97 . While their study in nuclear physics remains an experimental challenge, the recent experimental progress in atomic Bose-Einstein condensates (BECs) with internal spin degrees of freedom MAT99 ; LEA02  has raised the possibility of the existence of well-localized topological Skyrmions in atomic gases RUO01 ; Stoof ; BAT02 . In this Letter we identify, and show how to overcome, the specific instabilities of Skyrmions in trapped two-species atomic BECs, and hence demonstrate their energetic stability under realistic experimental conditions.

Battye et al. BAT02  recently considered an infinite homogeneous two-species BEC, with constant total atom density. They showed that an energetically stable Skyrmion may exist as a result of phase separation of the two species, which suppresses the decay. These calculations were extended to non-constant total atom density, and to the trapping of one component BAT02b . In this paper we show that in a harmonically trapped system there are additional instabilities, not considered in Refs. BAT02 ; BAT02b , which will play a crucial role in the experimental realization of Skyrmions in atomic BECs. Additional physical mechanisms, such as rotation or optical potentials, will be required for stability. We also show how density fluctuations, associated, e.g., with phonon emission adams , are important in the Skyrmion decay process.

There has been an explosion of interest in vortex and soliton experiments in atomic BECs ANG , and we anticipate similar developments for other topological objects. Hence we identify the Skyrmion energetic stability criteria for the parameters of the JILA two-species 87Rb experiments MAT99 . We find a threshold frequency below 0.1$`\omega `$, when only one species is rotated, and a narrow window of rotation frequencies for the entire system around 0.085$`\omega `$. We numerically evaluate the stable configurations (Fig. 1) by minimizing the energy of the full 3D mean-field theory of coupled Gross-Pitaevskii equations (GPEs) and calculate the associated topological charges.

A Skyrmion is a topological particle-like soliton solution with a coreless 3D texture. Besides their intrinsic fundamental interest, Skyrmions have important applications in nuclear physics SKY61 ; WIT83 ; BAT01 , and analogous structures are postulated for early Universe cosmology DAV89 . Skyrmions are localized objects such that the order parameter field becomes uniform sufficiently far from the particle. As a result, we may enclose the Skyrmion by a sphere on which the field configuration has a constant asymptotic value. This sphere can be mapped to a single point in the order parameter space. Topologically, it indicates that the 3D physical region inside the sphere can be represented by a compact 3-sphere $`S^3`$ (a sphere in 4D) com2 . Here we consider a SU(2) order parameter space which topologically also corresponds to $`S^3`$. Then the mappings from the compactified 3D physical region to the order parameter space are represented by the field $`U(𝐫)`$: $`S^3SU(2)`$. The crucial point is that having associated the enclosing sphere with a single point in $`S^3`$, we place the elements of the physical space into a correspondence with the elements of the compact order parameter space SU(2). The mappings from $`S^3`$ to SU(2) fall into topological classes, each characterized by an integer-valued, topologically invariant winding number $`W`$ SKY61 ; RUO01 ; BAT02 :

$$W=\frac{ϵ_{\alpha \beta \nu }}{24\pi ^2}d^3r\mathrm{Tr}[U(_\alpha U^{})U(_\beta U^{})U(_\nu U^{})].$$
(1)
Here $`ϵ_{\alpha \beta \nu }`$ denotes the completely antisymmetric tensor, the repeated indices are summed over, and $`_i`$ represents the derivative with respect to the spatial coordinate $`x_i`$. The topological charge describes how many times the set of points in physical space is ‘wrapped’ over the order parameter space of SU(2) for a given element $`U(𝐫)`$.

With cold atoms an SU(2) order parameter space is afforded most simply by a two-component BEC, whose interactions effectively fix the local value of the total density $`|\psi _+(𝐫)|^2+|\psi _{}(𝐫)|^2=\rho (𝐫)`$ of the two complex macroscopic wave functions $`\psi _+`$ and $`\psi _{}`$ RUO01 . We may then express the two-component BEC as

$$\left(\begin{array}{c}\psi _+(𝐫)\\ \psi _{}(𝐫)\end{array}\right)=\sqrt{\rho (𝐫)}U^{}(𝐫)\left(\begin{array}{c}1\\ 0\end{array}\right).$$
(2)
We search for a topologically nontrivial solution for $`U(𝐫)`$ with a nonvanishing winding number, determined by Eq. (1). By minimizing the energy of the corresponding trapped BECs, we may investigate the energetic stability of the Skyrmion and find its equilibrium configuration. As an initial state for the numerical simulations we use $`U(𝐫)=\mathrm{exp}[i\lambda (𝐫)\stackrel{}{\sigma }\widehat{𝐫}]`$ com2b , or a similar state displaced from the trap center, where $`\sigma _i`$ denote the Pauli spin matrices and $`\widehat{𝐫}`$ represents the unit radial vector. A direct substitution into Eq. (1) yields $`W=1`$ for a monotonic function $`\lambda `$ with $`\lambda (\mathrm{𝟎})=0`$ and $`\lambda =\pi `$ at the gas boundary. The corresponding BEC wave functions read:

$$\left(\begin{array}{c}\psi _+(𝐫)\\ \psi _{}(𝐫)\end{array}\right)=\sqrt{\rho }\left(\begin{array}{c}\mathrm{cos}[\lambda (𝐫)]i\mathrm{sin}[\lambda (𝐫)]\mathrm{cos}\theta \\ i\mathrm{sin}[\lambda (𝐫)]\mathrm{sin}\theta \mathrm{exp}(i\varphi )\end{array}\right).$$
(3)
Here $`(\lambda ,\theta ,\varphi )`$ can also be understood as the spherical angles of the 3-sphere. Then the boundary of the atom cloud corresponds to the pole $`\lambda =\pi `$ of the 3-sphere.

Due to the topological stability of the Skyrmion, any continuous deformation of Eq. (3), without altering the asymptotic boundary values, is still a Skyrmion with $`W=1`$. Here, $`\psi _{}`$ (the “line component”) forms a vortex line with one unit of circulation around the core oriented along the $`z`$ axis (Fig. 1). Moreover, $`\psi _+`$ (the “ring component”) forms a vortex ring with one unit of circulation around its core on the circle $`z=0,r=\lambda ^1(\frac{\pi }{2})`$ RUO01 . The line component is confined in a toroidal region inside the ring component in the neighborhood of the ring core.

Despite its topological stability, the Skyrmion solution in several physical systems is energetically unstable against shrinking to zero size without additional stabilizing features. This can be understood by means of simple scaling arguments: if the kinetic energy density, or the order parameter ‘bending energy’, is quadratic in the gradient of the order parameter, this energy density scales as $`1/R^2`$ with respect to the size $`R`$ of the Skyrmion. Because the Skyrmion occupies a volume proportional to $`R^3`$, the energy $`ER`$ monotonically decreases with the size of the Skyrmion. In the Skyrme model, the stability is provided by an additional interaction term, with the energy density scaling as $`1/R^4`$ SKY61 .

It was proposed that, in a two-component BEC, the separate conservation of the species can effectively stabilize the Skyrmion in a homogeneous space against collapse to zero size BAT02 . In the regime of phase separation, with the scattering lengths satisfying $`a_{++}a_{}a_+^2`$, the two species can strongly repel each other and the toroidal filling due to the line component can prevent the vortex ring from shrinking to zero radius RUO01 ; BAT02 . Consequently, a filled vortex ring, as opposed to an empty vortex ring, can be more stable against collapse. Moreover, in the Skyrmion the filling has one unit of angular momentum about the $`z`$ axis resulting in a $`1/r^2`$ centrifugal barrier that further prevents the shrinking. This should be contrasted to the Skyrmions in a ferromagnetic spin-1 BEC Stoof , which are closely related to the Skyrmion excitations proposed in superfluid liquid 3He VOL77 , and can freely mix the atoms between the different spin components. As a result, no stabilizing features due to the atom number conservation exist in that case.

We next turn to numerical studies of the energetic stability of the Skyrmion using the full 3D mean-field theory of the coupled GPEs without additional simplifications:

$$i\mathrm{}\dot{\psi }_i=\left(\frac{\mathrm{}^2}{2m}^2+V(𝐫)+\underset{k}{}\kappa _{ik}|\psi _k|^2\right)\psi _i.$$
(4)
Our simulations are performed for the parameters of the JILA two-species vortex experiment using perfectly overlapping isotropic trapping potentials for both states $`V(𝐫)=m\omega ^2r^2/2`$, with $`\omega =2\pi \times 7.8`$ Hz and harmonic oscillator length $`x_{ho}(\mathrm{}/m\omega )^{1/2}3.86\mu `$m MAT99 . In the interaction coefficients, $`\kappa _{ij}4\pi \mathrm{}^2a_{ij}N_i/m`$, the number of atoms in each species is represented by $`N_i`$, $`a_{ii}`$ denotes the intraspecies, and $`a_{ij}`$ ($`ij`$) the interspecies, scattering length. For the $`|F=2,m_f=2`$ and $`|1,1`$ spin states of 87Rb we have HALL98 : $`a_{++}=5.67`$nm, $`a_{}=5.34`$nm, $`a_+=5.50`$nm. The GPEs depend on the dimensionless interaction coefficients $`\kappa _{}^{}{}_{ij}{}^{}4\pi N_ia_{ij}/x_{ho}`$. Hence, the dynamics is unchanged in any scaling of length and time, which does not change $`N_i^2\omega `$.

We found the ground states by imaginary time evolution of the GPEs (4), using Runge-Kutta RK4IP  and split-step algorithms. At every time step, we separately normalized both wave functions to fix the atom number in each component. The numerical simulations were fully 3D, allowing for cylindrically asymmetric dynamics. The most demanding numerics was performed on a parallel multiprocessor supercomputer, using up to 32 processors. Spatial grids of $`128^3`$, $`256^3`$, and $`512^3`$ were used. For a typical spatial range of about $`30x_{ho}`$ the corresponding grid spacings are respectively $`0.23x_{ho}`$, $`0.12x_{ho}`$, and $`0.06x_{ho}`$. A typical imaginary time step was $`0.0025/\omega `$.

As a cylindrically symmetric initial state we used the Skyrmion (3) with the Thomas-Fermi density profile com2b , resulting in the winding number $`W=1`$. The winding number was calculated during the imaginary time evolution to determine the stability of the Skyrmion. For small BECs, with $`N_T=10^4`$ atoms ($`N_TN_++N_{}`$), the nonlinear repulsion between the two species was too weak to prevent the Skyrmion from shrinking for any $`N_{}/N_T`$. Also the line component $`\psi _{}`$ rapidly diffused to the boundary of the finite-size atomic cloud, altering the boundary conditions and resulting in a decreasing winding number $`W<1`$. Because the topological stability of the Skyrmion necessitates a well-defined asymptotic boundary condition with no phase variation, the Skyrmion was clearly lost. This decay mechanism is characteristic of trapped BECs and cannot occur in homogeneous systems, which implicitly assume large $`N_+/N_T`$.

For large BECs with $`N_T>10^6`$, represented in Fig. 2, the same instability mechanism was observed with a very large fraction of line component $`N_{}/N_T0.75`$. In other cases displayed in Fig. 2, the winding number typically evaluated to either 0 or 1 to better than 1%, with $`\lambda =\pi `$ well preserved at the boundary. A drop from $`W=1`$ to $`W=0`$ indicated instability of the Skyrmion against shrinking, see inset to Fig. 3. The Skyrmion collapse via shrinking occurred when the high density central part of the ring component, which passes through the line component vortex core, pinched off. The total density varies significantly at the trap center when the vortex ring core collapses, emphasizing the importance of density fluctuations in the decay process, see Fig. 3 (top). With sufficiently large BECs, even for the collapsed Skyrmions with no vortex ring, the local energetic minimum can correspond to line component trapped inside a toroidal region with $`\psi _+`$ forming an atom-optical confinement.

For sufficiently large total and line component atom numbers, the nonlinear repulsion between the two species was strong enough to inhibit the collapse of the vortex ring, stabilizing the Skyrmion against shrinking for cylindrically symmetric initial states (Fig. 2). However, due to the inhomogeneous potential, the Skyrmion was still unstable with respect to drift towards the edge of the BEC where its energy is lower. This was found using cylindrically asymmetric initial states, emphasizing the importance of the full 3D simulation. The drift occurred with the vortex line moving towards the boundary. Once the vortex line drifted to a low-density region, the nonlinear repulsion was no longer strong enough to prevent the shrinking of the Skyrmion and the collapse occurred as described above, see Fig. 3 (bottom). The drift reduces the total angular momentum of the atoms indicating an energetic instability. Physically, this results from dissipation, e.g., due to thermal atoms.

We have investigated several mechanisms by which Skyrmions might be stabilized against the vortex drift instability. Perhaps the simplest is to rotate the line vortex component about the $`z`$ axis, above its critical angular velocity, thus stabilizing the vortex. This should always work for cylindrically symmetric traps in which the symmetric density prevents the rotation coupling into the other BEC component. For the case of Fig. 1 a line component angular velocity of $`0.1\omega `$ stabilised the Skyrmion.

Rotation of both species introduced a new instability mechanism. For sufficiently rapid rotations the line component again reached the boundary, resulting in $`W<1`$. At the same time the ring vortex accomodated the rotation by twisting and finally breaking at the gas boundary. Nevertheless, for the case of Fig. 1, we found a small range of angular velocities around $`0.085\omega `$ for which the Skyrmion was stable against this, and against line vortex drift. It stabilized off the rotation axis, indicating an equilibrium between the lower energy of the line vortex and the higher energy of the ring component threading through it with increased distance from the trap center. The breaking of cylindrical symmetry implies a family of degenerate Skyrmions parametrized by the polar angle $`\varphi `$.

Another method for stabilizing the Skyrmion is to inhibit the drift towards low-density regions by creating a positive density gradient around the vortex line. This can be implemented with a blue-detuned Gaussian laser beam along the $`z`$ axis blue beam , providing a cylindrically symmetric repulsive Gaussian dipole potential perpendicular to the $`z`$ axis. The total potential then has the “Mexican hat” form: $`V=m\omega ^2r^2/2+V_0\mathrm{exp}[2(x^2+y^2)/w^2]`$. Setting $`w=7x_{ho}`$ and $`V_0=25\mathrm{}\omega `$, corresponding to about twice the width of the vortex line core, successfully stabilized a Skyrmion with $`N_+=N_{}=8\times 10^6`$.

Proving numerically that a Skyrmion is stable is difficult. The essential requirement is that it be stationary under imaginary time evolution. This was determined by plotting the density on 1D sections through the system, and the phase variation on 2D slice planes. In unstable cases these would evolve until the Skyrmion decayed. When they became asymptotically constant as a function of imaginary time the Skyrmion was judged to be stable. Convergence was particularly slow, and hence difficult to establish, when both components were rotated.

The Skyrmion can be created using electromagnetic fields to imprint topological phase singularities on the matter field while changing the internal state of the atoms, as proposed in Ref. RUO01 . In particular, a vortex ring can be engineered in a controlled way with an appropriate phase-coherent superposition of three orthogonal standing waves. Due to dissipation, the prepared Skyrmions relax to the ground states calculated here.

Different spin textures are often referred to as “Skyrmions”, even when they are not characterized by the topological invariant (1). It is important to emphasize that the Skyrmions studied in this paper are fundamentally different from the nonsingular Anderson-Toulouse or Mermin-Ho vortices AND77 , frequently also called Skyrmions Ho . Although these vortices may be energetically stable in a rotating trap MIZ02 , they are topologically trivial in the 3D ferromagnetic spin-1 BEC and, analogously to the vortex lines with two units of circulation, they can be continuously deformed to a topologically trivial uniform spin distribution.

###### Acknowledgements.

This research was supported by the EPSRC, the Australian Research Council, and the Australian Partnership for Advanced Computing.


format: text

text:
# 1 . Introduction

## 1 . Introduction

The main purpose of this paper is to bring to reader’s attention the current status of microlensing searches. It is likely that within a year or two these searches will lead to firm discoveries of stellar mass black holes and planets. The paper is centered on current developments at the Optical Gravitational Lensing Experiment (OGLE, Udalski et al. 1997, 2002). It is motivated by the fact that OGLE team is small and hence not in a position to do the diverse follow-up observations, which are needed to make conclusive case for black hole lenses, or to measure planetary masses. This is one of the reasons the OGLE data are made public domain with no strings attached, so other teams or individuals may use them as the basis for their own follow-up projects.

The search for gravitational microlensing, and follow-up observations, are a well developed ‘industry’, with dozens of events reported every year in real time by the MOA collaboration (Bond et al. 2001):

http://www.roe.ac.uk/%7Eiab/alert/alert.html

and hundreds reported by the OGLE collaboration (Udalski et al. 1997, 2002):

http://www.astrouw.edu.pl/~ogle/ogle3/ews/ews.html

A major increase in the discovery rate is to be expected when the upgrade of MOA capability is completed, with the first light from the new 1.8 meter telescope expected in 2004 or 2005:

http://www3.vuw.ac.nz/scps/moa/

While there is a very large diversity of research topics related to microlensing, as well as very many other results of the surveys (cf. Paczyński 1996a, Gould 2001, and references therein), this paper is focused on just two extreme cases: the search for black holes and the search for planets, i.e. very large and very small mass lenses. A definite detection of stellar mass black holes and planets is an interesting goal, almost certainly achievable within a year or two. Only microlensing is capable of detecting Earth mass planets and isolated black holes with the technology which already exists. This paper outlines OGLE contribution to the task. However, to be successful it requires diverse follow-up observations to be carried out by astronomers not necessarily associated with OGLE.

## 2 . Black Holes

The discovery of long duration microlensing events prompted suggestions that some of the lenses may be due to stellar mass black holes (Bennett et al. 2002a,b; Mao et al. 2002; Smith 2003). The argument is simple: the time scale of a microlensing event, i.e. the time $`t_E`$ it takes a lens to move with respect to the source by the Einstein radius $`\phi _E`$ is given as (cf. Paczyński 1996a, and references therein)

$$t_E=\frac{\phi _E}{\dot{\phi }}=1.01yr\left(\frac{M}{M_{}}\right)^{1/2}\left(\frac{8kpc}{D}\right)^{1/2}\left(\frac{1masyr^1}{\dot{\phi }}\right),$$
$`(1)`$
where $`\dot{\phi }`$ is a relative proper motion of a lens with respect to the source, $`M`$ is the lens mass, and $`D`$ is the effective distance defined as

$$D=\frac{D_sD_d}{D_sD_d},$$
$`(2)`$
where $`D_s`$ is the distance to the source, and $`D_d`$ is the distance to the lensing object (deflector). The angular Einstein radius is given as

$$\phi _E=1.01\mathrm{mas}\left(\frac{M}{M_{}}\right)^{1/2}\left(\frac{8kpc}{D}\right)^{1/2}.$$
$`(3)`$
In most cases the distance to the source is either known or may be estimated, but the distance to the lens, its mass, and its proper motion with respect to the source, are not known. The time scale of a microlensing event, $`t_E`$, is directly measured, and the longer it is the larger the lens mass, other things being equal. Unfortunately, other things ($`D_d`$ and $`\dot{\phi }`$) are not equal and may very a lot from one event to another. Therefore, a long duration is only an indication of a large mass, but only in a statistical sense.

For long events the effect of Earth’s motion around the sun makes the relative motion in the observer – lens – source system not linear (Refsdal 1966, Gould 1992), and leads to the so called photometric parallax effect, which is always measured for long events. The most spectacular case known so far is that of OGLE-1999-BUL-19 (Smith et al. 2002), in which several maxima were observed as a consequence of the Earth’s orbital motion. The amplitude of the parallax effect provides additional information, and partly removes the degeneracy in the eq. (1). The smaller the relative proper motion $`\dot{\phi }`$, the stronger the effect. In the case of OGLE-1999-BUL-19 it was clear that the relative proper motion $`\dot{\phi }`$ was very small, and that small velocity was responsible for the long time scale, $`t_E=1.02`$ years, while the lens mass was sub-solar.

More interesting was OGLE-1999-BUL-32 = MACHO-99-BLG-22 (Mao et al. 2002, Bennett et al. 2002b), the longest event so far, with $`t_E=1.75`$ years, and clearly measured, but small parallax effect, indicative of a large relative proper motion, and hence a large lens mass. This is currently the best candidate for gravitational microlensing due to stellar mass black hole, but the mass estimates requires statistical analysis, which is not very reliable. In fact the probability estimate that the lens is a black hole changed from 20% to 76% between v1 and v4 edition of astro-ph/0203257 (Agol et al. 2002). It is not possible to make a credible probability estimate as the stellar mass function, all the way from brown dwarfs to intermediate mass black holes, is not known.

The second strong black hole candidate is OGLE-SC5\_2859 (Woźniak et al. 2001). With $`t_E=1.5`$ years this is the second longest event known. The parallax effect is clearly detected, but its amplitude is small. However, there may be other complications (Smith 2003),

Unfortunately, in all cases so far the third quantity which is in principle observable: $`\phi _E`$, remains unknown. In fact there is not a single microlensing event for which Einstein radius was directly measured. The reason is simple: the expected value of $`\phi _E`$ is of the order of a milli arc second (Eq. 2), well below the resolution of existing optical instruments. There is only one case, EROS BLG-2000-5, a binary lens, for which it was possible to infer $`\phi _E=1.4`$ mas as the proper motion was measured with caustic crossing, and photometric parallax effect was also detected (An et al. 2002). The total mass of the lens, and the mass ratio, were found to be $`0.61M_{}`$, and $`M_2/M_1=0.748`$, respectively.

### 2.1 Resolving Microlensing - VLTI

A direct method to measure $`\phi _E`$ was proposed by Delplancke et al. (2001): to resolve the double image generated by microlensing with an optical interferometer, specifically VLTI. The first papers with the first results from VLTI are already published (Segransan et al. 2003), or posted (Kervella et al. 2003), and perhaps VLTI will soon be used to measure microlensed stars.

For an event well covered photometrically the dimensionless angular separation between the source and the lens is known as a function of time:

$$u=\frac{\phi }{\phi _E}=\left[u_{min}^2+\left(\frac{tt_{max}}{t_E}\right)^2\right]^{1/2},$$
$`(4)`$
where $`\phi `$ is the angular separation between the lens and the source. The combined brightness of the two images is larger than the source by the so called magnification factor A:

$$A=\frac{u^2+2}{u\sqrt{u^2+4}}.$$
$`(5)`$
The magnification is the largest when the separation is the smallest, i.e. when $`u=u_{min}`$, which happens at the time $`t=t_{max}`$. Photometric model of a well observed microlensing event provides accurate values of $`u_{min}`$, $`t_{max}`$, and $`t_E`$.

The two images have opposite parity: one positive, with magnification $`A_+`$, always located outside of Einstein circle, another negative, with magnification $`A_{}`$, always located within Einstein circle. We have

$$A_++A_{}=A,A_+A_{}=1.$$
$`(6)`$
The two images are separated by an angle $`\mathrm{\Delta }\phi _{+,}`$ given as

$$\frac{\mathrm{\Delta }\phi _{+,}}{\phi _E}=\sqrt{u^2+4}.$$
$`(7)`$
A single measurement of the angular separation $`\mathrm{\Delta }\phi _{+,}`$ between the two micro-images provides a direct determination of angular Einstein radius, as $`u(t)`$ is provided by the photometric model of the event. While such a measurement is very robust and elegant, the VLTI (or any other optical interferometer) has to be able to observe stars as faint as those which are microlensed. Several microlensed stars are bright, with $`I16`$ mag, and $`K14`$ mag. At peak magnification some reach $`I14`$ mag, i.e. $`K12`$ mag. Bright, long duration OGLE events are described in the Section: Long Current Events.

### 2.2 Astrometric Shift - HST

An alternative astrometric approach is not to resolve the double image, but to measure the motion of the light centroid (Hog et al. 1995, Miyamoto and Yoshi 1995, Walker 1995, Boden et al. 1998, Paczyński 1998). In all these papers it was suggested that a space probe, GAIA or SIM:

http://astro.estec.esa.nl/GAIA/

http://sim.jpl.nasa.gov/whatis/

will be able to do micro arc second precision astrometry, and therefore will be able to determine angular Einstein radii for virtually any microlensing event.

The light centroid shifts with respect to the source position by

$$\delta \phi =\frac{u}{u^2+2}\phi _E,$$
$`(8)`$
The displacement reaches maximum value $`\delta \phi _{max}`$

$$\delta \phi _{max}=2^{3/2}\phi _E0.354\phi _E\mathrm{for}u_{max}=2^{1/2}1.414.$$
$`(9)`$
This angular separation corresponds to a magnification $`A=1.155`$, i.e. $`0.156`$ mag above the baseline. The full range of displacements, from the moment when $`A=1.155`$ on the rising branch, to the moment when $`A=1.155`$ on the descending branch, is $`t_E\times 2\times 2^{1/2}2.8t_E`$. This is inconvenient, as we are targeting long duration events, with $`t_E>1`$ year.

While SIM and GAIA will be able to do a superb job in measuring the shifts in the light centroids of almost all selected microlensing events, they are scheduled for launch in 2009 and in $`20102012`$, respectively. This is a long time to wait. In this paper I am most interested in a possibility of identifying black hole lenses, and in measuring their masses in the next year or two. The best candidates for such studies are the longest microlensing events, which are likely to have the largest Einstein radii, probably in the range of several milli arc seconds. Small angle astrometry can be done with the HST with a precision of a fraction of a milliarcseconds (cf. Fritz et al. 2002, and references therein). Perhaps modern large telescope can also achieve such precision. The advantage of this approach is that the brightness of the lensed star is not a problem, but many images covering a substantial time interval have to be obtained. As we are focusing on long duration events, the time interval may be several years.

The following is an illustration of what may be expected. For simplicity I assume that the impact parameter of the lensing event is very small, i.e. $`u_{min}1`$, and hence the maximum magnification will be very large, $`A_{max}1`$ (cf. eqs. 4, 5). This implies that the source moves along a line in the sky which passes almost exactly behind the lens, the relative proper motion is one-dimensional, and so is the shift of the light centroid. This makes the presentation of motion very simple. The time dependence of centroid displacement with respect to the source is shown with solid line in Fig. 1. The maximum range is $`0.707\phi _E`$, but it takes time interval of $`2.83t_E`$, as given with the eqs. (9). Unfortunately, we do not know the asymptotic motion of the source, unperturbed by lensing. To learn this we would have to monitor it for a very long time indeed, as the solid line in Fig. 1 approaches the unperturbed source position very slowly, as $`1/t`$. What can be measured within a reasonable time interval is the non-linear motion of the light centroid. In most cases of interest the astrometric parallax effect due to Earth orbital motion is negligible (Figs. 1 and 2 in Paczyński 1998), so it is sufficient to calculate the departure from linear motion based on Fig. 1.

Let us assume that astrometric measurements are done over time interval $`\mathrm{\Delta }t_{1,2}`$, beginning at some time $`t_1`$, and continued till $`t_2`$. For a given value of $`\mathrm{\Delta }t_{1,2}`$ we want to select the beginning of astrometric observations $`t_1`$, to maximize the deviation of centroid position from a linear proper motion. For example, the two dots on the solid curve in Fig. 1, connected with a straight solid line represent an example for $`\mathrm{\Delta }t_{1,2}=2t_E`$. The vertical dashed line represents the maximum value of departure from linear motion, in this case it is $`0.18\phi _E`$.

More general results are shown in Fig. 2, where the beginning and the end of observations, $`t_1`$ and $`t_2`$ are shown with dashed lines as a function of $`\mathrm{\Delta }t_{1,2}`$, all in units of Einstein time scale $`t_E`$. The maximum displacement with respect to linear proper motion in the considered time interval is shown with thick solid line, with the value of $`\delta \phi _{max}`$ multiplied by a factor 10. It is interesting that $`\delta \phi _{max}`$ is approximately a linear function of $`\mathrm{\Delta }t_{1,2}`$.

Fig. 2 presents a practical problem with astrometry of centroid shift of long microlensing events: in order to detect a substantial non-linearity of centroid proper motion it is necessary to monitor it for a long time. For example, if $`t_E=1.5`$ years, and $`\phi _E=5`$ mas, then monitoring the motion over 3 year time interval provides us with the maximum displacement $`\delta \phi _{max}=0.18\times \phi _E=0.9`$ mas. It is convenient that astrometry has to begin not much earlier than the peak magnification. By that time all photometric parameters of the microlensing event are fairly well known,

### 2.3 Long Current Events

It is fortunate that to achieve the largest non-linear shift in the light centroid there is no need to begin astrometry while the microlensing event is not fully developed. The OGLE EWS alert is usually issued when the apparent magnitude of a candidate event is 0.2 – 0.3 mag above the baseline. At this early stage it is not always clear how long the event is going to be, what will be its peak magnification, and even will it be a microlensing event at all. As explained in previous section it is reasonable to begin accurate astrometry some time prior to the peak magnification.

In the past the longest microlensing events: OGLE-1999-BUL-32, OGLE-SC5\_2859, and OGLE-1999-BUL-19, were all recognized as very long well past their peak, analyzing archive data. There were two main reason for this time lag: improved photometric accuracy had to be developed, using new image subtraction software (DIA - Difference Image Analysis: Alard and Lupton 1998, Alard 2000, Woźniak 2002), and it had to be recognize that events with a time scale well over a year existed. The DIA is now used in real time data analysis, and it is known that very long events exist, so the task of early recognition of the very long events is now easy.

The following is the current status of OGLE EWS results. The unfolding OGLE-2003-BLG-192 has barely begun its rise, with the time scale estimate $`t_E235`$ days. The lensed star is bright: $`\mathrm{I}=16`$ mag, but the expected peak magnification is low. At this time it cannot be guaranteed that this is a microlensing event.

There are several other long events just beginning their rise, with $`I16`$ mag, or brighter: OGLE-2003-BLG-047 with $`t_E=156`$ days, OGLE-2003-BLG-188 with $`t=121`$ days. OGLE-2002-BLG-360, with $`t_E=270`$ days, is close to its peak at $`I=14.2`$. Two bright and long events are almost over: OGLE-2002-BLG-061 with $`t_E=305`$ days, and OGLE-2002-BLG-334 with $`t_E=160`$ days. Many more long events are faint at their baseline, heavily blended and therefore not very useful. However, several events every year are long and reasonably bright, with $`I16`$ mag. This may be time to propose follow-up astrometric observations for some of these events. Most lensed stars are in the Galactic Bulge and they are red, with $`(IK)2`$ mag, i.e. the bright events have $`K<14`$ mag.

## 3 . Planets

Another high profile discovery which is likely to result from microlensing is the first firm detection of a planetary signal, as proposed by Mao and Paczyński (1991). There is a rich literature on the subject, devoted to various aspects of planetary lensing, and how to recognize it. As planetary events are even less probable than stellar, a lot of theoretical effort was directed to increase the probability, at the cost of seeking low amplitude signals. This may be a good approach, but not at the beginning of a search. Right now we need credible events with no ambiguity about their interpretation. In other words, one needs an ‘overkill’ evidence to be persuasive the first time, and perhaps even the first several times. Such clear planetary signal will be present in caustic crossing events, like those shown in Fig. 10 of Paczyński (1996a). Even Earth mass planet can generate a complicated disturbance in a stellar microlensing light curve, with an amplitude of 20% or more. Such events are far less frequent than 5% or 2% amplitude planetary disturbances, but they offer unambiguous evidence of a planet, and an unambiguous mass ratio. Also, as any caustic crossing event, they will provide the information about the relative proper motion $`\dot{\phi }`$ (Gould 1994, Graff and Gould 2002).

While any strong planetary event is likely to provide the mass ratio and a relative proper motion, it is necessary to measure also photometric parallax effect to determine the mass. This may be possible in some cases, when the stellar event time scale $`t_E`$ is long, but in general this can be done only with spacecraft observations (Gould 1995, Gould et al. 2003), as they provide a baseline for the photometric parallax determination. Following the analysis of a possible planetary signal from OGLE-2002-BLG-055 (Jaroszyński and Paczyński 2002) the EWS alert system on OGLE has been upgraded to EEWS (Early Early Warning System), capable of recognizing and verifying in real time possible planetary disturbances in stellar microlensing events (A. Udalski, 2003, private communication). No structure is expected on a time scale shorter than $`1`$ hour, as it takes so long for a source star to move by its own radius. Therefore, instant follow-up observations at $`30`$ minute intervals are planned by the OGLE when a plausible planetary disturbance will be recognized. The new system should be much more sensitive to planet detection that the past attempts, which generated only upper limits (e.g. Gaudi et al. 2002).

The most natural follow-up observations are photometric measurements to be done from other location, to provide full time coverage. There are several other active microlensing projects: MOA (Mond et al. 2001), PLANET (Albrow et al. 2002), MPS (Rhie et al. 1999), GMAN (Becker et al. 1997, Bennett et al. 2002a), MicroFUN (A. Gould, 2003, private communication)

http://www.astronomy.ohio-state.edu/~microfun

This coverage should be adequate to obtain the mass ratio, a clear indicator of a planetary presence. Spacecraft observations are needed to determine the masses, as well as the linear planet - star separations (in Astronomical Units).

## 4 . Conclusions

It is very likely that stellar mass black holes and planets will be discovered in the next year or two. OGLE is generating several bright and long duration microlensing events every year (cf. Section 2.3), and these are promising candidates for stellar mass black holes. All these events have not only their times scales $`t_E`$ measured, but also the photometric parallax effect. To make a definite case for black holes it is necessary to measure their angular Einstein radii, $`\phi _E`$. This could be done best with an optical interferometer, possibly with VLTI (Delplancke et al. 2001). An alternative approach is to carry out accurate astrometric monitoring and to detect a non-linear motion of the light centroid with the HST, or perhaps even with new large ground based telescopes. A lens with $`10M_{}`$ is likely to have $`\phi _E3`$ mas. If intermediate mass black holes exist (Madau and Rees 2001, Chisholm et al. 2002, Wu et al. 2002, and references therein) the corresponding image splittings could be much larger, with $`\phi _E10`$ mas, and perhaps even more. Such events would make image centroid motion easily detectable not only with the HST, but also with ground based telescopes.

Once several firm planetary detections are made we shall be in a much better position to plan future upgrades of planetary microlensing searches. It is likely that the detection threshold could be lowered to improve the statistics without compromising credibility. But a major improvement will require a major increase in the photometric data rate, to monitor more stars more frequently and with a higher photometric accuracy. Note: vast majority of microlensing events are heavily blended and are relatively faint. It will take a lot of photons to reliably search them for small planetary disturbances.

While OGLE and other ground based telescopes will determine the mass ratios for star - planet systems, and the source - lens proper motion in caustic crossing events, only some of them will have large enough $`t_E`$ to allow photometric parallax to be measured from the ground (cf. Jaroszyński and Paczyński 2002), and to determine the masses as well. In most cases it will take space instrument to provide a large baseline to measure photometric parallax effect, and hence the masses for most of them (Gould 1995, Gould et al. 2003).

This paper is posted on astro-ph only. Its purpose is to outline prospects for various space and ground based projects which are needed to provide definite evidence for stellar mass black holes and planets, using OGLE as a search tool.

Acknowledgements. This research was supported by the NASA grant NAG5-12212 and the NSF grant AST-0204908.

## REFERENCES

* Agol, E. et al. 2002, Astrophys. J., 576, L131.
* Alard, C. 2000, Astron. Astrophys. Suppl. Ser., 133, 363.
* Alard, C., and Lupton, R. H. 1998, Astrophys. J., 503, 325.
* Albrow, M. D. et al. (PLANET) 2002, Astrophys. J., 572, 1031.
* An, J. H., et al. 2002, Astrophys. J., 572, 521.
* Becker, A. et al. (GMAN) 1997, Bull. AAS, 29, 1347.
* Bennett, D. P. et al. 2002a, Astrophys. J., 579, 639.
* Bennett, D. P. et al. 2002b, astro-ph/0207006, , .
* Boden, A. F., Shao, M., and van Buren, D. 1998, Astrophys. J., 502, 538.
* Chisholm, J. R., Dodelson, S., and Kolb, E. W. 2002, astro-ph/0205138, , .
* Delplancke, F., Górski, K., and Richichi, A. 2001, Astron. Astrophys., 375, 701.
* Bond, I. A. et al. (MOA) 2001, MNRAS, 327, 868.
* Fritz, B. G. et al. 2002, Astron. J., 124, 1695.
* Gaudi, B. S. et al. (PLANET) 2002, Astrophys. J., 566, 463.
* Gould, A. 1992, Astrophys. J., 392, 442.
* Gould, A. 1994, Astrophys. J. Letters, 421, L71.
* Gould, A. 1995, Astrophys. J. Letters, 441, L21.
* Gould, A., and Loeb, A. 1992, Astrophys. J. Letters, 396, 104.
* Gould, A., Gaudi, B. S., and Han, C. 2003, astro-ph/0304314, , .
* Gould, A. 2001, P.A.S.P., 113, 903.
* Graff, D. S., and Gould, A. 2002, Astrophys. J., 580, 253.
* Hog, E., Novikov, I. D., and Polnarev, A. G. 1995, Astron. Astrophys., 294, 287.
* Jaroszyński, M., and Paczyński, B 2002, Acta Astron., 52, 361.
* Kervella, P. et al. 2003, astro-ph/0303634, , .
* Madau, P., and Rees, M. 2001, Astrophys. J. Letters, 551, L27.
* Mao, S., and Paczyński, B. 1991, Astrophys. J. Letters, 374, L37.
* Mao, S., et al. (OGLE) 2002, MNRAS, 329, 349.
* Miyamoto, M., and Yoshi, Y. 1995, Astron. J., 110, 1427.
* Paczyński, B. 1996a, Ann. Rev. Astron. Ap., 34, 419.
* Paczyński, B. 1996b, Acta Astron., 46, 291.
* Paczyński, B. 1998, Astrophys. J. Letters, 494, L23.
* Refsdal, S. 1966, MNRAS, 134, 315.
* Rhie, S. H. et al. 1999, Astrophys. J., 522, 1037.
* Segransan, D. at al. 2003, Astron. Astrophys., 397, L5.
* Smith, M. C. 2003, astro-ph/0304442, , .
* Smith, M. C. et al. (OGLE) 2002, MNRAS, 336, 670.
* Udalski, A., Kubiak. M., and Szymański, M. 1997, Acta Astron., 47, 319.
* Udalski, A. et al. (OGLE) 2002, Acta Astron., 52, 1.
* Walker, M. A. 1995, Astrophys. J., 453, 37.
* Woźniak, P. R. 2000, Acta Astron., 50, 421.
* Woźniak, P. R. et al. (OGLE) 2001, Acta Astron., 51, 175.
* Wu, H. et al. 2002, Astrophys. J., 576, 738.


text:
# Reconstructing Three-dimensional Structure of Underlying Triaxial Dark Halos From X-ray and Sunyaev–Zel’dovich Effect Observations of Galaxy Clusters

## 1 INTRODUCTION

Understanding of statistical properties of dark matter halos has been significantly advanced in recent years, largely owing to the development of high-resolution numerical simulations (Navarro, Frenk, & White, 1996, 1997; Fukushige & Makino, 1997; Moore et al., 1998; Jing & Suto, 2000, 2002). Given those theoretical/empirical successes, a next natural question is how to apply them for the description of real galaxy clusters. While there exist a number of attempts along this line, they usually make the unrealistic assumption of the halo spherical symmetry (Makino, Sasaki & Suto, 1998; Suto, Sasaki & Makino, 1998; Yoshikawa & Suto, 1999). This paper describes a methodology to reconstruct the three-dimensional (3D) structure of dark halos from the two-dimensional (2D) surface brightness profiles of intra-cluster gas from the X-ray and/or Sunyaev-Zel’dovich (SZ) effect observations assuming that the halos are triaxial ellipsoids with a constant axis ratio.

Indeed it is a classical problem in astrophysics to determine the 3D properties of astronomical objects from the observed 2D counterparts (e.g., Lucy, 1974; Fabricant et al., 1984; Dehnen & Gerhard, 1993; Binney, Davies, & Illingworth, 1990; Gerhard & Binney, 1996). Deprojecting galaxy clusters is one of the most important applications. Fabricant et al. (1984) analyzed the X-ray surface brightness map of clusters, and showed that their mass distribution is far from being spherical. Zaroubi et al. (1998) developed a general method of deprojecting the 2D images of rich clusters based on the Fourier slice theorem to reconstruct the 3D cluster structures. Later their technique was tested against numerically simulated galaxy clusters (Zaroubi et al., 2001). Yoshikawa & Suto (1999) proposed a deprojection method for spherical clusters based on the Abel integral. Reblinsky (2000) provided a parameter-free Richardson-Lucy algorithm to reconstruct the 3D halo potential, and demonstrated its stability by applying it to gas-dynamical simulations. Dore et al. (2001) provided a perturbative approach to the cluster deprojection, taking into account the non-isothermality and asphericity of galaxy clusters. Recently, Fox & Pen (2002) also considered the problem of deprojecting aspherical clusters. They first constructed a parameterized 3D axisymmetric cluster model, and determined the 3D cluster shapes through the $`\chi ^2`$-fitting between the model predictions and the simulated cluster data.

All the previous approaches, however, were based on rather restrictive assumptions such as the cluster axial symmetry (e.g., oblate or prolate) about the line-of-sight, and/or the isothermality of galaxy clusters. Furthermore in their approaches, it was required to combine X-ray and/or SZ effect data with weak lensing (WL) map, which significantly limits the applicability of the previous methods. In fact, it is generally believed that it would be a difficult task to deproject the clusters from the 2D projected X-ray or SZ observables alone without having such restrictive assumptions and priors given the degeneracy of the parameters due to the the projection process itself.

Lee & Suto (2003; hereafter Paper I) suggested that one may understand the properties of the dark matter halos from the intra-cluster gas profiles. Assuming that the intra-cluster gas is in hydrostatic equilibrium in the triaxial dark matter halos, Paper I derived the 3D density and temperature profiles of the intra-cluster gas from the 1st principles using the perturbation theory, and found an analytic relation between the eccentricities of the intra-cluster gas and the underlying dark halo (see eq. in the Paper I). However, the perturbation result of the Paper I is valid only in the asymptotic limit of the low asphericity, and did not include the effect of the projection of the gas profiles on the plane of the sky. Here, we generalize the works of the Paper I to the case of highly aspherical clusters, and attempt to find a general relation between the dark halos and the intra-cluster gas taking into account the parameter-degeneracy caused by the projection. For the practical purpose of our theoretical modeling, we provide a series of empirical fitting formulae for the 2D and 3D gas profiles which can be used as templates in comparing with the observed profiles of clusters, and demonstrate the degree of the feasibility of the halo reconstruction using the template formulae and the simulated surface brightness maps of X-ray and SZ effect.

The organization of this paper is as follows. In §2 we construct a theoretical model for the intra-cluster gas in hydrostatic equilibrium within the gravitational potential of a triaxial dark matter halo with a constant axis ratio, and provide a series of empirical fitting formulae for the 3D gas profiles in a systematic manner. In §3 we consider the projection of the 3D gas profiles onto the plane of the sky, and provide another series of empirical fitting formulae for the 2D surface brightness density profiles. In §4 we describe a numerical algorithm to determine the 3D structures of the dark halos from the observed 2D cluster surface brightness profiles, and test the algorithm against a numerical toy model. Finally we discuss the final results, and draw our conclusions in §5.

## 2 MODELING INTRA-CLUSTER GAS PROFILES

### 2.1 Gravitational Potential of Dark Matter Halos

To predict theoretically the X-ray and SZ profiles of galaxy clusters, one first needs a good physical model for the intra-cluster gas. Here we assume that the intra-cluster gas is either isothermal or polytropic, and in hydrostatic equilibrium within the gravitational potential generated by a concentric and coaxial triaxial dark matter halo.

Consider a triaxial dark halo whose iso-density surfaces are given by the following equation:

$`R^2x^2+{\displaystyle \frac{y^2}{1e_b^2}}+{\displaystyle \frac{z^2}{1e_c^2}},`$ (1)
where the Cartesian system of coordinates $`(x,y,z)`$ is aligned with the halo principal axes, oriented such that the $`x`$-axis and $`z`$-axis run along the major and minor principal axes, respectively. The major axis length of the iso-density surface is denoted by $`R`$, while $`e_b`$ and $`e_c`$ ($`e_b<e_c`$) represent the two constant eccentricities of the ellipsoidal dark halos. Equation (1) implies that the density profile of a triaxial dark matter halo should be a function of the major axis length $`R`$.

We adopt the density profile of a triaxial dark halo proposed by Jing & Suto (2002):

$$\rho (R)=\frac{\delta _c\rho _{\mathrm{crit}}}{\left(R/R_0\right)^\alpha \left(1+R/R_0\right)^{3\alpha }},$$
(2)
where $`R_0`$ is the scale length, $`\delta _c`$ is the dimensionless characteristic density contrast with respect to the critical density $`\rho _{\mathrm{crit}}`$ of the universe at the present epoch, and $`\alpha `$ represents the inner slope of the density profile. Jing & Suto (2002) showed that $`\alpha 1`$ on the cluster scale and $`\alpha 3/2`$ on the galaxy scale. For simplicity and definiteness, we focus on the case of $`\alpha =1`$ throughout this paper.

The gravitational potential due to the ellipsoidal halos described by equation (2) is formally expressed as (Binney & Tremaine, 1987):

$`\mathrm{\Phi }(𝐫)`$ $`=`$ $`\pi G\sqrt{(1e_b^2)(1e_c^2)}{\displaystyle _0^{\mathrm{}}}{\displaystyle \frac{[\psi (\mathrm{})\psi (m)]}{\sqrt{(\tau +1)(\tau +1e_b^2)(\tau +1e_c^2)}}}𝑑\tau ,`$ (3)
$`\psi (m)`$ $`=`$ $`2{\displaystyle _0^m}\rho (R)R𝑑R,m^2={\displaystyle \frac{x^2}{\tau +1}}+{\displaystyle \frac{y^2}{\tau +1e_b^2}}+{\displaystyle \frac{z^2}{\tau +1e_c^2}}.`$ (4)
Equations (3) and (4) along with equation (2) allow one to compute the triaxial halo gravitational potential at least numerically. The halo gravitational potential $`\mathrm{\Phi }(𝐫)`$ depends on $`\delta _c`$, $`\rho _{\mathrm{crit}}`$,$`R_0`$ as well as $`e_b`$, $`e_c`$, $`\alpha `$. However, the dependence of $`\mathrm{\Phi }(𝐫)`$ on $`\delta _c`$, $`\rho _{\mathrm{crit}}`$, and $`R_0`$ can be separated out by introducing a dimensionless potential $`\stackrel{~}{\mathrm{\Phi }}`$ that depends only on $`e_b`$, $`e_c`$ and $`\alpha `$ such that

$$\stackrel{~}{\mathrm{\Phi }}(𝐫;e_b,e_c,\alpha )\frac{\mathrm{\Phi }(𝐫;e_b,e_c,\alpha ,\delta _c,\rho _{\mathrm{crit}},R_0)}{4\pi G\delta _c\rho _{\mathrm{crit}}R_0^2}.$$
(5)
### 2.2 Gas Density and Temperature Profiles in terms of Halo Potential

To determine the intra-cluster gas profiles in terms of the halo potential, it is necessary to specify the equation of state for the intra-cluster gas. We consider both the isothermal and the polytropic cases in order.

#### 2.2.1 Isothermal gas

The equation of state for the isothermal gas is given as

$$P_\mathrm{g}(𝐫)=P_{\mathrm{g0}}\frac{\rho _\mathrm{g}(𝐫)}{\rho _{\mathrm{g0}}},$$
(6)
where $`P_\mathrm{g}`$ and $`\rho _\mathrm{g}`$ represent the gas pressure and the gas density, respectively. In what follows, the subscript $`0`$ of some physical variable indicates the value of that physical variable at the center $`𝐫=0`$.

The density profile of the isothermal gas in hydrostatic equilibrium is given as

$$\rho _\mathrm{g}(𝐫)=\rho _{\mathrm{g0}}\mathrm{exp}\left[\kappa \{\stackrel{~}{\mathrm{\Phi }}(𝐫)\stackrel{~}{\mathrm{\Phi }}_0\}\right].$$
(7)
We define a dimensionless isothermal gas constant $`\kappa `$ as

$$\kappa \frac{4\pi G\mu _gm_p\delta _c\rho _{\mathrm{crit}}R_0^2}{k_BT_\mathrm{g}},$$
(8)
where $`m_p`$ is the proton mass, $`\mu _g`$ is the mean molecular weight of the intra-cluster gas, $`k_B`$ is the Boltzmann constant, and $`T_\mathrm{g}`$ is the (constant) gas temperature. Introducing $`F_\mathrm{\Phi }(𝐫)`$:

$$F_\mathrm{\Phi }(𝐫)\mathrm{exp}\left[\{\stackrel{~}{\mathrm{\Phi }}(𝐫)\stackrel{~}{\mathrm{\Phi }}_0\}\right],$$
(9)
one can rewrite equation (7) as

$$\rho _\mathrm{g}(𝐫)=\rho _{\mathrm{g0}}[F_\mathrm{\Phi }(𝐫)]^\kappa .$$
(10)
#### 2.2.2 Polytropic gas

The equation of state for the polytropic gas with the polytropic index of $`\gamma `$ $`(1)`$ is given as

$$P_\mathrm{g}(𝐫)=P_{\mathrm{g0}}\left[\frac{\rho _\mathrm{g}(𝐫)}{\rho _{\mathrm{g0}}}\right]^\gamma ,$$
(11)
and its density and temperature profiles in hydrostatic equilibrium are found as 111Note that the definition of $`\mathrm{\Phi }_0`$ in Paper I is different from that given here by a constant offset.

$$\rho _\mathrm{g}(𝐫)=\rho _{\mathrm{g0}}\left[1\kappa _p\{\stackrel{~}{\mathrm{\Phi }}(𝐫)\stackrel{~}{\mathrm{\Phi }}_0\}\right]^{1/(\gamma 1)},T_\mathrm{g}=T_{\mathrm{g0}}\left[1\kappa _p\{\stackrel{~}{\mathrm{\Phi }}(𝐫)\stackrel{~}{\mathrm{\Phi }}_0\}\right].$$
(12)
We define a dimensionless polytropic gas constant $`\kappa _p`$ as

$$\kappa _p\frac{\gamma 1}{\gamma }\frac{4\pi G\mu _gm_p\delta _c\rho _{\mathrm{crit}}R_0^2}{k_BT_{\mathrm{g0}}},$$
(13)
and introduce $`F_\mathrm{\Phi }(𝐫)`$:

$$F_\mathrm{\Phi }(𝐫)1[\stackrel{~}{\mathrm{\Phi }}(𝐫)\stackrel{~}{\mathrm{\Phi }}_0].$$
(14)
Then one can rewrite equation (12) as

$$\rho _\mathrm{g}(𝐫)=\rho _{\mathrm{g0}}\left[1\kappa _p+\kappa _pF_\mathrm{\Phi }(𝐫)\right]^{1/(\gamma 1)},T_\mathrm{g}(𝐫)=T_{\mathrm{g0}}\left[1\kappa _p+\kappa _pF_\mathrm{\Phi }(𝐫)\right].$$
(15)
### 2.3 Empirical Fitting Formulae for the Gravitational Potential

We have shown in $`\mathrm{\S }2.2`$ that the potential function $`F_\mathrm{\Phi }(𝐫)`$ defines the density and the temperature profiles of intra-cluster gas completely. In general, $`F_\mathrm{\Phi }(𝐫)`$ does not have a closed analytic form when the dark matter halos are triaxial ellipsoids. In Paper I, we computed $`F_\mathrm{\Phi }(𝐫)`$ analytically with the perturbation theory assuming $`e_b^2e_c^21`$, and related the iso-density surfaces of the intra-cluster gas to those of the dark halos; if the halo iso-density surfaces are triaxial ellipsoids with the eccentricities of $`e_\sigma `$ ($`\sigma =b,c`$), then the iso-density surfaces of the intra-cluster gas is also well approximated as triaxial ellipsoids with the eccentricities of $`ϵ_\sigma `$ that are related to $`e_\sigma `$ by

$$\frac{ϵ_\sigma ^2}{e_\sigma ^2}=\frac{6(1+u)\mathrm{ln}(1+u)+u^33u^26u}{2u^2[(1+u)\mathrm{ln}(1+u)u]},$$
(16)
where $`u|𝐫|/R_0`$. Note that the right-hand-side of equation (16) are written in terms of the rescaled spherical radius $`u`$ instead of the rescaled major axis length, which results from the fact that we neglected the higher-order terms $`O(e_\sigma ^4)`$ in deriving equation (16) with the perturbation theory. Although equation (16) seems complicated, it is a slowly varying function, and well approximated as a constant in the range of $`0<u<1`$ (see Fig. 3 of Paper I) such that

$$\frac{ϵ_\sigma ^2}{e_\sigma ^2}0.7^2.$$
(17)
Strictly speaking, equation (17) is valid only in the limit of $`e_b^2e_c^21`$ over the range of $`0<u<1`$. To proceed further in a tractable fashion, however, we extend the validity of equation (17) to $`e_b^2e_c^21`$ over the whole range of $`u`$. In other words, we assume that the iso-density surfaces of the intra-cluster gas are triaxial ellipsoids with the constant eccentricities of $`ϵ_b,ϵ_c`$ that are related to the halo eccentricities $`e_b,e_c`$ by equation (17). Actually this turns out to be a good approximation as will be shown below.

This assumption allows us to write the potential function $`F_\mathrm{\Phi }(𝐫)`$ in terms of the major axis length of the gas iso-density surfaces, $`\xi `$, defined as

$$\xi ^2u_x^2+\frac{u_y^2}{1ϵ_b^2}+\frac{u_z^2}{1ϵ_c^2}=\frac{1}{R_0^2}\left(x^2+\frac{y^2}{1ϵ_b^2}+\frac{z^2}{1ϵ_c^2}\right).$$
(18)
After some trials and errors, we find that the following empirical formula fits both equations (9) and (14) quite well:

$$F_\mathrm{\Phi }(\xi )=\left(\frac{1+\eta \xi ^p}{1+\beta \xi ^p}\right)^q,$$
(19)
where $`\beta `$, $`p`$, $`q`$, and $`\eta `$ are free parameters. The four free parameters are not constant but supposed to be functions of the halo eccentricities. We determine the functional form of each parameter by taking the following steps.

1. Compute $`F_\mathrm{\Phi }(𝐫)`$ numerically using equations (3) and (4) for various values of $`e_b`$ and $`e_c`$.
2. Fit those numerical data points to the empirical model (eq. ) and obtain the best-fit values of $`\beta `$, $`p`$, $`q`$, and $`\eta `$ using the Levenberg-Marquardt method (Press et al., 1992).
3. Finally model those best-fit parameters as functions of $`e_b`$ and $`e_c`$. In practice, we find that all the parameters can be written as functions of a single variable $`\mu e_b^3+e_c^3`$.

We find the following polynomials to provide good fits:

$`\beta `$ $`=`$ $`c_{\beta 0}+c_{\beta 1}\mu ,`$ (20)
$`p`$ $`=`$ $`const.,`$ (21)
$`q`$ $`=`$ $`c_{q0}+c_{q1}\mu +c_{q2}\mu ^2,`$ (22)
$`\eta `$ $`=`$ $`c_{\eta 0}+c_{\eta 1}\mu .`$ (23)
Table 1 lists the best-fit values of the polynomial coefficients. The best-fit values of the constant $`p`$ is determined to be unity for all cases. Figure 1 plots $`\beta `$, $`q`$, and $`\eta `$ as a function of $`\mu `$. The filled circles indicate the best-fit values of the parameters, while the solid lines represent their polynomial fits. Note that $`\eta `$ is one order-of-magnitude smaller than $`\beta `$ for all cases, suggesting that the term associated with $`\eta `$ can be safely neglected except at the outer part of clusters. The accuracy of the fits is illustrated in Figure 2 where $`F_\mathrm{\Phi }(\xi )`$ is plotted against $`\xi `$. The solid lines represent the numerical results while the dashed lines correspond to equation (19) with the best-fit coefficients of the polynomials plotted in Fig. 1 and listed in Table 1. Our empirical fitting formulae reproduce the numerical results very well for all cases over a wide range of $`\xi `$ within a fractional error less than $`0.5\%`$.

The characteristic feature of our fitting model is that all the parameters are expressed in terms of $`\mu =e_b^3+e_c^3`$. Given that the first-order perturbation approach of the Paper I indicated the dependence of the halo potential on $`e_b^2+e_c^2`$, one may have expected the fitting model to depend on $`e_b^2+e_c^2`$ rather than $`e_b^3+e_c^3`$. As a matter of fact, we indeed attempted first to model the parameters in terms of $`e_b^2+e_c^2`$. But, it turned out that the parameter values do not scale in terms of $`e_b^2+e_c^2`$ while they show very good scaling feature in terms of $`e_b^3+e_c^3`$ as shown in Figure 1. We have not yet completely understood the origin of this scaling, but the result is empirically quite robust.

In addition, we note that in the Paper I we derived the gas-halo eccentricity relation (eq. ) by using the first-order perturbation theory, and showed that the approximation error can be expressed as a function of $`e_b^3`$ and $`e_c^3`$ (see eqs  and  in the Paper I). The dependence of the approximation error on $`e_b^3`$ and $`e_c^3`$ implies that the higher order terms neglected in the first-order perturbation result may be scaled in terms of $`e_b^3`$ and $`e_c^3`$. Now that we have used equation (16) to construct the fitting model (eq. ), one may expect that equation (19) has similar scaling. In other words, we suspect that the dependence of equation (19) on $`e_b^3+e_c^3`$ rather than $`e_b^2+e_c^2`$ may be related to the scaling of the higher order terms neglected in the first-order perturbation result of the Paper I.

Note also that in our model the gas density/temperature profiles do not approach zero even at large $`u`$. This feature is ascribed to the hydrostatic equilibrium condition itself. In reality, the hydrostatic equilibrium condition may be satisfied only within some radius, e.g., the virial radius of the cluster. In fact, the halo density profiles at those regions are not well approximated by equation (2) because of the presence of substructures, and so on. This problem, however, is unlikely to limit the validity of our methodology in practice since the X-ray and SZ fluxes from regions beyond the cluster virial radius are usually negligible.

## 3 PROJECTED SURFACE BRIGHTNESS PROFILES OF INTRA-CLUSTER GAS

Now we are in a position to model the X-ray and SZ surface brightness profiles of triaxial galaxy clusters by incorporating the projection effect into the 3D models developed in $`\mathrm{\S }2`$. Just for comparison, let us recall that the surface brightness profiles of spherical clusters can be readily evaluated as

$$\mathrm{\Sigma }(\theta )=_{\mathrm{}}^{\mathrm{}}L(\sqrt{d_A^2\theta ^2+z^2})𝑑z,$$
(24)
where $`\theta `$ is the angular radius from the center of the cluster, $`d_A`$ is the angular diameter distance to the cluster, and $`L`$ is the emissivity given in terms of $`\rho _g`$ and $`T_g`$ such that $`L\rho _\mathrm{g}^2T_\mathrm{g}^{1/2}`$ for bolometric X-ray and $`L\rho _\mathrm{g}T_\mathrm{g}`$ for SZ observations, respectively. For the case of a triaxial cluster, however, the result becomes much more complicated because the 2D projection depends on the relative direction of the line-of-sight of an observer with respect to the principal axes of a triaxial halo. The projection effect of triaxial bodies on the plane of the sky is fully discussed by Stark (1977) and Binney (1985). In what follows we adopt the notation of Binney (1985).

Let $`(\theta ,\varphi )`$ be the polar angle of the line-of-sight in the halo principal coordinate system $`(x,y,z)`$, and let the observer’s coordinate system be defined by Cartesian axes $`(x^{},y^{},z^{})`$ with $`z^{}`$-axis aligned with the line-of-sight direction and $`x^{}`$-axis lying in the $`(x,y)`$ plane (see Fig. 1 in Oguri, Lee, & Suto, 2003). Then the halo principal coordinate system $`(x,y,z)`$ is related to the observer’s coordinate system $`(x^{},y^{},z^{})`$ by

$`x`$ $`=`$ $`\mathrm{sin}\varphi x^{}\mathrm{cos}\varphi \mathrm{cos}\theta y^{}+\mathrm{cos}\varphi \mathrm{sin}\theta z^{},`$ (25)
$`y`$ $`=`$ $`\mathrm{cos}\varphi x^{}\mathrm{sin}\varphi \mathrm{cos}\theta y^{}+\mathrm{sin}\varphi \mathrm{sin}\theta z^{},`$ (26)
$`z`$ $`=`$ $`\mathrm{sin}\theta y^{}+\mathrm{cos}\theta z^{}.`$ (27)
If we use the major axis length $`\xi `$ defined in equation (18) to characterize the iso-density surfaces, then the projection of $`L`$ onto the plane of the sky is written as

$$\mathrm{\Sigma }(x^{},y^{})_{\mathrm{}}^{\mathrm{}}L(\xi ^2)𝑑z^{}=\frac{2}{\sqrt{f}}_0^{\mathrm{}}L(z^{\prime \prime 2}+\lambda ^2)𝑑z^{\prime \prime },$$
(28)
where

$`z^{\prime \prime }`$ $`=`$ $`\sqrt{f}\left(z^{}+{\displaystyle \frac{g}{2f}}\right),`$ (29)
$`\lambda `$ $`=`$ $`{\displaystyle \frac{1}{\sqrt{f}}}(Ax^2+Bx^{}y^{}+Cy^2)^{1/2},`$ (30)
$`f`$ $`=`$ $`\mathrm{sin}^2\theta \left(\mathrm{cos}^2\varphi +{\displaystyle \frac{\mathrm{sin}^2\varphi }{1ϵ_b^2}}\right)+{\displaystyle \frac{\mathrm{cos}^2\theta }{1ϵ_c^2}},`$ (31)
$`g`$ $`=`$ $`\mathrm{sin}\theta \mathrm{sin}2\varphi \left({\displaystyle \frac{1}{1ϵ_b^2}}1\right)x^{}+\mathrm{sin}2\theta \left({\displaystyle \frac{1}{1ϵ_c^2}}\mathrm{cos}^2\varphi {\displaystyle \frac{\mathrm{sin}^2\varphi }{1ϵ_b^2}}\right)y^{},`$ (32)
$`A`$ $`=`$ $`{\displaystyle \frac{\mathrm{cos}^2\theta }{1ϵ_c^2}}\left(\mathrm{sin}^2\varphi +{\displaystyle \frac{\mathrm{cos}^2\varphi }{1ϵ_b^2}}\right)+{\displaystyle \frac{\mathrm{sin}^2\theta }{1ϵ_b^2}},`$ (33)
$`B`$ $`=`$ $`\mathrm{cos}\theta \mathrm{sin}2\varphi \left(1{\displaystyle \frac{1}{1ϵ_b^2}}\right){\displaystyle \frac{1}{1ϵ_c^2}},`$ (34)
$`C`$ $`=`$ $`\left({\displaystyle \frac{\mathrm{sin}^2\varphi }{1ϵ_b^2}}+\mathrm{cos}^2\varphi \right){\displaystyle \frac{1}{1ϵ_c^2}}.`$ (35)
Hence $`\mathrm{\Sigma }(x^{},y^{})=\mathrm{\Sigma }(\lambda )`$. In other words, the 2D isophotal curves of triaxial galaxy clusters are concentric and coaxial ellipses if their eccentricities $`ϵ_\sigma `$ of the gas iso-density surfaces are constants. As we noted before, however, it is not exactly the case for our model where the eccentricities $`e_\sigma `$ of the halo iso-density surfaces are constant. Nevertheless this holds approximately (see eq.).

Let us define a dimensionless surface brightness $`\stackrel{~}{\mathrm{\Sigma }}\mathrm{\Sigma }/\mathrm{\Sigma }_0`$ where $`\mathrm{\Sigma }_0\mathrm{\Sigma }(0)`$ and $`\stackrel{~}{\lambda }=\lambda /R_0`$. In the same spirit as for the 3D model (19), we propose the following empirical model for the 2D surface brightness profiles:

$$\stackrel{~}{\mathrm{\Sigma }}(\stackrel{~}{\lambda };e_b,e_c,\kappa ,\kappa _p,\gamma )=\left(\frac{1+\omega \stackrel{~}{\lambda }^s}{1+\zeta \stackrel{~}{\lambda }^s}\right)^t.$$
(36)
Note that $`\stackrel{~}{\mathrm{\Sigma }}`$ depends not only on $`\mu `$ but also on $`\kappa `$ (or $`\kappa _p`$), and even on the polytropic index $`\gamma `$ for the polytropic case. Thus, we consider the isothermal and polytropic cases separately, and obtain the following empirical fits.

For the isothermal gas, the parameters are fitted to the following polynomials of $`\mu `$ and $`\kappa `$:

$`\omega `$ $`=`$ $`{\displaystyle \underset{i,j=0}{\overset{2}{}}}\omega _{ij}\mu ^i\kappa ^j,`$ (37)
$`\zeta `$ $`=`$ $`{\displaystyle \underset{i,j=0}{\overset{2}{}}}\zeta _{ij}\mu ^i\kappa ^j,`$ (38)
$`s`$ $`=`$ $`const.,`$ (39)
$`t`$ $`=`$ $`{\displaystyle \underset{i,j=0}{\overset{2}{}}}t_{ij}\mu ^i\kappa ^j.`$ (40)
For the polytropic gas, we fix $`\gamma `$ to $`\gamma =1.15`$ following Komatsu & Seljak (2001), and fit the parameters to the following polynomials of $`\mu `$ and $`\kappa `$:

$`\omega `$ $`=`$ $`{\displaystyle \underset{i,j=0}{\overset{3}{}}}\omega _{ij}\mu ^i\kappa _p^j,`$ (41)
$`\zeta `$ $`=`$ $`{\displaystyle \underset{i,j=0}{\overset{3}{}}}\zeta _{ij}\mu ^i\kappa _p^j,`$ (42)
$`s`$ $`=`$ $`const.,`$ (43)
$`t`$ $`=`$ $`{\displaystyle \underset{i,j=0}{\overset{3}{}}}t_{ij}\mu ^i\kappa _p^j.`$ (44)
We determine the best-fit values of $`\omega _{ij}`$, $`\zeta _{ij}`$, and $`t_{ij}`$ using the same method as in $`\mathrm{\S }2`$. First we calculate $`\stackrel{~}{\mathrm{\Sigma }}`$ of equation (28) numerically, approximating the integration $`_{\mathrm{}}^{\mathrm{}}L𝑑z^{}`$ to $`_{z_c}^{z_c}L𝑑z^{}`$ where $`z_c=20R_0`$ which is roughly twice the virial radius of galaxy clusters (Makino, Sasaki & Suto, 1998). Then we compare the numerical data points of $`\stackrel{~}{\mathrm{\Sigma }}`$ with the fitting formula (36), and determine the best-fit values for each point using the Levenberg-Marquardt method. Finally we model the free parameters as polynomials of $`\mu `$ and $`\kappa `$ (or $`\kappa _p`$), and determine the best-fit polynomial coefficients.

The best-fit constant values of $`s`$ is determined to be unity for both the isothermal and polytropic cases just like the 3D model. Figures 3 and 4 plot the best-fit parameters as functions of $`\mu `$ for three different values of $`\kappa `$ (or $`\kappa _p`$). It is clear from these figures that the polynomial fitting works quite well for all cases. Tables 2 and 3 list the best-fit coefficients $`\zeta _{ij},t_{ij},\omega _{ij}`$ for the cases of bolometric X-ray and SZ observations, respectively.

Figure 5 illustrates the degree of the accuracy of our fitting formulae. The filled squares represent the numerical results while the solid lines represent our fitting models with the best-fit polynomials of $`\zeta `$, $`t`$, and $`\omega `$. We have found that our fits reproduce the numerical results within the fractional error of $`20\%`$ for all cases in a range of $`0\stackrel{~}{\lambda }10`$.

## 4 DARK HALO RECONSTRUCTION

In the triaxial dark halo model that we adopt here, the halo reconstruction from the X-ray and SZ cluster observations is basically to determine the shapes ($`e_b,e_c`$) and the orientations ($`\theta ,\varphi `$) of the dark halos from the observed surface brightness profiles.

In §3, we construct a parameterized model for the X-ray and SZ surface brightness profiles of galaxy clusters. The model is expressed as a function of the rescaled major axis length of the gas isophotes and characterized by three parameters. It turns out that the rescaled major axis length and the three parameters depend on the eccentricities and the orientation angles of the underlying dark halos, as well as the gas constant. Therefore, we find a way to link the observed 2D surface brightness profiles of galaxy clusters to the 3D structures of the underlying dark halo and the gas constants of the intra-cluster gas as well. In other words, one may expect to determine the values of $`(e_b,e_c,\theta ,\varphi )`$ and $`(\kappa ,\kappa _p)`$ by fitting our 2D parameterized models (eqs.  – ) to the X-ray and SZ cluster data.

To demonstrate how our modeling of cluster profiles can be used in the determination of the shapes and the orientations of dark matter halos, we apply the following reconstruction algorithm to the numerically projected profiles. We first numerically compute $`\stackrel{~}{\mathrm{\Sigma }}`$ directly using equations (3)-(28). Then we construct a pixeled map of surface brightness in $`N\times N`$ grids ($`N=32`$ in the present example) corresponding to the linear scale of $`10R_0x^{},y^{}10R_0`$ since $`10R_0`$ is roughly equal to the cluster virial radius (Jing & Suto, 2002). We create several realizations of X-ray and SZ profiles for both isothermal and polytropic cases using various different values of $`e_b,e_c,\theta ,\varphi `$.

Our reconstruction algorithm proceeds as follows;

* At each pixel point, say, $`(x_i^{},y_j^{})`$, build the model $`\stackrel{~}{\mathrm{\Sigma }}^{\mathrm{th}}(x_i^{},y_j^{})`$ using equation (36). The model is characterized by five free parameters $`e_b,e_c,\theta ,\varphi ,\kappa `$ (or $`\kappa _p`$) at each point.
* Fit $`\stackrel{~}{\mathrm{\Sigma }}^{\mathrm{th}}`$ to the observed surface brightness density profiles $`\stackrel{~}{\mathrm{\Sigma }}^{obs}`$, and calculate the $`\chi ^2`$:

$$\chi ^2=\underset{i,j=1}{\overset{N}{}}\frac{[\stackrel{~}{\mathrm{\Sigma }}^{\mathrm{th}}(x_i^{},y_j^{})\stackrel{~}{\mathrm{\Sigma }}^{obs}(x_i^{},y_j^{})]^2}{\sigma _{ij}^2},$$
(45)
where $`\sigma _{ij}`$ denotes the observational error at each pixel point $`(x_i^{},y_j^{})`$. Those systematic/random errors should depend on specific observation conditions and thus are not easy to estimate a priori. Thus in the following numerical tests, we simply assume that $`\sigma _{ij}`$ is unity, independent of $`i`$ and $`j`$.
* Determine the best-fit values of $`e_b,e_c,\theta ,\varphi ,\kappa `$ (or $`\kappa _p`$) through the $`\chi ^2`$ minimization.

The numerical testing results have revealed that the above algorithm works quite well in reconstructing the halo eccentricities within the fractional error of $`20\%`$ as long as the halo eccentricities are not so small ($`e_b,e_c>0.3`$), which can be understood considering that our reconstruction algorithm strongly relies on the non-spherical signature, and thus is expected to fail for the case of almost spherical halos with low eccentricities ($`e_b,e_c<0.3`$). For the low eccentricity cases, the perturbation results of the Paper I may be more useful. Since the low-eccentricity case is not the main interest of the present paper, we do not investigate those cases here. The algorithm also works in reconstructing the orientation angles but suffer relatively large fractional errors. We show the examples of our numerical reconstruction results in Figures 6 and 7 for the following cases of the halo eccentricities and orientation angles: $`(e_b,e_c)=(0.4,0.7)`$, $`(0.5,0.6)`$, and $`(0.6,0.8)`$; $`(\theta ,\varphi )=(15,75)`$,$`(45,45)`$, $`(75,15)`$ in units of degree.

Figure 6 plots the fractional error of the reconstructed eccentricities versus their input values for three different cases of the orientation angles. The left three panels are for $`e_b`$, and the right three panels for $`e_c`$. The top two panels correspond to the case of $`\theta =15^o`$ and $`\varphi =75^o`$, the middle two panels to the case of $`\theta =45^o`$ and $`\varphi =45^o`$, and the bottom three panels to the case of $`\theta =75^o`$ and $`\varphi =15^o`$. For almost all cases, the fractional errors fall within $`20\%`$.

Figure 7 plots the fractional error of the reconstructed orientation angles versus the original orientation angles for three different cases of the eccentricities. The left three panels are for $`\theta `$, and the right three panels for $`\varphi `$. The top two panels correspond to the case of $`e_b=0.4`$ and $`e_c=0.7`$, the middle two panels to the case of $`e_b=0.5`$ and $`e_c=0.6`$, and the bottom three panels to the case of $`e_b=0.6`$ and $`e_c=0.8`$. For most cases, the fractional errors fall within $`20\%`$. However, for a few cases, the fractional errors are larger than $`50\%`$, indicating that the angle reconstruction is not always stable compared with the eccentricity reconstruction.

## 5 DISCUSSIONS AND CONCLUSIONS

It was regarded as a difficult task to reconstruct the 3D structures of dark matter halos from the projected 2D surface-brightness maps of X-ray and/or SZ clusters given the parameter-degeneracy caused by the projection process itself. To break the degeneracy, previous approaches had to rely on such restrictive assumptions as the cluster axial symmetry about the line-of-sight and some priors like weak lensing informations. Here we have developed, for the first time, a theoretical framework within which the 3D halo reconstruction is possible in principle without such restrictive assumptions and additional priors.

We derived the density profiles of the intra-cluster gas from the 1st principles, assuming that the intra-cluster gas is either isothermal or polytropic, and in hydrostatic equilibrium within the gravitational potential of concentric and coaxial triaxial dark matter halos. In our theoretical modeling of the intra-cluster gas, the density profiles depend explicitly on the intrinsic properties of the underlying halos (the two halo eccentricities, $`e_b`$ and $`e_c`$). For the case of highly aspherical halos, however, we have found that the gas density profiles cannot be expressed in closed analytic forms. Therefore, we have attempted to find general fitting formulae for the gas density profiles as functions of the halo eccentricities which may be applicable even to highly aspherical clusters. We have found empirically a simple scaling relation with respect to $`\mu =e_b^3+e_c^3`$, and provided a set of fitting formulae for the 3D gas density profiles expressed in terms of $`\mu `$.

Then, we have incorporated the projection effect into the gas density profiles, and provided another set of fitting formulae for 2D X-ray and/or SZ surface brightness density profiles of galaxy clusters. Unlike the 3D formulae, the 2D formulae are expressed in terms of not only $`\mu `$ but also the orientation angles of the line-of-sight in the halo principal axes. In other words, we have found a link of the 2D X-ray and/or SZ observables to the shapes and the orientations of the dark matter halos.

We have proposed a numerical algorithm based on our model to determine the halo eccentricities and orientation angles from the observed X-ray and SZ surface brightness density profiles of galaxy clusters. Ideally we have to apply our algorithm to real observation data, but that is beyond the scope of the present paper since it should involve careful treatment of the data image analysis together with various observational systematic effects. Thus we have decided to test the reconstruction algorithm against simple numerical profiles. We found that the algorithm can reconstruct the halo eccentricities fairly accurately if the halo eccentricities are greater than $`0.3`$. On the other hand, it turns out that the reconstruction of the halo orientation angles are not always accurate, showing large scatters.

Finally, we conclude that our hydrostatic-equilibrium model for intra-cluster gas in triaxial dark halos makes it possible in principle to reconstruct the 3D structures of the dark halos from the X-ray and/or SZ cluster maps, as long as the intra-cluster gas and the dark halos are well approximated to be in hydrostatic equilibrium and the triaxial ellipsoids with a constant axis ratio, respectively. We plan to test our model against real observational data, and hope to report the results elsewhere in the future.

We thank the anonymous referee for his/her careful reading of the manuscript and several constructive criticisms. We are grateful to Masamune Oguri for useful discussions. J. L. acknowledges gratefully the research grant of the JSPS (Japan Society of Promotion of Science) fellowship. This research was supported in part by the Grant-in-Aid for Scientific Research of JSPS (12640231).