text:
# Quadrupolar interactions in Pr compounds: PrFe4P12 and PrBa2Cu3O6

## 1 Introduction

The highly degenerate $`f`$-shells of rare earth and actinide ions support a great variety of local degrees of freedom: magnetic (dipolar), quadrupolar, and octupolar. For localized $`f`$-electron systems, these (and some of their combinations) are potential order parameters, leading to rich phase diagrams, as for CeB6 . Field-induced octupoles are known to influence the shape of the phase boundaries for CeB6 or PrPb3 , but it is an outstanding question whether breaking time-reversal invariance with purely octupole order is possible; recent observations on NpO2 indicate that it may be so . The $`f`$-shell multipolar degrees of freedom may also interact and mix with those of other electrons: the unusual heavy fermion behavior of PrFe4P12 , and the exotic superconducting phases of PrOs4Sb12  have been tentatively ascribed to Pr3+ quadrupole fluctuations.

Detailed understanding has been achieved for CeB6 and some other Ce-based systems. The reason is that Ce3+ has only a single $`f`$-electron, and therefore microscopic analysis is here the easiest. Of particular interest is cubic CeB6 in which the fourfold degenerate $`\mathrm{\Gamma }_8`$ ground state supports the full array of 3 dipoles, 5 quadrupoles, and 7 octupoles . A microscopic analysis  indicates that octupole–octupole interactions should not be weaker than quadrupole–quadrupole interactions, which gives rise to a close competition between different kinds of order. In fact, if all interactions were exactly equal, the model would be SU(4)-symmetrical, which sparked early theoretical interest in the subject . The realistic multipolar model is not fully symmetrical, but the nearness of the SU(4) point is still felt in fluctuations . Whether the ground state of an SU(4) model is long range ordered, or it is rather liquid-like (the SU(4) version of RVB), depends on details like the lattice structure, and the range of the interaction . The approximate realization of high symmetries like SU(4) is of current interest for $`d`$\- and $`f`$-electron systems alike.

We expect that the result that various multipolar couplings are of the same order of magnitude, holds generally for light rare earths. The question arises what can we say about $`4f^2`$ systems.

Pr compounds show varied behavior. A few examples: PrPb3 undergoes an antiferroquadrupolar transition but it does not order magnetically, though both dipolar and octupolar interactions are thought to be important . PrB6 has two magnetic phase transitions but the magnetic anisotropy is explained by quadrupolar interactions . PrBa2Cu3O7-x is the only member of the group (RE)Ba2Cu3O7 (RE=rare earth), which is not a high-$`T_C`$ superconductor . It is suspected that this is related to the little-understood Pr ordering transition at 17K, but the connection is not clear. The parent Mott insulator PrBa2Cu3O6 has a similar phase transition at 11K .

## 2 Pr systems

In most compounds, Pr is trivalent. The existence of the dioxide PrO2 shows that it can also be tetravalent , though this is rare. The question of the valence state becomes crucial whenever the nature, or the driving force, of Pr ordering is in doubt. For Pr-filled skutterudites: the formation of a heavy Fermi sea, as seen in PrFe4P12 and in PrOs4Sb12, is most likely to involve the hybridization of Pr $`f`$-states with conduction band states, and consequently deviation from strictly integral valence.

These basic doubts notwithstanding, we propose to study the interplay of dipolar and quadrupolar interactions in localized $`f`$-electron models, and apply the results to PrFe4P12 and PrBa2Cu3O6. We assume the integer valence state Pr3+, and a Hund’s rule ground state with $`J=4`$. We note that a similar approach was successful in the case of CeB6 , though the system is known to possess the features of a dense Kondo system. We envisage a similar study of Pr systems.

Pr3+ is a non-Kramers ion. It follows that a unique ground state can be reached without breaking time-reversal invariance, either by having a singlet $`\mathrm{\Gamma }_1`$ crystal field ground state to begin with, or via the ordering of quadrupolar moments.

### 2.1 The $`\mathrm{\Gamma }_1`$$`\mathrm{\Gamma }_4`$ scenario for PrFe4P12

PrFe4P12 undergoes a second order transition at $`T_{\mathrm{tr}}=6.51\mathrm{K}`$, which is manifested in a susceptibility spike, and a $`\lambda `$-shaped anomaly of the specific heat . Though at first thought to be antiferromagnetic, the ordered phase is now known to be antiferroquadrupolar. Upon switching on a magnetic field, the transition temperature decreases. The character of the transition switches to first order at the tricritical point $`T_{\mathrm{tri}}5\mathrm{K}`$, $`H_{\mathrm{tri}}2\mathrm{T}`$, and the ordered phase is completely suppressed at $`4\mathrm{T}`$ . These parameters depend on the direction of the applied field. However, here we consider only the case of a field applied along the (100) direction.

The question whether the $`f`$-states are itinerant or localized can be posed again. PrFe4P12 is one of the few Pr-based heavy fermion systems . Though having an even number of $`f`$-electrons per site makes the distinction between large and small Fermi surfaces difficult, we take the view that the large Fermi sea of heavy $`f`$-electrons competes with a state in which the conduction electrons build a small Fermi sea, and multipolar inter-shell interactions are switched on . This is supported by inelastic neutron scattering finding that the crystal field excitations become sharply defined when the temperature is lowered below $`T_{\mathrm{tr}}`$. It follows that we may use a localized $`f`$-electron model to describe the ordered phase and the vicinity of the phase boundary. We use the crystal field states of the octahedral group $`O_h`$. We note that the icosahedron of the twelve P atoms surrounding a Pr site gives a novel tetrahedral component to the crystal field potential , but we assume that using cubic symmetry labelling is an acceptable approximation.

Since the local order parameter is one of the quadrupolar moments, it may look evident that the crystal field ground state should carry a (spontaneous) quadrupolar moment. However, an analysis of the high-$`T`$ behavior of the magnetization curves shows that the level scheme (ground state)–(first excited state) may be either of the following three: $`\mathrm{\Gamma }_1`$$`\mathrm{\Gamma }_4`$, $`\mathrm{\Gamma }_1`$$`\mathrm{\Gamma }_5`$, or $`\mathrm{\Gamma }_3`$$`\mathrm{\Gamma }_4`$ . Only the last, with the $`\mathrm{\Gamma }_3`$ ground state, carries an unfrozen quadrupolar moment which can be ordered by switching on (arbitrarily small) intersite interactions. It was also shown that this choice is consistent with a symmetry analysis of the structural distortion accompanying the antiferro-quadrupolar ordering . This latter argument relies only on the assumption of the $`\mathrm{\Gamma }_3`$ ground state, and does not refer to the nature of the first excited state.

Here we show that the alternative $`\mathrm{\Gamma }_1`$$`\mathrm{\Gamma }_4`$ scheme is also capable to account for most of the observed static properties of PrFe4P12. We list the crystal field states

$`\mathrm{\Gamma }_1`$ $`=`$ $`\sqrt{5/24}(|4+|4)+\sqrt{7/12}|0`$
$`\mathrm{\Gamma }_4^\pm `$ $`=`$ $`(1/4)(|3\pm |3+\sqrt{7}[\pm |1+|1])`$
$`\mathrm{\Gamma }_4^0`$ $`=`$ $`\sqrt{1/2}(|4|4).`$ (1)
where we use the quadrupolar eigenstates $`\mathrm{\Gamma }_4^\pm `$ in the representation $`\mathrm{\Gamma }_4`$. Since the ground state carries no kind of moment, the ordered quadrupolar moment must be interaction-induced, implying that the interaction has to exceed a threshold value of the order of the $`\mathrm{\Gamma }_1`$$`\mathrm{\Gamma }_4`$ splitting. Taking it for granted that intersite interactions mix the $`\mathrm{\Gamma }_4`$ states with $`\mathrm{\Gamma }_1`$, the possible local order parameters are those appearing in the decomposition

$$(\mathrm{\Gamma }_1\mathrm{\Gamma }_4)(\mathrm{\Gamma }_1\mathrm{\Gamma }_4)=2\mathrm{\Gamma }_1+\mathrm{\Gamma }_3+3\mathrm{\Gamma }_4+\mathrm{\Gamma }_5.$$
(2)
The system is capable of either $`\mathrm{\Gamma }_3`$ or $`\mathrm{\Gamma }_5`$ type quadrupolar ordering, or any (dipolar or octupolar) kind of $`\mathrm{\Gamma }_4`$ order. In principle, we have to allow for all the corresponding intersite interactions. However, we neglect the octupoles. We assume $`\mathrm{\Gamma }_3`$-type quadrupolar order and, recalling that the pair interaction has only tetragonal symmetry, and so $`𝒪_2^2`$ and $`𝒪_2^0`$ are not equivalent, we arbitrarily keep only the $`𝒪_2^2=J_x^2J_y^2`$ term. The following mean-field-decoupled hamiltonian includes also a dipole-dipole coupling:

$``$ $`=`$ $`_{CF}+_{\mathrm{Zeeman}}+_{\mathrm{quad}}+_{\mathrm{dipole}}`$ (3)
$`=`$ $`\mathrm{\Delta }{\displaystyle \underset{n=0,\pm }{}}|\mathrm{\Gamma }_4^n\mathrm{\Gamma }_4^n|g\mu _\mathrm{B}\mathrm{𝐇𝐉}zQ𝒪_2^2_{B(A)}𝒪_2^2zI𝐉_{B(A)}𝐉`$
where we allowed for two-sublattice (A and B) order on the bcc lattice ($`z=8`$). $`g=4/5`$ is the Landé factor, $`Q`$ the quadrupolar, and $`I`$ the dipolar coupling constant. We mention that a simplified interaction, with quadrupolar coupling only, was treated in a previous communication . However, the inclusion of $`_{\mathrm{dipole}}`$ is important for getting a better fit to the experimental results, in particular the susceptibility.

The diagonalization of (3) in the basis (1) is straightforward. We show only the results. We have found that the (ferromagnetic) dipolar coupling $`I=106\mathrm{m}\mathrm{K}`$, the antiferro-quadrupolar coupling $`Q=9.5\mathrm{mK}`$, and the crystal field splitting $`\mathrm{\Delta }=3\mathrm{K}`$ give good overall agreement with the observations (in case $`Q`$ appears implausibly small, check that $`𝒪_2^2`$ has big matrix elements). With these parameters, the only phase transition is the onset of antiferro-quadrupolar order at $`T_{\mathrm{tr}}=6.5\mathrm{K}`$. The transition is gradually suppressed by an external magnetic field (Fig. 1, left). We locate the tricritical point at $`T_{\mathrm{tri}}5\mathrm{K}`$, $`H_{\mathrm{tri}}2\mathrm{T}`$, in agreement with experiment. The same set of parameters gives also the observed Weiss temperature $`\mathrm{\Theta }=3.6\mathrm{K}`$ from the intercept in the inverse susceptibility (Fig. 1, middle).

We note that our estimate of the crystal field splitting $`\mathrm{\Delta }3\mathrm{K}`$ is substantially smaller than the one quoted in ; considering the scale of $`T_{\mathrm{tr}}`$, we have a quasi-quadruplet. However, this is not necessarily unreasonable: a calculation indicates a region where the gaps become very small (H. Harima, private communication).

To conclude our review of the results from the $`\mathrm{\Gamma }_1`$$`\mathrm{\Gamma }_4`$ scheme111We did not make a comparable study of the full range of models, so we could not claim that the present model gives best agreement with experiments. Rather, we wished to demonstrate that the $`\mathrm{\Gamma }_1`$$`\mathrm{\Gamma }_4`$ scheme is not obviously excluded, in spite of the absence of ground state quadrupolar moments., we show the magnetization curves for three representative temperatures (Fig. 2). In the region of first order transitions below the tricritical temperature, one finds a sharp metamagnetic transition which becomes smoother in the regime of continuous transitions. These results show a good overall resemblance to the measured curves .

### 2.2 The role of quadrupolar coupling in PrBa2Cu3O6

All the (RE)Ba2Cu3O7 systems, with the exception of RE=Pr, are $`T_C90\mathrm{K}`$ superconductors. Whenever the RE3+ ion has a magnetic moment, they also undergo RE antiferromagnetic ordering at $`T_\mathrm{N}=12\mathrm{K}`$. Analogous results hold for the parent Mott insulators (RE)Ba2Cu3O6, where most $`T_\mathrm{N}`$ follow de Gennes scaling , the peak belonging to RE=Gd with $`T_\mathrm{N}2.2\mathrm{K}`$ . This is well understood on the basis of the $`I_{\mathrm{Gd}\mathrm{Gd}}=156\mathrm{m}\mathrm{K}`$ meaured by ESR . We would expect that for RE=Pr,

$$T_\mathrm{N}(\mathrm{Pr})|_{\mathrm{deGennes}}=\frac{\left\{(g1)^2J(J+1)\right\}_{\mathrm{Pr}}}{\left\{(g1)^2J(J+1)\right\}_{\mathrm{Gd}}}T_{\mathrm{tr}}(\mathrm{Gd})0.1\mathrm{K}.$$

Instead, magnetic ordering and, at the same time, tetragonal-to-orthorombic distortion, is observed at $`T_\mathrm{N}17\mathrm{K}`$ for PrBa2Cu3O7, and $`T_\mathrm{N}11\mathrm{K}`$ for PrBa2Cu3O6 , exceeding the de Gennes estimate by two orders of magnitude. There must be a mechanism for producing substantially stronger exchange, or a different kind of ordering, or both.

ESR on Gd:PrBa2Cu3O6 gives the intersite exchange $`I_{\mathrm{Pr}\mathrm{Gd}}=140\mathrm{m}\mathrm{K}`$ . De Gennes scaling would lead us to expect this sign, but also an absolute value of about a factor of 4 smaller. The greater spatial extent of Pr $`4f`$ orbitals is likely to explain the enhancement of $`I`$; correspondingly, we would expect that $`I_{\mathrm{Pr}\mathrm{Pr}}`$ is even more enhanced, but would probably still fall short of accounting for the observed $`T_\mathrm{N}`$. Here we suggest that the presence of quadrupolar interactions may explain several features of the magnetic behavior of PrBa2Cu3O6, and also an additional enhancement of $`T_\mathrm{N}`$.

All previous works agree that in the level scheme of Pr3+ ions, a low-lying quasi-triplet consisting of the tetragonal doublet

$$\mathrm{\Gamma }_t^\pm =\alpha |\pm 3\beta |1$$
(4)
and the singlet

$$\mathrm{\Gamma }_t^0=\sqrt{1/2}(|2|2)$$
(5)
is well separated from the remaining six $`J=4`$ states, and therefore suffices for modelling (as far as the Pr sites are concerned) all low-$`T`$ phenomena . With $`\alpha =\sqrt{7/8}0.9354`$, (4) and (5) would constitute the cubic $`\mathrm{\Gamma }_5`$ triplet, but our fits yield $`\alpha 0.943`$, indicating a slight admixture from the doublet derived from $`\mathrm{\Gamma }_4`$. Let us observe that the doublet carries both $`J^z`$ dipole, and $`𝒪_2^2`$ (or alternatively $`𝒪_{xy}`$) type quadrupolar moment. If the doublet is the crystal field ground state, the degeneracy can be resolved either by magnetic, or by quadrupolar, ordering.

Using (4)–(5), and following the procedure of , we fitted the high-$`T`$ magnetization curves, and found that the dipole–dipole coupling has a strong planar anisotropy, and there is also a substantial ferroquadrupolar coupling of the $`𝒪_2^2`$ moments . The latter finding may lead us to ask whether the observed transition is perhaps purely of quadrupolar nature; however, neutron scattering shows that $`T_\mathrm{N}`$ is indeed a Néel temperature, with the $`T<T_\mathrm{N}`$ ordered moments strongly tilted out of the tetragonal $`c`$-direction . Within the crystal field model, it is rather mysterious why the system does not take advantage of the permanent $`c`$-axis moments of the doublet (4), and chooses instead $`ab`$-plane moments which have to be interaction-induced.

The mean field hamiltonian acting on the quasi-triplet is similar to (3)

$``$ $`=`$ $`_{CF}+_{\mathrm{Zeeman}}+_{\mathrm{quad}}+_{\mathrm{dipole}}`$ (6)
$`=`$ $`\mathrm{\Delta }|\mathrm{\Gamma }_t^0\mathrm{\Gamma }_t^0|g\mu _\mathrm{B}\mathrm{𝐇𝐉}zQ𝒪_2^2𝒪_2^2zIJ_x_{B(A)}J_x`$
where we allow for two-sublattice antiferromagnetism, but assume ferroquadrupolar coupling. In view of the bilayer structure of PrBa2Cu3O6, we take $`z=4`$. From our susceptibility fits, we estimate $`\mathrm{\Delta }20\mathrm{K}`$.

We note that a closely related three-state problem was considered by Líbero and Cox . They assume hexagonal symmetry, so their choice of the singlet and the doublet is quite different from ours, but nevertheless, when worked out fully, the results of their mean field theory should basically correspond to what we find. However, we calculate different quantities (such as the susceptibilities), and emphasize different aspects.

Fig. 3 (left) shows a representative phase diagram. At $`Q=0`$, the system gains energy from antiferromagnetic ordering, with the moments lying in the $`\pm x`$-direction. However, it is easy to convince ourselves that at the same time, the system develops quadrupolar order: $`J^x`$ connects $`\mathrm{\Gamma }_t^0`$ to $`\mathrm{\Gamma }_t^+\mathrm{\Gamma }_t^{}`$ (which is one of the quadrupolar eigenstates), but not to $`\mathrm{\Gamma }_t^++\mathrm{\Gamma }_t^{}`$ (which would be the other). Thus, for zero or weak $`Q`$, the primary order parameter is $`J_x`$, while $`𝒪_2^2`$ is a secondary order parameter (Fig. 3, middle). At sufficiently strong $`Q`$, both interactions are important, and the two orders appear at a first order transition (Fig. 3, right). The regime of first-order transitions is bounded by a lower, and an upper, tricritical point (big black dots in Fig. 3, left). Above a $`\mathrm{\Delta }`$-dependent threshold value of $`Q/I`$, pure quadrupolar order sets in at a higher critical temperature. From this point on, the quadrupolar splitting is added to the crystal field splitting, eventually suppressing magnetism. Note, however, that for a wide range of $`Q/I`$, magnetism is assisted by quadrupolar interactions, in the sense that the magnetic transition temperature increases with $`Q/I`$. This holds also for part of the regime where the two transitions are distinct.

A remarkable aspect of Pr ordering is its sensitivity to magnetic field, particularly for $`H(001)`$ . Quadrupolar interactions are known to have a strong signature in non-linear magnetic response . We show in Fig. 4 the results for $`I=0`$, $`Q=0.2`$, $`\mathrm{\Delta }=2`$, for the field in the $`ab`$ plane. The upward curvature of $`\chi `$, and the divergence of $`\chi _3`$ at the transition are reminiscent of the observed behavior. Fig. 5 shows the $`T`$-dependence of the magnetization for several fields $`H(100)`$. It is apparent that the sharp $`H=0`$ anomaly is quickly smeared out in higher fields, similar to the observed behavior of PrBa2Cu3O6.

To conclude, we have examined models where the local Hilbert space is spanned by the states from the lowest two crystal field levels. The near-degeneracy of two levels has interesting consequences when the splitting is comparable to both intersite intractions, and laboratory magnetic fields. We found that the induced quadrupolar moment scenario gives a good understanding of the static properties of PrFe4P12, and that the inclusion of quadrupolar interaction is helpful in understanding the non-linear magnetic behavior of PrBa2Cu3O6.

We have been supported by the Hungarian grants OTKA T 038162 and TS 040878. Support from a JSPS–HAS exchange program enabled P.F. to participate at ASR2002.

## References


format: text

text:
# AGNS AMONG ROSAT BLANK FIELD SOURCES

## 1 ROSAT blank field sources

We have found a class of ‘blank field sources’ or ‘blanks’ (Cagnoni et al. 2002). These are X-ray sources with $`f_X/f_V`$$`>10`$ and so are extremely X-ray loud. We selected our sample from the bright ($`f_X>10^{13}`$ erg cm-2 s-1) high latitude ($`|b|>20`$) ROSAT sources in WGACAT (1995 version, White Giommi & Angelini 1984) with no optical optical counterpart on the POSS (McMahon & Irwin 1992). While our sample contains only 16 blanks many more may exist, given our demanding search criteria (e.g. no optical counterpart within 39”, 99% radius). Blanks are a varied population. We have previously shown (Cagnoni et al. 2001) that this is an efficient method of finding high redshift ($`z>0.6`$) clusters of galaxies. Here we highlight the peculiar AGNs among the blanks. (We will use the term AGN in the text meaning unbeamed AGNs, i.e. excluding the blazar class).

## 2 AGNs expected among ROSAT blanks

Some “normal” BL Lacs are expected to contaminate our sample for $`f_X/f_V<35`$ (e.g. Maccacaro et al. 1988), but blanks reach $`f_X/f_V`$ up to a factor of 10 higher. An obvious possibility is that we are selecting the most extreme representatives of the BL Lacertae class, i.e. extremely variable (a factor of $`10`$) and/or BL Lacs with SED peaking at energies higher than the ROSAT band. This peak could be either the low energy (synchrotron, and some of these objects, which also proved to be extremely variable, were recently found by Costamante et al. 2001) or the high energy (inverse Compton, no such BL Lacs have been found so far) one shifted appropriately.
Normal type 1 AGNs with $`f_X>10^{13}`$ erg cm-2 s-1 have $`O18`$ and could be selected as blanks only if their optical counterpart is obscured in the O-band by $`A_O>3.5`$ mag, which corresponds to a local $`N_H>4\times 10^{21}`$ cm-2, or if their optical emission is intrinsically highly suppressed. Not all obscured AGNs could still be visible with $`f_X>10^{13}`$ erg cm-2 s-1 in the ROSAT PSPC. AGNs can reach such X-ray fluxes despite the obscuring material only if they are intrinsically bright, as quasars or QSOs, or if their spectrum is dominated by a conspicuous soft X-ray excess, as Narrow Line Seyfert galaxies (NLSy). QSO-2s and NLSy2 galaxies are the “missing links” of the unifying models of AGN (Urry & Padovani 1995) and their number density and importance for the X-ray background are still debated (e.g. Comastri et al. 2001). The other possibility for AGNs to be selected as blanks is when their optical emission, which is dominated by the big blue bump (BBB) component from the accretion disk (Elvis et al. 1994), is highly suppressed. This can happen in AGNs powered by radiatively inefficient mechanisms for accretion, e.g. advection dominated accretion flows (ADAFs, e.g. Narayan & Yi, 1994), where most of the energy is stored in the gas and advected toward the hole and only a small fraction is radiated.

## 3 AGNs and AGN candidates found among blanks

The main observational properties of the AGN and AGN-candidates in our sample are summarized in Table 1. While we were progressing our investigation one of the blanks, 1WGA J1340.1$`+`$2743, was classified as a BL Lac for the lack of features in its optical spectrum (Lamer, Brunner & Staubert (1997). The source is one of the most extreme in our sample and both its $`f_X/f_V`$ ($`250`$) and its broad-band radio indices ($`\alpha _{ro}=0.62`$ and $`\alpha _{ox}=0.61`$ for a NVSS counterpart of 4.6 mJy) are not compatible with normal BL Lacs. 1WGA J1340.1$`+`$2743 properties can be explained if this source is an extreme synchrotron dominated BL Lac such as those in Costamante et al. (2001) or if the source was flaring during the WGACAT observation. Chandra follow-up caught the source 10 times fainter than in WGACAT, suggesting that the source could indeed be an extreme and/or extremely variable BL Lac.
Two blanks (1WGA J1412.3$`+`$4355 and 1WGA J1415.2$`+`$4403) were spectroscopically identified in the RIXOS survey (Mason et al. 2000) as AGNs at $`z0.50.6`$. 1WGA J1415.2$`+`$4403 could be a normal type 1 AGN because its WGACAT flux was $`50`$% overestimated due to the steep nature of it spectrum. 1WGA J1412.3$`+`$4355 presents an unexpected red POSS counterpart while looks only moderately red in our optical-IR follow-up. Another 6 sources (see Tab. 1) are likely to be AGNs. 1WGA J0951.4$`+`$3916 shows signs of absorption in the X-ray band ($`N_H2_1^{+3}\times 10^{21}`$ cm-2) and has a red counterpart (POSS and our imaging). One more source, 1WGA J1420.0$`+`$0625, differs from the AGN-like ones for the lack of a red counterpart. Could we be seeing the reflected component of an intrinsically bright obscured AGN?

### 3.1 Perverse AGNs: Dusty, yet Unobscured in X-rays

The remaining 5 AGN-like sources seem to share the same peculiar combination of characteristics: type 1 AGNs X-ray spectra with low absorbing column densities (consistent with the Galactic values) but red counterparts ($`OE>2`$ in the POSS and $`RK>4`$ in our imaging)111With the exception of 1WGA J1535.0$`+`$2336, red in the POSS, but only moderately red in our imaging; of 1WGA J1220.6$`+`$3347 red in our imaging but not in the POSS. For a local AGN type 1, $`OE>2`$ implies $`A_V>4`$ (e.g. Risaliti et al. 2001) and thus strong obscuration ($`N_H>2\times 10^{21}`$ cm-2 for a dust to gas ratio typical of the ISM in our Galaxy), in contrast to the results of the X-ray data.
There are some possibilities to reconcile the optical and X-ray data: (1) the AGNs have $`z3.5`$ (e.g. Fig. 7 of Risaliti et al. 2001). This is excluded for 1WGA 1415.2$`+`$4403, but it is still possible for the other RIXOS AGN (1WGA J1412.3$`+`$4355) whose the redshift measurement is based on only 1 optical line; (2) a dust to gas ratio $`4060`$ times higher than the Galactic value; (3) a modified dust grain size distribution with respect to the Galactic ISM; (4) the presence of dust within a warm absorber (e.g. Komossa & Bade 1998); (5) a highly suppressed BBB emission (see § 2).
The planned optical spectroscopic follow-up is needed to investigate the nature of the AGN-like blanks.

## Acknowledgments

This work was supported by NASA ADP grant NAG5-9206 and by the Italian MIUR (IC and AC). IC acknowledges a CNAA fellowship.

## References


text:
# Conductance of single-atom platinum contacts: Voltage-dependence of the conductance histogram

thanks: Present address: Groupe de Physique des Solides, Campus Jussieu tour 23, 2 Place Jussieu, 75251 Paris cedex 05 Francethanks: Present address: Dept. of Physics, 510 Clark Hall, Cornell University, Ithaca, NY 14853

## Abstract

The conductance of a single-atom contact is sensitive to the coupling of this contact atom to the atoms in the leads. Notably for the transition metals this gives rise to a considerable spread in the observed conductance values. The mean conductance value and spread can be obtained from the first peak in conductance histograms recorded from a large set of contact-breaking cycles. In contrast to the monovalent metals, this mean value for Pt depends strongly on the applied voltage bias and other experimental conditions and values ranging from about 1 $`G_0`$ to 2.5 $`G_0`$ ($`G_0=2e^2/h`$) have been reported. We find that at low bias the first peak in the conductance histogram is centered around 1.5 $`G_0`$. However, as the bias increases past 300 mV the peak shifts to 1.8 $`G_0`$. Here we show that this bias dependence is due to a geometric effect where monatomic chains are replaced by single-atom contacts, where the former are destabilized by the electron current at high bias.

preprint: 7

The conductance of atomic-sized contacts (ASCs) is ballistic and can be written in terms of eigenchannels of the contact. In the limit of a single-atom contact the number of eigenchannels is governed by the number of valence orbitals of the atom. scheer98  The simplest description, with a single conductance channel, applies for monovalent metals (Au, Na, …), while for $`sp`$ metals three channels contribute, and for transition metals with open $`d`$ shells five channels contribute to the conductance. In general the transmission for each of the channels is smaller than 1, but one often finds that the single channel for the monovalent metals is nearly perfectly transmitted giving a conductance of almost one conductance quantum, $`G_0=2e^2/h`$, where $`e`$ is the electron charge and $`h`$ is Planck’s constant. The five channels for the transition metals add up to a total conductance that sensitively depends on the coupling of the atom to its neighboring atoms in the leads and usually ranges between 1.5 $`G_0`$ and 3.5 $`G_0`$. This can be judged from the position and width of the first peak in so-called conductance histograms. For a recent review we refer to Ref. Agrait02, .

The concept of conductance histograms was introduced to investigate possible conductance quantization in metallic contacts.Olesen95a ; Krans95b  The histograms are constructed from digitized traces of thousands of cycles of breaking (or making) ASCs. These traces are projected onto the conductance axis, and in this way preferred conductance values become visible as peaks. Only for monovalent metals, including the noble metalsBrandby95a ; Ludoph00a  and the alkali metals,Krans95b ; Ludoph00a  some form of quantization is observed as evidenced by prominent peaks in the histograms close to multiples of the conductance quantum. For most other metals there is no indication of quantization in the conductance histograms.Ludoph00a ; Agrait02  There is usually a single peak at low conductance in the histograms that indicates the preferred conductance of the single-atom contacts.Ludoph00a  For most metals this first peak, when measuring in a clean environment, is fairly reproducible and insensitive to the level of the applied bias voltage. For Au, which has been most studied, it has been shown that the conductance of the first peak is bias independent up to about 2 V, which can be seen in the statistical average represented by conductance histograms measured for various bias voltagesSakaiAu  and from individual curent-voltage ($`I`$-$`V`$) curves.Hansen00a ; Nielsen02a

Pt is a marked exception, and widely different results have been presented for this metal. The reported low bias position of the first peak varies between 1.0 $`G_0`$ and 2.5 $`G_0`$.Krans93a ; Olesen95a ; Sirvent96b ; Yuki00a ; Smit01a ; Smit02b  Bias dependences of Pt conductance histograms have also been reported,Yuki00a  where a low bias peak at 1.0 $`G_0`$ was replaced by a peak at 1.7 $`G_0`$ at higher biases. Recently it has been shown that at low bias the first peak is centered around 1.5 $`G_0`$ (as in Fig. 1), and that a peak at 1.0 $`G_0`$ can be caused by the presence of hydrogen molecules that act as the final bridge before the contact breaks.Smit02b

In this paper we investigate the bias dependence of Pt conductance histograms, compiled from thousands of ASCs formed using the mechanically controllable break-junction (MCBJ)Krans93a ; Agrait02  at liquid helium temperature (4.2 K) in a cryogenic vacuum environment. The conductance measured is the linear conductance $`G=I/V`$. A small but significant shift in the position of the first histogram peak from 1.5 $`G_0`$ to 1.8 $`G_0`$ is observed when the bias increases past 300 mV. Contrary to what would have been expected from measured $`I`$-$`V`$ curves of Pt,Nielsen02a  where the conductance decreases with voltage, the shift is towards higher conductance. This indicates that the shift is not caused by an electronic effect. We present evidence that the shift marks a geometric transition point where single-atom Pt contacts replace the monatomic chains known to exist for Pt.Smit01a

The bias dependence is illustrated by the conductance histograms shown in Fig. 1. For the curve recorded at a low bias of 10 mV the first peak is found close to 1.5 $`G_0`$.

The conductance histogram recorded at 100 mV is very similar to the one at 10 mV with a broad peak centered around 1.5 $`G_0`$. Increasing the bias further to 200 mV only causes a slight decrease in the intensity, whereas the position remains unchanged. However, when the bias increases to 300 mV, the main part of the histogram peak switches to be centered around 1.8 $`G_0`$ although a broad shoulder can still be seen around 1.5 $`G_0`$. At 400 mV, the shoulder has disappeared and the peak center has completely moved from 1.5 $`G_0`$ to 1.8 $`G_0`$. Also the peak intensity decreases considerably for the two high bias conductance histograms compared with those measured at lower bias. We have verified that the process is reversible, by lowering the bias back to zero observing the original histogram with a peak close to 1.5 $`G_0`$ reemerge.

The peak shift to higher $`G`$ is very puzzling when compared to the Pt $`I`$-$`V`$ curves presented previously by Nielsen et al. Nielsen02a  It was shown that $`I`$-$`V`$ curves of Pt are markedly non-linear but with a decreasing conductance as the voltage increases. The decrease in $`G`$ was found to be proportional to the voltage squared, resulting only in a limited decrease of the linear conductance at 300 mV which amounts to about $`0.1`$$`G_0`$ at 500 mV. The important observation is, however, that the conductance always decreases and never increases with voltage, the opposite effect of the one observed with the conductance histograms presented in Fig. 1. This clearly indicates that the peak shift is not caused by an electronic effect.

Investigating other possible causes we turned to the formation of monatomic chains, which was first discovered for Au.Ohnishi98a ; Yanson98a  Recently it was shown by Smit et al. that also Pt and Ir have this property.Smit01a  To demonstrate the formation of monatomic chains, Smit et al. recorded plateau-length histograms. It is important to emphasize the difference between plateau-length histogramsYanson98a ; Smit01a  and the conductance histograms used above. The peaks in the conductance histograms reveal the most probable conductance values obtained for breaking ASCs. A plateau-length histogram is constructed by measuring the lengths of the last conductance plateau, that corresponds to a conductor of a single atom in cross-section, from thousands of conductance traces of breaking ASCs. By plotting the number of times a given length occurs, a histogram is obtained with peaks revealing the typical disruption lengths, at which the monatomic chains tend to break. The first peak reflects the length of the single-atom contact and peaks for longer wires indicate the formation of monatomic chains that break at 2, 3, 4,… atoms in length.

In Fig. 2 we show two plateau-length histograms recorded at biases of (a) 200 mV and (b) 400 mV, respectively.

At 200 mV bias, the plateau-length histogram is very similar to the one presented by Smit et al. at 10 mV bias.Smit01a  Three clear peaks are visible, the first corresponding to a single-atom contact (a ‘one atom chain’), and the next peaks represent chains with two and three atoms, respectively. Thus, we conclude that monatomic chains still form at a bias as high as 200 mV. The three peaks in Fig. 2(a) are centered at 2.7, 4.8 and 7.0 Å, respectively. The distances between them are thus 2.1–2.2 Å, similar to the 2.3 Å found by Smit et al. at 10 mV bias,Smit01a  within the accuracy of $``$10% in our length calibration. The position of the first peak in Fig. 2(a) differs from this peak-to-peak distance as can be expected. It reflects the elastic stretching of the banks and the bonds to a single bridging atom at the verge of breaking.

The presence of only one single peak in the 400 mV plateau-length histogram of Fig. 2(b), centered at 1.9 Å, proves that chains no longer form at this high bias. The peak does not coincide with the corresponding first peak in the 200 mV plateau-length histogram. This difference is most likely due to the higher bias causing the single-atom contact to break at a smaller strain than at low bias (see below).

From the plateau-length histograms we find that the formation of atomic chains is inhibited above the bias voltage for which also the shift in the first peak in the conductance histogram is observed (300 mV). This suggests that the chain formation affects the peak position in the conductance histograms. To test this we recorded conductance histograms for curves measured while returning to contact from the vacuum tunneling regime, which we will refer to as return histograms. The conductance histograms presented previously have all been obtained from the conductance traces of breaking ASCs. When breaking the contacts, monatomic chain formation may occur. Instead, when we obtain the ASCs by returning the electrodes back into contact, chains cannot form. Measuring the conductance traces while forming ASCs, results in the return histograms presented in Fig. 3. From these return histograms it is clear that the first peak is centered close to 2 $`G_0`$, independent of bias. It thus seems that the peak shift occurs when the chain formation is inhibited, since the peak is located at an even higher conductance, independent of bias, when no chains can form.

For Au the conductance of a single-atom contact is indistinguishable from that of monatomic chains and we had expected a similar result for Pt, but the present results suggest otherwise. To get further insight into these findings, we have used the Transiesta programTransSiesta  to calculate the conductance and eigenchannel transmissions of single-atom Au and Pt contacts \[inset of Fig. 4(a)\] as the atom-electrode distance increases. The results are shown in Fig. 4. The method is based on density functional theory and takes the voltage bias and current explicitly into account in the self-consistent calculation of electronic density and potential. Our calculations show that, contrary to Au for which the conductance remains stable around 1 $`G_0`$, single-atom Pt contacts display a strong variation in the zero-bias conductance with a decrease from 2.1 $`G_0`$ to 1.1 $`G_0`$ as the distance increases from 2.65 to 3.5 Å. This variation is in accordance with the reported broad histogram peak centered around (1.5–2) $`G_0`$ for Pt.Sirvent96b ; Smit01a ; Smit02b  These results can help explain the behavior leading to the three types of histograms presented above.

The histogram with a peak at 2 $`G_0`$ (Fig. 3) is produced by single-atom contacts. In this case the electrodes are moved towards each other such that the atom-electrode distance is reduced to a minimum. The peak at 1.8 $`G_0`$ in the high-bias histograms (Fig. 1) also results from single-atom contacts, but in this case the breaking of the wire leads to a larger atom-electrode distance, and thus a lower conductance (Fig. 4). The peak at 1.5 $`G_0`$ in the low-bias histograms of Fig. 1 is due to the frequent occurrence of atomic chains, which are only stable for sufficiently low bias voltages, below 300 mV. The lower average conductance for atomic chains is suggested by the physics described in Fig. 4. As the single atom is positioned farther away from the banks the overlap of its orbitals with those in the electrodes is reduced leading to the observed decrease in the conductance. For monatomic chains the number of states overlapping with the central atoms in the chains is also reduced with respect to the single-atom case, leading to a similar reduction of the conductance. The suppression of chain formation at higher bias is partly due current-induced embrittlement, but mainly results from heating of the atomic degrees of freedom by the current Todorov01 ; Smit03 .

We can now attempt to classify the various results found in the literature. A first peak in the conductance histogram close to 1 $`G_0`$ was found in room temperature experiments under ultra-high vacuum,Olesen95a  or under an atmosphere of N2+5%H2 gas.Yuki00a  For the latter experiment the peak was seen to move to 1.7 $`G_0`$ at elevated bias. For all experiments performed under cryogenic vacuum and at a voltage bias below 100 mV (Refs. Krans93a, ; Sirvent96b, ; Smit01a, ; Smit02b,  and this work) the first peak is found at (1.5–1.6) $`G_0`$. However, when hydrogen is intentionally introduced into the vacuum system a peak near 1 $`G_0`$ appears, which has been interpreted as being due to the formation of a conducting hydrogen bridge forming the last contact.Smit02b  Assuming this mechanism is still effective at room temperature the results by Yuki et al.Yuki00a  find a natural explanation. Trace amounts of hydrogen in the UHV system at room temperature may similarly explain the observation by Olesen et al.Olesen95a

The remaining variation in the position of the first peak in the Pt conductance histograms can be attributed to a bias-dependence of the formation of monatomic chains. The shift of the first peak in the Pt conductance histogram from 1.5 $`G_0`$ to 1.8 $`G_0`$ when the bias increases past 300 mV, marks a geometric transition point where monatomic chains are replaced by single-atom contacts.NielsenPHD

We acknowledge financial support from The Center for Atomic-scale Materials Physics (CAMP) sponsored by the Danish National Research Foundation and from the EU network “Bottom up Nanomachines” (BUN). MB acknowledges support from the Danish Natural Science Research Council. YN has been supported by a Marie Curie fellowship of the European Community under contract number HPMF-CT-1999-00196.


text:
# 1 (A) HBS tiles and their decompositions to Penrose tiles. (b) Atomic decorations for AlCoCu. In (b), only TM and symmetry breaking Al atoms are shown. For AlCoNi both edge sites are decorated with Ni atoms.

## 1. Acknowledgments

The authors thank C.L. Henley and S. Naidu for useful discussions. IA wishes to thank King Abdul Aziz University (Saudi Arabia) for supporting his study and MW acknowledges support by the National Science Foundation under grant DMR-0111198. We thank the Pittsburgh Supercomputer Center for computer time used for this study.


text:
# Detection of TeV Gamma Rays from the BL Lac Object 1ES1959+650 with the Whipple 10m Telescope

## 1 Introduction

Blazars, consisting of flat spectrum radio quasars (FSRQs) and BL Lacertae objects (BL Lacs), are a relatively rare type of active galaxy characterized by extremely variable, non-thermal spectral energy distributions (SEDs). The emission extends from radio to $`\gamma `$-ray frequencies and is believed to be produced by a highly relativistic plasma jet closely orientated with the line of sight to the galaxy (Urry and Padovani, 1995; Blandford and Rees, 1978). In a $`\nu F_\nu `$ representation, the SEDs display two peaks. The lower energy peak is usually attributed to synchrotron radiation from relativistic electrons, while the higher energy peak is thought to be due to inverse Compton scattering by these same electrons. The photons involved in the scattering process may be either synchrotron photons produced in the jet, in the synchrotron self-Compton (SSC) model, or some external population, such as photons emitted from an accretion disc (Dermer, Schlickeiser and Mastichiadis, 1992) or reflected from emission line clouds (Sikora, Begelman and Rees, 1994). Alternative explanations for the high energy component invoke pair cascades triggered via pion and pair photoproduction from high energy protons in the jet (Mannheim, 1993, 1998), or proton synchrotron radiation (Aharonian, 2000; Mücke and Protheroe, 2001).

The EGRET detector on board the CGRO detected in excess of 66 blazars (Hartman et al., 1999) above $`100\mathrm{MeV}`$; however, only a few have been detected at higher energies ($`>300\mathrm{GeV}`$) by ground based atmospheric Cherenkov telescopes. Markarian 421 (Punch et al., 1992) and Markarian 501 (Quinn et al., 1996) were the first extragalactic TeV sources to be detected and the behaviour of these two nearby ($`z=0.031`$ and $`z=0.034`$ respectively) objects has been closely monitored and studied. Flux variability over two orders of magnitude has been observed, with doubling timescales as short as $`15\mathrm{minutes}`$ (Gaidos et al., 1996), placing strong constraints on the size of the emission region. The $`\gamma `$-ray flux is correlated with the X-ray flux (Buckley et al., 1996; Catanese et al., 1997a; Krawczynski et al., 1997) and, in the case of Markarian 421, with the spectral power law index in the $`\gamma `$-ray region (Krennrich et al., 2002; Aharonian et al., 2002a). Unconfirmed detections of 1ES 2344+514 (Catanese et al., 1998) and PKS 2155-304 (Chadwick et al., 1999) have also been reported and, most recently, the detection of very weak emission from H1426+428 has been confirmed by three groups (Horan et al., 2002; Aharonian et al., 2002b; Djannati-Atai et al., 2002). This last object is of particular interest since it is relatively distant ($`z=0.129`$) and hence the observed flux is expected to be strongly attenuated by absorption via pair production on the extragalactic infra-red background light (EBL) (Gould and Schrèder, 1967; Stecker, De Jager and Salamon, 1992).

All of the TeV blazars can be classified as high frequency peak BL Lacs (HBLs) (Padovani and Giommi, 1995) on the basis of the location of the lower energy peak in their SEDs, and predictions of TeV emission from other nearby sources of this type have driven observations by atmospheric Cherenkov detectors. 1ES1959+650 ($`z=0.048`$) was first suggested as a good TeV candidate of this type by Stecker, De Jager and Salamon (1996), who used simple scaling arguments to predict it as the third strongest TeV blazar, after Markarian 421 and Markarian 501, with a flux prediction of $`1.9\times 10^{11}\mathrm{cm}^2\mathrm{s}^1`$ above $`300\mathrm{GeV}`$. More recently Costamante and Ghisellini (2002) have predicted fluxes above $`300\mathrm{GeV}`$ of $`7.5\times 10^{11}\mathrm{cm}^2\mathrm{s}^1`$, based on a simple phenomenological parameterization of the SED adapted from Fossati et al. (1998), and $`0.03\times 10^{11}\mathrm{cm}^2\mathrm{s}^1`$ using a homogeneous, one-zone synchrotron self-Compton model (Ghisellini, Celotti and Costamante, 2002).

In fact, 1ES1959+650 is not a particularly extreme example of an HBL. Measurements with BeppoSAX in 1997 (Beckmann et al., 2002) placed the lower energy SED peak at $`10^{15}\mathrm{Hz}`$ ($`4\mathrm{eV}`$), as compared to typically $`1\mathrm{keV}`$ for Markarian 421 (Maraschi et al., 1999) and $`100\mathrm{keV}`$ for Markarian 501 during its 1997 flare state (Pian et al., 1998), although for both of these sources the position of the peak is known to vary widely depending upon the flux level. Giebels et al. (2002) report on X-ray observations of 1ES1959+650 obtained with the USA(1-16$`\mathrm{keV}`$) and RXTE (2-16$`\mathrm{keV}`$) missions during 2000, showing threefold increases in the X-ray flux on a timescale of a few days. The flux increase is correlated with spectral hardening, indicative of a shift of the lower energy peak of the SED towards higher frequencies. The source is also unusual among HBLs for its strong rapid optical variability. Observations by Villata et al. (2000) showed rapid flickering, including a decrease of 0.28 magnitudes in four days.

Prior to May 2002, only tentative evidence had been presented for TeV emission from 1ES1959+650. A detection with a statistical significance of $`3.9\sigma `$ was reported by Nishiyama et al. (1999) based on $`57\mathrm{hours}`$ of observations with the Utah Seven Telescope Array. More recently, Konopelko et al. (2002) reported a preliminary detection at $`5\sigma `$ for the HEGRA Cherenkov telescope system. Previous observations of 1ES1959+650 with the Whipple 10m telescope produced an upper limit at a flux level of $`1.3\times 10^{11}\mathrm{cm}^2\mathrm{s}^1`$ above $`350\mathrm{GeV}`$ (Catanese et al., 1997b).

Observations of 1ES1959+650 during May-July 2002 with the Whipple 10m telescope resulted in a clear detection of TeV $`\gamma `$-ray emission from this object (Dowdall, Moriarty and Kosack, 2002) which was rapidly confirmed by the HEGRA experiment (Horns and Konopelko, 2002). We report here the results of the Whipple 10m observations.

## 2 Observations and Data Analysis

The configuration of the Whipple 10m $`\gamma `$-ray telescope during these observations is described in detail in Finley et al. (2001). Briefly, the telescope consists of a $`10\mathrm{m}`$ reflector and a 490 pixel photomultiplier tube (PMT) camera. For the purposes of this analysis only the high resolution ($`0.12^{}`$ spacing) central 379 PMT pixels have been used. Cherenkov images are recorded and parameterized according to Hillas (1985) and then $`\gamma `$-ray like images are selected using the “supercuts” criteria (Reynolds et al., 1993) optimized on recent data from the Crab Nebula, the standard candle of TeV $`\gamma `$-ray astronomy (Weekes et al., 1989). Following a realignment of the optical system in February 2002, observations of the Crab Nebula showed that the telescope was at its most sensitive since the installation of the current camera, with one hour of Crab observations producing a $`6\sigma `$ detection; however, between February and July 2002 a decrease in the relative efficiency of the telescope of $`30\%`$ was measured by examining the response to the cosmic ray background. The cause of the effect is still under investigation - one explanation may be that it is due to increased atmospheric absorption caused by large forest fires in the region during this period.

Observations of 1ES1959+650 were made during moonless periods between May $`16^{\mathrm{th}}`$ and July $`8^{\mathrm{th}}`$ 2002. The total dataset consists of $`39.3\mathrm{hours}`$ of on-source data, together with $`7.6\mathrm{hours}`$ of off-source data for background comparison. For observations from Mt. Hopkins (latitude $`31^{}57.6^{}`$N) 1ES1959+650 culminates at a zenith angle ($`\theta _z`$) of $`33.5^{}`$ and so the data were necessarily taken at large $`\theta _z`$, between $`53.5^{}`$ and $`33.5^{}`$. After accounting for the zenith angle and reduced telescope efficiency, we estimate the energy threshold ($`E_{thresh}`$) for the majority of the observations to be $`600\mathrm{GeV}`$ . $`E_{thresh}`$ is defined here as the energy of the peak differential $`\gamma `$-ray flux for a source with the same spectrum as the Crab Nebula. Figure 1 shows the distribution of the Cherenkov image orientation angle alpha for the full on-source dataset, together with the distribution for the off-source dataset scaled such that the number of events in the region $`\alpha >30^{}`$ is the same. The average rate of the excess in the on-source $`\gamma `$-ray region at $`\alpha <15^{}`$ is $`1.08\pm 0.05\gamma \mathrm{min}^1`$, corresponding to a detection of greater than $`20\sigma `$.

## 3 Flux Variability

Figure 2 shows the daily averaged rates for 1ES1959+650 for all of the observations. The rates are expressed in multiples of the steady Crab Nebula flux and have been corrected in order to account for the varying zenith angle and changes in telescope efficiency according to the procedure of LeBohec and Holder (2002). This has been tested using Crab Nebula data taken over a wide range of zenith angles and atmospheric conditions. The correction factor varies from run to run and has a mean value for this dataset of $`2.0\pm 20\%`$. Strong night-to-night variability is evident; the largest change in rate, between MJD 52428 and MJD 52429, corresponds to a doubling time of $`7\mathrm{hours}`$ \- shorter than has ever been observed in other wavebands for this source. The mean flux over all observations was $`0.64\pm 0.03`$ Crab units.

Figure 3 shows the rate in $`5\mathrm{minute}`$ bins for two nights, May $`17^{\mathrm{th}}`$ and June $`4^{\mathrm{th}}`$, during which the source was most active. The statistical evidence for variability within these nights is given by the $`\chi ^2`$ probabilities of constant emission = 1% and 8% respectively. We conclude that there is no strong evidence for flux variability on this timescale.

## 4 Discussion

The detection of TeV $`\gamma `$-ray emission from 1ES1959+650 adds another member to the class of TeV blazars, all of which are close BL Lac objects having a low bolometric luminosity and the peaks in their SEDs at high frequency. The rapid flux variability, and the fact that 1ES1959+650 has not been detected during previous observations, indicates that the source was in an unusual flaring state during these observations. The flux level was at times orders of magnitude above the most recent model predictions. Throughout the period of the Whipple observations measurements in the $`212\mathrm{keV}`$ region by the All-Sky Monitor (ASM) on board RXTE have shown 1ES1959+650 to be active and variable, with daily average fluxes reaching $`20\mathrm{mCrab}`$ in May and July. Target of opportunity observations with the RXTE small field of view instruments were triggered following the $`\gamma `$-ray detection and will be reported on elsewhere (Krawczynski, private communication). The rapid flux variability observed at TeV energies implies a small emission region in the jet with a high Doppler factor (Mattox et al., 1993; Madejski, 1996; Buckley et al., 1996). The contemporaneous X-ray and TeV $`\gamma `$-ray data will allow us to constrain the jet parameters when modelling the emission processes.

Analysis of the Whipple data is ongoing, but attempts to reconstruct the source spectrum have been hampered by the effects of the decreased telescope efficiency during the observation period. The HEGRA collaboration measure a rather steep spectrum (spectral index $`\alpha =3.2\pm 0.3`$) for observations prior to 2002, while the spectrum during the flaring period exhibits pronounced curvature and deviates significantly from the spectrum seen during the quiescent state (Horns et al., 2002). The majority of models for the EBL lead to predictions of a cut-off in the $`\gamma `$-ray region beginning below $`10\mathrm{TeV}`$ for a source at $`z=0.048`$ (Primack (2001); but see Vassiliev (2000) for a discussion of an EBL model which does not produce a distinct feature in the observed spectrum). Deviations from a pure power law have now been resolved in the spectra of both Markarian 421 (Krennrich et al., 2001; Piron et al., 2001; Aharonian et al., 2002a) and Markarian 501 (Samuelson et al., 1998; Aharonian et al., 1999; Djannati-Atai et al., 1999) and the spectrum of the most distant TeV blazar, H1426+428, is measured to be very steep ($`\alpha =3.50\pm 0.35(\mathrm{stat})\pm 0.05(\mathrm{syst})`$) (Petry et al., 2002), but it is not yet clear whether these features are due to absorption on the EBL or are intrinsic to the sources. Clearly, further observations and spectral analysis of 1ES1959+650 may help to resolve this question.

The VERITAS Collaboration is supported by the U.S. Dept. of Energy, N.S.F., the Smithsonian Institution, PPARC (U.K.) and Enterprise Ireland.


text:
# Performance Characteristics of TreePM codes

## 1 Introduction

Cosmological N-Body simulations have played a crucial role in improving our understanding of formation of large scale structure. These have filled large gaps in a domain where analytical solutions do not exist. N-Body simulations have also played a useful role by testing scaling relations derived from physically motivated ansatze.

Limitations of computing resources and our ability to simulate physical processes numerically have meant that N-Body simulations give only approximate solutions. It is possible to test whether the approximations in evolution of the system affects estimation of physical quantities of interest so the results are generally more reliable than those based on approximate evolution of the system.

A large number of methods have been used for doing Cosmological N-Body simulations (Bertschinger, 1998). All of these use some approximation for calculation of force. Approximations are used because direct summation over all pairs of particles scales as $`O(N^2)`$, where $`N`$ is the number of particles, making the calculation very time consuming for large $`N`$. The use of these approximations reduces the number of calculations required to $`O(N\mathrm{ln}N)`$ or less. These approximations also introduce inaccuracies in the computed force.

The TreePM code (Bagla, 2002), that we study here, is a hybrid technique for carrying out large N-Body simulations to study formation and evolution of large scale structures in the universe. It is a combination of the Barnes and Hut (1986) Tree code and a Particle-Mesh code (Hockney and Eastwood, 1988). The TreePM method combines high resolution of tree codes with the ability of PM codes to compute the long range force with periodic boundary conditions.

In this paper we carry out a comprehensive study of the TreePM code. We analyse errors in estimation of force in both the tree and the PM components, and also study the distribution of errors for different distributions of particles. We also study the variation of CPU time required for the TreePM code as we vary parameters describing the mathematical model of this method.

## 2 TreePM method

A large number of methods have been suggested for improving the PM method by combining it with other methods of computing force at small separations (Efstathiou et al, 1985; Couchman, 1991; Xu, 1995) to reduce the softening length below the grid length cutoff imposed by the PM method. A common failing of nearly all such methods is that these continue to use the usual PM force, which is known to have large errors and anisotropies at scales comparable to the grid scale (Bouchet and Kandrup, 1985). These errors are present in the final force, and there is no natural way of reducing the errors in any of these methods.

The philosophy of TreePM method is to modify the PM force in order to have a better control over errors in the long range force. This is done by an explicit division of the potential and force into a long range and a short range component. PM method is then used to compute only the long range component, and the tree method is used to calculate the short range component. This can be expressed in terms of equations as:

$`\phi _𝐤`$ $`=`$ $`\phi _𝐤^s+\phi _𝐤^l`$ (1)
$`\phi _𝐤^l`$ $`=`$ $`\phi _𝐤F(kr_s,\mathrm{})`$ (2)
$`\phi _𝐤^s`$ $`=`$ $`\phi _𝐤\left(1F(kr_s,\mathrm{})\right)`$ (3)
Here, $`\phi _𝐤`$ is the Fourier transform of the gravitational potential, $`\phi _𝐤^s`$ and $`\phi _𝐤^l`$ are the short range and the long range components, respectively. $`F(kr_s,\mathrm{})`$ is a filter function that is unity for $`kr_s1`$ but decreases rapidly for $`kr_s>1`$. Thus $`r_s`$ is the scale at which the short range and long range forces are split. $`F`$ can depend on other parameters apart from $`r_s`$.

This manner of splitting force permits us to control errors in both the components of force. What we require is a function $`F(kr_s,\mathrm{})`$ that decreases sufficiently quickly at $`kr_s>1`$. Similarly in real space, the short range potential and force should become negligible beyond a few $`r_s`$, and the long range force should be negligible below $`r_s`$. Bagla (2002) suggested $`F=\mathrm{exp}\left[k^2r_s^2\right]`$. This function has an added feature that it is positive definite, as is its Fourier transform. Relevant equations for this case are:

$`\phi _𝐤^l`$ $`=`$ $`{\displaystyle \frac{4\pi G\varrho _𝐤}{k^2}}\mathrm{exp}\left[k^2r_s^2\right]`$ (4)
$`\phi ^l(r)`$ $`=`$ $`{\displaystyle \frac{Gm}{r}}\mathrm{erf}\left({\displaystyle \frac{r}{2r_s}}\right)`$ (5)
$`\phi ^s(r)`$ $`=`$ $`{\displaystyle \frac{Gm}{r}}\mathrm{erfc}\left({\displaystyle \frac{r}{2r_s}}\right)`$ (6)
$`𝐟^l(𝐫)`$ $`=`$ $`{\displaystyle \frac{Gm𝐫}{r^3}}\left(\mathrm{erf}\left({\displaystyle \frac{r}{2r_s}}\right){\displaystyle \frac{r}{r_s\sqrt{\pi }}}\mathrm{exp}\left({\displaystyle \frac{r^2}{4r_s^2}}\right)\right)`$ (7)
$`𝐟^s(𝐫)`$ $`=`$ $`{\displaystyle \frac{Gm𝐫}{r^3}}\left(\mathrm{erfc}\left({\displaystyle \frac{r}{2r_s}}\right)+{\displaystyle \frac{r}{r_s\sqrt{\pi }}}\mathrm{exp}\left({\displaystyle \frac{r^2}{4r_s^2}}\right)\right)`$ (8)
Here, $`\mathrm{erf}`$ is the error function and $`\mathrm{erfc}`$ is the complementary error function. Most of our study is restricted to the TreePM codes that use this filter.

We can generalise to more complex filter functions. We studied two families of functions:

$`F(kr_s,\alpha )`$ $`=`$ $`\mathrm{exp}\left[\left(k^2r_s^2\right)^{\alpha /2}\right]`$ (9)
$`F(kr_s,n)`$ $`=`$ $`\left(1+k^2r_s^2\right)^n`$ (10)
Here $`\alpha `$ and $`n`$ are additional parameters. For $`\alpha =2`$ in eqn.(9), we recover the one parameter filter used in Bagla (2002). Fig.1 shows the ratios $`|f^s/f^{tot}|`$ and $`|f^l/f^{tot}|`$ as functions of scale $`r`$ for these models. We have chosen $`r_s=1`$ for this figure. Left panel shows these ratios for models described by eqn.(9) and the right panel shows models described by eqn.(10). In the left panel, we see that the short range force for $`\alpha <2`$ decreases very slowly at large scales, making it an unsuitable choice. For $`\alpha >2`$, the short range force oscillates at large scales and the amplitude of peaks decreases very slowly with scale. Hence these models too are not useful as an alternative to the $`\alpha =2`$ model. Models described by eqn.(10) are shown in the right panel. In this case the behaviour is better than models with $`\alpha 2`$ in that the long and the short range forces both fall off much faster. However, the model with $`\alpha =2`$ is better than any of these models. This model is plotted in the right panel as the solid curve. We would like the overlap between the short and long range forces to be as small as possible. For one, a rapidly decreasing short range force implies that we need to sum this over a small region in space. Computation of the short range force is the most time consuming part and hence a sharply falling short range force is required for a well optimised code. The second reason which requires us to have a long range force that decreases rapidly below some scale is that the PM force typically has large errors at around grid scale. The main motivation of considering an explicit division of force in a long range and a short range component was to avoid carrying over the errors in the PM force by truncating the long range force at scales where PM contributes most significantly to errors. Thus we want to use that filter for which the range of scales at which both the long range and the short range components are relevant is the smallest, and the Gaussian filter is the optimum one in this regard and no other filter comes close to this. We shall not discuss any other filter from this point onwards as figure 1 suggests that the Gaussian filter is the most optimum one.

## 3 Error Analysis

We shall first study errors in the force field of one particle. As the tree approximation does not have any errors for a single particle, we shall focus exclusively on errors in the long range force.

### 3.1 Errors in the Long Range Force

The long range force (eqn.(7)) is computed by the Particle-Mesh method. We summarise the PM method here, and for details we refer the reader to other sources (Hockney and Eastwood, 1988; Bouchet and Kandrup, 1985; Bagla and Padmanabhan, 1997). Particle positions are used to compute the density field on a regular mesh.111We will use terms mesh and grid interchangeably in the context of PM codes. This is done using an interpolating function $`W(x,x_i)`$ that assigns mass of a particle at position $`x`$ to mesh points $`x_i`$. In three dimensions, the interpolating function is a product of three one dimensional interpolating functions. Density field on the mesh is Fourier transformed using Fast Fourier transform (FFT). In Fourier space, we solve eqn.(4) instead of the usual Poisson equation. We multiply this with $`\iota 𝐤`$ to get the negative gradient of the long range force in Fourier space. An inverse transform using FFT gives us the long range force on mesh points. This is interpolated back to particle positions using the same function $`W`$. In this scheme, mesh and the interpolating function $`W`$ are the main sources of anisotropy. We use the cloud-in-cell (CIC) interpolating function (Hockney and Eastwood, 1988).

$`W(x,x_i)`$ $`=`$ $`1|xx_i|/L|xx_i|L`$ (11)
$`=`$ $`0|xx_i|>L`$ (12)
Here $`L`$ is the separation between adjacent mesh points.

To estimate errors in the long range force, we follow the method described by Efstathiou et al (1985). We place a particle at a random position in a mesh cell and use our code to find the long range force. The force is evaluated at a large number of random points scattered within some distance from the particle. This process is repeated a number of times with a different position of the particle. We use the force calculated this way to compute the average long range force as a function of distance as well as the dispersion about this average.

Figure 2 shows a plot of the average long range force (solid curve) and eqn.(7) (dashed curve) as functions of scale. This curve was drawn for $`r_s=1`$. The average force clearly falls below the expected values at small scales.

This difference is caused mainly by the interpolating function. In the continuum limit we require the interpolating function to be like the Dirac delta function. However, we cannot use that and hence there is a serious shortfall in the force at small scales. This becomes clear in Fourier space analysis (Bouchet and Kandrup, 1985; Hockney and Eastwood, 1988). Same analysis also suggests that by de-convolving the interpolating function, we should be able to recover the expected long range force. We deconvolve222We divide the potential in k-space by the square of the Fourier transform of the interpolating function. Since the three dimensional interpolating function is product of the three one dimensional interpolating functions, we are essentially dividing $`\phi _k`$ by $`W_{k_x}^2W_{k_y}^2W_{k_z}^2`$. Square of each of these needs to be taken as we need to do the deconvolution two times. the interpolating function twice, as we use it twice in force calculation: once for calculating the density on the mesh and then for calculating force at the positions of particles. For reference, we mention that the Fourier transform of the one dimensional CIC function is $`4\mathrm{sin}^2(kL/2)/(kL)^2`$, where $`k`$ is the wave number and $`L`$ is the size of one mesh cell.

Such deconvolution of an arbitrary force will not work very well as it will only enhance anisotropies introduced by the mesh. However in this case contribution of wave modes near the Nyquist frequency is negligible due to the Gaussian filter. So we are essentially trying to correct for the amplitude at intermediate wave numbers. Partial correction can be obtained by de-convolving for the effect at small wave numbers by considering the Taylor expansion of $`W_k`$ about $`k=0`$ (White, 2000).

The effect of deconvolution on the long range force is shown in figure 3. We have plotted the difference of the average long range force and eqn.(7) divided by the long range force. We have also plotted the difference divided by the total force because it is this quantity that is relevant for errors in the total force. We have plotted these for the average force shown in figure 2 as well as the average force obtained after deconvolution described above. The effect of deconvolution is striking—the difference between the average long range force and eqn.(7) normalised by the total force drops by nearly one order of magnitude. Fractional error drops below $`1\%`$ and is below $`0.4\%`$ for the entire range of scales. One can lower this further by using a larger $`r_s`$, or use a lower $`r_s`$ if larger errors are acceptable. Figure 3 also shows the maximum of the difference between the average long range force and eqn.(7) normalised by the inverse square force as a function of $`r_s`$. As expected, this decreases as we increase $`r_s`$. Thus we can choose $`r_s`$ once we have fixed the error that is acceptable for our physical application. We will use long range force with deconvolution of the interpolating function $`W`$ in the following discussion.

From deviations of the average force from eqn.(7), we now turn to the dispersion in force about this average. Both these quantities are calculated in narrow bins in distance $`r`$. Figure 4 shows the dispersion in the long range force. We have plotted rms dispersion about the average force as a function of $`r`$. We have normalised the dispersion with the long range force (eqn.(7)), and the total force. The latter curve is relevant for total force. From this curve we can see that the rms dispersion about the average long range force reaches a peak of about $`1.4\%`$ near $`r=2r_s`$.

The dispersion, just like the deviation of the average long range force from the theoretical expectation eqn.(7), decreases as we increase $`r_s`$. Figure 4 shows the variation of the maximum of rms dispersion normalised by the inverse square force, as a function of $`r_s`$.

We can summarise our investigations of the errors in the long range force for a particle as follows: Maximum deviation of the long range force from the theoretical expectation is around $`0.4\%`$ and this deviation is much smaller at most scales. Dispersion about the average peaks at a value of $`1.4\%`$ and is well below $`1\%`$ for most scales. These values are for $`r_s=L`$ and the maxima of errors occur around $`2r_s`$, errors fall off rapidly on both sides. Errors decrease as we increase $`r_s`$.

We will not study errors in the long range force for a system of particles. We will, instead, study the error in the total force after studying the errors in the short range force.

## 4 Errors in the Short Range Force

The short range force is calculated using the tree approximation. We use the Barnes and Hut (1986) method. In this method, the simulation volume is taken to be a cube. Particles are arranged in a hierarchy of cells in an oct-tree, i.e., each cell can have up to eight daughter entities. These can either be cells or particles. The simulation volume is divided into smaller cells so that at the smallest level, no cell contains more than one particle. Once the tree structure is set up, we calculate the position of the centre of mass for each cell as well as the total mass contained within it.

Tree approximation is then used to calculate force on each particle. In this approximation, distant cells are treated as point masses with all the mass concentrated at the centre of mass. Error in force introduced by this approximation is proportional to the square of the ratio of size of the cell ($`d`$) to the distance ($`r`$) to the centre of mass from the point where force is being calculated. We define this ratio as $`\theta =d/r`$. The error in force is proportional to $`\theta ^2`$. In order to control errors, we would not like to use cells that are too close or too large. This is done by introducing a cell acceptance criterion, e.g. by computing the quantity $`\theta `$ and comparing it with a threshold $`\theta _c`$:

$$\theta =\frac{d}{r}\theta _c$$
(13)
The error in force will obviously increase with $`\theta _c`$. It can be shown that the maximum error for inverse square force is $`9\theta _c^2/4`$. For a detailed discussion of cell acceptance criteria and error in force, we refer the reader to several detailed studies (Hernquist, 1987; Salmon and Warren, 1994; Springel, Yoshida and White, 2001).

The tree approximation is used here for computing the short range force which differs from the usual inverse square force. We expect errors to be different in this case. We calculate the error at leading order in tree approximation and compare it with the corresponding error in the case of inverse square force.

We take a cubical cell of dimension $`d`$ and distribute point masses randomly inside it. Our aim is to calculate the error involved in computing the short range force in the tree approximation on a unit point mass situated at the origin due to the masses in the cubical cell described above. To calculate the error, we need to compare the short range force in the tree approximation ($`𝐟_{cm}^s`$) with the sum of short range forces due to all the particles in this cell ($`𝐟^s`$). We also have the inverse square force in the tree approximation $`𝐟_{cm}`$.

To quantify error in short range force we define $`ϵ^s`$ as a measure of the error in the computed short range force.

$$ϵ_s=\frac{\left|𝐟_{cm}^s𝐟^s\right|}{|𝐟_{cm}|}$$
(14)
We compute the error in case of the inverse square force for comparison. The worst case we can have is if two particles are at diagonally opposite ends of the cell, and the origin is along the same diagonal. The distance between the particles is $`d\sqrt{3}`$, and the fractional error in force upto leading order in $`\theta _c`$ is $`9d^2/(4r^2)=9\theta _c^2/4`$ where $`r`$ is the distance from the origin to the centre of mass. If the two particles are separated along an edge of the cell instead of the diagonal, the error is $`3\theta _c^2/4`$. A similar analysis for the short range force gives

$`ϵ_s`$ $`=`$ $`{\displaystyle \frac{3\theta _c^2}{4}}\left[\mathrm{erfc}\left({\displaystyle \frac{r}{2r_s}}\right)+\left({\displaystyle \frac{r}{r_s}}\right)\mathrm{exp}\left({\displaystyle \frac{r^2}{4r_s^2}}\right)\right]`$ (15)
$`+{\displaystyle \frac{\theta _c^2}{4}}\left[{\displaystyle \frac{2}{\sqrt{\pi }}}\left({\displaystyle \frac{r}{r_s}}\right)^3\mathrm{exp}\left({\displaystyle \frac{r^2}{4r_s^2}}\right)+{\displaystyle \frac{1}{8\sqrt{\pi }}}\left({\displaystyle \frac{r}{r_s}}\right)^5\mathrm{exp}\left({\displaystyle \frac{r^2}{4r_s^2}}\right)\right].`$
Here we have also assumed that the separation between the two particles is smaller than $`r_s`$. Error here has been calculated for the case when the two particles are separated along an edge of the cell and the origin (where the force is being calculated) and these two points are collinear. We have plotted the error as a function of distance $`r`$ in fig.5.

Now we turn to a more realistic estimate of error for random distribution of particles in the cell, instead of considering pathological situations which give rise to large error. Right panel of figure 5 shows the fractional error in short range force $`ϵ_s`$ as a function of distance $`r`$. Error is constant at small scales, increases around the scales where the contribution of short range force is decreasing rapidly, then plummets to very small values at larger scales. The range of scales where the error is large are those where the short range force contributes a significant fraction to the total force and its variation with distance is significantly different from inverse square. In the figure, error is plotted for two values of $`\theta _c`$. It is lower for smaller $`\theta _c`$: error for $`\theta _c=0.3`$ is a nearly a factor three smaller than for $`\theta _c=0.5`$. Here error has been averaged over $`100`$ different distributions of $`30`$ particles in the cell. Error is smaller for larger number of particles as chances of an extreme distribution of particles become smaller. Even for $`N_p=10`$, the error is considerably smaller than the extreme case we used for estimating errors above. Error for $`\theta _c=0.5`$ and $`N_p=30`$ is about $`2\%`$, slightly larger than the maximum error due to long range force ($`1.4\%`$ for $`r_s=L`$). Errors for $`\theta _c=0.3`$ are smaller and closer to the errors in the long range force. However, we should add that these are errors due to particles in a single cell. In a generic situation, these errors may cancel out or add up, depending upon the details of the particle distribution.

## 5 Errors in the TreePM force

We turn to the question of errors in the TreePM force calculated by adding the short range and the long range forces. Errors are calculated with respect to a reference force computed with the following configuration of the TreePM code. We chose $`\theta _c=0.01`$, hence the errors in short range force are negligible as these errors decrease in proportion to $`\theta _c^2`$ for the worst case scenario. We take $`r_s=4.0`$ thus error in the long range force is below $`0.2\%`$ if we extrapolate the maximum dispersion in long range force. With this setup, we can safely assume that the error in reference force (as compared to the inverse square force with periodic boundary conditions that we should use) is below $`0.2\%`$ at all scales, and this is sufficient for studying errors in the TreePM force calculated using more pragmatic values of $`r_s`$ and $`\theta _c`$ as we are aiming to keep errors below $`1\%`$ level for most particles in all situations.

We will calculate errors for two very different distributions of particles, one a random distribution of particles with uniform density (unclustered distribution) and the other a clustered distribution generated by an N-Body simulation. We will present the results for these distributions side by side.

We present the results of our analysis by plotting the cumulative distribution for fractional error. Errors were calculated for a distribution of $`2^{21}2\times 10^6`$ particles so that we have sufficient number of particles to study the distribution.

We have presented the variation of errors with $`r_{cut}`$ in fig. 6. For this figure we used $`r_s=L`$ and $`\theta _c=0.5`$. $`r_{cut}`$ is the radius within which we sum the short range force around each particle. We need to sum out to a distance where the contribution of the short range force to the inverse square force drops to an ignorable fraction. From fig. 1, we know that the short range force is less than $`1\%`$ of the inverse square force at $`5r_s`$, so we expect that $`r_{cut}`$ should be comparable to this value. We find that for the unclustered distribution, distribution of errors is very different for $`r_{cut}=3r_s`$ as compared to $`r_{cut}4r_s`$. For all $`r_{cut}4r_s`$, distribution of errors is very similar. This can result if particles at scales beyond $`4r_s`$ are not contributing to the force. This is to be expected in a uniform distribution of particles. Distribution of errors for the clustered distribution is shown in the right panel of fig. 6. For the clustered distribution, the particle distribution is not homogeneous and hence particles at larger distances contribute significantly to the force on each particle. As this contribution increases, force from $`r>4r_s`$ becomes more relevant. Error distribution for the clustered distribution shows small differences between $`r_{cut}=4r_s`$ and $`r_{cut}5r_s`$. From these figures we conclude that $`r_{cut}=5r_s`$ is a safe choice and that $`r_{cut}=4r_s`$ will suffice for most purposes.

Figure 7 shows the distribution of errors for different values of $`\theta _c`$. As before, we have shown results for the unclustered distribution in the left panel and that for the clustered distribution in the right panel. We used $`r_s=L`$ and $`r_{cut}=6r_s`$ for these plots. There is a marked difference in errors for different values of $`\theta _c`$ in case of unclustered distribution. Error decreases with $`\theta _c`$, the variation decreases as we get to smaller values of $`\theta _c`$. Situation is somewhat different for the clustered distribution as there is hardly any difference between errors for different values of $`\theta _c`$. This suggests that the errors are dominated by errors in the long range force for clustered distribution. For unclustered distribution, the total force on each particle is small, whereas force due to a cell with many particles is large and many large contributions cancel out to give a small net force. Numerical errors in adding and subtracting large numbers seem to dominate errors in case of the uniform distribution. Contribution of cells is large for larger $`\theta _c`$, hence we get a significant variation with $`\theta _c`$ for unclustered distribution.

Lastly, we look at variation of errors with $`r_s`$. Figure 8 shows the distribution of errors for difference values of $`r_s`$. We fixed the values of other parameters to $`r_{cut}=6r_s`$ and $`\theta _c=0.5`$. Left panel is for the unclustered distribution and the right panel is for clustered distribution of particles.

Errors for the clustered distribution increase monotonically as $`r_s`$ is decreased. This is to be expected from our analysis of errors in the long range force. The fractional error in the TreePM force for $`99\%`$ of the particles for $`r_s=L`$ is less than $`0.8\%`$. Variation of errors with $`r_s`$ for the unclustered distribution is not monotonic. Errors decrease till $`r_s`$ increases to unity, and then increase again. This happens because the small force on each particle for a uniform distribution again comes into play. For larger $`r_s`$, the force due to individual cells is large and numerical error in adding all the large contributions to get the small net force is also large. The fractional error in the TreePM force for $`99\%`$ of the particles for $`r_s=L`$ is less than $`1.5\%`$, thus the variation in error from unclustered to clustered situation is small compared to that in tree codes. From this figure, we conclude that $`r_s=L`$ is the optimum choice for the TreePM code.

We can summarise the results of this study of errors as follows: Errors for $`99\%`$ of particles are below $`1.5\%`$ for unclustered distribution and below $`0.8\%`$ for clustered distribution for $`r_s=L`$, $`r_{cut}5r_s`$ and $`\theta _c0.5`$. The best way to reduce errors is to use a larger $`r_s`$ as varying the other two parameters leads to large differences between the clustered and the unclustered cases.

## 6 CPU time requirements

In this section we will list the CPU time required for one time step of the TreePM code for a simulation of $`2^{11}`$ particles. We will also study variation of this as we vary the parameters of the TreePM method.

The computer on which these timings were obtained is powered by a $`1.6`$GHz Pentium 4 CPU. Programs were compiled using the Intel Fortran compiler (version 5.0) and double precision numbers were used throughout. Programs are written in Fortran 90/95.

The recommended configuration of the TreePM code uses $`r_s=L`$, $`\theta _c=0.5`$ and $`r_{cut}=5r_s`$. In this case one time step takes $`740`$ seconds. Times taken for building the tree, calculating the long range force, calculating the short range force and moving particles are, for one time step, equal to 6.38 seconds, 26.70 seconds, 711.51 seconds and 0.25 seconds respectively. This implies that the time taken for each particle per step is $`0.35`$ms. The quoted numbers are for an unclustered uniform distribution of $`128^3`$ particles.

Time required for one time step of the TreePM code varies with the choice of parameters $`r_{cut}`$, $`r_s`$ and $`\theta _c`$. We expect the time taken to increase with increasing $`r_s`$ and $`r_{cut}`$, and decrease with increasing $`\theta _c`$. The behaviour for these cases is shown in fig. 9. It is seen that the time taken falls sharply with increasing $`\theta _c`$, whereas errors are insensitive to the choice as long as $`\theta _c0.5`$. Therefore it makes sense to use $`\theta _c=0.5`$. Time taken rises sharply with increasing $`r_{cut}`$. Again, by comparing with the variation of errors, we recommend $`r_{cut}4.5r_s`$. Lastly, time taken increases rapidly as we increase $`r_s`$. We recommend $`r_s=L`$ as the error is minimum for this choice for the unclustered distribution.

## 7 Integrating the equation of motion

Our discussion so far has dealt only with evaluation of force. This is the main focus of this paper as the key difference between the TreePM and other methods is in the scheme used for evaluation of force. However for the sake of completeness, we give here details of integration of equations of motion. The usual form of the equation of motion is:

$`\ddot{𝐱}`$ $`+`$ $`2{\displaystyle \frac{\dot{a}}{a}}\dot{𝐱}={\displaystyle \frac{1}{a^2}}\varphi `$
$`^2\varphi `$ $`=`$ $`4\pi Ga^2\left(\rho \overline{\rho }\right)`$ (16)
Here $`𝐱`$ is the comoving coordinate, $`a`$ is the scale factor, $`\rho `$ is the density and $`\overline{\rho }`$ is the average density of the universe. Dot represents differentiation with respect to time. We can recast these equations in the following form:

$`𝐱^{\prime \prime }`$ $`+`$ $`{\displaystyle \frac{2}{a}}\left\{1+\left({\displaystyle \frac{1}{4}}\right){\displaystyle \frac{2\mathrm{\Omega }_\mathrm{\Lambda }a^3\mathrm{\Omega }_0}{\mathrm{\Omega }_\mathrm{\Lambda }a^3+\mathrm{\Omega }_0+a\left(1\mathrm{\Omega }_0\mathrm{\Omega }_\mathrm{\Lambda }\right)}}\right\}𝐱^{}`$
$`=`$ $`\left({\displaystyle \frac{3\mathrm{\Omega }_0}{2a^2}}\right){\displaystyle \frac{1}{\mathrm{\Omega }_\mathrm{\Lambda }a^3+\mathrm{\Omega }_0+a\left(1\mathrm{\Omega }_0\mathrm{\Omega }_\mathrm{\Lambda }\right)}}\psi `$
$`^2\psi `$ $`=`$ $`\delta ={\displaystyle \frac{\rho }{\overline{\rho }}}1`$ (17)
Here prime denotes differentiation with respect to the scale factor. $`\mathrm{\Omega }_0`$ is the density parameter of non-relativistic matter, $`\mathrm{\Omega }_\mathrm{\Lambda }`$ is the density parameter of the Cosmological constant and $`\delta `$ is the density contrast. We use a modified gravitational potential $`\psi `$. The equation of motion contains a velocity dependent term and hence we cannot use the usual leap-frog method. We recast the leap-frog method so that velocities and positions are defined at the same instant (Hut, Makino and McMillan, 1995). We solve the equation for velocity iteratively. Time step is chosen to be a small fraction of the smallest dynamical time in the system at any given stage. The fraction to be chosen is fixed by doing some tests.

One of the tests we did was to check for scale invariance in evolution of power law spectra. A simulation is repeated with different choices of time step until we find the largest time step for which we can reach the highly non-linear regime ($`\overline{\xi }>100`$) and retain scale invariance. We then use a time step that is half of this largest time step. Figure 10 shows $`\overline{\xi }`$ as a function of $`r/r_{nl}(t)`$ for several epochs. Here $`r_{nl}(t)`$ is the scale which is going non-linear at time $`t`$ and it varies in proportion with $`a^{2/\left(n+3\right)}`$ in the Einstein-deSitter model. $`n`$ is the index of the power spectrum. We have used $`n=1`$ for the left panel and $`n=1`$ for the right panel. We have only plotted $`\overline{\xi }`$ at scales more than four times larger than the artificial softening length used in the simulation. Scale invariance holds for both the spectra and we are able to probe the highly non-linear regime with this code.

## 8 Discussion

In this paper we have presented a detailed study of performance characteristics of the TreePM method. We have analysed different sources of error and suggested remedies for the main source of errors. The analysis of errors in realistic situations shows that the TreePM method performs very well and gives acceptably small errors.

This method compares favourably with other comparable methods such as implementations of the tree code like GADGET (Springel, Yoshida and White, 2001) and hybrid methods such as the TPM (Xu, 1995; Bode, Ostriker and Xu, 2002). From the numbers available in these papers, we find that the errors in the recommended configuration of the TreePM method are comparable with those in GADGET and lower than those in the TPM method. In terms of CPU time taken per step per particle, we again find that the numbers are comparable. Of course, it is not possible to make a detailed comparison of this quantity as the whole approach is different. E.g., we do not use multiple time steps whereas GADGET relies heavily on these to optimise the speed. On the other hand GADGET and TreePM have a uniform force resolution whereas TPM does not and hence the time taken is likely to vary more strongly with the amplitude of clustering for the TPM code as compared to the other two.

Splitting of force into a short range and a long range part has a useful implication in terms of parallel implementation. Unlike the tree methods, the TreePM method involves a considerably smaller overhead in terms of inter-process communication and hence is likely to score over other methods in terms of scaling for a very large number of CPUs. In summary we can state that TreePM is a competitive method for doing Cosmological N-Body simulations. Explicit use of three parameters gives users control over errors and CPU time required.

## Acknowledgements

A part of the work reported here was done using the Beowulf at the Harish-Chandra Research Institute (http://cluster.mri.ernet.in).


text:
# The Enigma of AA Dor

## References

Dorman B., Nelson L.A., Chau W.Y. 1989, ApJ 342, 1003

Kilkenny D., Harrop-Allin M., Marang F. 1991, Inf. Bull. on Variable Stars, No. 3569

Hilditch R.W., Harries T.J., Hill G. 1996, MNRAS 279, 1380

Iben I.Jr., Livio M. 1993, PASP 105, 1373

Rauch T. 2000, A&A 356, 665

Rossiter, R.A., 1924, ApJ 60, 15