text:
# 1 Introduction

## 1 Introduction

Recent progress in both experimental techniques and theoretical modeling of epitaxial growth allowed to improve our understanding of growth mechanisms on the atomic scale. In particular, a satisfactory understanding of growth in homoepitaxial systems is beginning to emerge. Many results were accumulated by basic material-oriented research as well as by general studies in the spirit of statistical mechanics. For example, the processes of nucleation and aggregation in single component submonolayer growth  are now well understood.

On the other hand, practically important materials prepared by growth contain several components or impurities. Our present understanding of the real microscopic processes going on in such systems is rather incomplete. The investigation of multicomponent growth is more complicated not only because we have to consider more material parameters than in homoepitaxy, but because there are additional possible atomistic processes which have to be considered (inlayer or interlayer exchange, floating of adsorbates on the surface or the island edges, additional energy barriers due to island or step edge decoration etc.).

Multicomponent growth is controlled by the interactions between species. The interactions are determined by the chemical nature of components, as well as by the geometrical positions of the atoms. There are two basic features of submonolayer multicomponent growth which affect the surface morphology: effects of the strain due to atomic mismatch, and effects of the different chemical bondings between components of the system. Although it is broadly accepted that the presence of strain plays a crucial role in heteroepitaxy, in particular in the formation of ordered nanostructures, the different chemical bondings by themselves can result in interesting effects. As an example we may mention Hwang’s STM study of growth of Co on Ru(0001) with predeposited Ag . He observed the formation of dendritic islands, despite the fact that neither of the two component metals exhibit this behavior. Hwang suggested that this can be explained by chemically induced step edge diffusion barriers between domains of different species, which should generally exist in multicomponent film growth. Another example are the effects of phase separation in multilayer molecular beam epitaxy .

In this paper, we concentrate on effects of the difference in chemical bondings between the atoms of a two-component system growing in the submonolayer regime. We discuss the effects on the resulting island morphology, and address the question of the modification of island density scaling.

We present some general considerations on different growth regimes in two component growth in Section 2. In the next Section 3, we describe a model for submonolayer two-component growth. Section 4 very briefly recalls the scaling law for the island density. In section 5, we discuss the case of a small concentration of one component. Section 6 contains results of kinetic Monte Carlo simulations of two component alloy growth. In particular, we show how the scaling of the island density is modified by the relative mobility of the components and interactions between them.

## 2 Different regimes of multicomponent growth

Growth in the presence of more than one type of particles is realized in various ways. The simplest scenario is the deposition of a certain type of particles on a substrate composed from a different material. This is usually called heteroepitaxy. Interesting effects like growth-induced surface alloying of species which are immiscible in the bulk have been observed in such systems.

There are many experimental studies of growth in which two or more different components were deposited on a substrate consisting in general of a different material (see e.g. ). In the case of multicomponent growth, we need to distinguish different regimes of growth. The species can be either codeposited simultaneously with generally different fluxes, or deposited sequentially, one of the species being predeposited. We restrict ourselves in this paper to the simplest situation of two component growth. We call the components A and B, and we consider the situation where an alloy consisting of A and B , or a material (A) in the presence of impurities (B) grows on a certain substrate.

Important material parameters in submonolayer growth are the strengths of the interactions, both laterally and perpendicular to the surface. In comparison to homoepitaxy, two-component growth contains many more material parameters. However, within a certain temperature window only some of them may be important. The crucial parameter is the interaction between particles of different type. It can be both attractive or repulsive. We consider here the former. Depending on the relative strengh of this interaction with respect to the interactions between particles of the same type, one can observe intermixing (possibly leading to antiferromagnetic-like ordering) or phase separation. The second important parameter is the relative mobility of the different species.

To simplify the classification, we fix one lateral interaction, let us say between A and A, and the mobility of one species, let us say of A. The interaction between A and B and the interaction between B and B can be stronger or weaker than the interaction between A and A, and the mobility of B can be higher or lower than the mobility of A. This yields 12 qualitatively different situations of submonolayer growth in a bicomponent system. (In the case of the multilayer growth, such classification cannot be easily done due to the presence of additional interlayer processes.) It is important to note that the precise values of parameters needed to decide in which regime growth proceeds are often not well known. In this situation, it seems to be useful to explore different possible effects on a general level. We are not going to analyze all the situation here but we discuss some partial results.

## 3 Two-component submonolayer growth model

To carry out an effective analysis we need a reasonably simple, yet flexible model. Kinetic Monte Carlo (KMC) simulations proved to be a useful tool for evaluating the importance of various microscopic processes. We use the variant of the solid-on-solid model with two surface species described in detail in . The basic microscopic processes are random deposition and migration. Two types of particles are deposited with fluxes $`F_A`$ and $`F_B`$. In the case of codeposition the concentration of the two components is changed by the variation of the ratio $`F_A/F_B`$.

Diffusion is modeled by hopping to nearest-neighbor sites with the rate $`R_D^X=k_0\mathrm{exp}(E_D^X/k_BT)`$ for a particle of the type $`X=A`$ or $`B`$. The energy barrier $`E_D^X=E_{\mathrm{sub}}^X+_{Y=A,B}n^{XY}E_\mathrm{n}^{XY}`$ depends on the interaction with the substrate (the first term) and on the lateral neighbors (the second term).

The interaction with the substrate is affected by the particle beneath the particle under consideration. It may be a substrate particle, which can be in general of a type different from A or B (let us call it C), or a previously deposited particle of type A or B. We have in general six possibilities: $`E_{\mathrm{sub}}^{AC}`$, $`E_{\mathrm{sub}}^{BC}`$, $`E_{\mathrm{sub}}^{AA}`$, $`E_{\mathrm{sub}}^{BA}`$ $`E_{\mathrm{sub}}^{AB}`$, $`E_{\mathrm{sub}}^{BB}`$. However, in the situation of a low coverage and fast descent to the lower terrace one can suppose that only $`E_{\mathrm{sub}}^A=E_{\mathrm{sub}}^{AC}`$ and $`E_{\mathrm{sub}}^B=E_{\mathrm{sub}}^{BC}`$ are the relevant two parameters controlling the mobility of both particle types.

The interaction with the nearest-neighbors is controlled by bond counting, i.e. $`n^{XY}`$ is the number of nearest-neighbor $`X`$-$`Y`$ pairs and $`E_\mathrm{n}^{XY}`$ the corresponding contribution to the barrier, with $`E_\mathrm{n}^{AB}=E_\mathrm{n}^{BA}`$. In all examples presented in this paper $`E_\mathrm{n}^{AA}=0.3`$ eV, $`E_{\mathrm{sub}}^A=0.8`$ eV and $`k_0=10^{13}`$ Hz.

## 4 Island density scaling

We call the island density $`N`$. In homoepitaxy, the scaling relation

$$N(F/D)^\chi ,\chi =\frac{i^{}}{i^{}+2}$$
(1)
between the deposition flux $`F`$ and the adatom diffusion coefficient $`D`$ has been well established theoretically \[6–8\] numerically \[8–12\] and experimentally . In (1), $`i^{}`$ denotes the size of the largest unstable cluster.

Little is known about the details of nucleation during submonolayer growth of two component systems. The presence of the second species affects the nucleation of the first one and vice versa. An obvious complication is that there are different types of islands. In principle, one can consider three island densities: $`N_A`$ ($`N_B`$) for islands composed only from A (B), and $`N`$ for islands composed from both particle types. In the case of mixing, there is no problem, since $`N`$ is the relevant quantity. The situation is less clear in the presence of phase separation. Here islands composed only from one type of particles may exist for weak or repulsive A-B interactions.

## 5 Effects of adsorbates

The situation is simpler if the coverage of one of the components is much smaller than the total coverage (let us say $`\theta _B\theta `$, then the B-component is an impurity), or if there is another mechanism like floating of one component (B) on the edges of the islands of the other component (A), since then $`N_AN`$.

Several microscopic mechanisms by which impurities could alter the relationship (1) were considered in the literature. First, impurities may act as nucleation centers, thus effectively decreasing $`i^{}`$ and therefore $`\chi `$; in the extreme case of immobile adatom traps, the limit of spontaneous nucleation with $`i^{}=\chi =0`$ would be realized . Recently, we have shown that the predeposition of only a low concentration of impurities ($`\theta _B`$ a few thousandth of ML) can lead to a severalfold increase of the island density . However, the increase significantly depends on the mobility of the adsorbates. In the case of essentially immobile adsorbates ($`E_{\mathrm{sub}}^B=5`$ eV), a new feature in the flux dependence of the island density appears. Instead of a single power law relationship, there is a plateau where $`N\theta _B`$, reflecting the dominance of heterogeneous nucleation in a certain flux  (temperature ) interval. When its origin is not recognized, the evaluation of such a plateau within the usual analysis of island density scaling  may lead to an anomalously low estimate of the adatom diffusion barrier, as well as to an anomalously small preexponential factor. This effect has been invoked  to explain the extremely low diffusion prefactors reported recently for several metal surfaces .

Second, impurities decorating the island edges may induce energy barriers to attachment. Kandel  predicted using a rate equation theory that, provided these barriers are sufficiently strong, the exponent $`\chi `$ in (1) is increased such that the above expression is replaced by $`\chi =2i^{}/(i^{}+3)`$. Also this mechanisms implies an increase of the island density compared to the case of pure homoepitaxy. We found in KMC simulation that the key kinetic process for growth of decorated islands is the activated exchange of impurities and adatoms . The rate of this process is taken as $`k_{\mathrm{ex}}=k_0\mathrm{exp}(E_{\mathrm{ex}}/k_BT)`$, where $`E_{\mathrm{ex}}`$ is the corresponding activation barrier. In this simulation the barriers for exchange $`E_{\mathrm{ex}}`$ and for the B-mobility $`E_{\mathrm{sub}}^B`$ are fixed at $`E_{\mathrm{ex}}=E_{\mathrm{sub}}^B=1`$ eV. We observed well decorated islands for both $`E_\mathrm{n}^{AB}<E_\mathrm{n}^{AA}`$ and $`E_\mathrm{n}^{AB}>E_\mathrm{n}^{AA}`$ (left panel of Fig. 1). We found that, in our case , the floating mobile adsorbates indeed strongly increase the island density but without appreciably changing its power-law dependence on the flux (right panel of Fig. 1). The increase is stronger for larger $`E_\mathrm{n}^{AB}`$, and only slightly higher for predeposition than for codeposition. We also observed a stronger coverage dependence of the island density in comparison with homoepitaxy (right panel of Fig. 1). This was interpreted as a delay of the saturation regime, where the island density becomes independent of coverage. The coverage dependence is again more pronounced for strongly interacting adsorbates, but it is weaker in the case of predeposition than for codeposition.

Third, it has been argued on the basis of KMC simulations  that *repulsive* impurities can alter the exponent $`\chi `$ in (1) by rendering the diffusion of the adatoms anomalous, in the sense that their mean square displacement behaves as $`x^2t^\alpha `$ with $`\alpha <1`$. A scaling analysis for $`i^{}=1`$ then predicts that $`\chi =\alpha /(1+2\alpha )<1/3`$, i.e. the island density scaling exponent is predicted to decrease in this case.

## 6 Island density in binary alloy

In this section, we study nucleation in a situation where the concentration of both species is comparable. We shall discuss effects of the different mobility of both components and of the variation of the A-B and B-B interactions; as before, the mobility of A-atoms and A-A interaction is fixed: $`E_{\mathrm{sub}}^A=0.8`$ eV, $`E_\mathrm{n}^{AA}=0.3`$ eV . We consider (i) the regime of phase separation, $`E_\mathrm{n}^{AB}=0.2`$ eV and $`E_\mathrm{n}^{BB}=0.4`$ eV, and (ii) the regime of the intermixing, $`E_\mathrm{n}^{AB}=0.4`$ eV and $`E_\mathrm{n}^{BB}=0.2`$ eV. In all simulations the particles are codeposited and the results are analyzed at the total coverage $`\theta =0.1`$ ML.

In Fig. 2, we show the island density $`N`$ as a function of the inverse temperature for different mobilities of the B-component in the regime of phase separation for two concentrations $`F_A=F_B=0.001`$ ML/s (left panel) and $`F_A=0.001`$ ML/s, $`F_B=0.0001`$ ML/s (right panel). The behaviour in homoepitaxy, the case when there is no B-component, is shown for comparison. One can see that there is a substantial increase of the island density. Similar as for the case of immobile impurities, there is the indication of a plateau for the case of the lowest mobility of the B-component for both concentrations.

Hence, the mobility of the atoms of one component strongly influences the island density of the growing alloy even its concentration is small. Note that if the mobility of B is higher than the mobility of A, then the island density is *lower* than in homoepitaxy for $`F_A=F_B=0.001`$ ML/s, but not for $`F_A=0.001`$ ML/s, $`F_B=0.0001`$ ML/s. If one analyzes separately the density $`N_A`$ of islands composed only from A (not shown) no clear picture appears.

The change of the island density with mobility is accompanied by a change of the island morphology. This is illustrated in the left column of Fig. 3 showing configurations at $`T=500`$ K for different mobilities of B-atoms ($`F_A=F_B=0.001`$ ML/s). We can see a transition from a situation in which A-islands are incorporated with comparably large B-islands into

one island, to a situation in which there are A-islands attached to small B-islands. The second column corresponds to the right panel in Fig. 2. The effect of B-mobility on the island morphology (as well as on the island dendity cf. Fig. 2) is now weaker. The right column shows the situation in the intermixing regime.

The left panel of Fig. 4 compares island densities in the regimes of phase separation and intermixing for low B mobility and two B-B interaction strengths. We can see that in the intermixing regime the presence of a small fraction of a less mobile component increases the island density more than in the phase separation regime (cf. the middle and the right columns in Fig. 3). The increase of the strength of B-B interaction has only a small effect, except for high temperatures in the phase separation regime.

In the right panel of Fig. 4 we display the coverage dependence of the island density and the density of free adatoms in the phase separation regime for different mobilities of the B-component. There is a delay of the saturation in the island density and of decrease of the free adatoms density. The effect is stronger for lower B-mobility. It can be explained by the continuous creation of centers for heterogeneous nucleation.

### Acknowledgements.

We thank P. Šmilauer for useful discussion. Support by the COST project P3.130 and by Volkswagenstiftung is gratefully acknowledged.


format: text

text:
# The properties of the X-ray holes in the intracluster medium of the Perseus cluster

## 1 Introduction

The Perseus cluster of galaxies, A 426, at a redshift $`z=0.0183`$ is the brightest cluster X-ray source in the sky. The X-ray surface brightness peaks around the radio source 3C84 in the central galaxy NGC 1275. The inner radio lobes of 3C84 clearly interact with the X-ray emitting intracluster gas, leading to X-ray holes North and South of the nucleus (Böhringer et al 1993, Churazov et al 2000; Fabian et al 2000; Fig. 1). Outer holes are also seen (Fabian et al 2000) which do not correspond to any structures on high frequency ($``$GHz) radio maps. Here we investigate the properties of the inner and outer X-ray holes and show that the outer holes lie at the ends of spurs of low-frequency radio emission seen at 74 MHz.

A comparison of the ROSAT High Resolution Imager data (Böhringer et al 1993) with 1.4 GHz radio data (Pedlar et al 1990) showed that equipartition radio plasma was likely underpressured with respect to the surrounding thermal intracluster medium. With the present Chandra data (Fabian et al 2000) we now have temperature information on the scale of the holes and can estimate the thermal pressure much better. We find that the radio plasma and the intracluster gas can only be in pressure equilibrium if the magnetic field has a pressure much lower than required for equipartition.

The combination of X-ray and radio data now allows improved estimates to be made of the total energy released by the radio source, and of $`k`$, the factor which accounts for the energy carried by the particles in the radio plasma in addition to that of the electrons radiating above a fiducial level of 10 MHz. In Section 2 we present the data and in Section 3 deduce the pressure within the radio lobes from standard equipartition arguments. We find a pressure mismatch with the surrounding thermal gas unless $`k/f`$ is large. Then assuming pressure equilibrium and a synchrotron lifetime argument we show that the field is below equipartition. Using dynamical constraints we then deduce the power needed to supply the lobes. The outer lobes are discussed in Section 4. Our conclusions are reported in Section 5.

## 2 The X-ray and radio data

The Chandra X-ray data have been presented by Fabian et al (2000). Three separate exposures were made of the Perseus cluster; in 1999 for 10 ks with ACIS-I, in 2000 Jan 10 for 10 ks with ACIS-S and in 2000 Jan 30 for 25 ks with ACIS-S. A study of the combined, cleaned data from the last two ACIS-S images, totaling 23.9 ks, will be presented by Schmidt et al (2001, in preparation); here we show an adaptively-smoothed version of this last image as the top panel of Fig. 1 and a Gaussian-smoothed image from the ACIS-I data covering the band 0.6-3 keV as the lower panel. A 2-7 keV image from ACIS-S is shown in Fig. 2. The adaptively-smoothed combined ACIS-S image is also used as the background to the contours in Figs. 3 and 4.

The two inner holes, which we term the N and S lobes, are clearly seen as is at least one outer hole to the NW (Figs. 1 and 5; also seen in Einstein Observatory images; Fabian et al 1981). A new hole is also seen about 1 arcmin S of the nucleus (Fabian et al 2000; see lower panel in Fig. 3). The inner holes have sharp bright rims in soft X-ray maps and contain the X-ray coolest gas in the images ($`kT3\mathrm{keV}`$; Fabian et al 2000) so the rims are not due to shocks. No additional hard X-ray structures are seen in the 2–7 keV image (Fig. 2).

Note that if the rims were due to strong shocks, the observed emissivity, which is proportional to $`n^2T^{0.5}`$ where $`n`$ and $`T`$ are the gas density and temperature, respectively, would change by a factor of $`16(T_{\mathrm{sh}}/T_0)^{0.5}`$. The ratio of the temperature of the shocked to unshocked gas, $`T_{\mathrm{sh}}/T_0`$ is most unlikely to be so high that the shocked gas is not readily detectable.

The present data are inadequate to determine clearly whether the holes are real, i.e. devoid of any thermal gas, or apparent, i.e. partially filled with hot gas at the cluster virial temperature of about 7 keV. At constant pressure the X-ray emissivity varies as the inverse square of the gas temperature, so deep soft X-ray holes can be made by a factor of two difference in temperature.

The three-dimensional geometry of the holes is unclear. A simple picture is that the jets (see Pedlar et al 1990) are close to the plane of the sky and that the holes are spherical. It is however possible that the jets and alignment of the lobes are more along our line of sight; e.g. the S radio jet may contribute to the N hole. Note from Fig. 2 that the S continuation of the NW rim avoids the nucleus. We adopt a spherical geometry for the N hole in our calculations.

The data have been deprojected (Fabian et al 1981) to obtain the gas pressure profile up to the rims around the inner X-ray holes. This confirmed that the simple (electron) pressure map shown in Fabian et al (2000) is accurate to within 20 per cent. In particular the total thermal pressure of the gas in the rim of the N radio lobe is $`0.5\mathrm{keV}\mathrm{cm}^3,`$ reducing by a factor of two at the position of the outer NW lobe.

3C84 was observed on 1998 March 7 simultaneously at 74 MHz and 330 MHz using all 27 antennae of the VLA in its A-configuration. The data were observed in spectral line mode with 32 channels at 74 MHz and 64 channels at 330 MHz: this facilitates RFI excision (which is largely narrow-band) and helps to avoid bandwidth smearing. The on-source time was 64 minutes which was spread throughout a 6 hour period to maximize hour-angle coverage. Bandpass, amplitude and initial phase calibration are derived from several short scans of Cygnus A made through out the observing run. Phase calibration utilized successive loops of self-calibration initiated with point source models at both frequencies. This was possible because 3C84 has enough short-spacing flux for self-calibration to work in its most elementary form and permit good phase solutions on time averages much shorter than the ionospheric incoherence time.

We show the 1.4 GHz (from G. Taylor, see Fabian et al 2000), 330 MHz and 74 MHz radio maps overlaid on an adaptively-smoothed X-ray image in Fig. 3. The spurs of the 74 MHz radio emission running NW and SSE both end at outer holes seen in the X-ray image (see also Fig. 5). The 74 MHz emission likely maps the regions containing the oldest relativistic electrons, plausibly indicating that the outer X-ray holes are old radio lobes, within which the most energetic electrons have lost their energy to synchrotron emission.

## 3 The inner radio lobes

### 3.1 Equipartition between relativistic particles and magnetic field

We now apply standard synchrotron radiation theory to the inner lobes, in order to quantify the properties of the relativistic (particles and field) component. In particular we use region ‘3’ of Pedlar et al (1990), which is the brightest region of the N inner lobe (see top panel in Fig. 3).

The total energy in relativistic electrons radiating from frequency $`\nu _1`$ to $`\nu _2`$ with spectral index $`\alpha `$ ($`S(\nu )\nu ^\alpha `$) in a magnetic field $`B`$, giving flux density $`S_\nu `$ at frequency $`\nu `$, is

$$E_\mathrm{e}=4\pi \times 10^{12}\left(\frac{cz}{H_0}\right)^2\left(1+\frac{z}{2}\right)^2\frac{S_\nu }{\nu ^\alpha }\frac{\nu _{2}^{}{}_{}{}^{0.5+\alpha }\nu _{1}^{}{}_{}{}^{0.5+\alpha }}{\alpha +0.5}B^{3/2}$$

$$aB^{3/2}\mathrm{erg},$$

where the estimate $`a=3.6\times 10^{48}`$ corresponds to the radio data reported by Pedlar et al (1990) at 332 MHz and $`H_0=75`$ km s-1 Mpc-1. It is also assumed that the particle energy distribution extends in the range responsible for the emission between $`\nu _1=10`$ MHz and $`\nu _2=1.4`$ GHz. These has been chosen to comprise the detected frequency interval 111Given the crucial role of $`\nu _1`$ we checked that our results are qualitatively consistent for $`\nu _1`$ within at least the range 1-50 MHz..

The total energy in particles and magnetic field is then

$$E_{\mathrm{tot}}=kE_\mathrm{e}+Vf\frac{B^2}{8\pi }=akB^{3/2}+bfB^2,$$

where $`b=8.9\times 10^{64}`$, $`V=4\pi R^3/3`$ is the volume and $`R=2.7\mathrm{kpc}`$ (Pedlar et al 1990). $`f`$ represents the filling factor of the relativistic plasma, while the factor $`k`$ accounts for the additional energy in relativistic particles accompanying the electrons which radiate above 10 MHz and any non-relativistic component including turbulence ($`k=1`$ for an electron-positron plasma emitting only in the above waveband; a typical value used in the literature is $`k=100`$). This implies an equipartition magnetic field strength

$$B_{\mathrm{eq}}=1.9\times 10^5\left(\frac{k}{f}\right)^{2/7}\mathrm{G},$$

and the corresponding equipartition pressure in particles and field of

$$P_\mathrm{p}+P_\mathrm{B}=\frac{ka}{3VfB^{3/2}}+\frac{B^2}{8\pi }1.3\times 10^2\left(\frac{k}{f}\right)^{4/7}\mathrm{keV}\mathrm{cm}^3.$$

### 3.2 Pressure equilibrium between the relativistic and thermal plasmas

We now consider the relation between the pressure of the relativistic component and that of the X-ray emitting gas. The thermal pressure in the rims immediately surrounding the lobes is $`P_{\mathrm{th}}0.5\mathrm{keV}\mathrm{cm}^3`$, and thus

$$\frac{P_{\mathrm{th}}}{P_\mathrm{p}+P_\mathrm{B}}40\left(\frac{k}{f}\right)^{4/7},$$

implying that in order to be in equilibrium at equipartition

$$\frac{k}{f}6\times 10^2.$$

This result can be interpreted in many ways. Perhaps $`k`$ is large indicating that the electron distribution extends to much lower energies than assumed by the fiducial synchrotron frequency of 10 MHz or the existence of heavy proton-rich jets. Alternatively $`f`$ may be small and the radio plasma occupies only a tiny fraction of the apparent volume of the radio lobes. The plasma may be in sheets and filaments with a large covering fraction which have swept much of the cooler gas to the rims.

Finally, the radio plasma may not be in equipartition, as will indeed be shown in the next Section. We plot in Fig. 6 the pressures against the magnetic field for various values of $`k/f`$. Pressure equilibrium between the thermal, magnetic and particle pressures and equipartition require a strong magnetic field and high value of $`k/f`$. From now on we do not assume that the plasma is in equipartition.

### 3.3 Radiative and dynamical constraints

### 3.4 Radiative losses

We now consider further constraints which can be inferred from the spatial distribution of the radio emission at different frequencies as well as dynamical constraints deriving from the conditions of the X-ray emitting gas. The only assumption we keep is that of equilibrium between the thermal and relativistic gas pressures, which constraints the relation between $`k/f`$ and $`B`$ such that

$$P_\mathrm{p}+P_\mathrm{B}=P_{\mathrm{th}}.$$

In the following we also consider constraints on the whole lobe, which is 7.5 kpc in radius (instead of the 2.7 kpc for region ‘3’), corresponding to an increase in volume by a factor of $``$ 21. The observed flux, and therefore particle energy, increases - using the maps of Pedlar et al (1990) and assuming the same spectral index - by about a factor of 5.5 in going from the modeled 2.7 kpc radius region to the rest of the hole. The overall effect is to change the ratio $`k/f`$ by less than a factor 2. We therefore plot only the allowed values of $`k/f`$ in Fig. 7 for the whole lobe, as this region will be considered in the following.

Under the assumption of pressure equilibrium let us then consider the radio emission. A limit on the maximum magnetic field is deduced by noting that GHz radio emission is seen throughout the N hole, so requiring that the synchrotron cooling time of the relativistic electrons there,

$$t_{\mathrm{sync}}=4\times 10^7B_5^{3/2}\nu _9^{1/2}\mathrm{yr},$$

be more than the age of the hole. A limit on this age is obtained by noting that the rims have not been shocked, so the hole has expanded at a velocity less than the sound speed $`c_\mathrm{s}`$ of the 3 keV gas in the rim, i.e.

$$t_{\mathrm{hole}}>r/c_\mathrm{s}=10^7\mathrm{yr},$$

where a radius of 7.5 kpc is used for the hole. This limits the magnetic field to $`B<2.5\times 10^5\mathrm{G}`$ and consequently rules out any equipartition solution (see Fig. 6). Note that if the electrons diffuse between higher and lower field regions within the lobe then the limit is strengthened. Note also that the above condition implicitly assumes that no reacceleration of electrons occurs in the lobe. This is supported by the distribution of radio spectral index (see Fig. 4) which shows a smooth behaviour through the hole and a gradual steepening of the index from the inner to the outer parts (see also Section 4).

This constraint thus imposes a strong upper limit on $`B`$ which is almost an order of magnitude below that required if $`P_\mathrm{p}=P_\mathrm{B}=P_{\mathrm{th}}.`$ We conclude here that $`P_\mathrm{p}=P_{\mathrm{th}}P_\mathrm{B}`$, and we obtain from Fig. 7 that $`k/f<500`$.

If the field is very weak then cooling of relativistic electrons by synchrotron emission is reduced and adiabatic and inverse Compton losses on the 3K radiation field might be important. In particular the latter becomes more efficient than synchrotron for $`B<10^6\mathrm{G}.`$

### 3.5 Dynamical constraints

We finally examine constraints from the dynamics of the system. $`PdV`$ work is of course done by the expanding bubble on the surrounding gas, with a power $`P_{\mathrm{th}}fV/t`$. Most of this power will however propagate to large radii in the cluster as a sound wave (Fabian et al 2000; Reynolds et al 2001) and deposit little energy locally. In any case the jet must be powerful enough to make the holes. In order to estimate such power we consider the evolution of an expanding bubble in a medium with constant pressure following Churazov et al (2000). We modify their expressions for the expansion of the hole to include the effect of gas clumping (through the filling factor $`f`$).

To make a hole with the observed radius of 7.5 kpc requires

$$\frac{L_{45}t_7}{P_{\mathrm{th}}f}=0.5,$$

when the expansion is already in the subsonic phase. We assume that the adiabatic index for the radio-emitting gas is $`4/3,`$ appropriate for a relativistic plasma, and parameterize the luminosity as $`10^{45}L_{45}\mathrm{erg}\mathrm{s}^1`$ and the age of the hole as $`10^7t_7\mathrm{yr}.`$ Secondly the expansion rate must not be so fast as to shock the rim gas, which means that

$$\frac{L_{45}}{P_{\mathrm{th}}t_7^2f}<14.$$

In Fig. 7 the oblique solid line representing a decrease in $`k/f`$ for increasing $`B`$ corresponds to the condition that this expansion timescale is shorter than the synchrotron cooling one. Finally the buoyancy of the hole must be considered (a detached low density region will rise at close to the local virial velocity; Churazov et al 2000) and the fact that it has not yet detached from the centre means (again modifying expressions from Churazov et al 2000) that

$$L_{45}>1.2f^{5/2}t_7^2.$$

These three expressions define the allowable region in the power vs time plane, $`L,t,`$ shown in Fig. 8.

Although it appears from Fig. 8 that very low powers are possible, work must be done in pushing the soft X-ray emitting gas out of the way to make the hole, even for low filling factors. The minimum luminosity required to move gas to make the hole is plausibly $`P_{\mathrm{th}}V/t`$, which is that needed to do work against gravity. This is shown as a dotted line in Fig. 8. Only if the soft X-ray emitting gas is removed by say mixing in with hotter gas, requiring little power from the relativistic plasma, could this limit be overcome. We assume therefore that the filling factor of the lobes $`f1`$.

Further independent support for the possibility that the relativistic gas may be clumped (i.e. $`f<1`$) comes from the requirement that relativistic electrons do not cool by Coulomb collisions in the thermal plasma. This would happen on timescales $`10^5`$ yr if the two phases were significantly mixed. Magnetic field presumably provides the required confinement and insulation. Highly filamentary structures extending over extremely large linear scales have been indeed detected in radio observations (see the case of Fornax A, Fomalont et al 1989). A similar situation can be envisaged for 3C 84. A direct estimate of the degree of clumping might be determined by the amount of inverse Compton scattered radiation in the X-ray band, as this depends directly on the (relativistic) particle density. However at the sensitivity level currently available such a constraint only translates into $`f^{2/3}B_5^{3/2}>10^3`$.

The total power is high compared with the radiated luminosity in the radio band, which is of order $`10^{40}10^{41}`$ erg s$`^1.`$ Most of the power either remains in the relativistic plasma in the lobes or is expended as $`PdV`$ work on the surroundings. Unless there is an efficient method of dissipating that energy locally, it is likely to be transported by buoyancy and low-frequency sound waves to the outer parts of the cluster.

Fig. 8 gives us an upper limit on the dynamical age of the lobes which, through the comparison with the synchrotron cooling time in Fig. 7, translates into lower limits on both $`k/f`$ and the magnetic field. Together with the constraint on $`f`$ from Fig. 8 (i.e. from the need to make the hole) we now have fairly strong constraints on $`k`$ (Fig. 9). It is restricted to be between 500 and 180, if $`f=1`$.

## 4 The outer lobes

We now consider the information which can be inferred from the multifrequency radio maps of the outer holes.

### 4.1 The NW hole

The dynamical time of the outer NW hole, based again on buoyancy, is $`6\times 10^7\mathrm{yr}`$ if its filling factor by (low-density) radio plasma, $`f`$, is high.

A possible weak spur in the 330 MHz map reaches towards it, but it is one of several spurs. A radio spur clearly points to it only in the 74 MHz map. This is borne out by the steep spectral index (down to $`1.7`$) seen on the S side of that spur (Fig. 4). The synchrotron cooling time of the radio spur is thus about $`6\times 10^7B_5^{3/2}\mathrm{yr}`$ if we assume that a spectral ageing break occurs around 100 MHz.

The good agreement between dynamical and synchrotron timescales if $`B_51`$ indicates that the NW hole is consistent with being an old detached radio lobe within which the magnetic field is in approximate pressure equilibrium with the surroundings. If this field $`B10^5\mathrm{G}`$ also pervades the N lobe, then $`k/f100`$.

Although the NW hole appears to be buoyant, suggesting that it is perhaps empty of thermal gas or at least has a high filling factor in relativistic plasma, a similar picture would emerge if the hole were filled with gas at the virial temperature ($`7\mathrm{keV}`$) and $`f`$ were small.

The edges of the NW hole appear sharp (Fig. 5). Presumably, large scale magnetic fields along the surface provide the tension necessary to maintain integrity of the hole against instabilities. It is possible too that magnetic fields in the thermal gas also restrict mixing, which might otherwise lead to fast cooling of the relativistic phase (note that very sharp boundaries between cooler and hotter gas have been seen in several clusters; e.g. Markevitch et al 2000).

### 4.2 The S hole

The Chandra hole to the S is also at the end of a spur of 74 MHz emission. This time there is also a stronger spur at 330 MHz. The spectral index map shows a steep index in and to the W of this spur. It is again, presumably, another old radio lobe for which similar conclusions to the NW hole apply.

## 5 Conclusions

We have conducted a detailed physical analysis of 3C84 in the Perseus cluster based on new observations of thermal X-ray and nonthermal low frequency radio emission. Our Chandra and 74 MHz radio images provide the highest resolution view of the core of the Perseus cluster yet obtained in either X-rays or at low radio frequencies ($`<100`$ MHz), respectively. The requirements of pressure equilibrium between the relativistic plasma and the hot gas, and of the synchrotron cooling timescale exceeding the dynamical one, restrict the ratio $`180<k/f<500`$ in the inner N hole and exclude the possibility that the relativistic plasma and mean magnetic field are in equipartition. Further indications from the outer lobes support this result. $`k340\pm 160`$ if $`f1`$, or if the filling factor is less than unity, tends toward $`k34\pm 16`$ when $`f0.1`$.

A value of $`k340`$ can be accounted for if the electron spectrum extends down to low energies corresponding to synchrotron emission at about 3 kHz. Otherwise energetic protons must accompany the electrons.

The value of $`f`$ greatly affects the estimate of the jet power: the dynamical constraints on the expansion and buoyancy of the hole require powers which greatly exceed (by up to a factor of 104) the minimum one, unless the filling factor is small. The most energetically efficient situation occurs when the filling factor is small. We estimate that $`f>0.1`$ unless the holes are created in a non-energetic manner.

Much deeper observations with Chandra are required to clearly detect hotter gas in the holes and so constrain the filling factor of the relativistic plasma. Future, higher resolution 74 MHz VLA observations would be useful for further delineating the structure of the low frequency spurs. More desirable, even lower frequency observations must await the proposed Low Frequency Array (LOFAR).

## Acknowledgments

We thank Robert Schmidt for help with the X-ray data. Basic research in radio astronomy at the Naval Research Laboratory is supported by the Office of Naval Research. The National Radio Astronomy Observatory is operated by Associated Universities Inc., under cooperative agreement with the National Science Foundation. AC acknowledges the Italian MURST for financial support. ACF and KMB thank the Royal Society for support.


text:
# The 𝜃-Vacua and the Leutwyler–Smilga Scaling Regime

## 1 Introduction

The vacuum angle, $`\theta `$, a fundamental parameter in QCD, is constrained to be zero to within one part in a trillion. Understanding why this parameter is so fine–tuned is a fascinating theoretical problem. Dashen’s phenomenon , the first–order phase transition at $`\theta =\pi `$ in which the discrete CP symmetry is spontaneously broken, is an example of novel physics at nonzero values of $`\theta `$. Since the physics at nonzero $`\theta `$ is inherently nonperturbative, the only two theoretical avenues to explore this region of the QCD parameter space are lattice QCD and effective theories. The former approach suffers from a sign–problem which is similar to that which plagues attempts to study QCD at nonzero baryonic chemical potential.

The approach employed here is to study the physics of nonzero $`\theta `$ using chiral perturbation theory in finite volume . Restricting the Euclidean four-volume, $`V=L^4`$, to the range $`\frac{1}{\mathrm{\Lambda }}L\frac{1}{m_\pi }`$ where $`\mathrm{\Lambda }`$ is the chiral symmetry breaking scale and $`m_\pi `$ is the mass of the Goldstone excitations, allows for an exact analytical treatment . This is possible since the lower limit implies that the partition function is dominated by Goldstone modes and the upper limit requires that these modes are constant. As a result, the partition function reduces to a finite–dimensional group integration.

## 2 Results

The full partition function for $`N_f=2`$ and $`N_c3`$ with quarks in the fundamental representation was calculated in Ref.  and the necessity of including the contributions from all topological sectors in the partition function was demonstrated in Ref. . For $`N_f3`$, the problem is technically more difficult. It was shown in Ref.  that the full partition function for arbitrary $`N_f`$, $`\theta `$ and quark masses, $`m_i`$, can be expressed as a $`N_f2`$ dimensional integral over single Bessel functions.

For simplicity, we focus on $`N_f=2`$ and $`N_c3`$ with quarks in the fundamental representation but the following results are generalizable . The partition function is $`𝒵(\theta ,\mu _1,\mu _2)=\frac{I_1\left(\mu _{12}\right)}{\mu _{12}}`$ where $`\mu _{12}=\sqrt{\mu _1^2+\mu _2^2+2\mu _1\mu _2\mathrm{cos}(\theta )}`$, $`\mu _i=\mathrm{\Sigma }Vm_i`$ and $`\mathrm{\Sigma }`$ is the infinite-volume chiral condensate at $`\theta =0`$ . Taking the two quark masses to be equal and $`\theta `$ to be an integer multiple of $`\pi `$ leads to the remarkable property that the partition function is independent of the scaling variable with $`𝒵^{(N_f=2)}(\theta =n\pi ,\mu )=\frac{1}{2}`$. For arbitrary $`\theta `$ and degenerate quark masses, the dependence on the volume, macroscopic chiral condensate, quark masses and $`\theta `$ is only through the dimensionless variable $`\mu |\mathrm{cos}(\theta /2)|`$, which as discussed in Ref. , has profound implications. On account of this very constrained dependence of the partition function on the parameters of the theory, the chiral and the topological properties are intimately intertwined. In particular, the chiral condensate and the topological density are proportional up to a simple $`\theta `$-dependent function .

Since the analytic form of the partition function is known, calculating the vacuum properties is straightforward. Taking the limit of degenerate quark masses, the chiral condensate is $`\mathrm{\Sigma }(\theta ,\mu )=\mathrm{\Sigma }|\mathrm{cos}(\theta /2)|I_2\left(2\mu |\mathrm{cos}(\theta /2)|\right)/\left[2I_1\left(2\mu |\mathrm{cos}(\theta /2)|\right)\right]`$. For fixed $`\mu `$, the chiral condensate decreases monotonically in the interval $`\theta [0,\pi )`$ and increases monotonically for $`\theta (\pi ,2\pi )`$. As $`\mu =m\mathrm{\Sigma }V\mathrm{}`$, the chiral condensate develops a cusp at $`\theta =\pi `$. For any value of $`\mu `$, the chiral condensate vanishes identically at $`\theta =\pi `$. In the limit of degenerate quark masses, the topological density, defined as the first derivative of the logarithm of the partition function with respect to $`\theta `$, is $`\sigma (\theta ,\mu )=m\mathrm{tan}(\theta /2)\mathrm{\Sigma }(\theta ,\mu )`$. The most interesting property of this relation is that in the limit of very large scaling variable $`\sigma (\theta ,\mu )`$ develops a discontinuity at $`\theta =\pi `$. This is the first–order phase transition proposed by Dashen . The topological susceptibility can also be calculated with the infinite volume result $`\chi (\theta ,m,V)=\frac{\mathrm{\Sigma }m}{2}|\mathrm{cos}(\theta /2)|`$. This is nothing but the flavor singlet Ward–Takahashi identity generalized to nonzero values of $`\theta `$ which predicts a linear rise in the topological susceptibility with the quark mass. All of the above results can be shown to be independent of both the pattern of chiral symmetry breaking and the representation of the matter fields .


text:
# Geophysical aspects of very long baseline neutrino experiments

## 1 Introduction

The Earth is a tectonically active planet. Large scale thermal convection, which is related to the motion of tectonic plates on the Earth’s surface, is taking place in the Earth’s crust and mantle (which collectively extend from the core-mantle boundary, $`r3480`$ km, to the Earth’s surface, $`r6371`$ km, where $`r`$ is the radius111The Earth is ellipsoidal, due primarily to its rotation; except where otherwise noted, depths are spherically averaged values.). The crust and mantle consist primarily of silicate rocks. The outermost layer is the crust, which has a thickness ranging from about 80 km under Tibet to about 5 km in beneath oceans. As discussed below, the physical properties of the crust are highly laterally heterogeneous. The mantle is divided into the upper mantle, which extends from the base of the crust to a depth of about 410 km ($`r5961`$ km); the transition zone, in the depth range 410 km to 660 km ($`r5711`$ km to $`r5961`$ km); and the lower mantle, in the depth range 660 km to 2891 km ($`r3480`$ km to $`r5710`$ km). The boundaries between the upper mantle and the transition zone, and between the transition zone and the lower mantle, are thought to be due to phase transitions in silicate minerals.

The Earth’s core consists primarily of iron and thus has a considerably greater density than the mantle. The outer core ($`r1222`$ km to $`r3480`$ km) is liquid; magnetohydrodynamic convection in the outer core is considered to be the cause of the Earth’s magnetic field. The inner core, which extends from the base of the outer core to the Earth’s center ($`r=0`$ to $`r1222`$ km), is solid.

For further general information on the structure of the Earth’s interior see recent textbooks (e.g., Lay and Wallace, 1995; Shearer, 1999) and the works cited therein.

Due to the increase of pressure with depth, the Earth’s density and elastic constants are vertically heterogeneous. However, because the Earth is tectonically active, its physical properties are also laterally heterogeneous. Let us denote the laterally averaged one-dimensional (1-D) density structure by $`\rho (r)`$, where $`\rho `$ is the density in units of gm/cm3 (or kg/m3), and denote the three dimensional (3-D) density distribution by $`\rho (r,\theta ,\varphi )`$, where $`\theta `$ and $`\varphi `$ are respectively the colatitude and longitude, in spherical polar coordinates.

The Earth’s average density can be determined from its total mass, $`m_e=5.97\times 10^{24}`$ kg, and its outer radius. If the Earth were a homogeneous sphere, its moment of inertia would be $`0.40Mr_e^2`$, where $`r_e`$ is the Earth’s outer radius. However, the observed moment of inertia has a much smaller value, approximately $`0.33Mr_e^2`$. This confirms that the Earth’s inner regions (i.e. the outer and inner core) are significantly denser than average.

Even if the Earth’s total mass and moment of inertia are combined with other geodetic data such as the spherical harmonic expansion of the Earth’s external gravity field (which is inferred from satellite data), these data provide integral constraints on the Earth’s density distribution but are insufficient to determine it uniquely. It therefore is necessary to use seismological data as the primary basis for inferring the Earth’s density distribution. However, for technical reasons that are not discussed in detail here, inferring the Earth’s density distribution directly from observed seismological data is not practically realizable (Bullen, 1975; Kennett, 1998). It thus is necessary to follow a two step inference process. First the spatial distribution of seismic wave velocities is inferred from seismological data; second, the density distribution is inferred from the seismic velocities, using the above integral constraints together with other empirical relations. Both of these steps introduce uncertainty into the density model.

### 1.1 Seismic velocities

For the purposes of very long baseline neutrino experiments, isotropic Earth models can probably be regarded as sufficiently accurate; the discussion in this paper is limited to such models. The most general anisotropic elastic solid has 21 independent elastic constants, but an isotropic elastic solid has only two independent elastic constants, the Lamé constants $`\lambda `$ and $`\mu `$.

In an isotropic elastic body the velocity of compressional elastic waves (P-waves), $`\alpha `$, and the velocity of transverse elastic waves (S-waves), $`\beta `$, are given respectively by

$`\alpha `$ $`=`$ $`\sqrt{(\lambda +2\mu )/\rho }`$ (1)
$`\beta `$ $`=`$ $`\sqrt{\mu /\rho }.`$ (2)
As a rough approximation, the ratio of P- and S-wave velocities in the solid Earth is given by

$$\alpha 1.7\beta ,$$
(3)
but the exact value of the proportionality constant varies with the chemical composition, pressure and temperature.

## 2 How Earth models are inferred

Inversion of observed seismic data for Earth structure is an underdetermined inverse problem, and all Earth models are subject to error and uncertainty. Regularization constraints of some type (e.g., smoothness, minimum variation from the starting model, etc.) must be applied to obtain a stable solution. Inverse theory allows formal error estimates to be made, but it is well known that systematic errors, which cannot be quantitatively estimated, may often be on the same order or larger. Systematic errors are due to factors such as the uneven distribution of seismic observatories on the Earth’s surface (in particular the lack of observatories on the ocean bottom) and the uneven spatial distribution of earthquakes, and thus cannot be easily reduced. The approximations (e.g., ray-theoretic, linearized perturbation with respect to a spherically symmetric model, etc.) used to model seismic wave propagation are another significant source of systematic errors; progress in forward modeling and inversion techniques is leading to reduction of such errors. Anelastic attenuation (absorption) of seismic waves also places inherent limits on resolving power, especially of deeper and shorter wavelength structure.

A 1-D model seismic velocity specifies $`\alpha (r)`$ and $`\beta (r)`$, while a 3-D model specifies $`\alpha (r,\theta ,\varphi )`$ and $`\beta (r,\theta ,\varphi )`$. A 1-D model may either be a globally averaged model or a model of the depth dependence under some region; similarly, a 3-D model may either be a global model or may be limited to some particular region. The main focus of seismological research on Earth structure has shifted to the quest to infer 3-D Earth models. In this context, the role of 1-D models is to provide the starting point for defining a 3-D model as a perturbation to the 1-D starting model.

The primary data used to obtain seismic velocity models are the arrival times of seismic body waves (P- and S-waves that travel through the Earth’s interior). The arrival time data are then analyzed to determine the location $`(r_0,\theta _0,\varphi _0)`$ and origin time $`t_0`$ of each earthquake and can then be converted to the travel time from the source to the receiver. A large dataset of travel time data for many earthquakes is then inverted to obtain a new Earth model, and the earthquake location process is then updated. This process is iterated several times until convergence is obtained. Travel time data are in some cases supplemented by surface wave dispersion data (the frequency dependence of the phase and group velocities of seismic surface waves) or free oscillation data (the frequencies of several hundred of the longest period modes, which are basically equivalent to surface waves). A recent trend is to use the seismic waveforms themselves (the recorded displacement of the ground as a function of time), rather than secondary data such as the travel times, as the data in the inversion.

Improvements in data and in inversion methodology over the past 20 years have led to steady improvement in seismic velocity models. Two well known 1-D models are the “Preliminary Reference Earth Model” (PREM) of Dziewonski and Anderson (1981) and model ak135 (Kennett et al., 1995). The latter is based on a more extensive dataset than the former, and is therefore more accurate. Research on 3-D Earth structure is a highly active field; recent reviews by Garnero (2000) and Nataf (2000) provide a useful starting point.

Lateral variation of elastic properties and density is greatest in the crust and uppermost mantle, but the density of broadband seismic observatories used for global seismology is far too small (especially in view of the non-uniform geographical distribution) to determine the lateral heterogeneity of the “crustal structure” (where this term includes both the crust and the uppermost mantle). Geophysicists must therefore use data collected from various local and regional surveys to correct for the effect of crustal structure so that their data can then be analyzed to determine 3-D Earth structure on a global scale. Two widely used models for this purpose are CRUST 5.1 (Mooney et al., 1998), which has a resolution of 5o (i.e., about 500 km), and its successor, Crust 2.0 (http://mahi.ucsd.edu/Gabi/rem.dir/crust/crust2.html), with a resolution of 2o. These models are not intended as accurate models of the crust, but rather are intended as “pretty good” models, for the purpose of removing crustal effects. Physicists planning neutrino beam experiments should exercise appropriate caution when using these models.

### 2.1 Density models

Both global and regional density models are subject to considerable uncertainty. Global scale density models are typically derived by applying an equation of state, which is an empirical approximation, to seismic velocity models (e.g., Bullen, 1975). Crustal density models are derived using a variety of empirical relations between seismic velocities and densities (see Mooney et al., 1998). It is striking that, especially for the case of sedimentary rocks, many of these empirical relations were published in the 1970s, which suggests that there has not recently been a high level of activity in this field.

It is difficult to quantify the uncertainty of published density models. One interesting approach is that of Kennett (1998). He exploited the fact that the frequencies of the longest period modes of the Earth’s free oscillations depend separately on the elastic constants and the density to a marginally resolvable extent to conduct the following numerical experiment. He fixed the seismic velocities and density of his Earth model to the values of the PREM model, and constructed a random ensemble of density models centered around the PREM model. He then calculated the free oscillation eigenfrequencies for each model and compared them to the observed eigenfrequencies to construct a set of the 50 best fitting models. These density models, all of which can be said to fit the free oscillation data acceptably, have a range of about $`\pm 23`$ per cent in the upper mantle. This should not be regarded as a conclusive error estimate, but it is one reasonable indication of the general level of uncertainty of present 1-D global density models.

## 3 Density models for neutrino beam experiments

Let us consider a hypothetical neutrino beam experiment (Fig. 1) with a neutrino source in Tokyo and a detector in Shanghai. Note that the neutrino beam follows a straight line, but a seismic wave from Tokyo to Shanghai (or vice versa) follows a curved path (the path of minimum travel time). Thus it is not possible to infer the physical properties of the neutrino beam path based only on observations of seismic waves traveling from Tokyo to Shanghai. Published 3-D Earth models, which were obtained by analyzing a large dataset using many sources and receivers, can be used to obtain a seismic velocity profile along the neutrino beam path, which can then be empirically converted to density. If the accuracy of the density profile obtained using the above procedure is deemed insufficient, further information could in principle be obtained by conducting a seismic observation campaign with receivers along the entire great circle from Tokyo to Shanghai. However, the fact that much of the beam path lies under the oceans would greatly complicate such a campaign.

Figure 2 shows the various density profiles under the hypothetical Tokyo-Shanghai path, taken from Model Crust 2.0. As shown in Fig. 2, the variation between the various density profiles is $`\pm 4`$ per cent in the depth range from 10–20 km and $`\pm 7`$ per cent in the depth range from 20–30 km. The variations in density are due to the differences in the physical properties of the various types of geological units, but can also be regarded as a crude indicator of the general level of uncertainty of the density. As, generally speaking, the amplitude of the Earth’s lateral hetereogeneity decreases with increasing depth, the variability of $`\pm 7`$ per cent in Fig. 2 can reasonably be regarded as as an upper bound on the uncertainty. Note that the density in the depth range 20–30 km in the rightmost column of Fig. 2 (3.35 g/cm3) is the value for the uppermost mantle, and is about 10 per cent higher than the density of the lowermost crust.

## 4 Discussion

Neutrino beam physicists should be aware of the various uncertainties and limitations of present geophysical knowledge of the Earth’s density distribution, as discussed in this paper. The planning of neutrino beam experiments should include simulation of the data reduction process, including a propagation of error analysis, to study the effect of this uncertainty. Three possible scenarios can be envisioned. (1) The uncertainty of present density models poses no significant problems; (2) moderate reduction of the uncertainty, through more detailed analysis of existing data, is required: (3) significant reduction of this uncertainty, by conducting a large scale campaign of geophysical observations, is required. Obviously, scenario (1) would be most desirable, while scenario (3) would be discouraging. This issue should be resolved at an early stage of the planning of neutrino beam experiments.

## References

Bullen, K., The Earth’s Density, Chapman & Hall (London, 1975).

Dziewonski, A. M., & Anderson, D. L., Preliminary reference Earth model, 1981, Phys. Earth Planet. Int., 25, 297-356.

Garnero, E. J., Heterogeneity of the lowermost mantle, 2000, Ann. Rev. Earth Planet. Sci., 28, 509-537.

Kennett, B. L. N., On the density distribution within the Earth, 1998, Geophys. J. Int., 132, 374-382.

Kennett, B. L. N., Engdahl, E. R., & Buland, R., Constraints on seismic velocities in the Earth from travel-times, 1995, Geophys. J. Int., 122, 108-124.

Lay, T., and T. C. Wallace, Modern Global Seismology, Academic (San Diego, 1995).

Mooney, W. D., Laske, G., & Masters, T. G., Crust 5.1: A global crustal model at $`5^o\times 5^o`$, 1998, J. Geophys. Res., 103, 727-747.

Nataf, H.-C., Seismic imaging of mantle plumes, 2000, Ann. Rev. Earth Planet. Sci., 28, 391-417.

Shearer, P. M., Introduction to Seismology, Cambridge U. Press (Cambridge, 1999).


text:
# High Resolution Polarimetry of the Inner Galaxy

## Introduction

The Milky Way was the first celestial radio source discovered, and was subsequently one of the first sources to be detected in linear polarization. There are two sources of this polarization: discrete objects such as supernova remnants (SNRs), and a diffuse polarized background produced by the relativistic component of the interstellar medium (ISM). All of this emission undergoes Faraday rotation as it propagates towards us, either in the source itself or in intervening material. With sufficiently high angular and frequency resolution, we can use the properties of this polarized emission to map out the distribution of ionized gas and magnetic fields in individual sources and in the ambient ISM. Only recently have instruments and techniques advanced to a point where such studies are feasible gld+99 ; hkd00 ; wdj+93 .

Motivated by the spectacular single-dish polarization surveys of Duncan et al dhjs97 ; drrf99 , we have made polarimetric images of the entire fourth quadrant of the Galaxy with the Australia Telescope Compact Array (ATCA). These data have been taken as part of the Southern Galactic Plane Survey (SGPS; mgd+01 ). While the primary focus of the SGPS is to study the Galactic distribution of H i, the ATCA simultaneously receives full polarimetric continuum data, which have allowed us to map out the distribution of linearly polarized emission in the survey region.

While the full survey has now been completed, a detailed analysis has only been carried out on a 28-deg2 test region, covering the range $`325.^{}5<l<332.^{}5`$, $`0.^{}5<b<+3.^{}5`$. We here summarize the main results of this analysis; this work is described in more detail by Gaensler et al gdm+00 .

## Observations and Reduction

The ATCA is a 6-element synthesis telescope, located near Narrabri, NSW, Australia. Observations of the test region of the SGPS were carried out in nine observing runs in 1997 and 1998, and comprised 190 separate telescope pointings (see mgd+01  for details). Data were recorded in nine spectral channels spread across 96 MHz of bandwidth and centered at a frequency of 1384 MHz. The two sources MRC B1438–481 and MRC B1613–586 were observed over a wide range in parallactic angle in order to solve for the instrumental polarization characteristics of each antenna shb96 . For each spectral channel, images of the field in Stokes $`I`$, $`Q`$, $`U`$ and $`V`$ were deconvolved jointly using the maximum entropy algorithm PMOSMEM sbd99  and then smoothed to a resolution of $``$1 arcmin. The final sensitivity in each image is $`\stackrel{<}{_{}}0.5`$ mJy beam-1.

Images of linearly polarized intensity, $`L=(Q^2+U^2)^{1/2}`$, linearly polarized position angle, $`\mathrm{\Theta }=\frac{1}{2}\mathrm{tan}^1(U/Q)`$, and uncertainty in position angle, $`\mathrm{\Delta }\mathrm{\Theta }=\sigma _{Q,U}/2L`$, were then formed from each pair of $`Q`$ and $`U`$ images. The nine $`L`$ maps (one per spectral channel) were then averaged together to make a final image of $`L`$ for the entire test region, while the nine $`\mathrm{\Theta }`$ and $`\mathrm{\Delta }\mathrm{\Theta }`$ maps were used to derive an image of the rotation measure (RM) over the field.

## Results

In Figure 1 we show images of $`I`$ and $`L`$ from both our 1.4-GHz ATCA observations and from the 2.4-GHz Parkes survey of Duncan et al dhjs97  (resolution $`10.^{}4`$). The total intensity images show the presence of SNRs and H ii regions (see mgd+01  for further discussion). Although the interferometric ATCA observations are not sensitive to the diffuse emission seen by Parkes, it is clear that the same features are present in both data-sets.

At first glance it seems that the $`L`$ images have very little in common with the Stokes $`I`$ emission. In particular, the ATCA $`L`$ image is dominated by diffuse polarization spread all over the field of view, composed of discrete patches separated by narrow “canals” of reduced polarization. While none of this emission is correlated with total intensity, there does seem to be a good match between the brightest polarized regions of the ATCA and Parkes data, despite the differing frequencies and resolutions of these data-sets. Using the images of $`\mathrm{\Theta }`$ and $`\mathrm{\Delta }\mathrm{\Theta }`$, we can determine the variation of polarization position angle with frequency wherever we detect polarized emission. The resulting RMs are generally small and negative, with a mean RM for the entire field of $`12.9\pm 0.1`$ rad m-2; 50% of the RMs have magnitudes smaller than $`\pm `$25 rad m-2 and 98% are smaller than $`\pm `$100 rad m-2.

The ATCA $`L`$ image reveals two large voids of reduced polarization, each elliptical and several degrees in extent. One void is centered on $`(l,b)=(332.^{}4,+1.^{}4)`$ (“void 1”) and the other on $`(328.^{}2,0.^{}5)`$ (“void 2”); both voids are also seen in the 2.4-GHz Parkes polarization map. The RMs around the edges of these voids range up to $`\pm 400`$ rad m-2, in distinction to the low RMs seen over the rest of the field.

A careful examination shows one marked correspondence between the ATCA Stokes $`I`$ and $`L`$ images: at $`(326.3,+0.8)`$, the bright H ii region RCW 94 shows reduced polarization towards its interior, and is further surrounded by a halo in which no polarization at all is seen. This is shown in more detail in Figure 2.

Finally, of the numerous unresolved sources distributed across the field, 21 of these sources show detectable linear polarization. The RMs for these sources fall in the range –1400 rad m-2 to +200 rad m-2.

## Discussion

### Diffuse Emission

We first note that the incomplete $`uv`$ coverage of an interferometer affects images of polarization in complicated ways. While it is physically required that $`LI`$, and we generally expect that structures seen in $`L`$ might correspond to similar structures in $`I`$, neither situation will be generally observed in interferometric data. This is because an interferometer can not detect structures larger than a certain size, corresponding to the closest spacings between its antenna elements (in the case of the ATCA, this maximum scale of $`35^{}`$). A source larger than this maximum scale will not be seen in Stokes I; if it is also a uniformly polarized source, it will not be detected in polarization either. However, magnetic field structure within the source, plus variations in the Faraday rotation along different lines-of-sight, can introduce power in Stokes $`Q`$ and $`U`$ on smaller scales, to which the interferometer is sensitive. We thus can often observe complicated structures in polarization which have no counterpart in total intensity gldt98 ; gld+99 ; hkd00 ; wdj+93 .

Clearly such an effect is occurring here, and is producing virtually all the linear polarization seen in Figure 1. We can crudely divide up the diffuse polarization we see into two components.

The brightest polarization seen with the ATCA matches well the bright polarized structures seen with Parkes. Since the amount of Faraday-induced polarization is very strongly dependent on both resolution and frequency, the fact that two such disparate data-sets show similar structures implies that these bright polarized structures are intrinsic to the emitting regions. By comparing the RMs observed for this emission to those observed for pulsars in this part of the sky, we can conclude that the distance to this emission is in the range 1.3–4.5 kpc. The depolarizing effects of RCW 94 (discussed further below) imply that the polarized emission is $`>`$3 kpc distant, while the lack of depolarization against other H ii regions gives an upper limit of $`6.5`$ kpc. Dickey dic97  has made H i absorption measurements against this emission to derive a lower limit on its distance of 2 kpc. Taking into account all these constraints, we argue that the mean distance to the source of polarized emission is $`3.5\pm 1.0`$ kpc, corresponding to the Crux spiral arm of our Galaxy.

The rest of the ATCA field is filled with fainter diffuse polarization, which does not have any counterpart in the Parkes data. This emission is best explained as being due to Faraday rotation in foreground material. The RMs measured for this emission imply that they are caused by foreground clouds of RM $`5`$ rad m-2, consistent with the conclusions made by Wieringa et al wdj+93 .

### Voids in Polarization

To the best of our knowledge, voids in polarization such as those described here have not been previously reported. There are two possible explanations to account for these structures: either they represent regions where the level of intrinsic polarization is low, or they are the result of propagation through a foreground object, whose properties have depolarized the emission at both 1.4 and 2.4 GHz.

If the voids are intrinsic to the emitting regions, then the distance of 3.5 kpc inferred above implies that they are hundreds of parsecs across — it is hard to see what could produce such uniformly low polarized intensity across such large regions. We thus think it unlikely that the voids are intrinsic to the emitting regions.

We thus favor the possibility that the voids are caused by depolarizing effects in foreground material. We have considered in detail the various ways in which foreground Faraday rotation can produce the observed structure, and can rule out bandwidth and gradient depolarization as possible mechanisms (see gdm+00  for details).

The only remaining possibility is that depolarization in the voids is due to beam depolarization, in which the RM varies randomly on small scales. We have developed a detailed model for “void 1” to confirm this. We consider void 1 to be a caused by a sphere of uniform electron density $`n_e`$ cm-3, centered on $`(332.^{}5,+1.^{}2)`$ with a radius of $`1.^{}4`$ and at a distance to us of $`d`$ kpc. Within the sphere, we suppose that there are random and ordered components to the magnetic field, and that these two components have identical amplitudes $`B`$ $`\mu `$G. The ordered component is uniformly oriented at an angle $`\theta `$ to the line of sight. We assume that the random component is coherent within individual cells of size $`l`$ pc, but that the orientation from cell to cell is random. Uniformly polarized rays which propagate through a different series of cells will experience differing levels of Faraday rotation, resulting in beam depolarization when averaged over many different paths.

By calculating the properties of the polarized signal which emerges after propagating through this source, we find that we can account for the observed properties of void 1 if $`n_e20`$ cm-3, $`B5`$ $`\mu `$G, $`\theta \stackrel{>}{_{}}80^{}`$, $`d300`$ pc and $`l0.2`$ pc (see gdm+00  for details). These properties are consistent with those of an H ii region of comparatively low emission measure. Indeed Figure 3 demonstrates that H$`\alpha `$ emission fills void 1, its morphology and perimeter matching exactly to that of the void. It is interesting to note that the O9V star HD 144695 is very close to the projected center of void 1, and is at a distance of $`300\pm 160`$ pc. The radius of the Strömgren sphere which this star would produce is consistent with the extent of the void. It is thus reasonable to propose that the star is embedded in and powers the surrounding ionized bubble.

Two properties of the voids which our simple model cannot account for are the requirement that the uniform component of the magnetic field be largely oriented in the plane of the sky, but that we generally observe coherent regions of large RM (of the order of a few hundred rad m-2) around the edges of the voids. We suggest that both these results can be explained by the field geometry which arises during the expansion phase of an H ii region as it interacts with surrounding material. This produces a magnetic field perpendicular to the line of sight over most of the void, but which is parallel to the line of sight (and can thus potentially produce high RMs) around the perimeter.

### Depolarization seen towards RCW 94

The reduced polarization seen coincident with RCW 94 in Figure 2 presumably results from beam depolarization, just as for the H ii region argued to produce void 1. However, the effects of beam depolarization are expected to be weakest around the edges of the source, and thus cannot account for the halo of complete depolarization surrounding RCW 94. We rather account for this depolarization halo by requiring the electron density to be approximately constant across RCW 94, but to fall off rapidly beyond the boundaries of the source. This produces a sharp gradient in RM around the edges of the source, resulting in complete depolarization.

The presence of significant CO emission at the same position and systemic velocity as for RCW 94 bact89  suggests that the H ii region is interacting with a molecular cloud. This possibility is supported by H i observations of the region, which show that RCW 94 is embedded in a shell of H i emission, which is further surrounded by a ring of decreased H i emission mdg+00b . McClure-Griffiths et al mdg+00b ; mgd+01  argue that this structure in H i confirms that RCW 94 is embedded in a molecular cloud, the shell of emission resulting from H2 molecules dissociated by the H ii region, and the surrounding region of reduced H i corresponding to regions of undisturbed molecular material. Simulations of H ii regions evolving within molecular clouds (rtf95  and references therein) show that for certain forms of the density profile within the parent cloud, the shock driven into the cloud by the embedded expanding H ii region can produce a halo of partially ionized material around the latter’s perimeter, which would produce the fall-off in $`n_e`$ required to produce the depolarization halo observed.

### Point Sources

With the exception of one source known to be a pulsar, the polarized point sources in our field are presumably extragalactic, and their RMs thus probe the entire line-of-sight through the Galaxy. When combined with information from pulsar RMs, we can use these data to constrain the geometry of the overall Galactic magnetic field. So far we have compared the RMs in our test region to those expected for a bisymmetric spiral configuration, and have found that pitch angles in the lower end of the range allowed by pulsars ($`p4.5^{}`$) are favored dom+01 . We are in the process of carrying out a more detailed study using the RMs of 163 background sources from the entire SGPS, in which we are comparing these measurements to the distributions expected for a wider variety of geometries and model parameters (Bizunok et al, in preparation).

## Conclusions

The ATCA’s sensitivity, spatial resolution and spectral flexibility have allowed us to study linear polarization and Faraday rotation from the inner Galaxy in an unprecedented detail. Even though the test region we have considered covers less than 7% of the full survey, we have been able to identify a variety of distinct polarimetric phenomena, and have used these to map out both global and turbulent structures in the magneto-ionized ISM. We anticipate that our analysis of the full SGPS will result in a comprehensive study of magnetic fields and turbulence in the inner Galaxy.

###### Acknowledgements.

The Australia Telescope is funded by the Commonwealth of Australia for operation as a National Facility managed by CSIRO. H$`\alpha `$ data were taken from the Southern H-Alpha Sky Survey Atlas (SHASSA), which is supported by the National Science Foundation. B.M.G. is supported by a Clay Fellowship awarded by the Harvard-Smithsonian Center for Astrophysics, while J.M.D. acknowledges the support of NSF grant AST-9732695 to the University of Minnesota.


text:
# o-Positronium scattering off H and He

## Abstract

Exploiting an approach similar to the R-matrix theory, the diffusion Monte Carlo method is employed to compute phase shifts and threshold cross sections for the elastic scattering of o-positronium off light atoms. Results are obtained for Ps-H and Ps-He as representative cases of open and closed shell targets. The method allows for an exact treatment of both correlation and exchange interactions, and represents the most promising approach to deal with these effects in more complicated targets. In particular the Ps-He threshold cross section, computed in a many body framework for the first time, represents a standard by which past and future numerical and experimental estimates can be judged.

Positronium (Ps) scattering off atomic and molecular targets has an overwhelming importance if an understanding of the interaction mechanism between an overthermal Ps and a condensed matter environment is required moge . For instance, by means of elastic and inelastic cross sections, it may be possible to model energy transfers from Ps to the surroundings or to describe the trapping of Ps in a free volume cavity. Despite its long history Drach1 ; Drach2 ; Page ; Drach3 ; Fraser ; Canter , and even in the case of light atoms, some quantitative aspects of the process still remain controversial and have recently been addressed by a number of authors, both experimentally Colema ; Nagash ; Garmer ; Skalsey  and theoretically Campb ; Adhik ; Biswas ; Blackw ; Jim . From the computational point of view, the difficulties with which almost every method is faced are related to the composite nature of both target and projectile. As a consequence, sensible results can be obtained only if correlation and exchange effects are properly treated. Moreover, an accurate description of the correlation effects between the target electrons and the positron in the Ps is important in computing ”pick off” annihilation rates. These effects have only recently been computed consistently for the case of positronium scattering off hydrogen and positronium atoms Jim ; Shumw . However, the full ab initio treatment (i.e. without the use of exchange or correlation model potentials) of systems with more than two electrons still represents a formidable task. A glance to the recent literature on bound systems containing a positron reveals a practically identical situation with only small numbers of electrons treated explicitly. In this context, it has been shown by several authors jiang ; maxcor  that flexible and accurate computational techniques for small and medium size systems are provided by the family of quantum Monte Carlo (QMC) methods. Among them, the diffusion Monte Carlo (DMC) scheme represents the most powerful approach to study strongly correlated systems thanks to its ability to sample a distribution proportional to the exact ground state wave function of a given Hamiltonian. Where fermions are concerned, the antisymmetric nature of the wave function and its consequent non-positiveness imply the introduction of a bias known as nodal error. For a given energy $`ϵ`$, the nodal error $`\mathrm{\Delta }ϵ`$, which disappears if the nodal surfaces of the exact wave function are known, has a value which commonly spans the range $`\mathrm{\Delta }ϵ/ϵ[10^5,10^4]`$ arnerev . Unless otherwise specified, the following results have been computed in the DMC framework.

The application of QMC methods to scattering problems was independently proposed in the eighties in two pioneering papers by two groups in the field of nuclear physics koonin ; Carlson . Their ideas have been recently applied to the exciton-exciton scattering problem Shumw , thus providing the first accurate calculation for the Ps-Ps system. This approach, which closely resembles the original idea behind the R-matrix theory of Wigner wigner , is briefly summarized in what follows for the case of an elastic collision. We define $`𝐫_{AB}=𝐑_A𝐑_B`$ as the relative position of the centers of mass of $`A`$ and $`B`$ and $`p`$ and $`\mu `$ as their relative momentum and reduced mass. The choice of a boundary $`r_{AB}=`$ under the condition $`V_{int}()\frac{p^2}{2\mu }`$, allows one to approximate the wave function in the region $`r_{AB}>`$ with the asymptotically exact form

$$\mathrm{\Psi }=𝒜\left[\mathrm{\Psi }_A(𝐬_A)\mathrm{\Psi }_B(𝐬_B)\frac{\mathrm{\Phi }_l(r_{AB})}{r_{AB}}Y_{lm}\right]$$
(1)
$`𝒜`$ is the antisymmetrization operator, $`𝐬_A`$ and $`𝐬_B`$ the internal coordinates of the two separate fragments, $`\mathrm{\Psi }_A`$ and $`\mathrm{\Psi }_B`$ their wave functions and $`Y_{lm}`$ and $`\mathrm{\Phi }_l`$ the angular and radial functions describing the dynamics of the relative motion of the two centers of mass.

The stationary form of $`\mathrm{\Phi }_l`$ can be expressed as

$$\mathrm{\Phi }_l(𝐫_{AB})=_l(pr_{AB})+S_l(p)𝒪_l(pr_{AB})$$
(2)
where $`_l`$ and $`𝒪_l`$ are Hankel functions, and $`S_l(p)`$ is the scattering matrix. Here, $`p`$ is connected to the total energy by $`E_{tot}=\frac{p^2}{2\mu }+E_A+E_B`$ where the computation of $`E_A`$ and $`E_B`$, both being bound state energies, can be exactly performed employing the DMC method. The condition $`\frac{\mathrm{\Phi }^{}()}{\mathrm{\Phi }()}=`$ and a corresponding value of $`p`$ (i.e. $`E_{tot}`$) are enough to precisely state the value of $`S_l(p)`$ in Eq.(2). The link between $``$ and $`p`$ is contained in the dynamics of the interior region and therefore the computation of $`E_{tot}`$ in this region, under the same boundary condition on $`r_{AB}=`$ (i.e. imposing the continuity of the logarithmic derivative), provides the value of $`S_l`$. Eventually the scattering matrix assumes the form

$$S_l(p)=\frac{_l^{}(p)_l(p)}{𝒪_l^{}(p)𝒪_l(p)}$$
(3)
Hence, as long as the boundary condition is exactly controlled, all the relevant information is contained in the interior region. In this subspace the wave function is normalizable and the energy, now parametrically dependent on $``$, can be studied by means of one of the QMC techniques e.g. VMC, DMC or one of their variants. In order to easily fulfill the above boundary condition in the DMC framework, one can choose the value of $``$ to be infinite. From a physical point of view, this choice corresponds to add a rigid wall located at $``$ along the distance between the target and projectile centers of mass. Results will be presented in terms of phase shift $`\delta _l(p)`$ which is defined by $`S_l=e^{2i\delta _l}`$ and can be expressed, starting from Eq.(3) and taking the $`\mathrm{}`$ limit, as

$$\mathrm{tan}\delta _l(p)=\frac{j_l(p)}{n_l(p)}$$
(4)
where $`j_l`$ and $`n_l`$ are respectively the spherical Bessel and Neumann functions.

Before going on we feel that it is worth stressing two important points. First, one has to satisfy the condition $`V\frac{p^2}{2\mu }`$, so an upper limit to the sphere radius $``$ does not exist, while it cannot be chosen smaller than some unfortunately not well specified threshold value. This imposes an upper limit to the relative kinetic energy. For Ps scattering off neutral atoms, the interaction potential between the target and the projectile dies off as $`1/r_{AB}^6`$. This short range potential allows the use of fairly small radii, a possibility not necessarily available for different colliding fragments. Secondly, since DMC samples the lowest energy state, it cannot be applied to scattering problems in presence of a bound state.

Both these two points highlight the importance that studying excited states could have within this approach. The possibility of raising the energy whilst keeping the surface constrain fixed can, in principle, allow the study of every system at any energy.

In this work, we applied the presented technique to the S-wave scattering of positronium off hydrogen and helium dynamically described by the full Hamiltonian

$$H=\frac{1}{2}\underset{i=1}{\overset{N_e}{}}_i^2\frac{1}{2}_p^2\underset{i=1}{\overset{N_e}{}}\frac{Z}{r_i}+\frac{Z}{r_p}+\underset{i>j}{}\frac{1}{r_{ij}}\underset{i}{\overset{N_e}{}}\frac{1}{r_{ip}}$$
(5)
where $`i`$ and $`j`$ refer to electrons, $`p`$ to the positron and $`Z`$ to the nuclear charge of the atom. The spatial part of the scattering wave function has been chosen to be of the form

$$\mathrm{\Psi }=𝒪[\mathrm{\Psi }_A(s_A)\mathrm{\Psi }_{Ps}(r_{1p})\frac{\mathrm{\Phi }(r_{PsA})}{r_{PsA}}\varphi _J(s_I)]$$
(6)
where $`\mathrm{\Psi }_A`$, $`\mathrm{\Psi }_{Ps}`$, and $`\mathrm{\Phi }`$ have the same meaning as in Eq. (1). $`\varphi _J`$ is a Jastrow factor for all the pairs of particles belonging to different fragments, $`s_I`$ is the set of distances for these pairs and $`𝒪`$ is the appropriate symmetry operator built according to Young diagrams. In the Ps-H case, the exact internal wave function of both fragments is known and $`𝒪`$ has the form

$$𝒪=1+(1)^SP_{12}$$
(7)
where $`S`$ (0 or 1) is the spin momentum of the state and $`P_{12}`$ the permutation operator between the two electrons. The $`S=0`$ space part of the wave function is everywhere positive, while the nodal surface for the $`S=1`$ state is exactly provided by the action of $`𝒪`$. Under this condition the energy can be computed by DMC without any approximation. The singlet state supports a bound state and, as said above, it is thus necessary to exploit an excited state technique. Evidently, in order to get the required scattering information, the boundary condition at the surface for every state must be controlled. As shown in Ref Shumw  the Correlation Function DMC methodCepBernu  with our choice of $``$ accomplishes automatically this requirement and it will therefore be used in the following. Detailed descriptions of this method are out of the scope of the present work and can be extensively found in the literature CepBernu ; CFQMC . In this respect, we limit ourselves to show the convergence for one calculation in Figure 1. All the simulations for the triplet state of Ps-H were carried out using a time step of 0.01 hartree-1, 2000 walkers, and a total of 100 blocks of 10000 steps each. Simulations for the singlet states were performed employing 2000 configurations, a time step of 0.01 hartree-1, and a grand total of 10000 decorrelated Euclidean time evolutions. Low energy phase shifts for both $`S=0`$ and $`S=1`$ systems are shown in Figure 2. Scattering lengths are reported in Table 1, together with three independent estimates of the same quantities very recently calculated Campb ; Adhik ; Jim .

Whereas all of them agree in assigning the value for the triplet state, there appears to be some controversy where the singlet state is concerned. Our results, which we believe to be statistically exact, are very close to the values proposed in Ref Jim  suggesting these could be safely considered as a definitive estimate. This fact can be also taken as a strong proof of the reliability of the method we are employing, as well as of the full-electron Stochastic Variational Method approach used in Ref Jim .

With this premise, we now address the more debated problem of positronium scattering off helium. Before discussing our computed quantities for this process, it is worth noting that the experimental measurements of the threshold value of the cross section span almost an entire order of magnitude Skalsey ; Nagash . The most recent theoretical estimates, obtained by different computational schemes, do not single out one of these as the correct one. The primary reason for this failure is the small size of the cross section and a consequent large fractional error associated with any approximation. To make a more direct comparison with experiments, Table 2 shows the most recently calculated and measured threshold cross sections.

In the present study, the system is treated with a genuinely many-body technique and no physical approximation have been made prior to the numerical simulation.

The absolute freedom one has in choosing the analytical form of the wave function in QMC methods allows us to employ the following explicitly correlated form for $`\mathrm{\Psi }_{He}`$

$$\mathrm{\Psi }_{He}=\mathrm{exp}\left(\frac{\alpha _1r_1+\alpha _2r_1^2}{1+\alpha _3r_1}+\frac{\beta _1r_2+\beta _2r_2^2}{1+\beta _3r_2}+\frac{\gamma _1r_{12}}{1+\gamma _2r_{12}}\right)$$
(8)
which has a DMC energy statistically exact. Moreover, the only Young diagram compatible with the choice of an helium atom in its ground state ($`S=0`$) gives the following form for $`𝒪`$:

$$𝒪=(1+P_{12})(1P_{13})$$
(9)
Simulations for this system were characterized by a time step of 0.005 hartree-1, 4000 walkers, and a total of 130 blocks of 25000 steps each. The numerical results for the phase shift are shown in Figure 3.

The value of the computed scattering length is 1.4046(6) bohr with a corresponding threshold cross section of $`7.892(2)\pi a_0^2`$. Comparison with numerical estimates and experimental results curiously shows this value to lie in a ’neutral’ zone, intermediate amongst the most recently proposed values (Table 2). The nodal error, being the only approximation introduced, deserves some comments. General considerations Chiesa  show this bias on the phase shift to be always negative and proportional to $`V^{1/3}`$ where $`V`$ is the sphere volume. As a result of this, our scattering length could be slightly lower than the exact one. More quantitatively, one can observe that in the interaction region the employed function closely resembles the functional form used in bound state calculations on similar systems, for which the nodal error roughly equals $`10^5`$ hartree maxcor . In the rest of the simulation volume the nodes of the trial wave function are practically exact thanks to the validity of Eq. 1. For such a reason we expect a bias on the energy of the same order of magnitude as the one in bound state calculations. If so, the nodal error would turn out to be smaller than the statistical fluctuations of our energy values, therefore warranting the statistical exactness of our results.

Among the numerous applications directly derivable from this method, we would like to emphasize that the possibility to sample the exact particle distributions in configurational space may allow one to obtain an effective interaction potential between Ps and a given atom or molecule. This potential, where all the physical effects are correctly accounted for, could be successively used to simulate Ps in condensed phases as molecular crystals and liquids. Moreover, it could also help in defining the preferential spatial location where the Ps positron would annihilate during a ”pick off” annihilation event, so that the interplay between the theoretical and the experimental results may enhance the diagnostic role played by Ps in condensed matter science. The DMC method is also suitable, as demonstrated in a slightly different context Mella , for the computation of $`Z_{eff}`$, and its formal extension to reactive processes (inelastic are still excluded) has been known since the seminal work of Alhassid and Koonin koonin .

The authors are in debt to Dr. Jim Mitroy for many helpful comments and discussions on methodological issues and positronium physics. Financial support from the Universitá of Milano is also acknowledged.