text:
# XMM-Newton Observation of the Merging Galaxy Cluster A1644

## 1 Introduction

Clusters of galaxies are formed in violent events, cluster mergers. These mergers can dramatically affect the observed properties of clusters, particularly in the X-ray band (e.g., Markevitch et al., 2002). Signatures of mergers include “cold fronts” (Markevitch et al. 2000; Vikhlinin, Markevitch, & Murray 2001b), merger shocks (e.g., Markevitch et al., 2002), and multiple cool cores (Markevitch, Sarazin, & Vikhlinin, 1999). The application of simple hydrodynamical analysis to these merger signatures allows the merger Mach number and other aspects of the merger kinematics to be determined (Vikhlinin et al. 2001b; Kempner, Sarazin, & Ricker 2002; Sarazin 2002). In general, merger Mach numbers are modest, $`2`$.

Mergers also provide interesting tests of the role of various physics processes in clusters. For example, very strong limits have been placed on the rate of thermal conduction and diffusion (Ettori & Fabian 2000; Vikhlinin, Markevitch, & Murray 2001a; Markevitch et al. 2003b) and Kelvin-Helmholtz instabilities and viscosity (Vikhlinin et al., 2001a). They also provide limits on the relative roles of thermal versus nonthermal pressure sources in the intracluster medium (Markevitch et al., 1999). The role of merger shocks in accelerating relativistic particles can also be studied (Markevitch & Vikhlinin, 2001). Recently, the relative distributions of gas, galaxies, and dark matter in merging clusters have been used to limit the self-collisional cross-section of the dark matter (Markevitch et al., 2003a); collisional dark matter has been suggested to explain the mass profiles in dwarf galaxies and to resolve other problems (e.g., Spergel & Steinhardt, 2000).

On the other hand, merger effects can complicate the use of clusters of galaxies as cosmological probes. Clusters provide a number of important cosmological tests, mainly because they are the only objects which are large enough to represent a fair sample of materials in the Universe, yet small enough to be relaxed (e.g., Bahcall, 2000). Most analyses of clusters for cosmology assume that they are equilibrium systems; for example, the masses of clusters are often determined from hydrostatic equilibrium, or more simply by assuming an equilibrium relation between the cluster X-ray temperature or luminosity and the cluster mass. During major mergers, clusters are not in hydrostatic equilibrium. Also, compression and heating during a merger can produce a large transient increase in the X-ray temperature and luminosity (Ricker & Sarazin, 2001). If equilibrium relations are applied to merging clusters, their masses can be greatly over-estimated (Randall, Sarazin, & Ricker, 2002). This can be particularly important for high redshift clusters; in flux limited samples, merging clusters can be over represented due to their boosted X-ray luminosities. It is important to understand these effects by studying relatively nearby clusters where more detailed dynamical information can be determined.

In this paper, we present new XMM-Newton observations of the Abell 1644 cluster. Our analysis of the Einstein observation suggested that this cluster was undergoing a major merger. Jones & Forman (1999) classified A1644 as a double cluster with roughly equal components based on the Einstein image. From a deprojection analysis of the Einstein data, White, Jones, & Forman (1997) found that A1644 harbors a small cooling flow of about 12 solar masses per year. Girardi et al. (1997) noted the possible presence of two components in the galaxy velocity distribution in the cluster, and listed two subclusters as part of this cluster. They also found that the cD galaxy had a large peculiar velocity, which suggests the possible presence of minor substructure. However, in the end they classified Abell 1644 as a regular cluster. No evidence for substructure has been found in recent optical and near-infrared observations by Tustin et al. (2001). A1644 is located close to a region of high galaxy cluster density. It is about 15 degrees north (in Galactic coordinates) of the core of the Shapley supercluster (Fig. 1) at about the same redshift, $`z(\mathrm{A1644})=0.0474`$ (Zabludoff et al., 1993).

## 2 Data Reduction and Analysis

XMM-Newton, with its combination of high throughput, broad spectral band, good spectral resolution, and fairly large field of view is ideal for studying this merging system. A1644 was observed by XMM-Newton for a total duration of 23 ks during revolution 199 (observation ID: 0010420201) on January 8, 2001. The EPIC-MOS cameras were operated in Full Frame (FF) mode and the EPIC-pn camera in Extended Full Frame (EFF) mode. The Thin1 filter was used in all EPIC observations for protection against optical light contamination. Calibrated event files have been created with SAS V5.3. The calibration index file has been created on May 7, 2002 based on the most up to date calibration files available at that time. Most further reduction and analysis steps have been performed using SAS V5.3 and the package HEAsoft V5.2.

### 2.1 Flare Rejection

XMM-Newton observations in general suffer from periods of very high background similar to Chandra observations (e.g., Lumb et al., 2002). There are several possible ways to ‘clean’ a data set of these flares. Our method is based on the fact that we will use independent blank-sky observations (§ 2.3) to determine the main non-flare hard background component, the ‘Particle Induced Background’ (PIB). (The PIB is distinct from the flares, which are believed to be produced by ‘soft’ protons.) Thus, we must be able to clean the blank-sky background observations in the same way as the source observation. When determining the lightcurve this requires ensuring a negligible ‘contamination’ by emission from astrophysical sources because otherwise the sensitivity to flare detection varies between source and background observation. Here we choose for the MOS cameras the energy band 10–12 keV and for pn the band 12–14 keV similar to previous works (e.g., Majerowicz, Neumann, & Reiprich, 2002), which are dominated by PIB (background induced directly or through fluorescent instrumental lines by high energy particles) and soft protons (if present). Note that this approach has the drawback of low statistics and furthermore assumes that the flare characteristics in this energy band are representative of those at lower energies as well. In § 2.3 we test the effect of increased statistics and the latter assumption.

Times were binned in 100 s intervals, which appears to be a reasonable compromise between good temporal resolution and low noise level. This time interval also is large compared to the read out cycles of the cameras. Further only events with pattern $``$ 12 and flag = 0 have been used to ensure a high rejection of potentially distorting events such as caused by hot/flickering pixels/columns or features close to CCD edges especially at the read out nodes.

Inspection of the resulting lightcurves showed no obvious flares in the pn data (Fig. 2). MOS count rates were higher by up to a factor $``$ 4 in the beginning of the MOS observations compared to the remaining times (Fig. 3). MOS exposures started earlier than those for the pn camera; the lack of obvious flares in pn is consistent with the MOS lightcurves for the overlapping time interval. The obviously contaminated times in the MOS observations were removed by taking into account only events with arrival times $``$ 95360200 s (XMM-Newton time) for all further steps.

Before calculating mean count rates in order to test for background flares, conservative ‘generic’ cuts have been applied on the count rates, $`C_\mathrm{X}`$. For pn (MOS) only times with $`2(1)C_\mathrm{X}45(25)`$ counts per 100 s have been used. The lower threshold is especially important for pn data in order not to bias the calculation of the mean count rate by inclusion of counting mode times, times during which transmission of information for individual events of one or more (most likely all four) quadrants to the ground station does not occur.

For the final cuts, upper *and* lower thresholds are required in order to treat statistical fluctuations in both directions symmetrically. Here we choose to define the good time intervals (GTIs) by including only times where the count rates lie within $`\pm `$ 2 $`\sigma `$ of the mean count rates. For the individual detectors the count rate limits for GTIs are found to be 8–25 (pn), 3–15 (MOS1), and 2–15 (MOS2) counts per 100 s. The effective exposure time (live time) after this procedure results in 11 ks for pn and 13 ks for MOS1 and MOS2.

### 2.2 Vignetting Correction

Depending on the nature of a detected event one should take into account vignetting caused by the baffles, mirrors, Reflection Grating Array (for the MOS cameras only), etc. in the analysis. For instance the X-ray photons from the target source (A1644) and the Cosmic X-ray Background (CXB) are vignetted in the same way. The soft protons show a different vignetting characteristic (but hopefully we got rid of most of those in § 2.1) and the PIB is not vignetted at all if the spatial variation of the fluorescent lines caused by the camera material itself is not interpreted as a vignetting characteristic. The noise in the CCDs, most easily visible as low energy events close to the CCD edges, is not vignetted, either.

In general, however, it’s not easily possible to determine the nature of an individual event. The combination of methods for flare rejection, vignetting correction and background subtraction employed here aims at a proper treatment of all components for the spectral analysis. This strategy involves correcting *all* events for the vignetting expected for X-ray photons regardless of their nature at this point.

The vignetting depends on detector position and energy. The vignetting correction has been accomplished with the SAS task evigweight which adds a weight column to the events file. The value of the weight is determined for each event by taking into account the detector position and energy of the event. Weights are assigned relative to the on-axis effective area so the weights have values $``$ 1. To accumulate images and spectra these weights are summed up, and errors are propagated including the vignetting correction.

### 2.3 Background Subtraction

In general, proper background subtraction is vital for an accurate spectral analysis of extended sources. In our case, the cluster emission fills most of the fields-of-view of the detectors, so it is difficult to determine the background from our observations. Thus, we determine the background from so called blank-sky observations. The blank-sky observations are mostly high galactic latitude observations of regions with very few bright sources. Several such observations have been combined for the EPIC cameras mostly for the Thin filter in (mostly) FF mode (Lumb et al., 2002) to achieve good statistics and reduce the effects of cosmic/sample variance. We have made sure that for the energy range and signal-to-noise ratio in spectral fits here, the differences in Charge Transfer Inefficiency (CTI) corrections between the pn FF and EFF mode are not significant. The blank-sky observations are available at the XMM-Newton calibration web page111http://xmm.vilspa.esa.es/external/xmm\_sw\_cal/calib/epic\_files.shtml.

The background observations have been recast on the sky in order to match the source observations using the SAS task attcalc222See Andy Read’s scripts at http://www.sr.bham.ac.uk/xmm3/ .. This allows using sky coordinates for regional selections and at the same time estimating the background from the same detector regions as the source, which in turn ensures that the significant variation of the fluorescent line emission with detector position is properly taken into account. By treating the blank-sky backgrounds identically to the data, the effect of incorrectly applying a vignetting correction to the PIB (in the source *and* background observations) is removed.

The same flare rejection method (§ 2.1) has been applied to the background observations. The count rate limits for GTIs here are found to be 8–31 (pn) and 2–16 (MOS1 and MOS2) counts per 100 s.

The PIB shows variation of about 8% (1 $`\sigma `$) as well as a long term decrease of about 20% between the beginning of XMM-Newton observations and January 2001 (Katayama et al., 2002). In order to account for these variations, total field-of-view count rates of source and background observations have been calculated using the same energy, pattern, and flag selection criteria as for the GTI determination (§ 2.1). The source to background count rate ratios $`0.99\pm 0.03`$ (MOS1), $`0.92\pm 0.03`$ (MOS2), and $`0.88\pm 0.02`$ (pn) have been found (1 $`\sigma `$ errors), consistent with the expected variation. For the calculation of these PIB scaling factors the events have not been corrected for vignetting.

As mentioned in § 2.1, cleaning of flares in the high energy bands generally used does not necessarily imply that all flares have been removed. Actually, we did find some weak remnant flares in the blank-sky observations when plotting lightcurves (GTI screened as described in § 2.1) in broad energy bands (see also Lumb et al., 2002). We therefore tested the variation in (10–12 keV) source to background count rate ratios when calculating GTIs in the energy bands 7–12 and 0.3–12 keV for MOS1 and MOS2. We found that the variation is less than the (1 $`\sigma `$) statistical uncertainty and conclude that the influence of remnant flares is negligible for calculation of background scaling factors here.

To check whether the limitation to high energy bands for the calculation of background scaling factors affects the results significantly we have calculated them also using the energy ranges 0.3–14 keV (pn) and 0.3–12 keV (MOS). In this case, of course, only the CCD parts shielded from X-ray photons (and soft protons) are usable. We found results consistent with the total field-of-view results for pn and MOS1 and a slightly larger value for MOS2 ($`1.00\pm 0.02`$). In order to exclude the possibility that an imperfect correlation between the normalization of the continuum and the fluorescence lines (dependent on detector position) for the PIB caused this slightly larger than expected variation, we also used the energy range 2–7 keV to calculate the background scaling factors. The results were consistent with the broader energy range used previously.

In summary, following these tests the background scaling factors 0.96 (MOS) and 0.88 (pn) have been assigned to the observations. Most of the tests above yield results consistent with the expected statistical uncertainties. To account for the difference between MOS and pn scaling factors a conservative 8% systematic error has been added to the full energy range of all background spectra. Note also that it is implicitly assumed here (and in most other works that follow similar procedures) that the spectrum of the PIB does not vary when its normalization changes. Currently, there seems to be no strong evidence for significant spectral variability (e.g., Katayama et al., 2002) apart from the fluorescence lines.

The main drawback of this method is that in general one gets the background contribution from the CXB wrong. The PIB dominates the hard part of the total background but the CXB is important for the lower energies. Ideally one would estimate the difference in CXB contribution to source and background observations from a region with negligible source emission in the source observation or another observation very close by (e.g., Markevitch & Vikhlinin 2001; Pratt, Arnaud, & Aghanim 2001; Majerowicz et al. 2002; Arnaud et al. 2002). This is not feasible in the case of A1644 since cluster emission fills the entire field of view and the only other XMM-Newton observation within 5 degrees to date is likewise filled entirely with NGC 5044’s emission. Here, we model a possible CXB difference in the spectral fitting procedure.

In order to estimate the spectral shape of a possible CXB difference between the source and background spectra we have extracted spectra from an irregular region with decreased cluster emission in the north-western part of the observation and fitted a two temperature model. A high temperature component which is supposed to account for the cluster emission and a low temperature, zero redshift, solar abundance component which should represent the CXB excess relative to the background observations. The model fits indicate that in this region the cluster emission still dominates. The ratio of the number of cluster photons to CXB photons equals about 4 in the energy range 0.3–10 keV and about 2 in the energy range 0.3–0.5 keV. A best fit value $`0.2`$ keV is found for the low temperature component. This estimate is consistent with values from studies of soft excess emission in other clusters (e.g., Kaastra et al., 2003). The origin of this component is uncertain; Kaastra et al. (2003) argue that it is due to intergalactic gas located around clusters. Our aim here is not to try to discuss the source of this soft CXB excess, but rather to determine its contribution to the cluster emission and account for it in the model fits (§ 2.4).

### 2.4 Spectral Fitting

Spectra from MOS1, MOS2, pn have been fit simultaneously. The same pattern and flag selections have been applied as for the generation of GTIs, except that only single and double events have been used for the pn camera. The response matrices were m11\_r7\_im\_all\_2000-11-09.rmf, m21\_r7\_im\_all\_2000-11-09.rmf, and epn\_ef2
0\_sdY9.rmf for MOS1, MOS2, and pn respectively, in addition to on-axis ancillary response files created with the SAS task arfgen.

A comparison of individual pn single and double event spectral fits gives consistent results. The response of the pn camera depends on detector position. However, it has been verified that for the accuracy needed here, use of the epn\_ef20\_sdY0.rmf response matrix does not result in significant changes of best fit parameter values. Individual fit results from the three cameras are consistent with each other. Point sources and hot pixels have been excluded.

XMM-Newton’s point spread function (PSF) is much improved compared to the PSFs of previous satellites capable of spatially resolved spectroscopy, e.g., ASCA. Nevertheless, for objects with strong temperature and surface brightness gradients its effect can become important. However, the main emphasis in this work is on the comparison of different segments of spectra at the same radii from cluster centers. In this case radial PSF effects should not cause any artificial differences.

Spectral fits have been performed within XSPEC using a model (wabs\*mekal)+mekal. Parameters for the second thermal component have been fixed at a redshift of zero, a metal abundance of $`A=1`$ solar (Anders & Grevesse, 1989), and temperature of $`kT_\mathrm{X}=0.2`$ keV, as discussed in § 2.3. The normalization is determined from the fits.

Source spectra have been regrouped to have at least 50 counts per bin. Spectral fits have been performed nominally in the energy range 0.3–10.0 keV for MOS and 0.3–12.0 keV for pn. All errors are quoted at the 90% confidence level for one interesting parameter unless noted otherwise.

## 3 Results

### 3.1 Imaging

Figure 4 shows the combined MOS1-MOS2-pn image of nearly the full XMM-Newton field-of-view of A1644. The cluster has a very complicated surface brightness distribution on all scales. A main cluster to the southwest and a smaller sub cluster to the northeast are easily identified. The emission surrounding both core regions is highly non-spherical. There is stronger emission to the south of the sub clump than to the north. The core region of the main clump itself contains a displaced core-within-a-core.

The Digitized Sky Survey (DSS) image in Fig. 5 shows that the X-ray peaks of main cluster and sub cluster coincide with two Bright Cluster Galaxies (BCGs). The relative line-of-sight velocity component of the two central galaxies is about 600 km s-1. The X-ray contours to the north of the sub clump appear to be slightly compressed. Radio observations of the two BCGs show that they are both radio sources (e.g., Owen & Ledlow, 1997).

The regions around the two BCGs are shown in more detail in Figs. 6 and 7. Included are X-ray contours of a 19 ks Chandra observation taken on the same day as the XMM-Newton observation and available from the public archive. These contours illustrate the complexity of the emission in the very cluster cores. Note that also in the very center of the sub clump the X-ray contours to the north appear slightly compressed (Fig. 7). This apparent compression is on a much smaller scale than the compression already noted in Fig. 5. Optical, radio, and X-ray centers coincide almost exactly for both sub clusters.

### 3.2 Spectroscopy

The intracluster gas temperature is a fundamental observable for a gravitational cluster mass determination based on the hydrostatic assumption. The overall best fit temperature based on a large elliptical region centered at R.A. $`=194.356\mathrm{deg}`$ and Dec. $`=17.363\mathrm{deg}`$ encompassing both sub clusters (semi major axis $`=733\mathrm{}`$, semi minor axis $`=528\mathrm{}`$, position angle $`=45\mathrm{deg}`$) is $`kT_\mathrm{X}=3.83\pm 0.06`$ keV. The best fit metal abundance is $`0.32_{0.03}^{+0.02}`$ solar for the same region.

Note that the column density of neutral Galactic hydrogen, $`N_\mathrm{H}`$, has been fixed at the value inferred from 21 cm radio measurements for our Galaxy ($`5.33\times 10^{20}`$ cm-2; Dickey & Lockman, 1990) for all results. In general this value is marginally larger but consistent with best fit values from direct spectral fits. Best fit values for temperatures and abundances change only well within their uncertainties if the column density is treated as a free parameter.333Note that neglecting the CXB difference between source and background observations, i.e., fitting a model wabs\*mekal to the spectra here, results in $`N_\mathrm{H}`$ values significantly lower than the Galactic value, or, if $`N_\mathrm{H}`$ is fixed to the Galactic value, to artificially decreased temperatures and increased $`\chi ^2`$.

Figures 8 and 9 show the radial temperature profiles for the main and sub clump, respectively. Note that the center for the innermost bin of the main clump is offset from the center of the other bins in order to account for the core-within-a-core structure. The reduced $`\chi ^2`$ values (and the best fit $`N_\mathrm{H}`$ values if left as a free parameter) for the innermost bin of both sub clusters are too large (Tab. 1 and 2). This is likely caused by higher temperature gas in the line of sight. For these two regions we therefore included a third thermal component, i.e., the model (wabs\*(mekal+mekal))+mekal was fitted. The limited statistics require freezing the temperatures and abundances of the third components to the values found further out. We therefore fixed the abundances at 0.5 solar and the temperatures at 3.5 and 3 keV for the main and sub clump, respectively. As a result the best fit central temperatures decrease further for both main and sub clump (open diamonds).

Each of the two temperature profiles appears surprisingly similar to temperature profiles of relaxed, apparently undisturbed clusters: a drop in the center to about 1/3 of the ambient gas temperature, an isothermal structure in the outer parts, and weak indications for a slight temperature drop in the very outermost regions accessible.

Note that the addition of a third component is only a crude approximation to account for effects caused by projection, PSF (see, e.g., Markevitch, 2002; Majerowicz et al., 2002), and a possible real multiphase temperature structure. The reduced $`\chi ^2`$ values of the central bins become significantly smaller but are still large even after addition of the third component. Also, these effects need not be limited to the central bins alone. A1644 is an extremely irregular object (see Figs. 4 and 12). This fact makes A1644 so interesting but at the same time complicates deprojection attempts, for which spherical symmetry usually is assumed. In this work, we therefore focus on projected temperature estimates and qualitative comparisons to temperature maps of simulated clusters.

The corresponding metal abundance profiles are shown in Figs. 10 and 11. In general, the abundances increase towards the centers of the two clumps, although there seems to be a drop in abundance in the innermost region of the sub clump. However, this abundance drop vanishes when performing the fit with three thermal components (open diamond). This drop can therefore be regarded as artificial, even though it should be noted that the precise values depend on the temperature and abundance values adopted for the third component.

It is worth noting that weak indications for higher normalizations of the CXB excess emission towards the center of the main clump are found which may support the findings of Kaastra et al. (2003), although a discussion of the nature of the CXB excess is beyond the scope of this paper (see § 2.3).

Figure 12 shows a hardness ratio (HR) map. To create this map, first, source photon images from MOS1, MOS2, and pn (using only single events) in the energy band 2.0–10.0 keV have been created and combined. Then, the corresponding exposure maps have been created and combined, which also account for vignetting. Using the combined image and exposure map, an adaptively smoothed and exposure and vignetting corrected source image has been created with the SAS task asmooth (signal-to-noise ratio $`=15`$). Subsequently, a similarly created (though simply smoothed with $`\sigma =7.5`$), exposure and vignetting corrected, combined background image has been subtracted from the source image. The same procedure has been followed to create an image in the 0.3–2.0 keV energy band, using the same smoothing templates as created for the 2.0–10.0 keV image. The ratio of the latter two images is shown in Fig. 12. Soft emission is shown in dark and hard emission in bright. We have verified that using the source smoothing template also for smoothing the background image does not result in significant changes in the relevant region. Note that the CXB excess emission has not been subtracted in the HR map. This may bias the HR map in the outermost regions of the observation. However, we will find below that spectral fits that take the CXB excess into account correspond well to the features seen in the HR map. Therefore, it is not expected that the region shown is biased significantly.

In Fig. 13 an adaptively binned temperature map is shown. The map has been created using all three detectors and making use of the Interactive Spectral Interpretation System (ISIS; Houck & Denicola, 2000). The pixel size is 14.4″. Pixels are nearly independent at the center of the main clump but are correlated in other regions. The actual spectral extraction region used for the fits depends on the surface brightness; the largest extraction region is 1′on a side. Most pixels shown black had fewer than 800 source counts using the maximum smoothing scale. A few more pixels with negative error $`>`$ 30% were masked out. The abundance of the main mekal component has been fixed here at the average value 0.32 solar (§ 3.2). The features seen in the HR map (to be discussed in detail below) seem to correspond to features in the temperature map. This indicates the consistency of both approaches and we’ll use the HR map, which is more closely related to the data, as guidance below.

The hardness ratio map and the X-ray image suggest a number of interesting regions in the cluster which may have differing temperatures and/or abundances. In general, clusters as irregular as A1644 often show more complex temperature structure. Therefore, we also extracted spectra in segments of annuli for the main clump. The radial boundaries for the nine segments of the main cluster correspond to the ranges for the three outermost annuli given in Tab. 1. Figure 12 shows these and other selected regions as well as surface brightness contours overlaid onto the HR map. (Note that the cluster emission extends beyond the outermost shown surface brightness contour.)

Best fit temperature and metal abundance values for the segments are shown in Figs. 14 and 15 (and Tab. 3), respectively. It is apparent that there is significant nonradial temperature structure, and that the relatively regular radial temperature profile in Figure 8 is the result of averaging over annuli with large azimuthal temperature variations. In particular, the temperatures of the segments to the west of the main clump (R3, R6, and R9) are all significantly lower than all regions to the east (R1, R2, R4, R5, and R7) except for R8. These findings correspond well to those in the HR map (Fig. 12).

We also determined temperatures and abundances for more complicated regions selected as being of interest based on the XMM-Newton image and the HR map. This allows us to search for cool/metal-rich trails (dark in the HR map) or hot spots (bright in the HR map). Whether brightness fluctuations in the HR map correspond to significant temperature/abundance changes is evaluated by direct spectral fits of selected regions. Note that not all artifacts, e.g., inexact exposure correction close to CCD chip boundaries, have been removed in the HR map in Fig. 12. The spectral analysis, however, is not affected by this.

The region to the south of the sub clump (R3) appears fairly soft, whereas a small region to the southeast (R4) appears hard. The spectral fit results (Fig. 16 and Tab. 4) reveal that indeed the temperature determination for R3 is significantly lower than the value for the rest of this annulus (R2, which excludes R1) and especially than that for R4. Furthermore, the metal abundance of R3 appears enhanced (Fig. 17). The best fit temperature for R5 (the small dark tail to the south of R3 in Fig. 12) is similar to that for R3. On the other hand, R1, the region north-northeast of the sub clump, has a high temperature similar to R4.

## 4 Discussion

We first consider some implications of these observations for the use of clusters of galaxies as cosmological probes. How much do the X-ray flux estimates based on resolved observations differ from simple estimates based on global measurements, such as those for clusters only observed with the ROSAT All-Sky Survey (RASS)? How much do mass estimates differ when only a broad beam overall gas temperature estimate is available?

The overall flux from A1644, which extends beyond XMM-Newton’s field-of-view, is $`f_\mathrm{X}(0.12.4\mathrm{keV})=4.03\times 10^{11}\mathrm{ergs}\mathrm{s}^1\mathrm{cm}^2`$ (Reiprich & Böhringer, 2002) based on RASS data. (A1644 was not observed in pointed mode with the PSPC.) It is worth noting that A1644 has been treated as a single cluster in the HIFLUGCS (Reiprich & Böhringer, 2002) and REFLEX (Böhringer et al., 2001) catalogs, primarily due to the limited exposure time ($`<`$ 300 s) of the RASS at A1644’s position and the relatively small projected separation of the two cluster components. The XMM-Newton image indicates that the X-ray flux ($`f_\mathrm{X}`$) ratio between the main and sub clump is about 3:1. This means instead of one cluster with $`f_\mathrm{X}4\times 10^{11}\mathrm{ergs}\mathrm{s}^1\mathrm{cm}^2`$ one actually has two clusters with about $`f_\mathrm{X}3`$ and $`1\times 10^{11}\mathrm{ergs}\mathrm{s}^1\mathrm{cm}^2`$ each. This could affect luminosity and mass functions significantly if such major mergers were common. The effect of an incorrect treatment would be to artificially increase the number of high mass clusters and decrease the number of low mass clusters. Using the shape of the mass function to constrain the cosmic mass density, $`\mathrm{\Omega }_\mathrm{m}`$, and the amplitude of mass fluctuations, expressed as $`\sigma _8`$, individually would then result in artificially low values for $`\mathrm{\Omega }_\mathrm{m}`$ and high values for $`\sigma _8`$. With the ongoing analysis of the complete sample of HIFLUGCS clusters with Chandra and XMM-Newton, we will be able to quantify the fraction of such clusters and their influence on derived cosmological parameters, as well as to study the difficulty of estimating survey selection functions in the vicinity of bright clusters (not only due to double clusters but also line of sight projections) in general.

A simple broad beam temperature estimate may be biased low compared to the ambient gas temperature due to cool emission in the dense (high emissivity) cores. For instance, the temperature in a large annulus around the main clump in A1644 (201–325$`\mathrm{}`$) is 4.40 keV. If we take this as the ambient temperature which is characteristic of the cluster gravitational potential, it is a factor of 1.15 higher than a broad beam temperature estimate including both clumps (3.83 keV). Since $`MT^{1.52.0}`$ (e.g., Finoguenov, Reiprich, & Böhringer, 2001), this factor translates into an underestimate of a factor of 1.23–1.32 in the cluster mass.

The surface brightness structure of both sub clumps is obviously messy. However, the radial temperature profiles (Figs. 8 and 9) are smooth and appear rather similar to profiles of relaxed clusters (e.g., compare to Fig. 1 in Allen, Schmidt, & Fabian, 2001). This might suggest that A1644 is moderately relaxed, and that the temperature structure may not be strongly affected by the interaction of the two sub clumps. However, a more detailed examination of the thermal structure shows that the cluster is not at all relaxed. As the HR map (Fig. 12) indicates and Figs. 14 and 16 confirm, the temperature structure is also quite complex. The region between the two clumps is significantly hotter than the region to the west of the main clump (the region which appears least disturbed in X-ray surface brightness). This may indicate that the gas in the region between the clumps has been heated up by adiabatic compression or shocks. Also, the core-within-a-core structure of the main clump (Fig. 4) may be caused by core oscillations induced by the passage of the sub clump (e.g., Tittley & Henriksen, 2003) or possibly by a remnant of a previous merger.

Furthermore, both the surface brightness and temperature structure around the sub clump are peculiar. Based on the X-ray image, the gas appears to be slightly compressed to the north and clearly elongated to the south and west. This suggests that this sub clump is moving through the intracluster medium (ICM) of the main clump, and losing some fraction of its gas in a trail behind the core of the sub clump. Given this dynamical picture what temperature structure might be expected? One may expect heated gas in front of the sub clump (projected roughly to the north) due to adiabatic compression or shocks. One might also expect the trail to the south to consist of somewhat cooler, more metal-rich gas stripped from the sub clump. This appears to be consistent both with the HR map and with the extracted spectra in the sub clump regions R1, R3, and R5.

Figs. 5 and 7 show that the X-ray peak of the sub clump and the optical center of the brightest galaxy coincide. This suggests that the sub clump has not passed through the core of the main clump; if it had, the very high ram pressure would have stripped the sub clump core, which would lag behind the collisionless dark matter and galaxies (e.g., Markevitch et al., 2002). The angle of the cool, metal-rich gas trail also suggests an off-center collision. One concern with this picture is that the dense gas in the core of the sub clump would then be interacting with lower density gas in the outer parts of the main clump. Would one expect the dense sub cluster to lose an observable amount of gas during its travel through the ICM?

We have compared A1644 to a hydrodynamical simulation of a small cluster undergoing an off-axis merger with a larger cluster. This simulation was done completely independently; no effort was made to adjust the initial conditions to match the observation of A1644. The temperature map of this simulation at a time when the smaller cluster has passed by the core of the main cluster for the first time is shown in Fig. 18. Comparison to the HR map of A1644 (Fig. 12) shows an interesting similarity. There is hotter gas in front of the smaller cluster and cooler, ram-pressure stripped gas behind. This trail of stripped gas cools adiabatically as it expands. There is a hotter region to the southwest, which is similar to region R4 seen in A1644. The larger scale $`T_\mathrm{X}`$ distribution around the main cluster in the simulation is also similar to that observed in A1644. There appears to be slightly cooler gas to the southeast in the simulation due to stripped gas from the smaller cluster, which agrees with the R8 region in the HR map and spectra.444 An animation of this simulation is available at http://A1644.dark-energy.net. Note that cool trails from sub clusters are commonly seen in simulations. Compare, e.g., to the simulations by Motl et al. (2003, see also the movie at http://www.astro.virginia.edu/coolflow/abs.php?regID=184).

The indication for an increased metal abundance in the cool trail of A1644’s sub clump makes other scenarios like cooling of intracluster gas onto the moving sub clump (as seen in A1795 on much smaller scales, Fabian et al., 2001) more difficult to reconcile with the data. Note that the radio source in the center of the sub clump (Fig. 7) may provide additional energy for gas removal from central denser regions.

## 5 Summary

XMM-Newton imaging spectroscopy of the galaxy cluster A1644 has shown two sub clusters undergoing an off-axis merger with remarkable detail. The X-ray analysis is also supported by optical data on the cluster, and the resulting dynamical picture agrees very well with an independent hydrodynamical simulation of such an off-axis merger. The findings imply that we see gas that has been removed from the smaller sub cluster possibly by ram pressure from the motion through intracluster gas of the main clump, possibly augmented by the effect of the central radio source of this sub clump.

The code to create the X-ray temperature map made use of the Interactive Spectral Interpretation System (ISIS; Houck & Denicola 2000). M. Ledlow made available the radio images. This work was supported by NASA XMM-Newton Grants NAG5-10075, NAG5-13088, and NAG5-13737. T. H. R. acknowledges support by the Celerity Foundation through a Post-Doctoral Fellowship. The XMM-Newton project is an ESA Science Mission with instruments and contributions directly funded by ESA Member States and the USA (NASA). The Digitized Sky Surveys were produced at the Space Telescope Science Institute under U.S. Government grant NAG W-2166. The images of these surveys are based on photographic data obtained using the Oschin Schmidt Telescope on Palomar Mountain and the UK Schmidt Telescope. This publication makes use of data products from 2MASS, which is a joint project of the University of Massachusetts and IPAC/Caltech, funded by NASA and NSF.


format: text

text:
# Stochastic description of the dynamics of a random-exchange Heisenberg chain

## I Introduction

In the last decades, a considerable number of dynamical systems has been studied and a great deal of attention has been paid to the analysis of its transport properties. In particular, the study of diffusion and transport properties of physical systems with short or long-range correlations in the disorder distribution has attracted a renewed interest Morgado ; flores1 ; dunlap ; sen ; datta ; adame1 ; adame2 ; flores2 ; ojeda ; chico ; izrailev ; evangelou1 . For instance, the unexpected high conductance of several doped quasi-one-dimensional polymers was explained by assuming pairwise correlations in the disorder distribution dunlap . Similarly, the suppression of Anderson localization was recently confirmed experimentally in semiconductor superlattices with intentional correlated disorder bellani . Further, it was demonstrated that long-range correlations in site also act towards the delocalization of 1D quasiparticle states chico ; izrailev . The 1D Anderson model with long-range correlated diagonal disorder displays a finite phase of extended states in the middle of the band of allowed energies, with two mobility edges separating localized and extended states chico . This result was experimentally validated by microwave transmission spectra through a single-mode waveguide with inserted correlated scatterers apl2 . The above results have motivated the study of other model systems that can be mapped onto the Anderson model, such as magnetic Fid1  and harmonic chains Fid2 .

In the context of stochastic process, Morgado et al. Morgado , studied diffusion in systems with long-range time correlation. First, they establish a direct connection between the noise density of states, $`\rho _n(\omega )`$, and the diffusive process. Second, they conjecture that the dynamics of a Hamiltonian system with space correlated disorder could be described by the same formalism if one supposes that the fluctuation in the density of states of the quasi-particle or elementary excitation, $`\rho _F(\omega )`$, plays the same role as $`\rho _n`$ in the stochastic description. In this work, we present a numerical analysis of the validity of this conjecture. We study it for the one-dimensional quantum Heisenberg model exhibiting long-range correlations in the random exchange couplings. For continuous ranges of the degree of correlation, this system presents super-diffusive and ballistic motions Fid1 . Here, we provide numerical evidence that the Hamiltonian description and the stochastic one can be unified through the referred conjecture, thus confirming early expectations Morgado .

## II Stochastic and Hamiltonian descriptions

Let us suppose that the equation of motion for an operator $`A`$ can be cast in the form Kubo ; Kubo2 ; Mori ; Lee2

$$\frac{dA(t)}{dt}=_0^t\mathrm{\Gamma }(tt^{})A(t^{})𝑑t^{}+h(t),$$
(1)
where $`h(t)`$ is a stochastic noise subject to the conditions $`h(t)=0`$, $`h(t)A(0)=0`$, and to the fluctuation-dissipation theorem Kubo :

$$C_h(t)=<h(t)h(0)>=<A^2>\mathrm{\Gamma }(t).$$
(2)
Here, $`\mathrm{}`$ indicates an ensemble average in thermal equilibrium. In principle, the presence of the kernel $`\mathrm{\Gamma }(t)`$ allows us to study a large number of correlated processes. For example, in analogy with the usual Langevin’s equation, we can study the asymptotic behavior of the second moment of the variable,

$$\sigma (t)=_0^tA(s)𝑑s,$$
(3)
namely

$$\underset{t\mathrm{}}{lim}\frac{<\sigma ^2(t)>}{t^\alpha }=K,$$
(4)
where $`K`$ is a constant. In Eq. (4), we have $`\alpha =1`$ for normal diffusion; for sub- and super-diffusion, $`\alpha <1`$ and $`\alpha >1`$, respectively.

The generalized field, $`h(t)`$, in Eq. (1) can be modeled by a thermal bath composed of harmonic oscillators; consequently, according to Eq. (2), the memory function can be written as

$$\mathrm{\Gamma }(t)=\rho _n(\omega )\mathrm{cos}(\omega t)𝑑\omega ,$$
(5)
where $`\rho _n(\omega )`$ is the noise density of states. It has been proved Morgado  that *If the Laplace transform of the memory function of a unidimensional system behaves as*

$$\stackrel{~}{\mathrm{\Gamma }}(z0)z^\nu ,$$
(6)
*then the diffusion exponent is*

$$\alpha =\nu +1.$$
(7)
In disordered Hamiltonian systems the diffusion process can be studied through direct integration of the equations of motion Fid1 ; Fid2 . Now then, how can we assure that the two approaches are compatible and lead to the same results? For those systems the density of states of the quasi-particle or of the elementary excitation, $`𝐃(E)`$, plays the most significant role. However, it displays fluctuations, which are intrinsically connected to the diffusive behavior. Here, we conjecture that, for the relaxation processes, the fluctuation in the density of states, $`\rho _F(E)`$, plays the same role as the noise density of states in the stochastic process, thus

$$\rho _n(E)\rho _F(E).$$
(8)
Consequently, if we calculate $`\rho _F(E)`$ we can obtain the memory function through Eq. (5) and the diffusion exponent using Eqs. (6) and (7). Below, we explore these ideas in the framework of the referred magnetic Hamiltonian system Fid1 .

The one-dimensional Heisenberg chain with long-range correlated random exchange couplings can be described by the Hamiltonian Fid1

$$H=\underset{l=1}{\overset{N}{}}J_l𝐒_l𝐒_{l+1},$$
(9)
where $`S=1/2`$, and the exchange integral is defined by

$$J_l=\underset{k=1}{\overset{N/2}{}}\mathrm{\Delta }(k)^{1/2}\mathrm{cos}(\pi kl/N+\varphi _k).$$
(10)
Here, $`\varphi _k`$ are random phases, and $`\mathrm{\Delta }(k)`$ is a power law spectrum in the k space given by

$$\mathrm{\Delta }(k)k^\beta .$$
(11)
By using a numerical renormalization group technique, it was predicted that this Hamiltonian supports a phase of low-energy extended spin waves in the strongly correlated regime $`\beta >1`$. Associated with these non-scattered modes, the spread of an initially localized wave packet displays a ballistic behavior in the long time limit Fid1 ; i.e. Eq. (4) with $`\alpha =2`$. In the weakly correlated regime, $`0<\beta <1`$, a super-diffusive behavior is obtained with $`\alpha =1.5`$. It is worth mentioning that the wave packet dynamics was investigated by means of a direct integration of the equations of motion using a fourth order Runge-Kutta method evangelou1 ; Fid1 . In this case the site spin operator $`𝐒_l`$ is here identified with operator $`A`$ in the stochastic formalism. Therefore if conjecture (8) is true the exponent of Eq. (7) should match the one calculated using the Hamiltonian dynamics Fid1 .

## III Numerical analysis

Now we describe the numerical method used to obtain the dynamical behavior, which consists in calculating the Laplace transform $`\stackrel{~}{\mathrm{\Gamma }}(z)`$ of the memory function $`\mathrm{\Gamma }(t)`$ defined by Eq. (5). However, we first compute the function $`\rho _F(E)`$, defined by

$$\rho _F(E)=<𝐃(E)^2>_C<𝐃(E)>_C^2,$$
(12)
where $`\mathrm{}_C`$ is an average over chains with different random sets of exchange integrals. To that end, we use Dean’s numerical method dean  to obtain the spin-wave density of states, $`𝐃(E)`$, and the corresponding $`\rho _F(E)`$, shown in Fig. 1, for $`\beta =0`$ and $`\beta =1.5`$. It is clear from the data that, for low energies, the fluctuations $`\rho _F(E)`$ in the correlated case, $`\beta =1.5`$ in this figure, are smaller than the uncorrelated case ($`\beta =0`$). For $`\beta >1`$, it was shown in ref Fid1  that the low energy states are extended, and that the density of states has a similar behavior as that of the pure chain. Therefore, we concluded that the finite value for $`\rho _F(E)`$ obtained in the correlated case ($`\beta =1.5`$) was not stable in the thermodynamic limit. In order to obtain the most probable value of the fluctuations in the correlated case, we perform a scaling analysis of the function $`\rho _F(E)`$ for chains with distinct lengths.

Hence, we consider the following scaling function

$$\mathrm{\Theta }(E,N_1,N_2)=\mathrm{exp}([\rho _F(E,N_1)N_1\rho _F(E,N_2)N_2]),$$
(13)
with $`N_1>N_2`$. This method has been successfully used to obtain the behavior of a harmonic chain in the thermodynamic limit Fid2 . In Fig. 2, we show the scaling as a function of the energy $`E`$. If $`\mathrm{\Theta }(E,N_1,N_2)`$ vanishes when we increase the chain size, the fluctuations are finite in the thermodynamic limit. However, if $`\mathrm{\Theta }(E,N_1,N_2)1`$, the fluctuations $`\rho _F(E)`$ decrease to zero when we increase the chain size. In order to clarify this figure, we interpolated the data using a Bezier interpolation function. In the insets we show the original data without interpolations. We use $`N_2=1.010^4`$, and three distinct values for $`N_1`$. Fig. 2(a) displays the results for $`\beta =0`$. From this figure we can see that $`\mathrm{\Theta }(E,N_1,N_2)1`$ in all energy ranges, except at $`E=0`$. This indicates the existence of an extended state at this energy, which is determinant in the super-diffusive behavior. Figure 2(b) shows $`\mathrm{\Theta }(E,N_1,N_2)`$ for $`\beta =1.5`$. In the range of low extended states, $`0<E<E_c`$, $`\mathrm{\Theta }(E,N_1,N_2)1`$, meaning that $`\rho _F(E)0`$ in the thermodynamic limit. Above $`E_c`$, $`\mathrm{\Theta }(E,N_1,N_2)1`$. The presence of a finite range of energy with extended states, $`\rho _F(E)=0`$ for $`0<E<E_c`$ ($`E_c4`$), is responsible for the ballistic behavior Morgado .

Now we proceed by numerically integrating $`\rho _F(E)`$ in order to obtain the limiting behavior of $`\stackrel{~}{\mathrm{\Gamma }}(z)`$ as $`z0`$, using $`E_C=4`$. In Fig. 3, we show $`\stackrel{~}{\mathrm{\Gamma }}(z)`$ as a function of $`z`$ for $`\beta =1.5`$. Observe the verification of the conjecture: $`\nu 1.0`$, implying $`\alpha 2.0`$. We find super-diffusive motion for $`\beta <1`$, and ballistic motion for $`\beta >1`$. Superdiffusive propagation of a wave packet was also obtained in the one-dimensional Fibonacci and Thue-Morse lattices, which present deterministic aperiodic order. Since the numerical method has greater errors for smaller values of $`E`$, we find it difficult to obtain the density of states for small $`E`$ (large times). This sort of error is expected to have an influence on small values of $`z`$ in the Laplace transform. Another numerical problem is due to the finite size of the chain, which has important implications for $`\beta >1`$, as seen by the presence of the plateau in Fig. 1.

Finally, Fig. 4 shows $`\alpha `$ as a function of $`\beta `$. We select $`0<\beta <2`$, covering both regimes. The numerical simulations are quite compatible with the stepwise behavior. The maximum deviation observed is at $`\beta =0.8`$, with an error of about $`8\%`$, which results in $`\alpha =1.39`$.

## IV Discussions and Conclusions

The main results reported in this work are the following: first, for the Heisenberg system we have calculated the superdiffusive and ballistic exponents with reasonable precision. Second, we have been able to associate a memory function to a Hamiltonian, through which we unite two powerful formalisms: a stochastic description and the quantum equations of motion. This result may have implications far beyond the specific example outlined in this work. Evidently, we need more simulations and experiments connecting the exponents $`\alpha `$ and $`\nu `$ in several diffusive processes in order to have a more complete picture of the validity of conjecture (8). It is important to notice that there are many phenomena associated with those described here. For example, in the ballistic regime ($`\alpha =2`$), there remains some open issues concerning the use of the fluctuation-dissipation theorem Costa . Since the early experiment of Kauzmann Kauzmann , situations have been found where the systems do not thermalize, usually associated with their slow dynamics Costa ; Kauzmann ; Grigera ; Parisi ; Rubia ; Ricci . This slow dynamics is due to the softening of the lower fluctuating modes, here identified by the presence of a finite range of energy with extended states: $`\rho _F(E)=0`$ for $`0<E<E_c`$. This effect has been used recently Bao  to obtain ballistic diffusion in ratchet devices. Another related phenomenon is chaos synchronization Oliveira3 , which presents a reduction in the allowed phase space. A large time to achieve equilibrium, or the absence of some regions of the phase space, leads to the same effect. Some of its important consequences, such as violation of ergodicity Costa ; Lee3  and of the fluctuation-dissipation theorem Costa , have been observed.

###### Acknowledgements.

We thank M.L. Lyra for useful discussions. Work supported by CAPES, CNPq, FINEP, FAPDF, FACEP and FAPEAL (Brazilian Federal and States research agencies).


text:
# References

Comment on “Conservative Quantum Computing”

In M. Ozawa (2002) Ozawa considers the limitations imposed by conservation laws on the possibility of accurately generating quantum logic gates. In particular, he argues that conservation of angular momentum imposes a fundamental limit on how accurately the controlled-not (CNOT) gate can be performed. His argument runs roughly as follows. Since CNOT does not commute with total angular momentum, it cannot be implemented directly in a system in which angular momentum is conserved. Therefore one must enlarge the two-qubit Hilbert space by attaching ancilla qubits, perform a unitary operation on the resulting enlarged Hilbert space, then trace out the ancillae. This defines a completely positive dynamical quantum map, which can only yield the CNOT gate on the original two qubits *probabilistically*. Ozawa calculates upper bounds on this probability of success, and concludes that the accuracy of quantum logic gates is fundamentally limited in the presence of conservation laws.

The purpose of this comment is to point out that this is an overly restrictive conclusion. Indeed, as Ozawa notes at the end of his paper: “The present investigation suggests that the current choice of the computational basis should be modified so that the computational basis commutes with the conserved quantity.” “… we may find such a computational basis comprised of orthogonal entangled states over a multiple-qubit system.” Such a modification of the computational basis is already well known. This is possible since while it is true that the CNOT gate sometimes does not commute with the total angular momentum operator, it does *not* follow that CNOT cannot be generated from an interaction that is rotationally invariant. It can, by using an encoding into multi-qubit states that define an *invariant subspace* with respect to the interaction at hand. As first shown in D. Bacon, J. Kempe, D.A. Lidar and K.B. Whaley (2000), using an encoding of a logical qubit into four physical qubits, *the isotropic Heisenberg exchange interaction* $`H_{\mathrm{Heis}}=J(X_iX_j+Y_iY_j+Z_iZ_j)`$, which commutes with total angular momentum, *is universal for quantum computing*. Here $`X_i`$ is the Pauli $`\sigma _x`$ matrix acting on qubit $`i`$, etc. This means that Heisenberg exchange can be used to generate the group $`U(2^K)`$ on $`K`$ encoded qubits, and, in particular, it can be used to generate an encoded CNOT gate, which treats only two out of the 16 dimensions of the 4-physical-qubit Hilbert space as computational basis states of an encoded qubit. This result was soon followed by a general proof that the Heisenberg exchange interaction is by itself universal over encodings into any number of qubits $`n3`$ J. Kempe, D. Bacon, D.A. Lidar and K.B. Whaley (2001). Specific gate sequences were then proposed for the $`n=3`$ qubit encoding D.P. DiVincenzo, D. Bacon, J. Kempe, G. Burkard, and K. B. Whaley (2000). Furthermore, by supplementing Heisenberg exchange with a Zeeman splitting $`\epsilon (Z_iZ_j)`$ (another interaction that commutes with the total angular momentum operator), a simple encoding into $`n=2`$ qubits, $`|0_L=|0_i|1_j`$ and $`|1_L=|1_i|0_j`$, already suffices for universal quantum computation J. Levy ; D.A. Lidar and L.-A. Wu (2002), and thus to generate an encoded CNOT gate. In fact, other exchange interactions that admit conserved quantities (e.g., the Hamiltonians $`H_{\mathrm{XY}}=J(X_iX_j+Y_iY_j)`$ and $`H_{\mathrm{XXZ}}=J(X_iX_j+Y_iY_j)+J_zZ_iZ_j`$, that have axial symmetry) are also universal for quantum computation, with D.A. Lidar and L.-A. Wu (2002) or without Zeeman splitting J. Kempe, D. Bacon, D.P. DiVincenzo and K.B. Whaley (2001).

The reason that encoding helps is that it enables the quantum logic gates to be executed on an invariant subspace of the full Hilbert space of the original qubits plus ancillae. The exchange Hamiltonians, and hence the encoded logic gates, preserve this subspace and, as is evident from their universality, commute with all the symmetries of the system. Thus, Ozawa’s result can be seen as a confirmation of the advantage that may be had by allowing a flexible definition of the computational basis: conservation laws do not impose fundamental limitations on the accuracy of quantum logic operations, as long as one uses Hamiltonians that act on an encoding of logical qubits into subspaces that are invariant under the symmetries generating the conservation laws.

This work was sponsored by the DARPA-QuIST program (managed by AFOSR under agreement No. F49620-01-1-0468). I am grateful to Prof. M. Ozawa for helpful correspondence on some the issues raised here.

| Da | niel A. Lidar |
| --- | --- |
|  | Chemical Physics Theory Group |
|  | Chemistry Department |
|  | 80 St. George St. |
|  | University of Toronto |
|  | Toronto, Ontario, Canad M5S 3H6 |


text:
# Nanolithography with metastable helium atoms in a high-power standing-wave light field

## I Introduction

Since the beginning of the previous decade, atom-optical manipulation of neutral atoms has been used to create nanoscale structures (for a recent review of atom nanofabrication, see Meschede and Metcalf Meschede and Metcalf (2003)). Two different processes can be distinguished for the fabrication of such structures. In a deposition process nanoscale structures are directly grown on a surface. Direct deposition experiments have been performed with sodium Timp et al. (1992), chromium McClelland et al. (1993) and aluminum McGowan et al. (1995). In a lithography procedure, on the other hand, the atoms first locally damage an organic resist layer and next the pattern is transferred to an underlying metal film in a wet-etching process Xia et al. (1995). Although this last technique is a two-step process, it has the advantage that in principle background-free structures can be created, and feature broadening caused by the growth process is avoided. To create a pattern on the sample, a mechanical mask (e.g. an electron microscope grid) or an optical mask (standing light wave) can be used. A mechanical mask partly blocks the atomic beam, and the grid structure is imaged onto the sample. Using grids in a lithography process, microscale structures have been created with metastable argon Berggren et al. (1995) and metastable helium Nowak et al. (1996); Lu et al. (1998), showing an edge resolution of 100 nm and 40 nm respectively. Atom lithography with an optical mask has been demonstrated with cesium Lison et al. (1997), metastable argon Johnson et al. (1998), and metastable neon Engels et al. (1999). In the cesium and argon experiments, structures with a feature size of respectively 120 nm and 65 nm were realized. In the latter experiment the atoms were optically quenched by the light field and an exposure time of eight hours was required. In this paper we report lithography measurements with a mechanical and an optical mask for metastable helium atoms. Helium has the highest internal energy (20 eV) of all metastable atoms. Therefore, the dose of atoms required to damage a resist molecule can be small, compared to other metastable species. This benefits the exposure time of the sample to the atomic beam.

An optical mask, formed by a standing-wave light field, acts as an array of cylindrical lenses, attracting the atoms to the potential minima of the light field. Two different regimes can be distinguished. A low-power regime, where the atoms are focused at the center of the Gaussian light beam, and a high-power regime, where the atoms oscillate around the potential minima of the light field. The lithography experiments mentioned above were all performed in the focusing regime using low-power optical masks. Theoretically, this regime can be described straightforwardly, since an analytical expression for the dipole force can be derived from the optical potential of the light field Ashkin (1978).

We have performed lithography experiments with a high-power standing-wave light field, where the simple expression for the dipole force is not valid. This is due to the fact that the dipole force becomes velocity dependent, because the transverse velocity of the atoms, channeling through the standing-wave light field, can no longer be neglected. In this high-power regime, the dipole force can be described by a numerical model, first presented by Minogin and Serimaa Minogin and Serimaa (1979). In a previous paper we have compared this model with the model that describes the conventional dipole force Petra et al. (2003). In the high-power regime, the background due to non-focused atoms is significantly reduced compared to the focusing regime. Also the alignment of the standing-wave laser beam with respect to the sample is far less critical. The sample can be placed behind the center of the laser beam, which reduces disturbances of the standing wave pattern due to diffraction of the laser beam on the sample edge significantly Anderson et al. (1999). Numerically calculated trajectories of atoms channeling through a standing-wave light field show an increased confinement of the atoms in the optical well along the atomic beam axis, due to a dissipative force that reduces the transverse velocity of the atoms. The Full Width at Half Maximum (FWHM) of the atom distribution on the sample is calculated to be 40 nm for parameters as used in our experiment Petra et al. (2003).

## II Experimental setup

The experimental setup is shown in figure 1. The metastable helium atoms are produced in a water-cooled DC discharge Rooijakkers et al. (1996). In the discharge helium atoms are transferred from the $`1{}_{}{}^{1}\text{S}_{0}^{}`$ ground state to the $`2{}_{}{}^{3}\text{S}_{1}^{}`$ metastable state, which has a lifetime of 8000 s. The atoms have a mean longitudinal velocity of 2000 m/s, which corresponds to an energy of 0.08 eV. The atomic beam is transversely cooled in two dimensions using curved-wavefront Doppler cooling Rooijakkers et al. (1996); Aspect et al. (1990). The collimation of the atomic beam reduces the transverse velocity spread of the atoms to about 3 m/s, corresponding to a beam divergence of 9 mrad. In the high-power regime, the guiding of the atoms through the standing wave is not very sensitive to the transverse velocity of the atoms, and calculations have shown that the divergence is sufficiently small Petra et al. (2003). The collimation of the beam enhances the flux of the atomic beam by a factor of 6 to $`1.1\times 10^{10}`$ atoms s-1 mm-2. The light for transverse cooling (20 mW) is provided by a DBR (Distributed Bragg Reflector) diode laser. This laser is locked to the $`2{}_{}{}^{3}\text{S}_{1}^{}2{}_{}{}^{3}\text{P}_{2}^{}`$ transition of the helium atom using Doppler-free saturation spectroscopy with a lock-in technique. The wavelength of the light corresponding to this transition is 1083 nm. The bandwidth of the laser is experimentally determined to be 2 MHz, which is comparable to the 1.6 MHz linewidth of the $`2{}_{}{}^{3}\text{P}_{2}^{}`$ excited state.

The standing-wave light field is obtained by retro-reflecting an 800 mW light beam, provided by a fiber laser (IPG, model YLD-1BC). This laser also operates at a wavelength of 1083 nm and is 375 MHz blue detuned from the $`2{}_{}{}^{3}\text{S}_{1}^{}2{}_{}{}^{3}\text{P}_{2}^{}`$ transition to reduce spontaneous emission of the atoms in the standing-wave light field. Higher detunings are unfavorable, because the atoms can then also interact with the light via the $`2{}_{}{}^{3}\text{S}_{1}^{}2{}_{}{}^{3}\text{P}_{1}^{}`$ transition. The frequency of this transition is only 2.3 GHz higher than the frequency of the $`2{}_{}{}^{3}\text{S}_{1}^{}2{}_{}{}^{3}\text{P}_{2}^{}`$ transition.

The fiber laser is unlocked but its frequency with respect to the locked diode laser is monitored with a Fabry-Perot interferometer, whose cavity length is scanned with a piezoelectric transducer. The fiber laser has a measured short-term linewidth of 1 MHz and a good long term stability. Only drifts on the order of 10 MHz are visible with respect to the locked diode laser during the exposure times. The laser beam is focused onto a retro-reflecting mirror to a size ($`1/\text{e}^2`$ radius) of 0.33 mm, which translates to a saturation parameter (defined as the intensity of the light divided by the saturation intensity of the atomic transition) of $`10^7`$. The retro-reflecting mirror and the sample are mechanically fixed in their relative position onto the sample holder, which is made of a solid aluminum block. For each experiment a substrate is mounted on the sample holder in free air. The sample holder is then introduced into the vacuum chamber via an air-lock. The exposure dose is determined by measuring the current through a fixed Faraday cup due to the metastable atoms, before and after the exposure Rooijakkers et al. (1996).

The laser beam is clipped by the sample, located $`100\mu `$m from the center of the laser beam, downstream the atomic beam. Numerical calculations that take into account the longitudinal and transverse velocity spread of the atomic beam have shown that this is the best position for the sample Petra et al. (2003).

For proper alignment of the standing wave, the laser beam is first adjusted such that it is not clipped by the sample and thus 100% of the light is back-reflected by the standing-wave mirror onto the power meter (see figure 1). Next, the beam is moved by a micrometer translation stage (not indicated in figure 1) and clipped by the sample such that only 50% of the light is back-reflected, which indicates that the sample position is at the center of the laser beam. Finally, additional fine tuning is done with the translation stage to move the laser beam to the desired position with respect to the sample.

The samples used in the experiment are silicon substrates covered with a 1 nm chromium layer and a 30 nm gold layer. These samples are cleaned with pure ethanol and with a solution of concentrated sulfuric acid (98% pure) and hydrogen peroxide (30% pure). In the experiments two different resist layers are used: nonanethiol \[CH3(CH2)8SH\] and dodecanethiol \[CH3(CH2)11SH\]. A Self-Assembled Monolayer (SAM) of either of these alkanethiols is deposited onto the substrates by leaving the samples overnight in a 1 mM solution of the alkanethiol in ethanol. The SAM forms a hydrofobic resist layer that protects the gold layer in the wet etching process.

For nonanethiol and dodecanethiol the samples are exposed to the atomic beam for eight minutes and ten minutes respectively. This corresponds to a dose of about one metastable atom per resist molecule. During the exposure, the SAM is locally damaged and loses its hydrofobic character. After being exposed, the samples are wet etched in a cyanide etching solution Nowak et al. (1996). The etching times for nonanethiol and dodecanethiol are also eight and ten minutes respectively. In the etching process, the SAM molecules that are damaged due to the impact of a metastable helium atom are removed, as well as the gold film underneath these damaged molecules. Longer etching times can affect areas untouched by the metastable helium atoms Xia et al. (1995). The samples are finally analyzed with an Atomic Force Microscope (AFM).

## III Results and discussion

A typical result of a sample image produced with a standing-wave light field is shown in figure 2. The image shows a $`10\mu \text{m}\times 10\mu \text{m}`$ scan of the substrate, taken with an AFM. The black lines indicate the positions where the metastable helium atoms have hit the sample and where the gold layer is removed in the etching process. The average distance between two successive lines is 542 nm, which is exactly half of the wavelength of the standing-wave light field (1083 nm). A more quantitative analysis of the pattern is made in figure 3. This figure shows that the average FWHM of the lines in figure 2 is about 100 nm. The results are the same for both the nonanethiol and dodecanethiol SAM, although the latter requires a 25% longer exposure and etching time. The irregular structure of the black lines is probably caused by the graininess of the gold layer, which is also visible on AFM images of untreated samples. For samples that are under-exposed and/or under-etched, the black lines are broken, while samples that are over-exposed and/or over-etched show connections between two successive lines. The chemistry involved in the experiments makes the reproducibility of the measurements difficult. In some cases no line structures were found on the sample although all physical parameters of the experiment were kept identical. When a sample is visibly inspected by eye, diffracted light indicates that the sample contains well written nanostructures.

To compare the measured structures with theory, numerical simulations using the Minogin and Serimaa model are shown in figure 4. The input parameters of these calculations correspond to the values used in the experiment. Based on these calculations, a FWHM of a single line of 40 nm is expected, which is a factor of 2.5 smaller than observed experimentally.

The discrepancy between the measured and the calculated width is attributed to both the graininess of the gold layer and the wet-etching process, which is an isotropic process that also removes the gold layer from inside a line structure during etching. This conclusion is in agreement with experiments performed with a mechanical mask, shown in figure 5. In these measurements a grid of 1000 mesh (1000 lines per inch) is used instead of an optical mask. From analysis with the AFM a $`40^{}`$ inclination in the gold layer is found, in agreement with an earlier observation Nowak et al. (1996). The AFM tip has a cone angle of $`10^{}`$ and does therefore not limit the measured angle.

The $`40^{}`$ side wall inclination translates to a step-edge of 40 nm and increases the expected FWHM linewidth of a single line with 40 nm. Nanolithography with a wet-etching technique incorporating a SAM as a resist layer therefore seems to be limited to structure sizes larger than 50 nm.

Comparing metastable helium nanolithography in the high-power optical mask regime with results employing other elements shows some clear advantages. First, the high internal energy of metastable helium atoms allows a relatively short exposure time of eight minutes. Second, no pedestal due to non-focused atoms is visible. Third, the high power available nowadays from fiber lasers allows an extremely large saturation parameter ($`10^7`$) and deep penetration in the channeling regime. In this regime a feature size of 1/11th of the optical-mask wavelength is observed, which is hard to realize in the focusing (low-power) regime, due to a large sensitivity to the longitudinal and transversal velocity spread in the atomic beam in this regime, and because the sample position with respect to the standing wave is very critical.

###### Acknowledgements.

The authors thank K.A.H. van Leeuwen for fruitful discussions. Financial support from the Foundation for Fundamental Research on Matter (FOM) is gratefully acknowledged.